黄佳轮,安凯旋,陈汉林,吴磊
(浙江大学 地球科学学院, 浙江 杭州310027)
砾石沉积作为地质沉积的重要组成部分,其粒径变化常被用于研究气候变化、构造演化和沉积盆地的关系[1-9],评估水坝等人类基础设施对自然环境的影响[10-11],描述生物栖息环境[12-13],建立河流搬运侵蚀理论模型[14-15],探索砾石粒径顺流变化的自然现象[16-21]等。因此,准确获取砾石的粒径数据是一项十分重要的基础工作。本文详细介绍了2 类砾石统计方法的适用条件和操作过程,以便研究者在选择砾石统计方法时参考。
根据砾石统计工作的研究目的,可将砾石统计方法分为地表砾石统计法和地层砾石统计法。
在地表的二维尺度上,根据所研究砾石的粒径范围和适用条件可分为沃尔曼法与面积样本法。沃尔曼法强调随机抽样过程,在抽样区以直线等间距或网格节点的形式完成抽样过程[22-23],多适用于粒径>2 mm 的砾质河床;面积样本法可适用的粒径范围更广,可达细砂级沉积物,适用于生物栖息环境和环境监测等领域的研究[24]。随着高分辨率相机和无人机航拍技术的出现,从拍摄的照片中获取砾石粒径信息的方法在地表砾石统计中应用越来越广。由沃尔曼法所衍生的照片网格节点法[1,4]和基于计算机程序的图像自提取砾石粒径方法[25-27]是现今砾石统计的主要方式。
地层砾石统计法更关注砾石的成层性,在研究沉积物粒度对气候与构造的响应和泥沙运移规律等方面有重要应用[14,28]。应用时可根据研究目的设置相应的筛网组合确定各粒度区间的砾石占比,进而根据插值法确定粒径参数。
本文从操作流程、抽样过程、样本容量、误差分析和适用条件等方面对砾石统计方法和应用进行了详细总结, 以方便研究者根据研究目的灵活选取砾石统计方法。
粒度是指碎屑颗粒的大小,通常有2 种表示方式:粒度体积值和粒度线性值。粒度体积值是基于等效球体原理所得到的标准粒径[29];粒度线性值是指碎屑颗粒长轴、中轴和短轴(a、b和c轴)的长度。粒度体积值多应用于砂级沉积物的粒度分析,粒度线性值多应用于砾岩研究。
国际上通常用伍登- 温特华斯 (Udden-Wentworth) 的几何级粒级进行划分,且用公式Φ=-log2D将其转换为Φ值,D为颗粒直径,单位为mm(见图1)。用此方法划分粒级的优点是各粒径区间构成等间距分布。在统计砾石粒径的平均值、分选度、偏度及峰态等粒径参数时较简单方便;此外,得到的Φ值呈正态分布,在图解粒径时,能很好地展现同等级粒径的尾部分布[30]。
图1 沉积物粒级划分(Udden-Wentworth scale)Fig 1 Size gradation for sediment(Udden-Wentworth scale)
沃尔曼法强调用随机抽样方式测量地表砾石,适用于粒径为2~256 mm 的砾石[23]。主要步骤:(1)首先在抽样区以步长或者测绳建立网格区域;(2)采样人员沿着设计的路线以1~2 步长(0.3~0.6 m)为单位,随机选择脚尖的砾石,并测量其砾石中轴;(3)重复步骤(1)和(2),直至样本砾石达到100 颗。
由于该方法操作简单,无需专业设备,在早期应用广泛[22,31-33]。由传统的沃尔曼法演变而来的网格节点法,强调在选定的区域建立等间距的网格节点,统计每个节点下的砾石,随着摄影技术的发展,现今多采用拍摄照片的方式获取每个节点下的砾石数据[1-2,4-5,7,14-15,34-38]。该方法可对地层和地表进行无损采样,且野外工作量少、精度高,广受青睐。
2.1.1 随机过程
沃尔曼法强调用随机抽样方式获取砾石粒径数据,因此,随机选择的过程尤为重要。第1,研究区各点位因受地形和水力环境等因素的影响,砾石分布具有明显的各向异性,因此所选择的点位应有代表性。为此,建议在点位周围10 m2区域内,选取视觉上最细、中等和最粗的各1 m2区域进行砾石统计,将综合结果作为此点位的粒径分布,一般每个点位拍照2~10 张[1-2,5,7,35]。第2,对网格节点法而言,网格的间距也是随机抽样的重要设置参数。理想情况下,合适的网格节点应满足一个节点下只有一颗砾石。因此,为了避免粗砾石被重复统计,网格的最小间距需根据采样区的最大粒径确定,一般为砾石最大粒径的1~2 倍[39]。但在特殊情况下,很难避免大砾石(D>128 mm)被重复计数[15]。因此,为评估粗砾石对D50和D84的影响,采用无缝立方体模型[40],即对大砾石覆盖n个节点统计n次,并且将去除或添加最大砾石后重新计算的D50和D84值作为相对应的粒径区间范围[15]。粗砾石重复计数的结果表明,对D50的影响较小,对D84的影响较大,但其上限误差不会超过15%[2]。基于可接受的误差范围,网格节点法的设置通常是在1 m×1 m 的区域内均匀地建立100 个节点,统计每个节点的砾石粒径,其中对覆盖n个节点的大砾石重复统计n次。
2.1.2 样本容量
砾石统计的核心内容是统计推断,即能够通过具有代表性的样本得到该点位砾石粒径的参数,而样本容量又关乎统计工作量与精度之间的关系,因此,在样本容量选取时,既要考虑满足精度要求,又要考虑合理的统计工作量。
样本容量主要取决于砾石的分选性、可接受误差范围和要求的置信区间。前人对砾石可重复随机抽样过程的样本容量与误差有系统的描述,并建立了样本容量n与误差e的经验关系[39]:
当误差表现形式为绝对误差且粒径单位为Φ时,
式(1)经相应变换,同样适用于mm 粒径下的百分误差计算:
粒径的Φ值更符合正态分布,当考虑粒径分布特征的影响时,
其中,t 表示学生分布的t 值,砾石统计中一般取置信度95% 且自由度(n-1)趋于无穷时的t 值(t =1.96);sI表示英曼分选系数,sI= |Φ84-Φ16|2;s表示粒径分布的几何标准差,s= (D84D16;e±Φm表示Φ均值粒径的绝对误差;e%Dm表示mm 均值粒径的百分误差。
目前,大多数研究普遍将砾石统计的可接受误差定为10%~15%,沃尔曼法的样本容量≥400(见图2);此外,基于二项分布的基本理论,样本容量与置信水平的相关研究也表明,当样本容量增至500时,D50与D84的误差分别 降至11%和9%[41]。 综上所述,在可接受的误差范围内(<15%),沃尔曼砾石统计法的样本容量应大于400。
2.1.3 误差分析
传统的沃尔曼法受人的主观影响较大[39],且无误差比对基准,本节主要针对目前广泛使用的照片网格节点式砾石统计法进行误差分析。照片网格节点法的误差来源主要有:粒径轴的判断问题、砾石沉积结构带来的影响以及最小截断粒径带来的误差。
2.1.3.1 粒径轴判断对误差的影响
在砾石统计中,通常将砾石b轴作为砾石粒径参数,由于受水动力的影响,野外砾石沉积出露的多为其最大扁平面,即a-b轴面。但在二维照片中,并不能完全排除获取的为砾石的b-c轴面[36],因此,在照片中不能对砾石的b轴进行明确的标记。前人对华盛顿Elwha River 31 个点位的512 个样本开展了野外砾石统计和照片网格节点法砾石统计,结果表明,照片中得到的主轴(PM)和次轴(Pm)与野外实测砾石的长轴(FL)和中轴(FI)具有高相关性,相关系数分别为0.99 和0.98,而和实测短轴(FS)的匹配度较低[27](见图3),表明砾石在垂直方向上以短轴为主。因此在二维平面上若要获取粒径的中轴(b轴),只需测量平面上的次轴(Pm),其强相关性也表明,粒径轴的判断在照片砾石统计工作中对误差影响较小。
2.1.3.2 砾石沉积结构对误差的影响
如果照片上的所有砾石b轴均完全可见,并且平行于拍摄平面,那么照片所得到的砾石粒径分布结果与用卷尺对砾石b轴进行测量得到的粒径分布是极为相似的。但由于受半掩埋、叠瓦构造、投影缩短等的影响(见图4),大多数砾石都不是理想的拍摄对象。关于对这类误差的改进,前人已有相关研究,建议通过在网格节点法的结果中增加5 mm[40]或乘以校正因子1.07 进行误差校正[38]。但由于缺乏误差的比对基础,且这并不是造成误差的主要原因,因此前人都将其归为可忽略的误差范畴。更具体的定量化分析尚待进一步研究[42]。
2.1.3.3 最小截断粒径对误差的影响
由于受照片精度和砾石统计目的的约束,一般将1~8 mm 作为统计下限[1,4,9],即最小截断粒径。由于最小截断粒径造成的误差对高百分位粒径(如D84,D50)影响较小,但对低百分位粒径(如D16,D5)影响较大[43]。文献[42]通过对加拿大境内3 条河流的74 个样本,共计约22 200 颗砾石的统计结果表明,合适的最小截断粒径取决于研究目的,对含沙量<5%的沉积物而言,当截断粒径为8 mm 时,D50误差为 -0.06Φ;而当截断粒径为32 mm 时,D50误差则增至-0.20Φ,当含沙量增加时(含沙量>5%),截断粒径带来的误差明显增大(见图5)。因此当砾石统计工作主要关注于粒径>2 mm 的高百分位(D50,D84)砾石,如研究冲积扇顺流粒径变化趋势[4,7,35]、探索古水流坡度[36]等时,最小截断粒径所带来的误差几乎可以忽略不计;而当统计工作更关注于低百分位粒径(D16,D5)砾石,如研究生物栖息环境[12]、山间河流的砾-砂转换带[44]等时,最小截断粒径所带来的误差则不可忽视,此时建议采用面积样本法和筛分法[31,45]。
图2 沃尔曼法样本容量与误差的关系Fig.2 The relation between Wolman sample size and error
面积样本法,在抽样面积范围内,统计地表所有砾石b轴,一般选定的区域面积为0.1~1 m2。对于面积样本法,其操作过程主要分为枚举法[46-48]、黏合剂法[49]和图像自提取法[25-27,42,50-60]3 种。
图3 野外实测粒径数据(F)与照片获取粒径数据(P)对比Fig.3 Comparisons of field(F)and photograph(P)measurements of gravel size
图4 砾石沉积结构导致砾石粒径偏小的三个潜在因素(据文献[42]修改)Fig.4 Three potential sources of grain size underestimation resulting from sediment structure and the use of photographic sampling methods(modified after[42])
2.2.1 枚举法与黏合剂法
枚举法,对野外抽样面积内的所有砾石手动拾取测量,适用于粗砾石(>8 mm)的统计。其最大的问题是难以确定所统计的砾石属于地表还是近地表。同时,在大砾石被移除后会出露更多小砾石,因此,该方法一般是从小砾石开始统计。前人对该方法做了一定的改进,其中一种改进方法是在抽样面积内对地表砾石进行喷漆,然后对沾染油漆的砾石进行测量,但此方法不能完全识别地表砾石,因为油漆可能沿着岩石的一侧流下,渗入地下沉积物[39]。另一种改进方法是将地表砾石喷上铁磁矿粉涂料,然后用强磁性的手持磁铁将所有具有磁性的砾石提取出来[48]。
图5 当截断粒径分别为2,8 和32 mm 时每5%砾石的粒径百分位均值误差Fig.5 The mean error in every 5th percentile when the truncated particle size are 2,8,and 32 mm,respectively
黏合剂法是一种适用于含砂砾石(通常适用于砾石粒径<32 mm)的统计方法。黏合剂的取样过程是将一块涂有黏合剂的纸板压在砾石表面,当黏合剂完全固结时,将纸板抬离地表,地表的砾石就会黏附在黏合剂上,之后将样本砾石与黏合剂分离,即可得到砾石粒径数据。常用的黏合剂有万能胶、环氧树脂、泥浆、黏土、油脂、蜡等[49,61]。但黏合剂取样技术一直饱受争议,因为黏合剂不能提供一致的穿透深度,会导致取样深度不一致,在操作过程中也难以确定黏合取样的砾石一定来自地表[39,48]。
2.2.2 图像自提取法
无论是沃尔曼法、枚举法还是黏合剂法,都需要繁重的野外工作及大样本才能准确估计每个点位的粒径数据。为了克服样本的约束,基于照片的粒度自测定方法在过去几十年里得到了长足的发展和广泛应用[26,56,60-62]。从图像中以面积样本法获取粒径数据的最简单方法是手动在Photoshop 中测量所有可见砾石的长轴或短轴。此方法与沃尔曼法的照片网格节点法类似,不再赘述。
图像自提取法,是基于程序算法自动识别图像中全部砾石b轴的粒径分布特征的方法。早期为了从图像中确定砾石大小,须对每颗砾石进行手工数字化处理[63],随着计算机科学的迅速发展,现今的图像数字化处理方法采用的是纹理识别法和图像分割法。
纹理识别法依赖于纹理之间的关系,其原理是图像上相邻像素的相关程度随颗粒大小的变化而变化,主要识别方式有半方差[64-65]、灰度共生矩阵[66]、自相关算法[26-27,53]等。这些方法只能得到砾石粒径的均值(D50),而且需要通过特定点位的砾石数据来校正,多适用于砂级粒度的测量。
相较纹理识别法,图像分割法的重点是对每个可见砾石进行完整的描述和测量[42,50,58-59,67-68]。该方法的局限性在于当图像较为复杂,如砾石重叠、砾石形状不规则、砾石颜色变化等时,误差很大;其优势在于不需要野外点位砾石数据的校正,且得到的是一个完整的砾石粒径分布曲线,可以得到所需的粒径百分数(D5,D16,D25,D50,D75,D84,D95)。
基于图像自提取砾石粒径的方法如表1 所示。
表1 基于图像自提取砾石粒径方法汇总Table 1 Summary of self-extracting gravel particle size methods based on image
2.2.3 样本容量
相较线性或网格统计法,面积样本法统计了采样区域所有表面砾石的粒径。因此,面积样本法的样本容量取决于区域面积而非砾石数量[42],而抽样区域的大小一般为区域内最大砾石的面积函数。前人的研究表明,面积样本法的抽样面积应该不小于抽样区内最大砾石面积的100 倍[45],此外,为了将误差降至10%以内,抽样面积应增大至最大砾石面积的400 倍[31];而威尔士Afon Ystwyth 地区的砾石统计结果表明,要得到所有粒径的完整分布,且误差控制在5%以内(粒径单位为mm),则抽样面积至少应为最大砾石面积的200~400 倍,若误差只需控制在10%以内(粒径单位为mm),抽样面积可以减半[42](见图6)。
2.2.4 误差分析
传统的面积样本法(枚举法和黏合剂法)的误差主要在于无法准确判断所统计的砾石来自于地表还是近地表。枚举法的误差主要取决于人的主观因素,一是因为操作者在砾石统计过程中更倾向于选取粗砾石,二是因为即使强调先统计选定区域的细砾石,但依然无法确定所统计的砾石来自于地表还是近地表;黏合剂法则强调根据砾石的粒径特性选择黏合剂的类型,黏合剂的渗透深度是关键影响因素。综上,传统的面积样本法其误差无比对基准,且误差大小因人而异。
图像自提取砾石粒径数据方法与照片网格节点法一样,存在粒径轴的判断问题、砾石沉积结构带来的影响及最小截断粒径所带来的误差(见图4)。其中,粒径轴的判断问题与照片网格节点法相同,此处不再赘述。自首次使用数字化方式获取砾石粒径以来,砾石沉积结构这一误差源得以确认,但其误差大小一直没有确切的结论[63]。文献[42]对法国Ain River 的10 个0.6 m2区域的枚举法(喷漆式)和手动数字化的粒径结果分析表明,在砾石沉积结构影响下,手动数字化比枚举法(喷漆式)得到的粒径结果要大(见图7),导致这一结果的原因可能是喷漆式枚举法更倾向于高估细砾的数量,而手动数字化则会低估细砾的数量。因此,相比于其他影响因素,砾石的沉积结构所带来的误差相对较小。
图6 面积样本法D50和D90 百分位误差和抽样面积与最大砾石面积比值的关系(据文献[42]修改)Fig.6 Plot of percentile errors versus the ratio between sampled area and the area of the largest grain in the population(Dmax)(modified after[42])
图7 手动数字化与枚举法(喷漆式)百分位粒径结果对比(据文献[42]修改)Fig.7 Comparison between manual digitization and enumeration method(paint-and-pick sampling)of percentile particle size(modified after[42])
图像自提取砾石粒径中的纹理识别方法的最小截断粒径主要取决于所采用的编程算法和图像分辨率,此外,拍摄时的灯光及沉积物颜色等也会对砾石统计的精度产生重要影响[50]。当照片精度为0.7 mm/像素时,将23 个像素(大约16 mm)作为截断粒径是比较合适的[58];若将20 个像素作为截断粒径,研究结果表明,当照片精度为1.16 mm/像素时,随机误差为0.15Φ;当照片精度为0.32 mm/像素时,随机误差则降为0.05Φ[60]。
综上所述,面积样本法的误差源主要来自编程算法和像素精度,因此,如何更精确地识别砾石轮廓是编程算法应关注的方向。
地层砾石主要强调砾石的成层性分布,多适用于河床监测、泥沙运移分析以及沉积地层中的砾岩对构造/气候响应的研究[14-15,28]。层状砾石多采用体积样本法获取粒径数据,在野外工作中,一般用网格筛得到粒径>1 mm 的砾石数据,而粒径<1 mm 的砾石通常需在实验室条件下进行更细的网格筛分,或借助粒度分析仪获取。
对于缺乏泥沙供应量的砾石河床,地表的砾石比近地表的砾石粗,因此,在地表形成具有一定厚度的“铠甲层”(armor layer)[69-70](见图8) 。而在河床监测以及泥沙运移等方面,地表砾石“铠甲层”的砾石粒径是判断水力环境的重要参数。由于“铠甲层”是具有一定厚度的砾石层,通常采用体积样本法对其抽样统计,即从预选区域选取一定体积的砾石,将其
筛分得到砾石粒径数据,因此,体积抽样的深度很关键。前人建议将最大砾石粒径(Dmax)的1~2 倍作为“铠甲层”与近地表砾石的界限[35,44-46],但通常将最大砾石粒径(Dmax)作为“铠甲层”砾石的深度[35,44]。
图8 砾质河床砾石的垂向分层Fig.8 Vertical stratification of gravel size in gravel bed river
沉积砾岩的粒径变化记录了对应地质时期的气候变化和构造活动信息[28,71],因此,研究其粒径变化特征具有重要意义。在沉积砾岩粒径研究中,对于松散砾岩,体积样本法是获取粒径数据的最好方法[28],在某些特殊情况下(如砾岩固结很好),也可采用面积样本法或网格节点法[1,72]。
服从正态分布的样本表明,准确描述尾部分布需要比描述中心分布更大的样本。因此,一个足够描述尾部分布的样本也足以描述粒度的正态分布特征。砾质河床粒径分布的粗尾只由少量大砾石组成,但它们在样品总重量中的占比却很大。大砾石的存在与缺失不仅影响粗尾的百分位粒径,而且也影响中值粒径。因此,样本体积需要足够大,以便将大砾石包含在具有代表性的样品中。样本容量由最大砾石的质量决定[39,45]。由于砾石质量是粒径的三次方函数,因此,样本容量也定义为最大砾石粒径的三次方函数。前人在砾石统计工作基础上,总结出样本容量与最大砾石粒径的经验关系[39],这一抽样标准在砾石研究中得到了广泛应用[15,35]:
(a)当最大砾石粒径<32 mm,其质量占比应不超过0.1%,
(b)当最大砾石粒径为32~128 mm,其质量占比应不超过1%,
(c)当最大砾石粒径>128 mm,其质量占比应不超过5%,
其中,Ms为样本容量的质量(kg),Dmax为最大砾石粒径(m)。
体积样本法是将抽样样本经过一系列方格筛网筛选后确定样本各粒度区间所含砾石的质量(筛网大小及间距由研究目标而定),从而根据插值法确定各粒径的百分位大小。方格筛网法直接将筛网的边长s作为对应砾石的b轴长度(即b=s,b为砾石b轴的长度,s为筛网边长),实际上只有当砾石的b-c截面是圆形时才成立,而自然界中存在大量b-c截面为类椭圆形的砾石,其真实b值往往比筛网边长s要大(见图9)。因此,方格筛网方法往往会低估砾石的b值。为了矫正这一误差,需要测定对应的校正因子(k值)。
图9 方格筛网与砾石b-c 截面关系示意图Fig.9 Schematic diagram of the relationship between square-hole sieve and the b-c axis of gravel
在方格筛网方法中,直接将通过筛网的最大b值作为筛网边长s,即b=s;而实际情况是通过筛网的最大b值与s满足函数关系:
假设砾石的b-c截面为椭圆,参数方程为
椭圆离心率k=当b值最大时,筛网方程与椭圆方程相切。由式(8)和式(9),可得解方程(10),得
由此可得原方格筛网法测得的b轴长度的误差为
由式(12)可得,砾石b轴的测量误差取决于离心率的k值,而k值是基于研究区域的经验值确定的。文献[14] 对喜马拉雅山脉南缘的马斯岩迪河(Marsyandi River) 的砾石统计研究表明,k的均值为1.9;帕米尔东北缘的乌泊尔冲积扇与青藏高原北缘的疏勒河冲积扇k的均值分别为1.57 和2.20(待发表数据)(见表2)。由式(12)可得,以上3 个例子中,b轴的统计结果误差分别为-20%,-16% 和-22%(负值表示测量值小于真实值)。
以上结果表明,在沉积物磨圆度不高的沉积环境中,如冲积扇扇根、冲积河流等,筛分法产生的误差较大。因此,在筛分法使用中,需配合沃尔曼砾石统计法辅助计算研究区的k值,以校正体积样本法带来的误差[14]。
粒度分析方法与野外工作紧密结合,具有很强的实践性和针对性,且所用方法由研究目标的粒径范围及野外工作条件决定(见表3)。如在探索河床砾石顺流变化的自然现象时,研究的重点在于粒径>2 mm 的砾石,且河床砾石的顺流变化不仅体现在粒径上,还体现在磨圆度和砾石成分等上,据此可判断在粒径顺流中占主导地位的是磨蚀作用还是选择性沉积[18,73-74],而沃尔曼砾石统计法则能很好地满足这些条件;在环境监测、泥沙运移和描述生物栖息环境等领域,本着无损取样原则,使用面积样本法是十分必要的[24,75];在研究冲积河流的砾-砂转换带时,自然界河流在顺流方向普遍存在由砾(~5 mm)向砂(<2 mm)的突变[76],因此,对上游的砾质河床多采用沃尔曼法,对下游的砂质河床则多采用体积样本法[44,77]。
此外,由于砾石统计的随机性和多样性,若条件允许,采用不同的砾石统计方法对同一位置的统计结果进行对比是非常必要的,以便对不同的砾石统计方法进行互检,进而改进方法。
表2 不同区域砾石k 的均值Table 2 The mean of k of gravel in different regions
表3 基于研究目的的砾石统计方法总结Table 3 Summary of gravel statistical methods based on research purposes
系统总结了地表砾石统计法和地层砾石统计法。地表砾石统计法又分为沃尔曼法与面积样本法。沃尔曼法强调抽样的随机过程,通常以网格节点的形式选择目标砾石,适用于对粒径2 mm 以上砾石的统计,在可接受误差区间(<15%),样本容量需大于400,多用于探索砾石粒径的顺流变化、研究河流搬运侵蚀规律和冲积河流砾-砂转换带等领域;面积样本法较沃尔曼法的优势在于适用粒径更广,可达细砂级别,而且可以做到无损采样,主要方法为图像自提取法,样本容量需达最大砾石面积的200~400 倍,多用于环境监测和生物栖息地划分等,但因受计算机程序算法、沉积物特性等的影响,在部分领域得到的粒径数据误差较大。
地层砾石统计法多使用体积样本法,样本容量需根据最大砾石的粒径和质量确定,适用于砂砾级的沉积物粒度分析,多用于河床监测、泥沙运移变化以及地层中沉积物的粒度变化研究等。在沉积物磨圆度不高的地区,需通过沃尔曼法测定校正因子k,以校正粒径数据。