长江中下游地区历来都是我国洪涝灾害最严重的地区，对其洪水规律的研究意义重大。然而历史或现代仪器所记录的洪水时间尺度较短，不足以认识和把握洪水发生的规律，因此有必要利用地质资料对古洪水进行研究^[1-5]。湖泊沉积记录因其沉积连续、记录介质丰富、分辨率高等特点，近年来被广泛的应用于古洪水的研究中^[6-8]。长江中游地区多数湖泊主要分布在长江天然堤后，当湖泊形成后便与长江没有直接的水道联系，只有在洪水期长江水位上涨时才会漫过天然堤，将外源物质带入湖泊中沉积。较之正常沉积，洪水期的沉积物在粒度、磁化率指标上存在差异，这为长江中游地区的古洪水研究提供了良好的载体。

沉积物的粒度特征可以用来指示碎屑颗粒物质沉积时的水动力强弱以及颗粒物质的来源^[9-12]。相比传统的粒度频率分布曲线，粒度端元分析法^[13-17]的优势在于：1)更有效的利用高精度粒度数据；2)从多组分混合物中分离出每一单组分含量和众数值，从而更深入的对粒度组分对应的沉积动力过程进行研究^[18]。端元模型分析法在各个领域均有着广泛的应用^[19-22]。传统的古洪水识别方法往往直接对各指标曲线划分古洪水期，容易受到数据中的误差以及人的主观因素影响。对此，笔者借助端元分析法，从整体数据的协方差矩阵计算端元的数量及分布，结合湖泊演化及区域地质条件，获得了指示古洪水的端元模型，并且对古洪水期进行了有效识别。

武汉市(29°58′~31°22′ N、113°41′~115°05′ E)位于我国长江中游地区，地处江汉平原东部，境内河道纵横，湖泊密布(图 1)。主要的河流有长江和汉江，其中长江自武汉西南的小军山流入境内，主要流向为北东方向，经金口斜穿市区，北纳东荆河、汉水、府河和倒水，南接巡司河、青山港等水域，流过谌家矶后转向东南，于营盘山出境，境内流长约64 km；汉江自新沟进入市区，自西向东在汉阳龟山的集家咀汇入长江，境内流长约36 km。武汉境内共有湖泊166个，形态各异，大多以天然堤后湖泊为主，湖岸线相对曲折，常见湖滨沼泽湿地，区内湖泊多经历复杂的历史变迁。本研究以较为典型的ZK145钻孔为载体，该钻孔位于长江南岸青山船厂(30°39′41″ N、114°26′56″ E)，钻孔位置原来为湖泊，后期湖泊范围不断萎缩，近代因人工填湖而成为陆地。

图  1  ZK145钻孔位置及研究区地貌略图

Figure 1.  Borehole location of ZK145 and overview of research area

ZK145钻孔第四纪沉积物厚度约40 m左右，平均取芯率为79.88%，除地表 1.3 m为人工填土外，其沉积环境大致可分为：下部为粉砂、细砂和砾石构成的典型二元结构，属河流相沉积环境；上部为粉砂、黏土，属湖相沉积环境(图 2)。具体分层情况如下：1.3~3.0 m(第1层)整体为褐黄色黏土质粉砂，稍湿，可塑；3.0~17.0 m(第2层)整体为粉砂质黏土层，自上而下颜色依次为褐黄色—灰绿色—深灰色，含两层浅灰色细砂(深度分别为13.2~13.8 m和15.0~15.5 m)；17.0~28.5 m(第3层)整体为深灰色粉砂，含植物残骸，底部有一层黄褐色细砂；28.5~31 m(第4层)整体为深灰色黏土层；31.0~34.0 m(第5层)整体为黄褐色细砂—中砂层，底部为中砂，含水率较高；34.0~37.9 m(第6层)整体为砾石层，34~35.3 m含褐黄色粗砂，砾石成分以石英、燧石为主，粒径在2~

图  2  ZK145钻孔剖面及年代学框架

Figure 2.  Lithology profile and chronological framework in ZK145 borehole

5 cm之间，最大粒径为15 cm，分选较差，磨圆以次棱—次圆为主，37.9 m以下为基岩。

该钻孔综合样品(粒度、磁化率)为连续采样，样品间隔为10 cm，共采集257组。粒度样品前处理以及所有测试工作均在南京师范大学地理科学学院粒度分析实验室进行。测试前称取适量自然风干的松散样品放入烧杯，依次加入10%的稀双氧水溶液和10%的稀盐酸溶液，以去除有机质和碳酸盐。然后用蒸馏水清洗样品并静置24 h后，加入10%的六偏磷酸钠溶液，并在超声波清洗仪中震荡以使颗粒充分分散。粒度样品的测量采用仪器为Malvern公司生产的Mastersize 2000型激光粒度仪，检测范围为0.02~2 000 μm，重复测量误差小于2%。磁化率样品的测试工作由中国地质大学(武汉)地球物理与空间信息学院岩石磁学实验室完成，采用的仪器为捷克生产的KLY-3S型卡帕桥仪器，灵敏度高达2×10^-8(SI)。

此外，ZK145孔共采集5组AMS¹⁴C测年样品，分别位于钻孔深度为3.2 m、14.0 m、19.6 m、24.9 m和30.1 m处，所采样品均为黑色黏土。AMS¹⁴C测年样品在位于波兰的Poznań Radiocarbon Laboratory进行测试，并对测试结果运用OxCal 4.3软件进行树轮校正，从而得到样品的测年结果。根据测年结果，采用分段线性内插和外推的方法，获得了钻孔各深度的年龄值(图 2)。按此方法推测该钻孔岩芯40 m处的年龄约12.45 cal. ka B.P.。地质时代属于晚更新世末期。

沉积物携带、传输和沉积是一个分选的过程，因此与特定沉积过程有关的沉积物通常具有特定的粒度分布，即粒度端元^[15]。端元的提取有参数化与非参数两种方式。对于参数化方法，常通过指定函数类型、数量及位置参数来对数据进行分解，其中，Weibull函数是广泛应用的拟合函数之一^[23]。但是参数化方法存在一些问题^[14]：1)只能使用单个粒度样本导致无法充分提取整体沉积信息，而且指定的函数与沉积地质过程并无明显关联；2)确定端元数量存在一定难度。对于非参数化方法，Weltje^[13]首次提出了端元建模分析，基于主成分分析、因子旋转、非负最小二乘法等提取端元，并由Dietze et al. ^[16]为其实现了稳定的算法。非参数化方法无需规定端元主峰作为先验条件，而是通过整体数据协方差反映端元数量及组成，反映的是整体样品的粒度组成特征，更适用于钻孔沉积物样品的粒度分析。该算法主要包含两种途径：1)选取不同的端元数量q进行迭代，根据其方差贡献率大小及与原始数据的相关性来综合确定最佳端元数量；2)将确定的端元数目值带入计算得到端元模型，并对端元模型进行评估。

算法流程如下：

(1) 为了避免每列数据因尺度差异而带来的弱相关性，需要对粒度矩阵进行缩放处理。采用的方式为选取分位数l，使得每一列数据下界为P_l，上界为P_100-l：

式中：W_m×n为进行缩放处理后得到矩阵；X_m×n为粒度矩阵；X_Pl为分位数取l时，该列上的取值大小，X_P100-l为分位数取100-l时，该列上取值大小，更一般的，l=0则反映了每一列的最小值与最大值(对l取不同值的原因在于，如果在某些粒级上存在大量数据的值为0，则需要通过取不同的分位数来确保模型的稳定性)。

(2) 对W_m×n做内积处理，然后再做特征值分解。

式中：A_n×n为W_m×n与自身求内积后的方阵；λ_n×n为特征值，除了对角线元素，其余元素均为0。v_n×n为特征向量。

(3) 通过特征值方差的累计贡献观察潜在的端元数量。类似于因子分析，选取不同的端元数量q对特征空间采用方差最大化因子旋转法(Varimax)进行正交旋转(在Matlab中调用rotatefactors函数)，于是计算得到因子载荷矩阵V_n×q。

(4) 将因子载荷矩阵V_n×q进行放缩，使其值范围在(0，1)之间，于是得到B_n×q。再由：

通过线性非负最小二乘法(在Matlab中调用lsqnonneg函数)，经过计算可得到特征向量的得分S_m×q。

(5) 再对最后计算得到的因子载荷矩阵和因子得分矩阵进行放缩处理，目的是使其与原数据单位一致，于是得到端元载荷矩阵B_n×q^*和端元得分矩阵S_m×q^*。

对模型误差进行评估。

式中：X_m×n^*为模型拟合得到的数据；E_m×n为模型拟合数据与原始数据的误差矩阵；r_n为每列求出的相关性，对r_n求平方和得到变量间相关性r_变量²；r_m^*为每行求出的相关性，对r_m^*求平方和得到样本间相关性r_样本²；r²表示总的相关性，选取不同的q值进行计算，绘制q与r²关系图，可以用于选取合适的端元数量。

ZK145钻孔混合样品的频率分布曲线如图 3所示。频率分布平均值曲线的局部极大值分别在3.5 ϕ和6.5 ϕ附近。该钻孔的频率分布曲线大体上有4种分布特征：1)峰值在3~4 ϕ和6~8 ϕ的双峰分布(如：样品ZK145-004、ZK145-008和ZK145-010)；2)峰值在7 ϕ附近的单峰分布(如：样品ZK145-030、ZK145-090和ZK145-120)；3)峰值在1~1.5 ϕ和4~6 ϕ的双峰分布(如：样品ZK-138、ZK-152和ZK-213)；4)峰值在1~2 ϕ和4~6 ϕ的双峰分布(如：样品ZK145-241、ZK145-244和ZK145-250)。

图  3  ZK145钻孔粒度频率分布曲线

Figure 3.  Grain size distribution curves of representative samples from the ZK145 borehole

为了从混合组分中分离出与特定的沉积过程有关的粒度分布，从而提取出指示古洪水的粒度端元指标，首先通过Mastersizer 2000软件从粒度数据中获得了大小为257 × 47的粒度矩阵(范围由0到11.5 ϕ，间隔为0.25 ϕ)，以运用特征向量旋转算法对其进行端元建模分析。

对ZK145钻孔粒度矩阵由(1)~(3)式计算得到不同分位数l下的特征值与特征向量，由此得到不同端元数量下的累计方差贡献(ECV)(图 4)。

图  4  ZK145钻孔粒度累计方差贡献率(ECV)

Figure 4.  Explained cumulative variance (ECV) of grain size data for the ZK145 borehole

可以发现，当端元数量达到4时，ECV的值基本超过95%，即可以认为端元数量q ≥ 4。同时，为确保端元数目的准确性，还需对模型拟合得到的数据与原始数据之间进行相关性分析：选取不同的q值，由(4)~(11)式计算X_m×n与X_m×n^*的平均相关系数r²，得到端元数量q与r²的关系图(图 5)。

图  5  ZK145钻孔粒度数据平均相关系数

Figure 5.  Mean coefficient of determination r² of grain size data for the ZK145 borehole

在q-r²图上可以得到端元数量q为4时其相关性r² =0.823(针对该钻孔数据，l取值为0~2较为合适)，综合以上分析，确定端元数为4。当端元数目q =4时(l =0)，通过①~⑤步骤计算可以得到EM1~EM4端元在不同粒径上的载荷矩阵B_n×q^*，由此可做出端元模型的分布曲线图(图 6a)。为了进一步对模型拟合好坏进行评估，通过(9)、(10)式计算模型与原始数据之间的相关性，得到变量相关系数图(衡量每列变量的平均相关性)及样本相关系数图(衡量每行样本的平均相关性)(图 6b，c)。

图  6  ZK145钻孔粒度端元模型及其评价

从模型的拟合程度来看，与原始数据X_m×n相比较，两者的总相关性r²达到了82.71%(图 5)；在不同粒级上的平均相关系数r_变量²达到了82.28%，其中在3.5 ϕ、7 ϕ的附近粒级范围为最高(图 6b)；样本间的平均相关性r_样本²达到了83.13%，大部分样品的相关系数均在75%以上(图 6c)。从方差的解释上来看，4个端元分别解释了总方差的4.5%、36.1%、19.3%、40.1%，其中端元2和端元4占到了总解释方差的76.2%。

经分析计算，得到了ZK145钻孔的端元模型EM1~EM4(图 6a)。总体来看，分解出的4个端元均显示非对称、多峰分布的特点，显示多种营力的综合或同一营力不同强度沉积的混合作用。各粒度端元的特征以及与河流水动力之间的联系如下：

EM1的主峰峰值在1.5 ϕ(353.6 μm，属中砂范围)附近，在6 ϕ(15.6 μm，属粉砂范围)处有一个次峰。河流沉积物一般具有多峰的特征，而最典型的特征是粗粒组分的峰值在200~400 μm，甚至更粗的范围^[24]。粒径介于250~600 μm的中—粗砂端元常见于河流、沉积物分选较好的湖泊沉积以及冲积扇沉积，这种端元可视为河流冲积砂^[16]。EM1在钻孔深度28.3~22.9 m、33.8~31.5 m的含量为最高，平均百分含量分别为：10.3%、15.9%，而EM1含量在其他深度均接近于0，通过EM1可用来指示河流冲积过程。

EM2的主峰峰值在3 ϕ(108.8 μm，属细砂范围)附近，在8 ϕ(3.9 μm，属粉砂范围)附近有一小而平缓的次峰。河流沉积物跃移组分的众数粒级一般为200~400 μm，并且分选较好^[25]。EM2在形态上与EM1较为相似，显示其沉积动力过程与之相似。EM2主要是由外源水动力输送进入湖泊沉积而成，代表经流水搬运的细砂。EM2在钻孔深度34.0~22.9 m的含量为最高，平均百分含量为44.6%，在其他深度只在局部出现高值，通过EM2可用来指示外部强水流入湖过程。

EM3的主峰峰值在5 ϕ(31.3 μm，属粉砂范围)附近，在2~0 ϕ(250~1 000 μm，属中—粗砂范围)有一平缓的次峰。封闭湖泊碎屑沉积物分布主体的组分粒度一般在8~35 μm之间，都为粗偏态^[25]。EM3主要通过降水产生的地表径流搬运沉积，众数粒级的大小取决于降水量控制的地表径流的大小。EM3在钻孔中的分布较为均匀，平均百分含量为35.8%，通过EM3可以用来指示封闭湖泊沉积过程。

EM4的主峰峰值在7 ϕ左右(7.8 μm，属细粉砂范围)，在3.5~2 ϕ(88.4~250 μm，属细砂范围)处有一平缓的次峰。EM4与EM3形态相类似，同为地表径流搬运而成，但EM4的众数粒级较EM3小，显示地表径流的水动力条件较之更弱。EM4在钻孔深度22.9~3 m的含量为最高，平均百分含量为50.3%，同样的，可以通过EM4指示封闭湖泊沉积过程。

ZK145孔所在的青山船厂湖，属于长江武汉段堤后湖泊，可以作为长江中游古洪水记录的良好信息载体，主要有以下原因：1)长江武汉段位于江汉平原东部边缘，是长江中游江汉—洞庭平原洪水出口河段，对古洪水有较好的反映；2)不同于平原腹地河道长期振荡，长江武汉段受到两岸基岩丘陵(龟山、蛇山、小军山、白浒山等)、岗地夹持，其河道摆动较小^[26]，天然堤稳定且不断加积发展，堤后湖泊得以长期存在；3)船厂湖位于长江天然堤与中晚更新世网纹红土及下蜀土组成的岗地之间，无较大水系注入，其物源主要为湖沼内原生沉积、周缘岗地侵蚀(面流或沟谷汇入)沉积和长江洪水沉积。当长江发生较大洪水时，洪水将会携带泥沙越过天然堤，从而使湖泊沉积的粒度组成及其物质组分发生变化。届时，洪水带来的外源物质会使得沉积物磁化率数值出现较大波动。本研究结合所分离的端元组分(EM1、EM2、EM3及EM4)及磁化率指标对钻孔沉积阶段进行了划分，并且在湖相沉积阶段(22.9~3 m)成功识别出了9次古洪水事件。

第Ⅰ阶段：37.9~34.0 m，相应的年代约为12.09~11.41 cal. ka B.P.。该阶段岩性主要以砾石为主，为河床相沉积。

第Ⅱ阶段：34.0~22.9 m，相应的年代约为11.41~9.69 cal. ka B.P.。在该阶段粗粒端元组分(EM1、EM2)总体为最高，曲线波动剧烈，并且呈现出由高降低再升高的趋势。EM1在该阶段含量为最高。体积磁化率曲线在该阶段由波动剧烈逐渐变为平缓。该阶段岩性主要为中粗砂、细砂、粉砂，中部夹少量黏土层，可判断其为河湖交替相。此时，粗粒组分端元含量的增加不仅与长江古洪水有关，还与周缘岗地侵蚀(面流、沟谷汇入等)有关。

第Ⅲ阶段：22.9~3 m，相应的年代约为9.69~3.43 cal. ka B.P.。在该阶段EM1含量趋于0，EM2峰值呈现交替出现的态势，磁化率曲线波动幅度为最大，并且与EM2组分的峰值有较好的对应关系，主要原因为古洪水携带泥沙越过天然堤在湖泊中沉积。该阶段岩性主要以粉砂、黏土为主，属湖相沉积。此时，受周缘岗地侵蚀影响较小，粗粒组分端元(主要是EM2)含量增加的主要原因为古洪水携带泥沙越过天然堤沉积而造成。

第Ⅳ阶段：3.0~1.3 m，相应的年代约为3.43~1.48 cal. ka B.P.。在该阶段EM2含量曲线波动剧烈，其值呈现由低升高再降低的趋势，磁化率曲线也呈现由低升高的趋势。该阶段岩性以粉砂质黏土为主，属湖沼相沉积。

根据前人研究，体积磁化率值大小与沉积时的水动力强弱因素相关^[27]，古洪水一般携带上游流域更大范围泥沙进入中下游沉积区域，由于物源改变，使得体积磁化率值发生波动，因此可将同一时期体积磁化率的高值或低值作为古洪水判别的依据之一^[28]。第Ⅲ阶段(22.9~3 m)较好的记载了古洪水事件，为了对其进行识别，本文主要以EM2含量的高值作为识别标志，并以同一时期体积磁化率的高值或低值作为辅助判别，共识别出了9期古洪水沉积层(图 7)，深度分别为：21.5~20.9 m、20.0~19.4 m、18.6~18.1 m、17.5~17.0 m、15.8~15.0 m、14.0~13.2 m、8.9~8.7 m、7.6~7.3 m、4.1~3.8 m，通过AMS¹⁴C测年数据线性插值，确定了9期古洪水事件的年代分别为：9 692~9 667 cal.a B.P.、9 629~9 402 cal.a B.P.、8 767~8 343 cal.a B.P.、7 801~7 411 cal.a B.P.、6 361~5 649 cal.a B.P.、4 844~4 754 cal.a B.P.、4 270~4 245 cal.a B.P.、4 126~4 083 cal.a B.P.、3 732~3 696 cal.a B.P.。Zhu et al.^[29]对长江中游地区古洪水事件进行了比较，指出在9 000~8 400 a B.P.、7 500~7 200 a B.P.、4 200~4 000 a B.P.等时段均发生古洪水事件；朱诚等^[30]对长江三峡及江汉平原地区全新世洪涝灾害进行研究，指出8 000~5 500 a B.P.、4 700~3 500 a B.P.为两个古洪水频发期；葛正帅等^[31]对长江上游三峡河段古洪水记录研究发现，7 810±445 a B.P.、4 840±465 a B.P.、3 983±483 a B.P.均为特大洪水时期；张玉芬等^[28]基于磁组构特征对江陵剖面进行研究，在9 874~4 600 a B.P.左右共识别了17期古洪水事件。本文研究结果与长江流域全新世洪水研究的结果相一致。

图  7  粗粒组分端元含量与磁化率指标对古洪水的指示

Figure 7.  Indication of the endmember content for the coarse-grained components and the magnetic susceptibility index to the paleo-flood

(1) 对具有代表性的长江堤后湖泊ZK145钻孔晚第四纪沉积物进行端元建模分析，可以分离出4种端元模型：EM1代表河流冲积砂，EM2代表经流水作用的细砂，两者共同组成了ZK145钻孔的粗粒组分；EM3、EM4代表静水湖泊沉积，两者共同组成了ZK145钻孔的细粒组分。同时也对模型的拟合进行了评估：对原始数据的累计方差贡献超过了95%，与原始数据的总相关性达到82.71%，均表明模型能较好对数据进行拟合。此外，经与样品的频率分布曲线比较，端元建模所分离的4种端元模型与混合样品的粒度频率分布曲线具有较好的一致性。

(2) 利用端元含量结合磁化率指标将ZK145钻孔主要划分为四个阶段：河床相沉积阶段、河湖交替相沉积阶段、湖相沉积阶段、湖沼相沉积阶段。经过比较发现粗粒组分端元含量的波峰对应着磁化率呈现剧烈波动的态势，利用粗粒组分的端元模型(主要为EM2)有效识别出9期特大古洪水，分别为：9 692~9 667 cal.a B.P.、9 629~9 402 cal.a B.P.、8 767~8 343 cal.a B.P.、7 801~7 411 cal.a B.P.、6 361~5 649 cal.a B.P.、4 844~4 754 cal.a B.P.、4 270~4 245 cal.a B.P.、4 126~4 083 cal.a B.P.、3 732~3 696 cal.a B.P.，与长江流域全新世洪水研究的结果基本一致。

武汉地区古洪水的研究意义重大，本研究利用端元建模分析为长江堤后湖泊的古洪水研究提供了新方法、新思路。同时，分离出了代表古洪水的端元模型，也为反演古洪水创造了条件。