Application of Environmentally Sensitive Factors in Bay Sedimentary Environments

海湾是被陆地环绕且面积不小于口门宽度为直径的半圆面积的海域^[1],海湾内往往沉积带有环境变化信息的泥质体,它保存了不同历史时期的气候变化、沉积环境和海陆变化信息,同时保存了人类活动的信号^[2⁃4]。因此,对海湾中的沉积记录信息的提取,在研究水陆交换、人文信息和气候历史等方面具有重要意义。

粒度是沉积物基本性质之一,是衡量环境变化的重要指标^[5]。它主要受水动力强弱、搬运方式、物源等因素的影响,因而沉积物粒度分布特征可以表征沉积物的物理性质及分类、沉积环境等方面的变化^[6]。在以往的研究中,学者们往往以沉积物的组分含量（黏土、粉砂、砂）以及粒度参数（分选系数、平均粒径、偏态、峰态）为标准揭示环境的变化^[7⁃9],但这些粒度指标只能近似地指代环境变化,为了更直观地研究受到多种物源和不同沉积环境影响的古环境,人们尝试采用近似的替代指标^[10]。近年来对于柱状样粒度资料进行数据处理,提取对沉积环境变化敏感的粒度组分或端元组分,以此重建地质历史时期的古环境、古气候已成为研究热点^[11⁃12]。为准确提取环境敏感因子,各国学者分别探讨了沉积物粒度的多组分分离法,如冷传旭等^[11]利用因子分析法提取东亚季风敏感粒级时将中心化对数比变换引入季风敏感粒级的提取过程是可行的,能有效地揭示更真实的中国东部地区东亚冬季风演化历史；Prins et al.^[13]使用Weibull分布的函数拟合法揭示了马克兰大陆坡沉积物是由浊积砂、浊积粉砂或风积砂及河流沉积物三个端元组成的混合物；滕珊等^[14]利用粒级—标准偏差法提取了沉积物的敏感性粒级,探讨了研究区沉积物的搬运机制及对东亚冬季风强度的响应；Yi et al.^[15]利用V-PCA揭示了粒度变化是东亚季风强度的指标,其强度主要受太阳日照和全球冰量的影响。

环境敏感因子在近海环境演化的研究中已取得丰硕成果^[13⁃15],但不同环境敏感因子提取方法在海湾沉积环境中的对比研究较少。不同环境敏感因子在海湾沉积环境中的适用性如何?不同方法反映东亚季风的方式是否一致或存在差异?海湾与相邻近海沉积提取的敏感粒级反映的东亚季风变化是否具有一致性?基于此,本文选取了BasEMMA、粒级—标准偏差法、V-PCA三种方法,提取威海湾WH-05岩心沉积物的环境敏感因子,研究沉积环境特征,探讨不同方法在海湾沉积环境中的适用性及对东亚季风的指示作用。研究成果有助于丰富环境敏感因子对季风事件响应研究成果,对于完善海湾沉积环境演化具有重要意义。

威海湾位于山东半岛东北部,濒临北黄海,赵北嘴和北山嘴之间,面积约54 km²,海岸线长约26 km,岸线稳定。水深一般介于6~9 m,最大水深13.5 m^[16]。研究区地处东亚季风区,冬季盛行西北风,夏季以东南风为主^[17]。威海湾东北为低温低盐的鲁北沿岸流,它在冬季风的驱动下从山东半岛北岸向东运动,之后越过成山头继续向南延伸^[18]。表层沉积物以粉砂、黏土质粉砂为主,中值粒径介于5.7~6.7 Φ^16]。注入威海湾的河流有李家夼河、庙琦河、望岛河、长峰河等,这些河流全是季节性河流,其径流水源主要来自大气降水,属于雨水补给型河流^[19]。威海湾得天独厚的自然条件和地理位置使湾内波浪和海流均较小,地形稳定。

2019年12月,中国海洋大学在威海湾通过两艘渔船（鲁V渔60611和鲁V渔60623）搭建的海上平台钻探获取了WH-05岩心（37°28′ N、122°10′ E）样品,取心率为60%,无沉积间断。水深6 m,岩心深度为18.2 m,钻孔位置见图1。

Figure 1.  Geographic location and borehole location of the study area (modified from Su^[20])

钻取的岩心按照2 cm间隔分样,粒度分析在中国海洋大学海底科学与探测技术教育部重点实验室完成。取0.5~1.0 g沉积物样品置于100 mL烧杯中,加5 mL 30%的过氧化氢溶液浸泡24 h,去除沉积物中的有机质；再加入10%的盐酸浸泡24 h,去除沉积物中的钙质胶结物等；后加入0.05 mol/L的六偏磷酸钠溶液以分散沉积物后上机测试。分析仪器采用Beckman Coulter LS 13 320激光衍射粒度分析仪（测量范围0.017~2 000 μm,分辨率0.1 Φ）进行测试,选取20%样品进行重复测试,重复测试误差≤1%。本文使用Folk et al.^[21]提出的方法进行沉积分类命名,采用McManus^[22]提出的矩法对粒度参数进行计算和分级。

BasEMMA^[23]使用VBA(Visual Basic for Applications)嵌入到Microsoft Excel程序中,使用Excel图形实时显示中间结果和最终结果。Zhang et al.^[23]建议BasEMMA应该执行两次。首先,利用“EM数从2到5”、“最大生成数=30”和“种群数=300”等参数进行粗略快速的分析。其次,在确定合适的EM数后,应使用“EM数从3到3”和“最大生成数=300”等参数进行全面而缓慢地分析。“粒度数”和“样本数”取决于所涉及的数据。

粒级—标准偏差法^[6,24]是通过计算每一粒级所对应含量的标准偏差以获得粒度组分的个数和分布范围,计算公式如下：

s=[∑i=1n(si-s)-2]/n

式中：s为偏差,s_i 为样本值,n为样本数。

V-PCA曾成功应用于中国渤海晚第四纪的沉积记录以重建亚洲季风的变化^[15],Darby et al.^[25]用此方法识别北极沉积物中的海冰与底流搬运作用。通过该方法可以获得粒度主成分,其可以反映对应于某种动力机制。粒级旋转分析原理是：将一系列具有相关性变量集（粒级）,通过正交分解,以获得不相关的、相对独立的变量集（粒级谱）。粒级旋转的目的是获得最大的方差贡献来代表整体样品^[26]。将柱状沉积物不同粒级的百分含量作为变量,通过SPSS 20对这些粒级进行粒级旋转主成分分析,以粒径为0.3~300 μm粒径谱的相关矩阵作为输入矩阵,获得沉积物粒度主控因子。该方法提取相对独立的粒度主成分,这些粒度主成分可以对应于一定的搬运机制。不同的搬运机制对应于不同的 V-PCA粒度谱。V-PCA粒级谱的众数为因子载荷最大的粒级,它最能代表粒级谱^[25]。

结合AMS¹⁴C测年结果^[27]（表1）,8.4 ka B.P.以来威海湾沉积物以粉砂为主,含量介于26.86%~97.54%,平均值为78.30%；黏土次之,含量介于2.45%~33.43%,平均值为19.27%；砂含量最低,含量介于0~65.68%,平均值为2.42%。平均粒径介于3.11~7.54 Φ；分选系数介于0.76~2.95,分选较差、差；偏态介于-1.28~2.64,为极负偏和极正偏；峰度值介于0.87~3.58,属尖锐（表2、图2）。根据粒度参数变化,将研究区8.4 ka B.P.以来地层可大致划分为5个沉积单元,具体分层情况如下。

Figure 2.  Grain size distribution for core WH⁃05 sedimentary formation (dating data from reference [27])

DU3层可以分为上下两个单元（DU3-1和DU3-2）。DU3-1层（16.37~17.76）,主要成分是粉砂（82.39%）,其次是黏土（17.09%）,砂的含量最少（0.52%）。平均粒径为6.32~7.30 Φ；分选介于0.77~1.66,分选中等、较差；偏态介于-0.97~1.14,为负偏到正偏；峰态介于0.88~2.09,为很窄到宽平（表2、图2）。粒度参数高频波动,指示水动力较强。频率分布曲线主要呈双峰分布,主峰介于7~8 Φ,次峰介于5~6 Φ；概率累积曲线主要呈三段式,以跳跃组分为主,含量介于10%~80%,其次是悬浮组分,含量在20%左右,滚动组分含量较少,不足10%。DU3-2层（17.76~18.2 m）,主要成分是粉砂（59.01%）,其次是砂（24.51%）,黏土含量最少（16.48%）。平均粒径为3.11~7.17；分选介于1.36~2.95,分选较差、差；偏态介于-1.23~2.64,偏态由负偏到极正偏；峰态介于1.82~3.58,为宽平到很宽平（表2、图2）。频率分布曲线主要为双峰曲线,主峰较宽平,次峰尖锐。主峰值介于6~8 Φ,次峰值介于4~5 Φ；概率累积曲线主要是三段式,以跳跃组分为主,含量介于20%~98%,其次是滚动组分,含量在20%左右,悬浮组分最少,不足2%。该层砂的含量逐渐增加,DU3-2层位砂的含量首次超过了黏土的含量,说明该层水动力较强,沉积环境动荡。

DU2层可以分上下两个单元（DU2-1和DU2-2）。DU2-1层（8.37~15.17 m）,主要成分是粉砂（81.37%）,其次是黏土（18.08%）,砂含量明显减少,只有0.54%。平均粒径为5.92~7.54 Φ；分选系数介于0.76~1.87,分选中等、较差；偏态介于-1.28~1.53,为负偏到正偏；峰度值介于0.87~2.49,为很窄到宽平（表2、图2）。此段粒度参数为高频波动,一方面是因为对沉积物样品进行高分辨率粒度测试；另一方面可能是因为柱状样位于海湾内,水深较浅,易受风暴沉积的影响。频率分布曲线（图3）呈单峰到多峰分布,主峰值介于5~6 Φ,次峰值介于8~9 Φ,主要成分为粉砂；概率累积曲线主要呈三段式,跳跃组分介于10%~70%,悬浮组分介于10%~30%,此外还有少量滚动组分,含量在10%左右,由顶到底沉积物悬浮组分逐渐增多,曲线斜率由顶到底逐渐增大,说明分选性越来越好,水动力自上而下水动力由弱变强。DU2-2层（15.17~16.37 m）,主要成分是粉砂（79.01%）,其次是黏土（20.38%）,砂的含量最少（0.61%）。平均粒径为6.13~7.29 Φ；分选介于0.83~1.66,分选中等、较差；偏态介于0.23~1.25,正偏；峰度介于0.96~2.14,为很窄到宽平（表2、图2）。概率分布曲线主要呈单峰特征,主峰尖锐,介于7~8 Φ,部分样品呈双峰特征,次峰的众数在4.2 Φ；概率累积曲线主要呈三段式,跳跃组分含量介于20%~80%,滚动组分含量在20%左右,悬浮组分含量在10%以下。该段整体水动力较强,沉积环境不稳定。该层整体属于水动力较强的环境,随着水深逐渐增大,粉砂和黏土含量增多,水动力相应地减弱。

Figure 3.  Frequency distributions and cumulative probabilities for typical samples in each sedimentary unit of core WH⁃05

DU1层（0~8.37m）沉积物以中细粉砂为主,平均含量为77.27%；其次为黏土,平均含量为21.90%,砂平均含量为0.83%。平均粒径为5.73~7.32 Φ；分选系数介于1.33~1.89,分选较差；偏态介于-1.08~1.55,为负偏到正偏；峰度值介于1.78~2.42,峰值宽平（表2、图2）。频率分布曲线主要为单峰曲线（图3）,曲线尖锐,峰值在5.5 Φ左右,主要是粉砂成分；概率累积曲线为典型的三段式,以跃移组分为主,含量介于30%~95%,滚动组分介于10%~20%,悬浮组分在5%之内。此段水动力较强,沉积环境不稳定。

对威海湾WH-05岩心按照2 cm间隔选取432个样品沉积物粒度进行端元模拟分析,获得了具有统计意义的五个粒度端元。所获得的五个端元组分,随其所解释的方差量的大小而变化。而拟合优度可以用决定系数（coefficient of determination,R²）所指示,它可以用于评估计算丰度与测试丰度之间的相似性,可以确定近似于整个数据变化的模拟最小端元组分个数^[26]。基本端元模型最终输出的是具有统计学意义并且符合现实地质状况的真实的粒度分布。而且,由于它和测量粒度分布具有一样的纲量,所以方便进行对比、分析以及端元主控因素判定^[28]。在取端元模拟组分个数2~5情况下（图4）,每一个粒级的决定系数（R²）的变化。当取端元数为2时,在26.7 μm处决定系数达到最低（图5a）；当端元数为3时,决定系数基本可以大于0.6。图5b是端元数取1~11时随着累计方差贡献的变化,当端元数取3时,累计方差贡献达到95.92%。Dietze et al.^[28]对端元模拟分析研究结果显示,在端元模拟分析过程中,所模拟端元的方差贡献和必须达到95%才是可靠的,决定系数到三个端元时达到0.6,此时可以确定三个粒度端元可以解释所有的样品（图4）。

Figure 4.  Distribution of the five end⁃members of the BasEMMA

Figure 5.  Determinant coefficients and variance contributions of end⁃member components

粒度端元为3时,体积含量峰值分别出现在7.4 μm、39.8 μm、493.6 μm,所对应的粒度组分范围是7.8~11.0 μm、31.3~44.2 μm、353.6~500.0 μm。EM1（7.8~11.0 μm）粒级沉积物含量最低,约为15%,粒级区间与概率累积曲线上悬浮组分对应,指示为悬浮运动；EM2（31.3~44.2 μm）含量较高,介于25%~30%,粒级区间与概率累积曲线跳跃组分第一段所对应的,其斜率较大,表现为较高速率的沉积,反映了该粒度组分的构成主要来自水动力较强的沉积环境中；EM3（353.6~500.0 μm）含量介于20%~25%,是最粗的沉积物,粒度区间与概率累积曲线图中跳跃、滚动组分均有对应,说明EM3是在水动力条件较为强烈的沉积环境中形成的。

通过粒级—标准偏差法对432个粒度样品提取敏感粒级,获得的每个粒级组分的标准偏差会随着粒级的变化而变化（图6）,图中较高标准偏差值反映了沉积物粒度含量在某一粒级范围内差异大,标准偏差的大小代表了对应粒级的百分含量随年代以及沉积环境改变的波动强度。对WH-05岩心而言,3个明显的标准偏差峰值分别出现在7.5 Φ（5.6 μm）、4.5 Φ（43.7 μm）、1.5 Φ（356.5 μm）,粒级组分区间是7.1~8.6 Φ（2.6~7.4 μm）、4.0~4.9 Φ（33.7~63.4 μm）、1.3~2.0 Φ（256.9~409.6 μm）。粒级—标准偏差变化曲线反映了不同层位的样品的粒度含量在某一粒级中。

Figure 6.  Grain⁃grade standard deviation curve

采用V-PCA法,从WH-05岩心中分离出4个粒度主成分：F1（6.2 μm,30.684%）、F2（18.9 μm,26.741%）、F3（0.8 μm,26.661%）、F4（14.3 μm,4.638%）,所解释的方差贡献和为88.723%（图7、表3）。

Figure 7.  V⁃PCA results for core WH⁃05

F1（6.2 μm,30.684%）粒度主成分解释30.684%的方差贡献,在6.2 μm处有一个宽广的正相关峰,在39.7 μm处有一个较狭窄的负相关峰。F2（18.9 μm,26.741%）粒度主成分解释26.741%的方差贡献,在18.9 μm处有一个狭窄的正相关峰,在194.2 μm处有一个起伏较低的负相关峰。F3（0.8 μm,26.661%）粒度主成分解释为26.661%的方差贡献,在0.8 μm处有一个宽广的正相关峰,在92.1 μm处有一个较狭窄的负相关峰。F4（14.3 μm,4.638%）粒度主成分解释为4.638%的方差贡献,在14.3 μm处有一个狭窄的正相关峰,在56.6 μm处有一个宽广的负相关峰。

根据F1（6.2 μm,30.684%）和F2（18.9 μm,26.741%）的方差贡献和粒度主成分众数,可以将其解释为WH-05岩心的主控因子,这两个主控因子基本反映了沉积物样品的总体粒度趋势。

粒级组分（GS）是通过V-PCA程序计算出的主成分F1-F4与其平方载荷之和求得,公式如下：

GS=（54.438*F1 + 18.437*F2 + 14.095*F3-3.965* F4）/100

GS序列表示东亚季风强度^[15]。当东亚冬季风增强时,水动力增强,粗粒物质增多,GS值减小；当东亚夏季风增强时,水动力减弱,细粒物质增多,GS值增大。

BasEMMA和粒级—标准偏差法在选用的粒度参数和计算过程上存在一定差异,但均在WH-05沉积地层中提取了3个敏感粒级且粒级范围基本一致,表明两种方法能反映相同的环境变化信息。下文以BasEMMA为代表,与V-PCA进行对比分析。

Yi et al.^[15]推断沉积物GS序列与东亚季风强度有关,东亚夏季风增强时,沉积物GS值相应增大,与本研究区分析结果较为一致。EM1/（EM1+EM2+EM3）表示的是BasEMMA算出的细粒物质在沉积物总体的含量,这两项指标虽然原理不同,但都可以表示沉积物中的细颗粒物质（图8b）。为探讨两种方法在同尺度下与东亚季风的强度对比,改绘了部分高质量的中国黄土纪录的堆叠归一化平均粒径（MGS）^[29]（图8a）,它通常被用作东亚季风强度替代指标。

Figure 8.  East Asian monsoon intensity surrogate index (MGS) trend and comparison between BasEMMA and V⁃PCA extraction results for sensitive fractionation components

BasEMMA与V-PCA提取的粒度端元随深度变化趋势基本一致,环境敏感组分具有较好的相关性,决定系数（R²）达0.781 2（图9）。两种方法在DU2-1和DU2-2相关性较高,决定系数分别是0.856 0和0.855 8；其次是DU3,决定系数达0.772 1；DU1相关性较低,只有0.670 8（图10）。

Figure 9.  Correlation between BasEMMA and V⁃PCA extraction of sensitive fractions

Figure 10.  Correlation comparison of environmentally sensitive factors extracted by BasEMMA and V⁃PCA in different sedimentary units

对比中国沿海其他岩心关于研究敏感粒级的成果,发现粒级—标准偏差法在渤海、黄海和东海沉积环境研究中的应用较为普遍^[14,30]（图11、表4）；V-PCA主要用于分析某种动力机制的粒度主成分^[15],应用相对较少；BasEMMA作为一种新方法^[23],尚未得到广泛应用。以上三种方法在研究区应用效果较好,说明这三种方法在海湾沉积环境演化研究中均具有适用性。

Figure 11.  Core locations in Bohai Sea, Yellow Sea and East China Sea (borehole locations from references [14⁃15,23,29])

为探讨WH-05在山东半岛泥质区内对东亚季风的指示作用以及两种方法在此研究区内的指示意义,参考了前人文献中在研究区附近的岩心所展现出的东亚季风强度变化^[31⁃33]并改绘部分高质量的中国黄土纪录的堆叠归一化平均粒径^[29]（图8a）。虽然BasEMMA和V-PCA是从不同的角度提取敏感粒级因子,两者在原理和方法上有巨大的差异,但在8.4 ka B.P.以来东亚冬季风的变化上,二者在总体上都呈先减小再增大的趋势。

通过对比发现,在8.4~7.5 ka B.P.阶段,对应沉积单元DU3。两种方法趋势均为增大（图8b）,说明在此阶段细颗粒物质增多,水动力减弱,较低的东亚冬季风强度与中国黄土记录的东亚冬季风^[29]减弱相对应（图8a）。

在7.5~7 ka B.P.阶段,对应沉积单元DU2-2。在7.2 ka B.P.左右两种方法均达到峰值（图8b）,细粒物质组分含量达到最高,这与前人研究中的7.2 ka弱季风事件^[34⁃35]一致。7.2 ka B.P.之后,两种数值均呈下降的趋势,细粒组分减少,但幅度不大,此时仍处于季风较弱阶段。

在7~4 ka B.P.阶段,对应沉积单元DU2-1。从约6.5 ka B.P.至今,两种方法呈逐渐减弱的趋势（图8b）,细颗粒成分减少,水动力增强,东亚冬季风增强。这与前人研究中全新世早期东亚冬季风呈降低趋势,全新世晚期（6.5 ka B.P.之后）东亚冬季风呈增强趋势相对应^[29,31⁃32]。

从4 ka B.P.至今,对应沉积单元DU1。由图8可见,EM1/（EM1+EM2+EM3）与东亚季风强度拟合比较好,说明BasEMMA在此阶段对沉积环境有着更好的指示作用。两种方法的趋势均为减小,在约4 ka B.P.到2 ka B.P.相对稳定,从2 ka B.P.至今相对减弱。这说明沿岸水动力呈持续增强趋势,东亚冬季风也在相应增强。

选取2 ka B.P.以来的冬季温度变化数据^[36]（图12a）,用以比较BasEMMA与V-PCA在反映季风事件上的差异性。由图12b可以看出,GS数值在“小事件”时间段有更明显的变化。2 ka B.P.内GS显示3个对环境反映敏感的点。a点和c点均为细粒物质含量较少的点,两点的冬季温度变化都为负值,冬季温度降低,此时东亚季风较强（图12a）,a点与15世纪初明清小冰期所对应^[37⁃38],此时东亚冬夏季风都强,二者势均力敌,两种气团的交界停留在同一区域,会形成较长时间的降水和降温天气^[35]；c点与公元前29年西汉小冰期所对应^[39⁃40],之后是持续十几年的寒冷和干旱气候,此时冬季风强度高且持续时间长。b点冬季温度变化为正值,气温上升（图12a）,与13世纪初中国第四个温暖期所对应^[41],此时细颗粒物质含量较高,沿岸流较弱,细颗粒组分被留在原地,季风较弱。两种方法出现这种区别是因为从6.5 ka B.P.至今东亚季风逐渐增强,在此增强趋势中有些许年份是减弱,这些不同的输运机制对应不同的V-PCA粒度谱^[25],在寒暖过渡期更能体现东亚季风的变化。而BasEMMA和粒级—标准偏差法是从整体数据考虑,所以在一定程度上忽略了“小事件”。而且粒级—标准偏差法对粒度变化过于敏感,些许误差也会被计算其中^[42]。不过BasEMMA和粒级—标准偏差法较准确提取柱状样的敏感粒级,这是V-PCA所不及的。

Figure 12.  2 ka B.P. winter temperature changes to date, correlated with BasEMMA and V⁃PCA

因此,三种方法均能准确反映山东半岛海湾沉积物中的环境敏感组分,其中BasEMMA和粒级—标准偏差法更适用于整体敏感粒级提取,V-PCA对“小事件”的季风变化较敏感。

本文以威海湾WH-05岩心为例,研究了8.4 ka B.P.以来沉积环境演化以及不同方法提取的环境敏感因子的差异性及应用,取得如下认识。

（1） WH-05岩心可划分为3个沉积单元（DU1、DU2和DU3）。8.4 ka B.P.以来威海湾沉积物以粉砂为主,黏土次之,砂最少。粉砂、黏土和砂含量变化差异均较大,以跳跃组分为主,悬浮组分次之。DU1水动力较强,沉积环境不稳定；DU2上层水动力较强,随着水深逐渐增大水动力相对减弱；DU3沉积动力环境整体较强。

（2） BasEMMA、粒级—标准偏差法、V-PCA提取的环境敏感因子均能较好地反映海湾沉积环境演化。但BasEMMA和粒级—标准偏差法更适用于整体敏感粒级提取,V-PCA对特殊年份的粒度变化较敏感。

（3） 8.4~6.5 ka B.P.时期,GS序列和EM1/（EM1+EM2+EM3）都呈缓慢增大的趋势,此阶段东亚季风减小；从6.5 ka B.P.至今,GS序列和EM1/（EM1+EM2+EM3）呈降低趋势,指示东亚季风加强。