半干旱草原型流域表层土壤饱和导水率传递函数及遥感反演研究*

黎明扬 刘廷玺，2† 罗艳云，2 段利民，2 张俊怡 周亚军 Scharaw Buren

（1 内蒙古农业大学水利与土木建筑工程学院，呼和浩特 010018）

（2 内蒙古自治区水资源保护与利用重点实验室，呼和浩特 010018）

（3 Application Center for System Technologies，Fraunhofer IOSB，llmenau 98693，Germany）

土壤是植被生活的基质，不仅为植物生长提供必需的养分，同时也在整个生态水文循环中起到传输和贮存水分、物质和能量的重要作用。土壤传递函数（Pedo-Transfer Functions，简称PTFs）是土壤科学中被广泛应用的一种方法，可利用有机质含量、容重、粒径组成等土壤基本特性指标，间接计算出土壤饱和导水率（KS）、水分特征曲线等常用的参数和曲线［1］。国内外学者通过大量室内和野外实验对KS和PTFs进行了广泛的研究［2-5］，Aimrun和Amin［6］在Sempadan水稻种植区域的研究表明通过以自然对数形式输入土壤干容重、颗粒百分比、有机质和几何平均粒径等因变量拟合出来的KS是最佳的PTFs形式；Xu等［7］将土壤光谱信息作为附加输入，使用人工神经网络模型（ANN）重新评估了PTFs在估算水力传导率方面的优越性。通过各种方法建立的PTFs虽然能够较为准确地预测试验点上的KS，但要估算一个流域或区域不同空间点的KS，如果去实地测试每个点上的土壤基本特性指标，再通过PTFs估算各点的KS，这无疑是费时、费力、不科学的。合成孔径雷达（SAR）问世于20世纪60年代，是目前反演和解译流域或区域土壤特征参数方面最为有效的技术［8-9］。Radarsat-2是一款搭载C波段传感器的全天时、全天候的空间微波遥感成像雷达，其作用距离远，抗干扰性能好，对云层和地面植被有一定穿透能力，该雷达被大量用于表层土壤体积含水率的反演中［10-12］，但目前利用其来反演土壤KS的研究还不多，具有重要的学术意义。

干旱半干旱地区草原生态系统十分脆弱，无论是单纯采用大面积土壤采样的PTFs还是直接利用全极化SAR数据反演的方法来预测，均无法在面尺度上准确地体现区域空间变化，因此研究一种通过PTFs与SAR相结合且对生态干扰破坏最小的新方法变得至关重要。本文采用野外原位采样与同时期Radarsat-2卫星数据相结合的方法，首先拟合了研究区表层土壤KS的PTFs，之后结合全极化SAR数据与表层土壤参数均值建立函数关系，最终预测出雷达数据研究区范围内表层土壤的KS。

1 材料与方法

1.1 研究区概况

选择内蒙古高原半干旱的锡林河流域为研究区，其位于内蒙古自治区中部锡林郭勒盟，发源于赤峰市克什克腾旗宝尔图山，流经锡林郭勒盟阿巴嘎旗，在贝力克牧场转向西北流经锡林浩特市，最终注入查干诺尔沼泽地自然消失［13］。开展野外采样的区域位于锡林河支流浩勒图郭勒与干流锡林高勒河交汇处以上流域，43°24″～44°4″N，116°17″～117°15″E，面积约为1 852 km2，区域地势三面环山，植被多为天然牧草，以茅针和羊草最为常见。研究区属中温带半干旱大陆性季风气候，多年平均降水量为266.8 mm，其中，6—8月降水量占年总降水量的50%以上［14-15］。

图1 研究区位置、采样点分布及雷达影像Fig. 1 Location and sampling points in and radar image of the study area

1.2 试验设计

首先，利用全国1∶100万土壤类型图（http://westdc.westgis.ac.cn）和全国DEM高程图（http://www.gscloud.cn）细致刻画出研究区五种土壤类型分区图和等高线分布图，在此基础上，按照每种土壤类型所占研究区的面积比例，选择多个具有代表性的样地，在雷达影像范围内的采样点共布设42个，其中分为32个建模采样点，10个验证采样点，土壤采样点垂直于等高线布设，包含研究区90%以上的海拔高度、坡度、坡向及不同地上生物生长状况。使用手持GPS记录采样点坐标，在保证选取土壤质地均一且完整的前提下，去除地表浮土及植被后，每个样点处挖50 cm深的剖面，从表层向下取样，每隔10 cm一层，每层采用100 cm3的环刀重复采集土样2个和自封袋重复样3个，分别用于测定水力参数、干容重、土壤粒径分布及其他物理化学特性。反演土壤参数和验证模型时均使用0～30 cm表层土壤数据。本次试验于2017年7月23日至30日共采集原状环刀土样420个，自封袋土样630个。

1.3 土壤理化指标测定

土壤粒径利用HELOS & RODOS激光粒度分析仪进行干法测定，分级标准采用美国制：黏粒（＜0.002 mm），粉粒（0.05～0.002 mm），砂粒（2～0.05 mm）；土壤干容重采用环刀法测定；土壤有机质采用浓硫酸-重铬酸钾外加热法测定；饱和含水率采用环刀底面缠纱布，由底部浸水测定，每隔6 h称重一次，待两次称重误差小于2%时，取平均值为其试验值；土壤饱和导水率利用定水头下的马氏瓶渗透仪进行测定。

1.4 雷达影像的处理

为保证野外采样时间与雷达成像时间对应一致，提前购置2017年7月26日RADARSAT-2精细四极化SLC（single look complex）格式雷达影像一景，相幅为25 km×25 km，分辨率为8 m。使用ENVI软件的SARscape模块对原始SLC影像进行多视（Multilooking）、Refined-Lee滤波（Refined-Lee Filtering）、地理编码和辐射定标（Geocoding and Radiometric Calibration）等处理，在Google影像上选取GCP（Ground Control Point）对影像进行几何配准后，得到标准四极化后向散射系数影像图，利用Output ROIs to ASCII功能将所需的采样点后向散射系数输出成文本文档以供使用。

由于试验前期极度干旱，连续4个月总有效降水不超过20 mm，导致试验期间研究区平均植被盖度仅为46%，平均植株高度低于10 cm，且植株长势差，各采样点平均地表生物量为58 g·m-2，故研究认为植被对雷达返回信息的影响可以忽略，未进行植被的去极化处理。具体研究区位置及高程信息、土壤类型分区、采样布设及处理后雷达影像如图1所示。研究采用同极化差SAR数据计算组合地表粗糙度参数ZS［8］，同极化差Δσo=σovv-σohh与ZS的关系可表示为

式中，A(θ)和B(θ)的表达式如式(2)及式(3)所示，θ为本地入射角，为了简化计算，统一取雷达入射角40.917°为本地入射角。

1.5 研究方法

将研究区32个建模采样点0～30 cm每隔10 cm一层的表层土壤数据分别代入Cosby［16］、Saxton［17］、Wosten［18］及非线性多元经验回归模型［19］，建立锡林河流域表层土壤KS的PTFs，通过32个建模采样点表层土壤参数均值对PTFs进行验证，并与雷达后向散射系数建立多元线性模型，利用10个未参与建模的验证采样点雷达后向散射系数反演表层土壤参数，代入PTFs中进行验证，最终选择模拟效果最佳的PTF进行区域表层土壤KS的预测。研究采用判定系数R2、均方根误差RMSE及 t 检验对模拟值进行精度分析。

2 结果与讨论

2.1 土壤参数的统计特征

研究区主要被草甸地为主的大粒径砂质土壤覆盖，地貌多变，根据各采样点所属的土壤类型，将42个参与模型计算的采样点实验数据划分为厚栗黄土、草甸沼泽土、荒漠风沙土、石灰性草甸砂土、淡黑土5种类型，并按照土层的不同深度进行分类统计，其结果如图2所示。

研究区土壤砂粒含量极高，几乎不存在黏粒，草甸沼泽土、厚栗黄土及石灰性草甸砂土砂粒含量最高，各层砂粒含量皆高于80%，荒漠风沙土最低，其平均含量仍达80.19%。在随深度变化方面，5种土壤的表层粒径变化均不明显，这说明长期的降水会让向下运移的水分不断携带体积小、粒径细的颗粒一同运动，最终导致土壤在运移较剧烈的表层达到粒径分配方面的平衡。土壤较高的砂粒含量导致其保水能力差，水分下渗速度快，这使得研究区植被稀松，有机质含量较低，这与Zhang等［20］和Yu等［21］对半干旱区植被与土壤理化性质关系的研究结果近似。在地形平坦的区域，土壤有机质含量随深度增加呈下降趋势，而主要分布在地势陡峭的厚栗黄土，其土壤有机质含量在20 cm土层有明显增高，这是因为植物为了生长在畜水困难的坡地上，深层表土根系更加发达，拥有大量毛细根的20 cm土层经过常年的积累有机质含量有明显提高。

图2 5种土壤粒径、容重及有机质含量随深度变化规律Fig. 2 Variation of soil particle size, bulk density and organic matter content with depth relative to soil type

2.2 饱和导水率土壤传递函数的适用性

土壤的结构、物理化学性质直接或间接影响土壤KS，将建模采样点表土各层数据输入到Cosby、Saxton、Wosten三种传统模型及多元经验回归模型（Empirical model）中建立研究区表层土壤KS的PTFs，代入建模采样点0～30 cm表层土壤参数的均值进行模型的检验。各PTFs形式、建模及检验的模拟实测值对比如表1和图3所示，图3中层值与均值的线性结果分别位于各分图的左上及右下。

表1 土壤饱和导水率的传递函数Table 1 Transfer functions of soil saturated hydraulic conductivity

注: 方程中，Ks为饱和导水率(m·d-1)，c1、c2、c3分别为土壤的黏粒、粉粒和砂粒含量(%)，ρ和c4分别为土壤容重(g·cm-3)和土壤有机质含量(g·kg-1)，θs为土壤饱和含水量(%)，dg为平均粒径(mm)，σg为粒径标准偏差(%) Note: In the equation, Ksstands for saturated hydraulic conductivity (m·d-1); c1for clay content (%); c2for silt content (%); c3for sand content (%); ρ for soil bulk density(g·cm-3); c4for soil organic matter content (g·kg-1); θsfor saturated water content (%); dgfor average particle size (mm); σgfor standard deviation in particle size (%)

从图3中4种PTFs的模拟值与实测值拟合关系可以发现，Cosby模型效果最差，这与Mermoud和Xu［22］的研究结果相同；其余3种模型的拟合及检验结果均高于0.958，Saxton模型拟合效果最佳，建模精度可达0.985，RMSE为0.262，Wosten与多元经验回归模型次之；各模型置信度均大于95%，说明PTFs模型在多种土壤类型模拟中不存在显著性差异并具有普遍适用性，这与前人对多种类型PTFs的研究结果相近［23-25］。

PTFs的不同形式侧面说明了模拟结果的差异性，由于Cosby模型只涉及到土壤的颗粒组成，所以拟合结果偏差较大，其主要问题表现为当实测土壤砂粒含量大于80%且土壤容重超过1.6 g·cm-3时，模型KS的模拟值会出现较大幅度的偏差，这说明Cosby模型在拟合砂粒含量和土壤容重同时较高的土壤时，具有一定的局限性。

2.3 合成孔径雷达对土壤参数的反演

图3 饱和导水率土壤传递函数模拟值与实测值对比Fig. 3 Comparison between measured and PTFs-estimated values of soil transfer function of saturated hydraulic conductivity

辐射校准后的SAR图像是目标特征的综合反映，SAR能同时实现对地面目标的距离向和方位向的高分辨率成像，全极化SAR系统交替发射两路相互正交的极化脉冲，并同时用两路相互正交的极化通道接收信号，其获得的类似光学照片的目标图像与传统光学遥感相比，具有一定的优越性［8］。使用位于雷达影像范围内建模采样点0～30 cm表层土壤参数的平均值与微波雷达RADARSAT-2四极化后向散射系数HH、HV、VH、VV及其两种组合HH/VV、HV/VH以及组合地表粗糙度ZS建立多元线性回归方程，以刻画雷达影像研究区范围内各土壤参数的分布情况。各土壤参数的结果及反演形式如表2及式(4)至式(9)所示。

表2 0～30 cm表层土壤参数反演结果Table 2 Inversion of the parameters of the 0～30 cm surface soil layer

通过表2可以看出，由于微波雷达对地物的高穿透性，各表层土壤参数与全极化SAR数据的拟合结果较为理想。土壤的粒径分布也可以较好地表示，C1作为粒径处于0.002～0.05 mm的粉粒含量，在反演中效果最差，本文通过反演效果较好的C1及C3推求获得，这样做既满足粒径分布的平衡要求，其结果也有一定的提升；θs与ρ有较好的线性关系，这与前人的研究结果一致［3,26-29］，可通过式(10)计算获得。如式(9)所示，利用SAR数据也可以直接推求表层土壤的KS，相较前文中的四种模型，拟合效果仅优于Cosby模型，这说明在结合SAR数据对KS进行面尺度推广时，使用PTFs可以获得更精准的结果。

2.4 土壤传递函数精度检验及饱和导水率的预测

利用未参与模型运算的10个验证采样点表层土壤数据的平均值进行模型的验证，结果如图4所示。研究表明，Saxton模型的雷达预测检验效果最好，模拟值最为收敛且误差最小，这说明Saxton模型更倾向于黏粒含量低的砂质土壤，这与Buccigrossi等［29］的结论相同，相较其他模型而言，Saxton模型中的参数并不是最多的，这说明参数多的模型不一定较简单的模型效果好。直接利用全极化SAR数据进行KS预测的方法精度有待提高，其模拟值收敛性不高，容重高于1.7 g·cm-1的模拟值偏低，甚至出现负值；Cosby和Wosten模型的检验结果较差，当利用SAR估计的表层土壤容重偏高时，原本对容重变化不敏感的Cosby模型模拟效果会更差，Wosten模型包含了土壤有机质含量一次方及二次方项，当表层土壤有机质含量的预测出现问题时，KS模拟值便会受到一定程度的影响，这导致在建模阶段表现尚好的Wosten模型在验证阶段失准；多元经验回归模型由于其计算了平均粒径及粒径的标准偏差，这可以在一定程度上修正SAR对表层土壤粒径分布方面的预测偏差。

选择效果最佳的Saxton模型，使用ENVI软件的波段运算工具（band math）对雷达影像研究区范围的表层土壤KS进行预测，并利用彩虹条从高到低表示。图5为雷达影像研究区范围内表层土壤KS的预测分布图。

图5反映了Saxton模型模拟的雷达影像研究区范围内表层土壤KS的分布情况，可以看出，人类活动对表层土壤的KS影响较大。圆形的滴灌区颜色深而纯，说明砂质土壤在经过农业改造后表层土壤的KS可以获得大幅度降低，最低不足1 m·d-1，表层土壤较高的持水保水能力将更有利于植物的生长；灌区附近分布着少量村庄，砖瓦和混凝土影响了水分的下渗，导致该区域土壤KS虽大于滴灌区但仍下渗较慢。河流两侧山地地势起伏大，为表面裸露、植物稀疏的岩浆岩山区，土壤中包含大量孔隙发达的火成岩，导致土壤下渗速度较快，最高可达25 m·d-1；靠近河流的河谷地区水分充足，水分下渗速度较缓，使得该区域生长了大量的水生植物，含水量较多的地表植物其叶片形状和植株倒伏特性是随机的，这会在一定程度上影响雷达对土壤参数反演的精确性，这与Millard等［30］和Liao等［31］在湿地中植被对RADARSAT-2反演土壤参数时产生较大影响的观点一致。通过预测可以发现研究区表层土壤KS的分布规律，即：山区裸地＞沙丘沙地＞河间湿地＞村镇建筑用地＞滴灌区。

图4 雷达反演饱和导水率精度分析Fig. 4 Accuracy analysis of the inversion of saturated hydraulic conductivity by radar

图5 雷达影像研究区范围内Saxton模型表层土壤饱和导水率的遥感预测Fig. 5 Remote-sensing based prediction of saturated hydraulic conductivity of the soil surface layer with Saxton model within the radar image of the study area

3 结 论

在利用SAR大面积预测土壤KS时，使用PTFs预测的结果明显优于直接利用SAR数据反演的结果；PTFs模型的拟合效果与其中包含参数的个数无关，Saxton模型更适用于砂粒含量较高、有机质含量较低的半干旱草原型流域。环境退化的草原表层土壤的KS普遍为4～8 m·d-1，自然条件下缺少植被、砂粒含量高是导水较快的主要原因，在灌区、村镇等经过开垦或建设开发的地区，表层土壤的KS大幅降低，说明人类活动是其变化的重要影响因素。而利用SAR进行区域表层土壤KS的大面积预测技术尚不成熟，有待进一步验证。