密排卵砾石粗糙剖面分形特征的试验研究

钟亮1,2，廖尚超2，潘云文3，张建梅1,2，许光祥1,2，倪志辉1,2

(1.重庆交通大学国家内河航道整治工程技术研究中心，重庆　400074；2.重庆交通大学水利水运工程教育部重点实验室，重庆　400074；3.四川大学水力学与山区河流开发保护国家重点实验室，四川　成都　610065)

【摘　要】　卵砾石粗糙剖面特征对水流泥沙运动规律具有重要影响，为克服传统方法描述粗糙剖面总体特征存在缺陷的问题，本文将分形几何理论引入密排非均匀卵砾石粗糙剖面特征研究，以床面激光扫描试验资料为基础，探讨了分形尺度区宽度与中值粒径d50的关系，分析了粗糙剖面分形特征的各向异性，研究了分形维数频率分布及其与颗粒组成的关系，获得了分形维数随剖面高程统计参数的变化规律。结果表明：粗糙剖面分形尺度区宽度与床面颗粒的中值粒径d50有关，上临界尺度可取为d50；粗糙剖面分形特征具有不同程度的各向异性，中值粒径相同时，级配组成越宽，各向异性指数总体越大；分形维数频率分布具有双峰特征，中值粒径越小，颗粒组成越均匀，双峰分布越明显；分形维数随床面高程极差、标准差和峰度的增大而增大，随偏度的增大而减小。

【关键词】　密排卵砾石；粗糙剖面；分形尺度区；分形维数

基金项目：国家自然科学基金资助项目(51509026)；重庆市基础科学与前沿技术研究项目(cstc2017jcyjAX0278)；国家内河航道整治工程技术研究中心科学技术研究计划专项项目(NERC2014A02)。

作者简介：钟亮(1980—　)，男，江西赣州人，副教授，博士，主要从事水力学及河流动力学研究。

E-mail：zlcqjtu@163.com

1　引言

卵砾石颗粒是山区河流床面结构的重要组成部分，在颗粒群体运动作用下，常形成凹凸不平、形态各异的密排粗糙剖面。卵砾石粗糙剖面形态与河流动力学及河床演变学中的许多问题密切相关，如紊流结构、水流阻力、泥沙输移等，但受颗粒大小、形状、方位、级配和聚集效应等因素的影响，粗糙剖面常具有明显的非线性，如何对该类形态复杂剖面进行定量综合描述一直备受关注。粗糙剖面特征量化方法包括糙率尺寸法、暴露度法、统计参数法和分形分析法等。糙率尺寸法采用某一糙率尺寸(如床沙粒径[1]、凸起高度等)来表征剖面粗糙，该方法对于均匀颗粒具有良好的表征作用，而用于非均匀颗粒时却较难体现颗粒级配组成及随机排列的影响。暴露度法[2]考虑了颗粒相对位置和隐蔽作用，但因实测资料较少和问题较为复杂，目前非均匀颗粒暴露度随机分布公式尚不成熟。统计参数法将粗糙剖面视为高程的随机分布场，采用高程统计参数(如频率分布、构造函数或谱参数等[3-5])描述粗糙特征，此类方法虽考虑了剖面总体特征，但其量化方式较为复杂。分形分析法是随着非线性数学理论的发展和对粗糙形态认识的深入而提出来的，采用该方法可较好体现剖面特征的整体复杂性和局部精细性，且量化参数较为简单，因而也得到了一些应用，如Robert[6]、金德生等[7]、郭碧云等[8]和段光磊等[9]探讨了河道纵、横剖面的分形维数(简称“分维”)，认为可用分形维数刻画河道剖面复杂程度；倪志辉等[10]的研究发现，纵、横剖面分维反映了河道地形特征且存在一定的响应关系；这些研究成果说明分形理论可用于量化剖面特征，且主要针对河道宏观形貌。在不考虑底床宏观形貌的情况下，卵砾石粗糙剖面多属颗粒尺度区的微观分形[11]，本文部分作者[12,13]前期从随机分形模拟和分形插值模拟等角度对该问题进行了初步探讨，但现有研究主要针对数值模拟剖面，还较缺乏高精度实测资料，基于实测资料的相关研究也很薄弱，许多问题仍有待探讨，如分形尺度区宽度与颗粒代表粒径的关系、剖面分形特征的各向异性、颗粒组成对剖面分形维数频率分布的影响、分形维数与高程统计粗糙参数的联系等，因此，本文基于非均匀卵砾石床面激光扫描试验资料开展上述问题研究，并试图建立剖面粗糙特征不同量化方法的结果联系，研究成果对全面认识卵砾石剖面的粗糙特征、科学量化复杂剖面轮廓形态具有重要意义。

2　资料与方法

2.1　试验资料

调查发现，天然卵砾石床面按颗粒排列方式的不同大致可分为密铺型、散铺型、瓦叠型和散叠型等类型，其中散叠型床面由多层松散颗粒随机积聚而成，颗粒排列较为凌乱，没有明显的趋向性，是最常见的天然河床结构。因此，本文采用非均匀天然卵砾石颗粒，通过人工铺制的方式形成散叠型床面开展相关问题的研究。我国水利部颁布的《河流泥沙颗粒分析规程》(SL 42—2010)规定，砾石和卵石颗粒的粒径范围分别为2～16mm和16～250mm，本文根据试验条件选用了粒径5～21mm的卵砾石颗粒，经定制圆孔筛盘(孔径范围5～21mm，17级)筛分分组和等量配置后，形成5组非均匀颗粒用于床面铺制，颗粒粒径分组及床面颗粒参数见表1；其中床面B1、B2、B3的各组颗粒配比相同(各组号颗粒重量均占25%)，中值粒径不同；床面B2、B4、B5的中值粒径相同，各组颗粒配比不同(各组号颗粒重量分别占25.00%、12.50%和6.25%)。

为形成多层河床结构，散叠型床面铺制要求盒体净高大于2倍最大粒径，结合激光扫描采样情况，确定采用边长300mm的正方形铺制盒，盒体内部净高45mm，用平板制作并置于水平操作台调平与固定。铺制前先将各组颗粒搅拌均匀，然后将颗粒放入盒内水平方向的左端，采用刮板从左向右将床面刮平，形成散叠型床面，并定义水平铺制方向为x轴(类似于水流作用方向)，垂直方向为y轴。卵砾石床面高程采用激光扫描仪测量(图1)，测点密度通过设定测桥传感器的行走速度控制，平面采样分辨率为1mm×1mm，采样精度为0.1mm。为减小盒体边壁对床面结构的影响，研究中未考虑近壁10mm范围，选用了内部280mm×280mm区域，该区域的床面数字高程模型(Digital Elevation Model，DEM)见图2，图中长度及高程单位均为mm，高程采用了将床面平均高程设为0后的相对高程，以便直观对比，同时将x轴统一为了水平铺制方向。

2.2　计算方法

分形维数是定量描述粗糙剖面复杂性的重要参数，分形维数越大，剖面结构越复杂。作者前期研究[14]显示，结构函数法更适于粗糙剖面分形维数计算，该方法获得的分形维数在数学上定义为关联维数。结构函数法将粗糙剖面视为一个高程函数序列z(x)，定义结构函数S(ε)为序列z(x)中相距为ε(ε也称为区间尺度)的任意两点高差平方的算术平均值，则具有分形特征的序列z(x)满足

S(ε)=E[z(x+ε)-z(x)]2=cε4-2D

(1)

选用不同的区间尺度ε，对高程函数序列计算出相应的S(ε)，然后在双对数坐标系中对lnε～lnS(ε)进行线性拟合，则粗糙剖面的分形维数D与直线斜率α的转换关系为D=2-α/2。

3　结果与讨论

3.1　分形尺度区宽度与d50的关系

粗糙剖面仅在一定尺度范围内存在分形特征，该尺度范围称为“分形尺度区”(或“无标度区”)，分形尺度区反映了分形满足自相似性的范围，并对分形维数的计算结果有较大影响。分形尺度区的判定方法包括人工判定法、相关系数检验法、强化系数法、拟合误差法和分维值误差法等，本文根据前期研究成果[14]，采用人工判定和相关系数检验相结合的方法确定分形尺度区，并在此基础上分析分形尺度区宽度与颗粒中值粒径d50的关系。

粗糙剖面分形维数计算采用的结构函数也称变异函数，现有研究[15，16]认为，变异函数的最大作用尺度范围为0.5L，L为粗糙剖面长度，因此，在0<ε≤0.5L范围内选择多个尺度ε，计算对应的S(ε)，并在双对数坐标系中点绘lnε-lnS(ε)关系。lnε-lnS(ε)线性拟合效果可用平均标准偏差s和相关系数R综合评判，s越小，R越大，拟合效果越好，实际计算时，取能通过相关系数检验(显著性水平0.03)且比值R/s最大的区段作为分形尺度区。研究区域x、y向的粗糙剖面均为281个，对于给定的尺度ε，可计算各粗糙剖面对应测度S(ε)i的算术平均值S(ε)，图3给出了x、y向平均标度律lnε-lnS(ε)关系。图中显示：①各床面x、y向粗糙剖面的lnε-lnS(ε)关系曲线形状相似，级配组成相同，S(ε)随中值粒径的增大而增大(床面B1、B2、B3)，中值粒径相同，S(ε)随颗粒非均匀性的增大而减小(床面B2、B4、B5)；②lnε-lnS(ε)变化关系存在较为明显的分区，当lnε≤2.3时，lnε-lnS(ε)线性关系显著，粗糙剖面具有较强的分形特征，当lnε>2.3时，lnε-lnS(ε)接近水平线，斜率约为0，分形维数约为2，已属欧氏几何中的绝对光滑平面，粗糙剖面不具有分形特征；③存在分形特征的上临界尺度εc=10mm(lnεc=2.3)，分形尺度区宽约(0.5～1.5)d50，平均宽约d50，d50为颗粒的中值粒径，与Qin等[17]的清水冲刷粗化床面研究结果较为吻合，该尺度区即为Robert[11]和Butler等[18]定义的颗粒尺度区，反映颗粒排列形式及组织结构，ε>εc区域已属床面形态尺度区，粗糙剖面在该区域不具有分形特征。

3.2　分形特征的各向异性

为探讨粗糙剖面分形特征的各向异性，选择了图4所示的8个代表方向剖面a～h，采用结构函数法计算各方向粗糙剖面线的分形维数D(i)，令Da和ΔDmax分别为各方向粗糙剖面分形维数的算术平均值和极差值，则各向异性指数V可定义为

(2)

式中：n为粗糙剖面的数量。

表2给出了床面B1～B5各方向粗糙剖面的分形维数，数据表明：①各床面代表剖面线的分形维数范围1.230～1.910，异性指数范围0.29～0.39，差异较大；②倾斜方向粗糙剖面的分形维数总体要大于x方向和y方向，高分形维数主要集中在剖面h、b、f、d；③颗粒组成对粗糙剖面分形特征的各向异性有影响，中值粒径相同时，级配组成越宽，各向异性指数总体越大；④同一床面不同方向粗糙剖面分形维数的差异说明，粗糙剖面分形特征具有不同程度的各向异性。

3.3　分形维数的频率分布

前文指出，x向为卵砾石颗粒的平铺方向，类似于河道中的水流作用方向(纵剖面，共281个)，y向相当于与水流正交的方向(横剖面，共281个)。采用结构函数法分别计算x、y向粗糙剖面的分形维数，图5给出了粗糙剖面分形维数的区间频率分布。图中显示：①全体粗糙剖面的分形维数呈明显的双峰分布，y向分形维数总体大于x向；②级配组成相同的情况下，中值粒径越小，分形维数值基本无交叉，双峰分布越明显(床面B1、B2、B3)；③中值粒径相同的情况下，颗粒越均匀，双峰分布越明显(床面B2、B4、B5)；④若将x、y向分别视为水流运动方向及其正交方向，纵、横粗糙剖面具有不同的分形特征，说明实际河道中水流塑造作用形成的不同颗粒排列方式，对粗糙剖面分形特征有明显影响，本文试验范围内横剖面的分形维数更大，复杂性更强。

3.4　分形维数与统计粗糙参数的关系

图6给出了本文5种试验床面纵、横方向共2810个粗糙剖面的分形维数D与该剖面高程极差Rz、标准差σz、偏度Sk和峰度Ku的变化关系，图中显示：①剖面分形维数D随极差Rz和标准差σz均呈递减变化趋势，Rz和σz反映高程数据的离散程度，因此剖面高程数据越离散，分形维数越小，剖面轮廓的复杂性将减弱；②各剖面高程偏度Sk范围在-1.2～0.6之间，其中Sk<0的点据占主导，说明剖面高程分布主要呈负偏态，分形维数D总体随偏度Sk的增大而增大，由于Sk描述高程频率分布的对称性，因此剖面高程分布对称性越好，分形维数越大，粗糙剖面的复杂性将增强；③峰度Ku描述高程频率分布的陡缓性，正态分布的Ku=3，以正态分布为参照，Ku<3为低阔峰，Ku>3为高狭峰，即Ku越大分布形态越陡峭；各剖面高程的Ku值在2～10之间变化，且Ku>3的点据占主导，主要呈现为高狭峰，分形维数D随Ku总体呈缓慢递减变化趋势，说明剖面高程频率分布越陡峭，高程值分布越不均匀，粗糙剖面的复杂性将减弱。

4　结论

(1)密排卵砾石粗糙剖面具有分形特征，分形尺度区宽度与床面颗粒的中值粒径d50有关，上临界尺度可取为d50，其特征属反映颗粒排列形式及组织结构的颗粒尺度区分形。

(2)不同方向粗糙剖面的分形维数差异较大，粗糙剖面分形特征具有不同程度的各向异性，异性指数范围0.29～0.39，中值粒径相同时，级配组成越宽，各向异性指数总体越大。

(3)粗糙剖面分形维数的频率分布具有双峰特征，中值粒径越小，颗粒组成越均匀，双峰分布越明显。

(4)粗糙剖面的分形维数与剖面高程统计参数有关，床面高程极差、标准差和峰度越大，分形维数总体越小，床面高程偏度越大，分形维数总体越大。

参考文献：

[1]　YEN B C.Open channel flow resistance[J].Journal of Hydraulic Engineering，ASCE，2002，128(1)：20-39.

[2]　韩其为，何明民.泥沙运动统计理论[M].北京：科学出版社，1984.

[3]　NIKORA V,WALSH J.Water-worked gravel surfaces:High-order structure functions at the particle scale[J].Water Resources Research,2004,40,W12601,doi:10.1029/2004WR003346.

[4]　ABERLE J,NIKORA V.Statistical properties of armored gravel bed surfaces[J].Water Resources Research,2006，42,W11414,doi:10.1029/2005WR004674.

[5]　ABERLE J,NIKORA V,HENNING M,et al.Statistical characterization of bed roughness due to bed forms:A field study in the Elbe River at Aken,Germany[J].Water Resources Research,2010，46,W03521,doi:10.1029/2008WR007406.

[6]　ROBERT A.Statistical properties of sediment bed profiles in alluvial channels[J].Mathematical Geology,1988,20(3):205-225.

[7]　金德生，陈浩，郭庆伍.河道纵剖面分形-非线性形态特征[J].地理学报，1997，52(2)：154-162.

[8]　郭碧云，傅旭东，张正峰.震后龙溪河流域河床纵剖面变化和沟谷分形[J].应用基础与工程科学学报，2014，22：118-128.

[9]　段光磊，彭严波，郭满姣.河道实测冲淤量不同计算方法结果比较分析[J].长江科学院院报，2014，31(2)：108-113，118.

[10]　倪志辉，李明伟，吴立春.长江上游河道剖面分维数与通航水力要素的关系[J].水运工程，2016(6)：116-124.

[11]　ROBERT A.Fractal properties of simulated bed profiles in coarse-grained channels[J].Mathematical Geology,1991,23(3):367-382.

[12]　钟亮，许光祥.河道形态特征的随机分形模拟[J].人民黄河，2012，34(4)：30-32.

[13]　许光祥，钟亮.河道剖面轮廓的一元分形插值模拟[J].人民长江，2012，43(3)：20-23.

[14]　ZHONG L,ZENG F,XU G X.Comparison of fractal dimension calculation methods for channel bed profiles[J].Procedia Engineering,2012,28:252-257.

[15]　WEBSTER R.Quantitative spatial analysis of soil in the field[J].Advance in Soil Science,1985,99(596):537-542.

[16]　ROSSI R E,MULLA D J,JOURNEL A G,et al.Geostatistics tools for modeling and interpreting ecological spatial dependence[J].Ecological Monographs,1992,62(2):277-314.

[17]　QIN J,NG S L.Estimation of effective roughness for water-worked gravel surfaces[J].Journal of Hydraulic Engineering,2012,138(11):923-934.

[18]　BUTLER J B,LANE S N,CHANDLER J H.Characterization of the structure of river-bed gravels using two-dimensional fractal analysis[J].Mathematical Geology,2001,33(3):301-330.

Experimental Study on Fractal Characteristics of Dense Rough Gravel Profiles

ZHONG Liang1,2，LIAO Shangchao2，PAN Yunwen3，ZHANG Jianmei1,2，XU Guangxiang1,2，NI Zhihui1,2

(1.National Engineering Research Center for Inland Waterway Regulation,Chongqing Jiaotong University，Chongqing　400074；2.Key Laboratory of Hydraulic and Waterway Engineering of the Ministry of Education,Chongqing Jiaotong University，Chongqing　400074；3.State Key Laboratory of Hydraulics and Mountain River Engineering，Sichuan University，Chengdu Sichuan Province　610065)

Abstract:The characteristics of dense rough gravel profiles(RGP)have an important influence on flow and sediment movement,in order to overcome the shortcomings of traditional methods to describe the general characteristics of RGP,this study presents a fractal geometry method to study the characterization of non-uniform RGP,based on the experimental data of laser scanning of non-uniform gravel bed surfaces,the relationship between the width of the fractal scale region and the mean particle size d50 was discussed,the anisotropy of fractal features of RGP was analyzed,the frequency distribution of fractal dimension and its relationship with particle composition were studied,and the variation regularity of the fractal dimension with the statistical parameters of the elevation was obtained.Results showed that the fractal scale width is related to the median size d50 of the bed sand particles,and the upper critical value can be taken as d50.The fractal feature of RGP has different degrees of anisotropy,with the same median particle diameter d50,the larger the gradation composition,the larger the anisotropy index.The frequency distribution of fractal dimension has the characteristics of double peak value,the smaller the median particle size,the more uniform particle composition,the more obvious of the double peak distribution.The fractal dimension of RGP increases with the increase of the value of bed elevation range,standard deviation and kurtosis,and decreases with the increase of skewness.

Key words:Dense Gravel；Rough Bed Profiles；Fractal Scale Region；Fractal Dimension