摘要：我国西部戈壁荒漠地区含有大量砾石，这些区域土壤含水量的精准反演对我国西北地区生态和气象环境监测、植树造林及水利工程建设等均具有重要意义。然而，现有的SAR土壤水分反演模型均假设土壤由细颗粒物组成，并未考虑砾石的影响。对此，本文提出了一种利用全极化SAR遥感数据

本文内容来源于《测绘学报》2024年第11期（审图号GS京(2024)2421号）

高砾石地表全极化SAR土壤水分反演方法

1.中国矿业大学江苏省资源环境时空信息工程高校重点实验室，江苏 徐州 2211162.中国矿业大学环境与测绘学院，江苏 徐州 221116

摘要：我国西部戈壁荒漠地区含有大量砾石，这些区域土壤含水量的精准反演对我国西北地区生态和气象环境监测、植树造林及水利工程建设等均具有重要意义。然而，现有的SAR土壤水分反演模型均假设土壤由细颗粒物组成，并未考虑砾石的影响。对此，本文提出了一种利用全极化SAR遥感数据反演高砾石地表土壤水分的方法。首先，将高砾石地表后向散射建模为地面造成的面散射和砾石造成的体散射两部分；对于面散射，使用高级积分方程模型（AIEM）及Oh模型进行土壤水分反演，对于体散射，则使用致密介质辐射传输模型（DMRT）进行反演；最后，将两部分的反演结果进行加权求和，作为最终的土壤水分反演结果。利用内蒙古自治区乌海市野外土壤水分实测数据及ALOS-2全极化SAR数据对本文方法进行精度评价，并与现有的土壤水分反演方法进行比较，结果表明，本文方法（R2=0.60）比现有方法（R2=0.35）的反演精度有显著的提高。

关键词：土壤水分, 土壤湿度, 高砾石地表, AIEM模型, DMRT模型, 全极化SAR

基金项目

国家自然科学基金(41977220);武汉大学测绘遥感信息工程国家重点实验室开放研究基金(22R05);浙江省微波目标特性测量与遥感重点实验室开放研究基金(2022-KFJJ-003);新疆维吾尔自治区重点研发计划项目(2022B03003-4)

本文引用格式

郎丰铠, 何苏颖, 邱奥深, 时洪涛, 郑南山. 高砾石地表全极化SAR土壤水分反演方法[J]. 测绘学报, 2024, 53(11): 2189-2200 doi:10.11947/j.AGCS.2024.20220719

LANG Fengkai, HE Suying, QIU Aoshen, SHI Hongtao, ZHENG Nanshan.

阅读全文

土壤水分是重要的陆表参数之一，在生态、环境、农业、气象等领域均具有重要的应用价值及科学意义。目前常用的获取地表土壤水分的方法主要有3种：基于站点的实测方法、基于气象数据的模拟与同化方法，以及基于遥感数据的反演方法[1]。其中，基于遥感数据的反演方法以其监测范围、监测精度、分辨率等方面的综合优势，成为目前监测土壤水分最有效的手段。微波遥感土壤水分反演方法具有明确的物理基础，再加上微波遥感独有的全天时、全天候观测能力，使该类方法成为土壤水分遥感反演的研究热点。此外，基于微波遥感原理的GNSS-R（global navigation satellite system reflectometry）土壤水分反演技术研究也越来越多[2-4]。

综上所述，目前主动微波土壤水分反演领域的主要研究内容为如何减小入射角、地表粗糙度及植被覆盖等因素对雷达后向散射系数的影响，并未考虑土壤中的砾石。然而，我国西北地区自然地表（荒地）土壤中含有大量砾石。如图1所示，农田与荒地土壤中的砾石含量有显著区别。有研究表明，天然戈壁中，粒径大于10 mm的砾石含量达到近60%[22]。在这些区域，雷达后向散射系数除受土壤水分、地表粗糙度、植被覆盖情况及SAR系统参数的影响外[23-26]，还受到土壤层中砾石颗粒大小、形态、密度等多种因素的影响。当地表砾石颗粒较大、含量较高时，砾石间的相干散射及砾石本身的散射会显著影响地表总的后向散射，因此，现有的土壤微波散射模型不再适用。随着土壤水分在水文学、气象学、生态学等研究中的应用，人们对于该类特殊区域土壤水分监测的要求明显提高，探索一种适用于高砾石地表的高精度土壤水分反演算法，实现基于主动微波遥感的高砾石地表土壤水分反演，是目前一个迫切需要解决的问题。本文研究旨在分析砾石区土壤散射机制以及介电特性的基础上，建立高砾石区土壤微波散射模型，根据该散射模型建立高砾石地表土壤水分反演算法，最终实现高砾石地表覆盖地区土壤水分反演。

图1 不同类型的土壤表面及剖面

Fig.1 Different soil surfaces and profiles

本文基于实地调研数据及L波段全极化SAR数据，对我国西北干旱地带高砾石地表的微波散射特性进行分析，构建了高砾石地表土壤水分反演模型，提出了适用于低矮植被或无植被覆盖的高砾石地表土壤水分主动微波遥感反演方法。主要工作包含以下两个部分：①根据高砾石地表的特点，结合典型地物散射类型，研究总结已有的理论模型及半经验模型，建立适用于高砾石地表的后向散射模型，对其微波散射特性进行模拟；②根据所建立的高砾石地表后向散射模型，推导高砾石地表土壤水分反演公式，建立相应的后向散射系数数据库，结合土壤介电模型，构建高砾石地表土壤水分反演算法，实现该类特殊地表的土壤水分反演。

图2 高砾石地表雷达后向散射模型

Fig.2 Radar backscattering model of high gravel surface

（1）面散射部分。对于地表面散射部分，本文采用目前应用较多且精度较好的AIEM模型[13]，该模型能在很宽的地表粗糙度范围内精确地模拟地表散射特性。在SAR系统中，来自空气与地表土壤交界粗糙表面的后向散射信息（σs），可表示为

(1)

单基情况下，根据互易关系有

(2)

在AIEM模型中，SAR观测系统所接收到的经过目标地物单次散射的后向散射信息可表示为

(3)

(4)

(5)

(8)

则土壤表面介电常数εs可表示为

(11)

由于高砾石区土壤中砾石之间的距离要比入射波长小，砾石间的相干散射对高砾石区土壤散射量有很大的影响，即产生了多散射效应（multiple scattering effect），本文用文献[29—30]建立的致密介质辐射传输理论（dense medium radiative transfer theory, DMRT）模型来对高砾石土壤层产生的体散射部分后向散射信息进行分析。

文献[31]的研究表明，雷达波接收到的致密介质后向散射系数是电容率、入射角、体散射反照率及表面粗糙度的函数，除了雷达参数会对高砾石地表的散射强度产生制约之外，它还会受到土壤层中砾石颗粒大小、形状、密度及土壤水分等因素的影响。可表示为

(12)

式中，下标pp指的是HH和VV，其中，

，

，分别表示水平和垂直极化的功率透射系数，ε为介质的介电常数；exp-2s2（k1cosθ-k2cosθ）2为损耗函数，表示地表粗糙度对功率透射系数的影响；ω为体散射反照率，取决于高砾石土壤层的土壤密度、土壤水分、砾石颗粒形状、大小等。

一般来说，土壤中的散射颗粒为随机分布，假设颗粒为球体或类球体，体散射反照率ω与SAR的极化方式无关，则可通过VV极化和HH极化的比值将其消除，从而将比值RT表示为土壤介电常数εs和入射角θ的函数

(14)

则土壤内部介电常数εs可表示为

(15)

1.1节中，本文构建了适用于高砾石地表的微波后向散射模型，根据该模型，可计算出该类地表所对应的土壤介电常数，结合土壤介电常数模型，即可实现高砾石地表土壤水分的反演。在过去的几十年中，研究者们提出了许多用于描述土壤介电特性与含水量的模型（包括理论和半经验模型），其中在遥感定量反演领域，使用最多的是半经验土壤介电常数模型[32]

(16)

根据以上模型可知，土壤介电常数与土壤水分之间存在一定函数关系

(17)

结合式（9）和式（15）可分别推导出基于面散射和体散射反演土壤水分的关系式为

(18)

(19)

选定高砾石地表研究区，基于实测的土壤水分及地表粗糙度数据，结合式（18）和式（19），可建立起面、体散射后向散射系数与土壤水分之间的关系，进而通过全极化SAR数据实现研究区土壤水分的反演。图3为本文所提出的基于高砾石地表散射模型的土壤水分反演算法流程图。该算法具体步骤如下。

图3 高砾石地表土壤水分反演流程

Fig.3 Flowchart soil moisture inversion process of high gravel surface

（1）对全极化SAR影像进行几何校正、辐射校正、滤波等预处理，对完成预处理之后的影像进行Freeman二分量分解[33]，将总的后向散射σt分解成面散射σs和体散射σv两部分。

（4）将上述反演得到的面散射、体散射部分土壤水分按照各散射分量在总散射中所占比例进行加权求和（权重分别为A=σs/（σs+σv），B=σv/（σs+σv）），即可获取最终总的土壤水分。

本文研究区位于内蒙古自治区乌海市，其中心位置经纬度约为106.817°E，39.651°N，海拔约为1174 m，地处自治区境的西南部，且处于鄂尔多斯高原的西部和乌兰布和沙漠的东南边缘，其北部和东部与伊克昭盟为邻，西部与阿拉善盟接壤，西南方与隶属于宁夏回族自治区的石嘴山市交界。乌海市融合了华北和西北的特点，是沿黄河经济带的中心及“宁蒙陕甘”经济发展区域的接合部。在气候方面，由于研究区位于亚欧大陆的中部，具有十分典型的大陆性气候特点。在该气候条件下，研究区降水量较少，降水季节和地区分布不均匀，历年平均降水量可达到159.8 mm，平均相对湿度约为42%，平均蒸发量约为3289 mm。其中采样区域以裸地为主，植被稀疏，部分采样区有少量植被覆盖，多为沙生类低矮植被如沙蒿、四合木等。采样区位置范围及采样点分布情况如图4所示。

图4 研究区位置及采样点分布（底图来自天地图）

Fig.4 Location of study area and distribution of sampling points (the base map comes from Tianditu)

为了验证本文方法的有效性，同步收集了研究区SAR数据、光学数据及野外实测数据，具体如下。

（1）SAR数据。ALOS-2卫星是日本先进陆地观测卫星ALOS的后继星，是目前在轨运行的L波段合成孔径雷达卫星，其频率为1.2 GHz。本文试验所用ALOS-2数据为保留了相位和幅度信息的单视复数（single look complex, SLC）数据，表1给出了该数据的具体信息。分别用ENVI软件中的SARScape模块和PolSARPro软件对ALOS-2数据进行单极化影像预处理和全极化预处理，其中单极化预处理主要包括多视、地理编码、滤波，最后获得单极化影像后向散射信息；全极化预处理包括地理编码、极化分解、极化滤波等步骤，最后得到各散射分量的后向散射信息。

表1 所用ALOS-2数据参数

Tab.1

新窗口打开| 下载CSV

（2）光学数据。本次试验下载的数据为与SAR卫星过境时间大致同步的无云或少云的Landsat-8 OLI数据（L1级），该数据主要是用于计算研究区归一化植被指数（normalized difference vegetation index, NDVI）、归一化水体指数（normalized difference water index, NDWI）及植被覆盖度，以此来判断研究区植被覆盖程度并去除植被覆盖的影响。由于所获取的L1级Landsat-8数据已完成地形校正和几何校正，在使用影像之前，仅需要对遥感影像进行辐射定标以及大气校正、重采样等预处理工作，预处理操作均在ENVI 5.3软件中进行。

（3）实测数据。试验所用验证数据来源于2021年7月内蒙古自治区乌海市调研的实测采样点数据，在SAR卫星过境时间同步进行实地采样。在地势较为平坦的区域选择了高砾石地表和沙地两类地表，共测量了119个采样点，实时采集高砾石地表土壤水分、地表粗糙度、采样点位置信息（经纬度坐标）、地表植被覆盖情况、土壤层中砾石粒径大小、砾石含量情况等地面参数。其中植被覆盖沙地地表有效采样点为69个，无植被覆盖高砾石地表有效采样点为50个，高砾石地表均选择了目测粒径大于5 mm的砾石地表覆盖面积超过30%的区域，有的采样区地表砾石覆盖超过了70%。采样点分类信息见表2。

表2 采样点分类信息

Tab.2

新窗口打开| 下载CSV

其中，地表粗糙度使用探针式粗糙度板进行量测，每个采样点的土壤水分分别使用TDR仪器法、FDR仪器法、采土烘干法3种不同的方法进行多次测量（砾石含量较高地表探针无法插入，仅采用采土烘干法）。野外实测土壤体积含水量、土壤容重数据及地表粗糙度的统计结果见表3。

表3 野外采样数据信息

Tab.3

新窗口打开| 下载CSV

为了对本文所提出的高砾石地表土壤水分反演模型的可行性及在该类特殊地表的适用性进行评价分析，本次试验使用反演精度较高且适用范围较宽的改进水云模型（MWCM）结合AIEM模型、Oh模型的主动微波遥感土壤水分反演方法[20]作为对照方法。本小节将结合ALOS-2全极化SAR数据、Landsat-8 OLI数据和采样区实测数据，分别基于本文模型和对比模型，进行土壤水分反演试验。为了证明对比模型在沙地土壤水分反演的有效性并检验不同模型在高砾石地表土壤水分反演的适用性，增强本次试验的可信度，讨论了4种土壤水分反演的方案。

方案1：基于本文模型进行高砾石地表土壤水分反演试验；

方案2：基于本文模型进行沙地土壤水分反演试验；

方案3：基于对比模型进行高砾石地表土壤水分反演试验；

方案4：基于对比模型进行沙地土壤水分反演试验。

研究表明，HH极化方式的SAR信号对于植被的主干和枝茎不敏感，比其他极化方式具有很好的穿透性，能很好地表现地表土壤水分的信息，在裸露地表，一般在小入射角、L波段（1.2 GHz）、HH极化时，SAR后向散射系数基本上不受地表粗糙度变化影响，仅为土壤水分的函数[34]。因此本文在进行土壤水分反演试验时除全极化数据外，均使用HH极化的ALOS-2数据。

本文研究中无法忽略沙地采样区植被覆盖的影响。对于经过预处理操作之后的SAR影像，为了获取裸露沙地的回波散射信息，需要先使用植被覆盖地表微波散射模型来校正植被对于回波信号的影响，分别使用水云模型（WCM）和改进水云模型（MWCM）进行植被去除试验，并对两个模型的去除效果做出评估，最后选择了校正精度较高的改进水云模型（MWCM）模型。

在获取了沙地裸露地表后向散射系数与高砾石地表后向散射系数之后，按照试验设计方案开展基于实测数据的土壤水分反演试验。地表实测有效采样点共计119个点，其中高砾石地表采样点为50个，沙地采样点为69个。进行高砾石地表土壤水分反演试验时，随机选取其中的20个点进行模型的测试训练，剩余的实测点用于模型的精度验证试验；同理进行沙地土壤水分反演时亦随机选取30个采样点用于模型测试，其余采样点用于模型精度验证。

图5为不同模型分别在高砾石地表和沙地进行土壤水分反演试验的反演值与实测值的散点图分布。由图5（d）可以看出，对比模型在沙地反演结果的决定系数（coefficient of determination, R2）到0.69，说明该对比模型对于沙地等砾石含量较少的地表土壤水分计算具有一定的有效性。而从对比模型在砾石地样本集的反演结果中（图5（c），R2=0.35）可以看出，对比模型，即传统的主动微波遥感土壤水分反演方法，在应用于高砾石含量地表土壤水分反演的时候具有一定的局限性，反演精度并不理想。同时，将去除了植被影响的沙地样本集数据放入本文模型中进行反演验证，发现土壤水分的反演值与实测值之间的相关程度仅为0.41（图5（b）），说明本文模型计算的结果对于沙地土壤水分敏感度还不够高，同时也说明在砾石含量可以忽略的土壤中，若土壤不是十分的湿润，则可将整个来自于地表的散射视为面散射源，使用常用的土壤水分反演算法进行反演研究即可，不需要引入体散射分量对其进行描述分析。而本文方法反演的高砾石地表土壤水分决定系数较高（图5（a）），达到R2=0.60，相较于对比模型，反演精度有了明显的提升。

图5 土壤水分反演值与实测值的散点分布

Fig.5 Dispersion distribution of soil moisture inversion values and measured values

由图5（a）可知，当使用本文模型进行高砾石地表土壤水分反演时，在30个验证样点中，模型反演值与实测值相差0.6%以内的有23个点，相差0.6%~1.13%的有7个，整体反演精度较高。由图5（b）可知，在运用对比模型反演高砾石地表土壤水分时，模型反演值与地表实测值之间的关系总体上呈现出模型反演结果变化大于实测数据的趋势，说明该对比模型对高砾石地表土壤水分的敏感性较弱。

除了上述散点图所表现出来的决定系数之外，常用的模型评价指标还有均方根误差（root mean square error, RMSE）、平均绝对误差（mean absolute error, MAE）等，这3个指标常作为定量化评价的标准，用于评估土壤水分反演模型在不同覆盖类型地表中的适用性。一般情况下RMSE可以很好地反映模型的反演值与实测值之间偏离程度，但在实际情况下，如果存在个别偏离程度非常大的离群点，即使离群点的数量非常的少，也会使得RMSE指标变得很差。此时结合MAE评价指标能更好地反映反演的土壤水分与实测值之间误差的实际情况。

表4展示的是本文模型和对比模型在不同地表类型进行土壤水分反演的精度对比情况。由表4可以发现，对比模型在高砾石含量地表进行土壤水分反演试验的RMSE与MAE均为4次试验中的最大值，精度不甚理想，而对比模型在沙地反演的RMSE与MAE分别为0.62%和0.48%，取得了较好的结果。本文方法在高砾石含量地表反演的土壤水分的值与实测值之间的RMSE为0.6%，MAE为0.47%，是4种试验方案中精度最高的，说明本文模型反演的土壤水分值与实测值最为接近，误差最小。

表4 不同反演方法精度对比

Tab.4

为了深入理解本文模型中各参数对土壤水分反演的影响，本文结合实测数据，对本文模型中各参数关于土壤水分的敏感度进行了分析。由于本文将SAR后向散射分解成了面散射和体散射两部分，因此将对这两部分分别进行分析。

图6 面散射部分及体散射部分模型敏感度分析

Fig.6 Model sensitivity analysis of surface scattering and volume scattering

综上可知，本文采用的面散射部分后向散射模型（AIEM+Oh模型）和体散射部分后向散射模型（DMRT模型）在本文试验采用的SAR数据系统参数及采样区地表参数条件下，面、体散射后向散射系数对土壤水分的响应是敏感的，因此可用于土壤水分反演。

本文所选择的试验区有一定植被覆盖，为了获取土壤表面的回波散射信息，需要先去除植被对于回波信号的影响。目前常用的去除植被影响的模型是水云模型（WCM）和改进水云模型（MWCM），为了选择更适合试验区的模型，本文进行了植被去除对比试验。针对本次试验的沙地采样点，分别拟合出进行植被校正试验前后HH、HV、VH、VV 4个极化后向散射系数的变化。为便于展示和分析，仅选择了其中10个点，如图7所示。

图7 采样点不同极化组合去除植被影响前后的后向散射系数

Fig.7 Backscattering coefficients of HH and HV before and after vegetation correction, and their corresponding NDVI

由图7可知，经过WCM校正之后，沙地采样点后向散射系数出现了不同程度的变化，多个点有较大程度的减小。而经过MWCM校正的SAR后向散射系数基本未变。为了进一步分析，将HH极化植被校正前后的后向散射系数及对应的NDVI值列出（表5)。可以看到，2号点和9号点NDVI值较大，而其他点NDVI值较小。

表5 采样点HH极化后向散射系数及对应的NDVI值

Tab.5

新窗口打开| 下载CSV

结合图7和表5可知，WCM校正后，采样点后向散射系数的变化情况与NDVI并没有明显的关系。6号、8号点NDVI并不大，但校正后不同极化方式的后向散射系数均显著减小，而9号点NDVI值较大，但校正后HH、HV极化均增大，这显然与理论不符。

由表5可知，MWCM校正后，仅2号点和9号点略有降低。这是因为MWCM加入了植被覆盖度参数，在计算植被覆盖度时，充分考虑了土壤对NDVI的贡献，当NDVI值较小时，认为该点为裸土，只有大于一定阈值，才认为该点有植被覆盖。参考图4，结合实地情况，本文沙地采样区的确仅有个别点有少量低矮植被覆盖。因此，MWCM校正结果更符合实际情况。

综上，可以看出，WCM植被校正方法具有较大的不确定性，总体上WCM模型高估了植被影响，会导致土壤水分反演结果变小；而MWCM加入了植被覆盖度参数，并且充分考虑了土壤对NDVI的影响，因此估计较为准确，会提高反演精度。

在高砾石含量地表，砾石本身的散射及砾石间的相干散射会对SAR系统所接收到的总回波散射造成较大的影响，而现有的土壤微波散射模型及土壤水分反演方法忽略了砾石散射的影响。本文针对高砾石区地表土壤水分反演问题，考虑了土壤中砾石的形状、颗粒大小、密度等对于土壤水分造成的影响，引入体散射分量对含砾石土壤层部分的散射进行描述，从而构建了适用于高砾石地表的微波散射模型，在此基础上，进一步提出了用于反演高砾石地表土壤水分的方法。基于L波段ALOS-2全极化SAR数据和Landsat-8光学数据，结合高砾石地表采样区实测数据，对比分析了本文方法与现有方法在不同采样区域内的土壤水分反演精度。结果表明，对比方法在沙地有较好的反演效果（R2=0.69, RMSE=0.62%, MAE=0.48%），而本文所提出的高砾石地表土壤水分反演方法在高砾石地表采样点的反演结果与实测数据之间具有更高的一致性（R2=0.60, RMSE=0.60%, MAE=0.47%)。

需要注意的是，本文方法在反演过程中需要实测数据作为参考，来进行各参数值的设定，在没有实测数据时，根据其他方法得出的各参数值具有较大的偶然性，容易造成反演精度的下降。且当研究区内砾石过多及粒径过大时，可能会导致研究区内SAR后向散射对于土壤水分的变化不再敏感，从而产生较大的反演误差。因此，对于砾石含量、尺寸的界定及模型使用条件界定等一系列问题，需要后续不断地进行模型优化和更深入的研究。此外，试验中使用的数据是L波段的SAR数据，不同的波段具有不同的穿透性，进而影响到地物目标微波散射特性的分析，因此本文方法对于其他波段（如C波段等）数据的适用性仍有待进一步验证。

来源：测绘学报