基于粒度分布曲线的邻近传播聚类算法在沉积环境识别中的应用<br/>——以白洋淀地区为例

基于粒度分布曲线的邻近传播聚类算法在沉积环境识别中的应用
——以白洋淀地区为例

2020-02-28刘祥奇宋磊吴奇龙李国民毛欣

海洋地质与第四纪地质 2020年1期

刘祥奇，宋磊，吴奇龙,3，李国民，毛欣

1. 中国地质科学院水文地质环境地质研究所，石家庄 050061

2. 中国地质大学（北京）工程技术学院，北京 100083

3. 中国地质大学（武汉）地球科学学院，武汉 430074

沉积物粒度组成及其分布模式是沉积物的基本特征，能够揭示沉积物形成的动力条件，是判别沉积环境和沉积物源的重要指标。不同沉积环境形成的沉积物受其搬运方式、动力条件等的约束，可分选出特定的粒级组合，因此粒度分布特征分析具有重要的沉积学意义[1]。沉积物粒度频率分布曲线可以定量反映某一粒径分布含量、次总体类型、分选程度等沉积过程的细节信息[2]。这些信息的提取有助于理解沉积物的来源、成因及形成过程，尤其是粒度分布子体分离法[3]、Weibull分布函数拟合法[4]，端元模型法[5]等统计方法的应用，对于研究不同粒度组成端元或者组分的气候环境意义具有重要的推动作用。上述统计学方法重点关注某一沉积类型沉积物形成的动力过程和沉积成因，对于沉积环境复杂、沉积类型多样的沉积区，其沉积物成因相对复杂，如何快速、准确判别沉积物的沉积相显得尤为重要。聚类分析作为数据挖掘领域重要的手段之一，直接比较数据之间的性质差异，将数据性质相同或相近的聚为一类，在沉积物粒度数据分析与处理中具有较好的适用性[6]。然而，常规的聚类分析方法需要预设聚类数量，不能真实地反映根据沉积物粒度分布曲线类型反演出研究区存在的沉积动力类型，而且以往研究应用聚类分析大多以粒度参数表征粒度频率分布曲线作为研究对象，不能精确地表示粒度分布曲线形态以及各粒径分布含量[7]。值得注意的是，AP聚类算法（Affinity Propagation Clustering Algorithm）作为新兴的聚类方法，具有快速、高效、不需要预设聚类数量的优点，对于高维数据聚类效果明显，在图像、文本聚类以及基因表达数据分析等领域应用广泛[8-10]，将AP聚类算法方法引入沉积学研究中，对于粒度数据的深度挖掘与分析具有重要的推动作用。

白洋淀地区分布深厚的第四纪沉积物，沉积环境较为复杂，以冲积、洪积和湖积作用为主，是开展不同沉积环境条件下沉积物粒度分布曲线判别与划分研究的理想区域。区内主要侧重于利用气候环境代用指标进行古环境演化方面的研究，缺乏沉积环境动力条件的研究[11-13]。通过采集研究区内沉积剖面中已知沉积相的沉积物，确定不同沉积环境的沉积物粒度特征，研究AP聚类算法在沉积学研究中的可行性与适用性，尝试建立沉积物粒度特征与沉积动力条件的内在联系。

1 自然地理概况

研究区位于华北平原白洋淀及其周边区域，是永定河冲积扇与滹沱河冲积扇的前缘低洼区，位置为 38.7°～39.2°N、115.6°～116.5°E。研究区属于温带大陆性半湿润半干旱气候，年平均温度为7.3～12.7 ℃，多年平均降水量564 mm，主要集中在6—8月份[14]。区内地势平坦，由西向东地势略有降低。白洋淀有拒马河、府河、唐河等九条河流汇入以及赵王新河、大清河等出淀河流。第四纪以来，该地区不断沉降，并普遍接受沉积，而且区内水文条件的变化对第四纪沉积物的形成也具有重要的影响，逐渐形成了一套以冲洪积和湖积为主的沉积地层[15]。

2 材料与方法

2.1 样品采集与测试

本研究对白洋淀周边地区露头剖面的沉积地层进行了详细考察，选取其中22个剖面作为研究对象（图1，剖面地层描述见表1），对各剖面典型沉积相层位进行样品采集，共获得85个沉积样品。现代河流漫滩相沉积样品来源于白沟引河，共采集2个沉积样品，沉积物为棕黄色黏土质粉砂，具有典型的二元结构。现代湖泊表层沉积样品来源于白洋淀湖泊，共采集12个沉积样品，沉积物为灰黑色粉砂质黏土、黑色黏土，有机质含量高（图2）。

粒度测试在中国地质科学院第四纪年代学与水文地质演变重点实验室完成，采用mastersizer 2000激光粒度仪进行测试，该仪器的测试范围0.2～2 000 μm。测试前加入足量的10%双氧水并加热，去除有机质。待烧杯中不再冒泡之后，加入足量的10%盐酸，直至不再冒气泡才停止加热，然后加入蒸馏水静止24 h。利用虹吸法抽取上层清液，加入蒸馏水，静止24 h，反复几次直至呈中性。加入10 mL 10%的六偏磷酸钠，用超声波振荡10 min后上机测试。

图1 研究区域采样点Fig.1 study area and sampling sites

2.2 邻近传播算法（AP聚类算法）

邻近传播聚类算法又称AP聚类算法，是B J Frey 和 D Dueck 于 2007 年首次提出的一种新的聚类算法[16]。它是一种不需要预设聚类数，通过数据点之间的“信息传递”来寻找自己的类代表点（exampler），将属于同一类代表性点的样品聚成一类，从而实现聚类过程[17]。

AP聚类算法将所有的样品点都看作为潜在的类代表性点，以数据之间的相似度系数基础，两两之间进行“信息传递”。在信息传递的过程中，利用因子图（图3）中The Max-Sμm algorithm准则进行信息的迭代更新，迭代收敛时，使得公式（1）在约束条件下网络相似度S（c）达到最大。

式中，c为由类代表性点构成的向量，i为样品点，s（i，ci）表示样品点与其对应的类样品点之间的相似度；δk（c）是约束条件，表示数据点 i选取 k作为它的类代表性点，但是数据点k并不是类代表性点，因此使其值为-∞，用于限制没有类代表性点的类簇产生。

相似性数据矩阵s（n，n）是由n个数据之间的相似度系数构成，对角线上的元素由人为设定的偏向参数所构成，它是用于评价各数据点作为类代表性点可能性的程度。数据之间传递的信息包括：从潜在的类代表性点传递到数据点的归属度和从数据点传递到类代表性点的吸引度[16]。在因子图中，基于The Max-Sμm algorithm准则信息更新的公式：

当信息迭代到更新收敛时，点i所对应的类代表点ci：

基于AP聚类算法的基本思想及其原理，可得算法的基本步骤：

第1步：根据样品点之间的相关距离，计算相似的矩阵S以及设置各个样品点对应的偏向参数。

第2步：设置各个样品点的初始归属度a (i,k)=0。

第 3 步：利用公式（2）、（3）进行信息更新，在更新的过程中，为了防止迭代震荡，引入阻尼因子λ，范围（0，1）。

表1 采样点描述Table 1 Description of sampling sites

第4步：当各个点对应的类代表性点ci在迭代过程中不再发生改变或者迭代次数达到设定的最大值时，迭代终止。将相同类代表性点的样品点归为一个类簇。

3 沉积物粒度分布曲线的AP聚类算法

沉积物粒度分布曲线特征是判断沉积作用类型的重要手段之一。通常，研究者采用粒度参数代替粒度频率分布曲线，近似地作为沉积作用类型的代用指标，常常会丢失曲线中蕴含的重要信息并且具有多解性。为了充分描述曲线的分布特征，在-1～11.3Φ（0.4～2 000 μm）范围内按间距 0.16Φ 划分 75个粒级。同时，考虑到某些粒级变量的含量较低并且在不同样品之间存在显著差异，使用对数比转换的数据，可能会掩盖住数据的结构信息，本文选用经过标准化的原始数据[18-20]。采用欧式距离作为沉积物粒度频率分布曲线之间的相似性的度量，则样品x、y之间的相似度：

图2 大步村剖面及其地层描述Fig.2 Dabu Village section and its stratigraphic description

图3 邻近传播因子图Fig.3 Factor GrAPh for affinity propagation

式中：xik，xjk分别表示沉积物样品i，j的第项变量含量。

经过上述数据的前处理程序，AP聚类算法将85个剖面典型沉积相沉积物的粒度频率曲线划分为 11 簇（图 4）。

11簇聚类结果显示，每一簇的粒度频率分布曲线分布范围比较集中，各粒级所对应的含量之间的差异很小，粒度频率分布曲线相似度较高。根据AP聚类算法的聚类结果采用silhouette函数绘制每一个点的轮廓值。轮廓值作为聚类分析有效性的常用内部指标，它被定义为：

式中，a（i）是第i个点与同簇的其他点之间的平均距离；b为一个向量，其元素是第i个点与属于不同类的各个点之间的平均距离。轮廓值的取值范围为[-1,1]，此值越大说明该点的分类越合理[21]。本文中11簇样品点的轮廓图如图5所示。

图4 划分的11类簇曲线Fig.4 11-group cluster curves

图5 各簇样品点的轮廓值Fig.5 Sihouette value of the sample points of each cluster

其中，各簇内样品点的轮廓值大部分大于0.3，仅有5个点的值小于零，最小为-0.08。考虑到划分聚类的数量以及各簇之间数据差异性相对较小，采用AP聚类算法对于沉积物粒度频率分布曲线的聚类划分是合理的，并且聚类效果显著。值得一提的是，AP聚类算法可以通过调整偏向参数能够增加或者减小聚类数量进而改变各簇内数据之间的差异大小，使得根据研究目的获得适合的聚类效果。本文选用数据点相似系数最大值作为偏向系数，以此保证各簇内数据差异很小，进而将动力条件差异不明显的沉积环境形成的粒度样品聚为一类。

4 分析与讨论

4.1 基于AP聚类算法的沉积物沉积相识别

不同沉积类型的粒度分布曲线各粒级对应的百分含量不同，沉积类型对应的动力条件、运输介质等差异越大则粒度分布曲线分布区间、形态差异越明显[22]。当沉积物所受的搬运介质、搬运方式一定，并且介质动力大小稳定时，沉积物粒度组分的总体服从于某种分布的单因子随机事件的集合，在频率分布曲线和概率累计曲线形态上表现为单峰和一段式，否则呈多峰和多段式[23-24]。如前文所述，AP聚类算法将研究区沉积物的粒度频率分布曲线划分为11类簇，各簇内曲线形态、分布范围相近。但是，考虑到沉积物颗粒对于动力条件变化十分灵敏，同一沉积环境条件下粒度频率曲线会随着流体动力条件的波动在曲线形态、分布范围上产生相应的变化。因此，需要对曲线形态相似、分布范围发生轻微变化的曲线进一步归纳分析。

从11类簇粒度分布曲线可知，每一簇曲线与该簇的类代表性点之间的差异很小，利用类代表性曲线作为每簇曲线的典型样品曲线，并按曲线峰态、分布范围以及搬运方式将11簇曲线进一步划分为4组（图6）。第一组包含Ⅳ、Ⅵ、Ⅸ三簇对应的沉积粒度频率分布曲线，曲线形式表现为单峰、一段式（图 4、6a），峰值分布于 5～20 μm，分布范围较为集中，仅存在悬浮搬运方式，表明沉积物形成于动力条件较弱的静水环境；第二组包含Ⅱ、Ⅴ、Ⅷ、Ⅹ、Ⅺ五簇沉积物粒度频率分布曲线，曲线形态为双峰、两段式（图 4、6b），主峰在 20～200 μm，次峰在5～20 μm，跳跃和悬浮总体的搬运截点在4Φ左右，各总体的含量有所差异，表明沉积环境较第一组动力强度有所增加，并且水动力相对变化。第三组包括Ⅰ、Ⅶ类簇粒度分布曲线，为双峰、两段式（图4、6c），主峰在 100～500 μm，次峰在 10～100 μm，主要为跳跃总体存在，搬运截点变化较大，表明沉积物形成于强动力的沉积环境；第四组为Ⅳ类簇，为双峰、两段式（图4、6d），跳跃与悬浮总体含量相当，表明沉积物形成于动力条件强并且变化明显的沉积环境。

根据上述分析，4类粒度频率分布曲线的组合可能代表了4种较为一致的沉积环境，将各组中代表性样品粒度频率分布曲线与稳定湖相（湖心相）、湖滨相、河流相等典型沉积相粒度频率分布曲线进行对比（图7）[23]。图中可以看出聚类结果得到的曲线与已知沉积环境的粒度曲线主次峰位置、曲线形态、分布范围较为一致。结果表明，第一组沉积物形成于稳定湖相，包括动力条件相似的湖心相、湖沼相；第二组沉积物形成于湖滨相或漫滩相；第三组沉积物形成于河流相；第四组为洪水沉积物。由此可见，AP聚类算法确定的4个组合分别代表了4种沉积环境。为了进一步验证上述分析的准确性与可行性，将上述结果与野外判别的4类85个沉积样品的沉积相进行对比发现，AP聚类算法确定的4个组合可以将野外识别的沉积相类别聚集到一起，两者具有较好的一致性。

前人对大量已知沉积相的沉积物粒度参数进行定量描述，建立了多种沉积环境之间沉积物粒度参数判别公式。利用张平等根据典型稳定湖相与典型风成黄土沉积物粒度参数建立了稳定湖相与风成黄土粒度判别式（式8）[25]。本文利用此判别式对采集的所有样品进行沉积环境的判别分析。

图6 四组代表性样品频率曲线、概率曲线Fig.6 Frequency distribution and probability cumulative curves of four groups of representative samples

图7 各簇代表性与典型沉积相粒度频率分布曲线对比Fig.7 Comparison of frequency distribution curves of typical sedimentary facies and representative sediments in each cluster

式中，F表示判别值，若为稳定湖相沉积物；Mz、Sd、SK以及Ku分别是利用矩值法计算以为单位的平均值、分选系数、偏态以及峰态。

其中第Ⅳ、Ⅵ、Ⅸ簇为正值，代表稳定湖相，判别结果与对比法得出的结果一致。

4.2 研究区不同沉积环境的沉积物粒度特征

基于粒度频率分布曲线的邻近传播聚类算法将全部85个粒度样品聚集成4类，并与已知沉积相粒度频率曲线进行对比，确定4类曲线对应的沉积环境分别为稳定湖相、湖滨相或漫滩相、河流相以及洪积相四种沉积相。各类沉积物之间在粒度频率分布曲线形态、分布范围以及组成成分上存在差异，通过定量描述这些差异能够为沉积环境、类型的区分提供依据[26]。粒度频率分布曲线利用粒度参数进行定量描述，粒度参数主要有粒度平均值、分选系数、偏态和峰态。为了更加充分揭示粒度参数对于沉积环境的指示意义，选用物理意义明确的Folk-Ward图解法进行计算（表2）[27]。表3给出了聚类后类代表性样品的粒度组成成分、粒度参数统计。

从表3中可以看出聚类结果的11簇样品中，第Ⅱ、Ⅴ、Ⅷ、Ⅹ、Ⅺ簇样品之间粒度参数和组成成分相近，第Ⅳ、Ⅵ、Ⅸ簇样品之间粒度参数和组成成分相近，第Ⅰ、Ⅶ簇样品与第Ⅲ簇样品粒度分布形成于高能环境。4类典型沉积相的粒度组成与粒度参数之间存在显著差异。

第一类包括Ⅳ、Ⅵ、Ⅸ簇沉积物，粒度组分以细粉砂为主，黏土、细粉砂和粗粉砂三者含量之和占95%以上，平均粒径6.227～7.537Φ，在4类沉积环境中粒度成分最细，分选系数为1.301～1.722，分选差，除Ⅸ簇中部分样品呈正偏外，其余全为近对称分布，峰态变化于1.066～1.158，第Ⅳ 簇呈常峰态，其余全为尖峰态。

第二类包括Ⅱ、Ⅴ、Ⅷ、Ⅹ、Ⅺ簇沉积物，其中Ⅱ、Ⅷ、Ⅹ、Ⅺ簇沉积物以粗粉砂为主，占40%以上，黏土、细粉砂、细砂含量均大于10%，四者之和占98%以上，平均粒径为5.255～5.841Φ ，第Ⅴ簇沉积物以细砂（52%）为主，粗粉砂（34%）次之，平均粒径为4.213Φ，较Ⅱ、Ⅷ、Ⅹ、Ⅺ簇沉积物颗粒稍粗，表明强度略大于Ⅱ、Ⅷ、Ⅹ、Ⅺ簇沉积物对应的动力环境。该组分选系数为1.343～1.799，分选差，偏态为0.209～0.467，属于正偏、极正偏，峰态为0.887～1.561，第Ⅱ簇沉积物粒度分布曲线呈宽峰态，第Ⅴ、Ⅷ簇粒度分布曲线呈尖峰态或极尖峰态，其余两簇为常峰态。

表2 Folk&Ward 图解法公式Table 2 Formulas of Folk-Ward graphic methods

表3 各簇代表性样品粒度组成成分及粒度参数Table 3 The size composition and size parameters of representative samples in each cluster

第三类包括Ⅰ、Ⅶ簇沉积物，该组粒度组分较粗并且含量差异相对较大。其中，第Ⅰ簇沉积物以中砂（45%）为主，细砂（18%）、粗砂（22%）次之，平均粒径为1.462，样品分选差，呈极正偏，为尖峰态，代表强动力沉积环境；第Ⅶ簇沉积物中细砂含量占70%以上，平均粒径为2.919Φ，颗粒较第Ⅰ簇沉积物偏细，表明介质动力条件强并且稳定，相比于Ⅰ簇沉积物对应的沉积环境动力条件较低。

第四类中第Ⅲ簇沉积物以细砂为主，粗粉砂、中砂、粗砂含量均大于15%，平均粒径为1.582Φ，分选系数1.976，样品分选差，极正偏，宽峰态，表明沉积动力复杂多变且动力较强。

4.3 AP聚类算法在沉积环境识别中应用的可行性

根据上述方法对粒度分布曲线直接进行聚类分析，避免人为干扰聚类效果，有效地揭示研究区内的沉积环境差异是此方法的优点，也是能够利用此方法对粒度分布曲线进行研究区所存在的沉积环境类型反演的出发点。本文以白洋淀地区作为研究对象，尝试利用AP聚类算法建立粒度分布曲线与沉积环境之间的内在联系，通过识别结果与野外识别进行对比，发现所得结果较为一致，85个样品中有6个不一致。如下给出了沉积环境识别过程中存在差异的沉积物样品信息（表4）以及样品对应的粒度分布曲线（图8）。

表4 沉积环境识别过程中存在差异的沉积物样品Table 4 Sediment samples with differences in the process of sedimentary environment identification

图8 沉积环境识别过程中存在差异的样品粒度曲线Fig.8 Sediment samples with differences in sedimentary environment identification

其中，SITE01-03、SITE14-01、SITE14-02三个样品的粒度曲线具有明显双峰，分布范围和曲线形态与典型湖滨相或漫滩相较为符合，野外观察结果存在识别错误（图 8）。SITE03-04、SITE05-01、SITE13-02样品峰值粒径比典型湖沼相粒径较粗，但是没有明显次峰（图8）。对比其他簇曲线形态及分布范围，此类样品更为可能形成于沉积环境相对稳定，但动力强度较强的湖沼相沉积环境。AP聚类算法进行聚类分析时，由于样品相对典型湖沼相样品粒径较粗，被划分到湖滨相沉积物一类，出现了偏差。因此，对于85个样品中有3个样品沉积环境识别出现误差，错误率为3.53%。

鉴于聚类方法是尽可能将相似的数据点聚集成一类，沉积物的形成受物源、介质等多个因素影响，同一沉积环境条件下，粒度分布曲线会在一定范围内发生波动，因此采用邻近传播聚类分析时需要合理选取偏向系数控制各簇之间数据的差异。本文为了充分挖掘粒度数据中蕴含的沉积环境类型，选取了相似系数最大值作为偏向系数，尽可能减小数据之间差异，导致出现同一沉积环境条件下形成的粒度曲线可能被划分到不同簇中的现象。