闽江口硅藻-盐度转换函数建立

徐 斌，陈 敏，吴诗勇，周国华，王承涛，张爱梅，方 琦

（1.自然资源部第三海洋研究所，福建 厦门361005；2.江西省地质工程（集团）公司，江西 南昌330029；3.安徽理工大学地球与环境学院，安徽 淮南232001；4.福建省地质调查研究院，福建 福州350013）

硅藻是一类单细胞生物，只要有水的地方均能发现其踪迹，具有种类丰富、环境适应性敏感以及壳体不易破坏等特点，这使得硅藻可以用于研究地质历史时期的古环境状况。随着研究程度的不断加深以及计算机技术的革新，定量古环境研究已逐渐替代定性分析。目前有效的定量分析方法之一就是建立硅藻与环境间的转换函数，将硅藻群落和特定的环境变量联系起来，利用加权平均计算属种的环境最适值和耐受值，用属种的最适值和耐受值去推导环境变量值。具体分析步骤可分为因子分析和回归模型建立，因子分析采用典型对应分析（Canonical Correspondence A-nalysis，CCA），这是目前相对最有效的直接梯度分析方法，主要作用是能将生物属种、环境因子和取样站位同时显示在一个低维的空间中，可直观揭示出属种排列分布和各环境指标的关系。回归模型包括非线性回归模型下的加权平均回归分析（Weighted Averaging Analysis，WA）和偏最小二乘法回归模型（Weighted Average-Partial Least Square，WA-PLS）。

目前转换函数在海洋、湖泊和河口均有应用。海洋中主要应用于海表温度和盐度，例如Sève（1999）在拉布拉多海建立了海表温度和盐度转换函数［1］。湖泊中的研究较为深入，国内外已建成一些完善的转换函数数据库。Charles（1985）对迪朗达克湖建立了pH、碱度、导电率、总磷等转换函数数据库［2］；董旭辉等（2006）建立了长江中下游湖泊硅藻-总磷转换函数数据库［3］；张慧清等（2013）、吴燕妮等（2016）以及李小平等（2012）分别利用已建立的长江中下游湖泊硅藻-总磷转换函数数据库研究了石塘湖、白马湖和淀山湖的古环境［4-6］。在河口区主要应用于盐度和海平面变化，在古盐度重建研究中发挥着重要的作用。Zong等（2010）在珠江口成功建立了硅藻-盐度转换函数［7］；庄陈程等（2014）建立了长江口现代潮滩硅藻海平面转换函数［8］。闽江是福建省最大的河流，重建器测时代之前闽江口水体盐度对研究闽江口海侵历史、海平面变化及淡水输入量变化等均具有重要的科学意义。目前关于闽江口古盐度重建的相关研究未见发表，本研究利用闽江口表层沉积硅藻数据建立了该河口区的硅藻-盐度转换函数。

1 材料与方法

1.1 研究区概述

闽江位于福建省东部、台湾海峡西岸，发源于武夷山脉，水系全长2 872 km，流经福建和浙江两省8个县、市，是福建省的母亲河。闽江正源为沙溪，由沙溪、建溪和富屯溪三大支流流经南平汇合为闽江干流，南平以上称之为闽江上游，建溪尤西段称为中游，水口以下称之为闽江下游。闽江河口大部分位于福州盆地，盆地属于典型断陷盆地，河流具有坡降较小，河面开阔，水流稳定等特点。闽江流经竹岐进入福州盆地后，到达淮安时受到南台岛阻隔分成南北支流，北支穿过福州市区至马尾称为北港，南支称为南港。南北港在马尾处汇合后折向北东，径直穿过闽安至亭江一线，到达亭江附近受到琅岐岛阻隔再度分为南北两汊，南汊途径梅花直奔东海称为梅花水道；北汊沿长门水道继续向北东方向前进，受粗芦岛、川石岛、壶江岛的阻断又分为乌猪水道、熨斗水道和川石水道汇入东海，其中川石水道是闽江主航道［9］。闽江口属于强潮型河口，据多年监测资料表明，河口最大潮差达到7.04 m，平均潮差为4.37 m，河口以内潮差迅速减弱，导致闽江感潮河段不长。一般认为枯季潮区界可达侯官，潮流界可至文山里附近；洪季潮区界至解放大桥附近，潮流界到达马尾附近。最近也有学者提出潮区界和潮流界均有上延，枯季可到达竹岐附近。

1.2 采样站位

本次研究共采集了23个表层沉积样品，具体采样位置见图1。23个表层沉积样品于2009年11月在对闽江口及其邻近海域进行多科学综合调查时采集，调查期间共设置了31个取样站位，利用沉积物抓斗进行表层沉积物采集，共在25个站位采集到沉积物样品，由于其中2个站位样品很少，被应用于其他研究，故选用剩余23个站位的表层沉积物进行硅藻分析研究，取样站位信息如表1所示。

图1 闽江口及其邻近海域采样站位分布Fig.1 Locations of sampling stations in Minjiang Estuary and its adjacent sea areas

表1 研究区站位采样情况Tab.1 Descriptions of sampling stations in the study area

续表1

1.3 硅藻分析方法

处理步骤主要包括前期处理、制片、鉴定和统计分析。大致流程为从样品中取10 g左右放入烘箱烘干4 h，称重后加入蒸馏水和双氧水浸润去除有机质，经超声波分散仪（频率80 Hz，时间2 min）分散后，倒入15μm孔径网筛过滤淘洗。经反复清洗后取上层富含硅藻悬浮液转移至15 cm3离心管中，静置24 h并离心，吸取20 mm3涂匀于盖玻片上，晾干后用加拿大树脂制成固定片。利用Olympus BX51光学显微镜进行硅藻薄片的观察和鉴定（物镜40×，目镜20×）。每个样品鉴定和统计硅藻300粒左右，不足300粒的样品以统计3张标准片为准。硅藻的分类鉴定依据参考金德祥（1965、1982、1992）、程兆第（1996）等所著的大量硅藻分类学书籍和文献［10-14］。

1.4 环境数据采集

因在近岸河口区硅藻分布往往与水深、盐度、沉积物粒径具有最大相关性［15］，故本研究选择水深、盐度、沉积物平均粒径（Mz）和标准偏差（σi）形成环境因子数据集。水深和盐度数据均取自于航次现场调查过程，其中水深为高潮时水深，采用水坨现场测定；盐度为日平均盐度，由WTW Vario Cond便携式盐度仪测定。样品粒度分析按照《海洋调查规范》［16］的规定在自然资源部第三海洋研究所完成。

1.5 数据处理

数据分析前，应先对硅藻数据和环境数据进行前期处理，我们选取相对百分含量大于2%的种类进行分析。本研究MJ1站位位于河口最外侧，沉积环境较为复杂，硅藻丰度极低，该站位仅含有少量多束圆筛藻（Coscinodiscus divisus），与最近的MJ2站位丰度差距巨大，与其他站位硅藻组成也相差甚远，MJ1站位又位于浙闽沿岸流影响范围，综合分析其硅藻分布可能受到沿岸流影响。为排除MJ1可能带来的误差，后期因子分析时将MJ1站位排除。由于在研究区所处的我国东部近岸水域及河流中硅藻季节性更替和变化非常大［17］，因此构建转换函数的盐度数据应选择当季的平均值（即2009年秋季）。而鉴于采样时日平均盐度值与该季平均值差异不会很大，因此选用日平均盐度来代替秋季平均盐度。环境数据虽量纲不同，但其值差异不大，在本次分析中不做标准化处理。

1.6 理论基础

1.6.1 因子分析 所有因子分析皆由Canoco 4.5软件完成，该软件包括4个核心模块。第一模块为数据转换模块Wcanolmp，其作用是将Excel表格中物种数据转换为软件可识别数据；第二模块为数据拼接及低频数据删除模块CanoMerge；第三模块为核心模块Canoco for window，所有的排序手段都集中于此模块，其作用是通过排序手段建立样方、物种和环境之间的关系；第四模块为成图模块CanoDraw for window，此模块的作用能将排序结果成图［18］。具体分析步骤如下：

①消拱对应分析（Detrended Correspondence A-nalysis，DCA）获取硅藻数据的单峰响应值，即4个轴梯度长度。若最大梯度长度大于2，则利用单峰模型下CCA分析硅藻与环境之间的关系。

②利用CCA分析环境数据的膨胀因子（Variance Inflation Factor，VIF），当某环境因子VIF＞20，说明该环境因子对硅藻分布影响很小，在因子分析过程中应将其删除。

③限制性典型对应分析（Constrained Canonical Correspondence Analysis，cCCA）可判别各环境变量对硅藻群落的影响程度。当某环境因子第一轴和第二轴特征值的比值λ1/λ2＞0.5时，则该环境因子可作为转换函数的目标因子。

1.6.2 模型建立 转换函数建立通过C2软件完成，软件由英国纽卡斯尔大学的Juggins（2007）研发［19］。目前运用最多的回归模型为非线性模型，非线性模型包括WA和WA-PLS。WA方法根据不同的还原方式和对属种忍耐值是否降权进行分类，可分为传统加权平均（WA-Cla）、反向加权平均（WAInv）、降权的传统加权平均（WAtol-Cla）和反向加权平均（WAtol-Inv）［20］。WA-PLS方法是最小二乘方法与反向还原回归方法的结合，考虑硅藻数据的残差结构，通过每一次残差结构组分的不断提取进而提高硅藻属种的环境适宜值计算的准确性。评价模型的推导能力，主要根据其均方根方差（RMSE）及相关系数（R2）进行判断。

2 结果与讨论

2.1 硅藻种属组成及总丰度

23个表层沉积样品中共鉴定出硅藻24属59种，其优势种（相对百分含量＞10%）主要以咸水或半咸水浮游种为主，包括弓束圆筛藻（Coscinodiscus curvatulus）、多束圆筛藻（Coscinodiscus divisus）、琼氏圆筛 藻 （Coscinodiscus jonesianus）、辐射圆筛藻（Coscinodiscus radiatus）等（表2）。次优势种（5%＜相对百分含量＜10%）主要包括星脐圆筛藻（Coscinodiscus asteromphalus）、洛氏圆筛藻（Coscinodiscus rothii）和细弱圆筛藻（Coscinodiscus subtilis）等。MJ3站位以淡水种为主，主要包括海氏窗纹藻（Epithemia hyndmanii）、黄埔水链藻（Hydrosera whampoensis）和软双菱藻（Surirella tenera）。

表2 研究区表层沉积硅藻优势种Tab.2 Surface sediment diatom dominant species in the study area

通过表层沉积硅藻总丰度分布图（图2）可知，河口段包括MJ1～5站位，其丰度值变化区间为13～49 906个/g，最大丰度值出现在MJ2站位，MJ1站位丰度值最小；长门至马尾段包括MJ6～8这3个站位，其丰度值变化范围为544～78 356个/g，最大丰度值位于MJ6站位，最小丰度值位于MJ8站位；闽江北港共6个站位，最大丰度位于MJ11站位，丰度值为114 886个/g，MJ12站位丰度值仅为2 801个/g；南港包括MJ17～23这7个站位，最大丰度值位于MJ21站位，丰度值为87 904个/g，最小值为239个/g；竹岐至侯官段MJ15和MJ16站位丰度值分别为143个/g和40 330个/g。

2.2 因子分析

DCA分析结果显示4个轴中最大梯度长度为2.180，值大于2，故选用单峰模型下的CCA分析。下文依次进行CCA、cCCA和DCCA分析。

表层站位、硅藻与环境因子的CCA排序结果如图3所示，盐度和水深两个环境变量线段较长，且盐度变量与横轴夹角呈锐角，因此横轴应与盐度相关性最好，水深变量与纵轴夹角呈锐角，故认为控制研究区硅藻分布的主要因素是盐度和水深。

图2 研究区表层沉积硅藻总丰度分布Fig.2 Distribution pattern of total surface sediment diatoms in the study area

图3 表层采样站位与环境因子间CCA排序图Fig.3 CCA ranking for stations and environmental factors

排序图中站位与环境变量的关系，根据沈林南等（2014）的相关研究［21］可知，站位在环境变量正延长线上投影长度表示其正相关性大小，负延长线上投影表示该站位与环境变量呈负相关性。由此可以将研究区划分为4个区。Ⅰ区包括MJ2、3、4等3个站位，与盐度有较好的正相关性，区内主要含有爱式辐环藻（Actinocyclus ehrenbergii）、波状辐裥藻（Actinoptychus undulatus）和离心列海链藻（Thalassiosira excentrica）等近岸咸水种，其中MJ3站位还含有海氏窗纹藻、黄埔水链藻、软双菱藻等淡水种，各淡水种相对百分含量分别为7.14%、7.14%、7.14%，且这几个站位分布在闽江口口门外，受近岸咸水影响强烈，代表受河口外沿岸水体影响的区域；Ⅱ区包括MJ5、6、7、15、17等5个站位，它们与水深及标准偏差成正相关关系，主要含有星脐圆筛藻、圆筛藻（Coscinodiscus sp.）等广温性远洋种，反映了这些原本应该存在于近海的外洋性种类是通过潮汐作用沿着河口潮汐通道和河道上溯而沉积在闽江口门内的，其中MJ15站位位于闽江竹岐至候官段，距离口门较远，但其硅藻组成与其他4个站位相似，因此可能是天文大潮及风暴作用时潮水影响达到闽江该处河段，从而将近海的硅藻携带沉积至此处，代表受潮汐上溯海水入侵影响的区域；Ⅲ区包括MJ9、10、11、12、13、14站位，与水深及盐度均呈负相关关系，主要出现小眼圆筛藻（Coscinodiscus oculatus）、辐射圆筛藻和细弱圆筛藻等近岸咸水和半咸水种，该区可能受到海水入侵与径流的共同影响；Ⅳ区包括MJ8、16、18、19、20、21、22、23站位，与水深及盐度均呈负相关，且与平均粒径具有较好的相关性。区内出现了披针桥弯藻（Cymbella lanceolata）和纯脆杆藻（Fragilaria capusina）等淡水种，其相对百分含量分别为0.26%和2.42%，说明该区受到闽江径流影响。

硅藻与环境因子的CCA排序结果如图4所示，排序图中环境因子箭头越长代表该因子与种类数据的相关性越强，由此可见盐度和水深是研究区控制硅藻种类分布的主要环境因子。由CCA分析结果可知（表3），轴1和轴2的特征值分别为0.114和0.035，轴1是最大解释变量轴，盐度与轴1最接近，故轴1代表与盐度接近的解释变量；轴2则代表与水深和标准偏差相关的解释变量。轴1和轴2累计解释了81.9%的硅藻与环境因子间的关系，这也解释了硅藻属种与环境因子相关性极强。

图4 表层沉积硅藻与环境因子间CCA排序图Fig.4 CCA ranking for surface sediment diatoms and environmental factors

经cCCA分析后，盐度、水深、标准偏差、平均粒径的特征值 λ1/λ2值依次分别为0.524、0.107、0.056、0.068，其中盐度特征值比值大于0.5，其他均小于0.5，故盐度可作为转换函数的目标因子。

表3 CCA分析结果Tab.3 Results of CCA analysis

2.3 转换函数建立

DCCA分析结果显示，最大梯度长度为2.568，大于2.4个标准离差单元，因此使用非线性模型下的WA和WA-PLS建立硅藻-盐度转换函数更为合适［22］。前文已指出MJ1站位沉积环境复杂，可能受沿岸流等的影响，因此本次模型建立仅使用剩余22个站位数据进行模型构建。

转换函数建立通过异常样点删除和模型比选等手段最终确定。首先利用WA模型建立硅藻-盐度转换函数，模型成果表显示（表4），WA模型下的WA_Cla具有最小的R2，并且具有较小的RMSE和Max_Bias，故此模型中WA_Cla模型推导能力最强。但根据残差图可知（图5），MJ3站位残差最大，说明MJ3站位对模型建立影响加大，通过分析可知MJ3站位位于河口区，主要种类为近岸咸水种，但与I区的其他两个站位不同，其中还含有淡水硅藻，并且其淡水种组成在其他站位均无相似情况。因此我们认为MJ3站位淡水种可能是南下的浙闵沿岸流将敖江的淡水种硅藻携带至此处沉积。这种情况在北仑河口地区也出现过，沿岸流携带来的外地种类改变了河口区原本的硅藻组成［23］。因此在闽江硅藻-盐度转换函数模型构建时需要对MJ3站位硅藻种类数据进行修正，即删除该站位中的淡水种。

重新运行软件建立WA模型，分析结果如表4所示，WA_Inv模型的R2达到0.908 15，相比未修正MJ3站位硅藻数据时有显著提高，其中RMSE为2.341 7，较之前变化不大，Max_Bias为5.430 9，相比之下有所提升。综合来看，剔除MJ3站位淡水种后，模型推导能力显著提升。利用WA-PLS模型再次运行软件，WA-PLS模型分析结果见表4，此时Component 5的R2高达0.988 51，RMSE为0.637 9，Max_Bias为2.063 1，相比WA_Inv模型，R2提高了8.85%，RMSE降低了72.76%，Max_Bias降低了60.02%。整体而言，Component5模型推导能力最强，盐度拟合结果如图6所示。

表4 模型成果Tab.4 Results from modeling

图5 WA_Cla模型残差图Fig.5 Residual graph of WA_CLa model

根据Component 5模型拟合结果可以得出转换函数关系式表示为S拟=0.25+1.007S实。该转换函数可通过柱状硅藻数据进行古盐度反演，柱状硅藻属种必须选取表层硅藻训练集数据中出现过的物种并且相对百分含量大于2%的导入模型，软件通过对柱状硅藻数据与表层数据拟合，得出最为接近表层硅藻数据，此表层数据对应的实测盐度值带入转换函数即可得出预测盐度。

图6 Component 5模型盐度拟合图Fig.6 Salinity fitting diagram of Component 5 model

2.4 讨论

近年来，随着数理统计方法的不断发展，转换函数的应用逐步被应用开来，有孔虫、硅藻等海洋生物群体均有利用转换函数反演古环境。转换函数的精度控制目标为获得最低的RMSE值与最高的R2值，目前为提高转换函数精度方式主要包括模型比较和数据过滤。根据前人研究成果显示，转换函数模型包括WA、WA-PLS、ML、MAT等，以上模型为我们提供了更多的选择，但仍以WA和WA-PLS最受青睐，其实各模型差异性不大，并且在许多情况下，各模型预测性能基本无差异，只是通过比较RMSE值与R2值选取最优模型。通过已建立转换函数残差图剔除异常数这是本研究建立转换函数所选取模型的原因。本研究利用WA模型剔除了MJ1站位，再使用WA模型和WA-PLS模型重建时发现MJ3站位中淡水种同样影响转换函数精度，最终得出RMSE及R2均较为理想，与Zong等在珠江口所建立转换函数精度相当［7］。文中通过剔除站位提升精度的同时也减少了样本量，缩短了环境梯度宽度，盐度梯度宽度从27.6缩短为23.7，从图5中可以看出目前预测盐度值大于23的仅为MJ2站位，导致转换函数预测宽度减小，影响转换函数精度。针对此问题本研究采取前人普遍接受以20%切割点作为误差补偿［24］，本区如果柱状硅藻数据与现代训练集本集一致，盐度值以2.76作为平均误差作为补偿提升精度。

3 结论

（1）对23个表层沉积样品进行分析，共鉴定出硅藻24属59种，优势种有弓束圆筛藻、多束圆筛藻、琼氏圆筛藻、辐射圆筛藻、星脐圆筛藻、洛氏圆筛藻和细弱圆筛藻，硅藻丰度范围为13～49 906个/g。

（2）通过CCA分析将研究区表层沉积硅藻划分为4个区，Ⅰ区主要含有爱式辐环藻、波状辐裥藻和离心列海链藻等近岸咸水种，代表受河口外沿岸水体影响的区域；Ⅱ区主要含有星脐圆筛藻、圆筛藻属等广温性外洋种，代表受潮汐上溯海水入侵影响的区域；Ⅲ区主要出现小眼圆筛藻、辐射圆筛藻和细弱圆筛藻等近岸咸水种和半咸水种，该区可能受海水入侵与径流共同影响；Ⅳ区内出现了披针桥弯藻和纯脆杆藻等淡水种，说明该区受到闽江径流影响。

（3）通过WA-PLS模型下Component 5模型建立了转换函数，转换函数关系式为S拟=0.25+1.007S实。

（4）本研究转换函数精度较高，虽因样本量减少引发盐度梯度减小，但以2.76作为平均误差作为补偿提升转换函数精度。未来的研究中可以以该转换函数为基础，反演闽江口地区古盐度变化情况。若想建立一个能适用于闽江口地区的完整环境数据库，需完善闽江口其他环境数据以及样品量。

致谢：感谢厦门大学高爱国教授提供23个表层沉积物样品，感谢沈林南博士鉴定部分薄片。