1、引　言

土壤水分是全球水循环和能量交换过程中一个关键变量，在水文、气候和农业生产中扮演着重要角色(Bai 等，2015；Kim和Van Zyl，2009；周鹏 等，2010)。传统的土壤水分监测方法虽然精度很高，但很难大范围、高效率以及全过程地获取水分信息(Bai 等，2015；Kim和Van Zyl，2009；周鹏 等，2010)。而遥感技术的发展则为不同时间、空间尺度的土壤水分变化监测提供了有效的途径(Bai 等，2016；Kornelsen和Coulibaly，2013)。遥感监测土壤水分的方法主要有：热惯量法(余涛和田国良，1997)、植被指数法(吴海龙 等，2014)、光谱特征空间法(Jackson，1983)以及微波遥感法等(赵昕 等，2016)。可见光和红外遥感易受到天气的限制，而微波遥感具有全天时、全天候的监测能力，受云、雨、大气等影响较小，因此相较其他遥感监测方法具有先天技术优势。

荒漠化生境地表植被对土壤水分变化极为敏感，而荒漠化采矿塌陷区，植被生态更加脆弱易损，目前这一领域的研究相当薄弱。仅有卞正富等人(2009)利用野外实测和TM影像反演方法从工作面尺度分析了采空区与非采空区上方的土壤含水量，发现实测微观尺度上两个地区的土壤含水量有着明显的差别，而遥感面域尺度上两个地区的土壤含水量并没有明显的差别。另一些研究则立足于矿区尺度，如刘英等人(2011)利用MODIS数据的NDVI产品以及地表温度数据评估了2010年10月8日当日神东矿区的旱情；马保东等人(2011)利用MODIS数据对神东矿区进行了9年的土壤水分监测，发现矿区与背景区相比，土壤并未发生干化；刘英等人(2016)进一步利用12期TM/ETM+/OLI和HJ-CCD影像反射率数据，监测分析了过去25年神东矿区土壤湿度时空变化特点。然而，现有的研究没有严格界定开采区与非开采区，且研究手段较为单一，多是利用光学数据反演植被冠层及土壤表层水分，并没有获得采矿造成土壤水分变化的规律性，更没有将开采沉陷与土壤水分变化建立起必然的联系。

本文以内蒙—陕西交界处的干旱—半干旱矿区为研究区，首先利用Sentinel-1A VH，VV双极化后向散射系数数据，结合Landsat 8 OLI地表反射率数据，采用水云模型WCM(Water-Cloud Model)分离出研究区植被散射对于雷达后向散射系数的贡献，得到裸土后向散射系数；然后根据李震等人(2006)、韩桂红等人(2013)依据AIEM模型(Chen 等，2003)和Oh模型(Oh，2004)所建立的经验模型得到研究区土壤水分空间分布图，将此结果与基于Landsat 8 OLI地表反射率数据计算的土壤水分结果进行采样对比，验证反演结果的合理性；最后基于Sentinel-1A的DInSAR结果，重点分析了沉陷区内、外的土壤浅层水分含量变化，探究了长时序下高强度开采对于地表土壤水分含量的扰动效应。

2、研究区与数据预处理     (2.1) 研究区概况

研究区位于内蒙古自治区和陕西省交界处，东经109°12′48″—110°46′05″，北纬38°32′24″—39°38′08″，是黄土高原和毛乌素沙漠的过渡带，属于典型的干旱—半干旱大陆性气候，年降水量在400 mm左右，土壤贫瘠，风蚀沙化比较严重。植被类型多是些低矮、稀疏的旱生、半旱生草本植物，生态环境十分脆弱(马保东 等，2011；Ma 等，2016)。研究区覆盖了整个神东主矿区，矿区大规模、高强度的矿产资源开采活动不可避免的会影响到当地生态系统，而土壤含水率对于植物生长至关重要，因此，研究干旱—半干旱区，尤其是荒漠化矿区的地表土壤水分对于当地生态保护、恢复与重建具有重要意义(刘英 等，2016)。

    (2.2) 数据及预处理          2.2.1. 雷达遥感数据(Sentinel-1A)处理

Sentinel-1A卫星是欧洲航天局哥白尼计划GMES(Global Monitoring for Environment and Security)中的一颗地球观测卫星，于2014年4月3日发射上空。卫星载有C波段SAR传感器，工作频率为5.405 GHz，重访周期12 d，共有4种观测模式，本文获取的研究区数据为干涉宽幅IW(Interferometric Wide Swath)模式的斜距单视复数SLC(Single Look Complex)产品，数据来源于欧空局(https://scihub.copernicus.eu/[2017-07-21])，分辨率为5 m×20 m，主要参数见表1。

(1)沉陷区提取。为了对比沉陷区内外的土壤含水量变化规律，可利用Sentinel-1A数据，基于DInSAR技术获取研究区内受开采扰动引起的沉陷区。ENVI5.3下的高级雷达处理扩展模块SARscape5.2.1提供了完整的SAR数据处理功能，可提取毫米级的地表形变。本文利用2016-02-10和2016-03-29两期SLC数据，在SARscape5.2.1下完成参数设置、数据导入、DEM数据下载和基线评估等前期准备，然后利用软件提供的DInSAR工作流工具经干涉图生成、滤波和相干性计算、相位解缠、控制点选取、轨道精炼和重去平、相位转形变以及地理编码等过程获得地表形变图(图2)。研究区在这一时期内最大下沉值为0.079 m，矿区内的沉陷区比较明显，但也要注意排除大气扰动、水体等因素造成的干扰。

(2)后向散射系数生成。Sentinel-1A的后向散射系数生成主要包括了多视、滤波、辐射定标和地理编码等过程。采用窗口大小为7×7的Refined Lee滤波去除斑点噪声(何连 等，2016)，通过地理编码和辐射定标(王娇 等，2017)生成具有地理坐标信息的后向散射系数图像(图3)。

         2.2.2. 光学遥感数据(Landsat 8)处理

为了去除植被对雷达后向散射系数的影响，进一步验证Sentinel-1A的土壤水分反演结果，结合空间分辨率、时间、云量等条件，研究中同时下载了Landsat 8 OLI地表反射率数据，数据来源于美国地质勘探局(https://earthexplorer.usgs.gov/[2017-07-21])，此数据不需要再做辐射定标和大气校正处理，主要参数见表2。

由于Landsat 8 OLI影像与Sentinel-1A影像的坐标系和空间分辨率不统一，因此需要对OLI影像进行重投影和重采样预处理。

3、研究方法     (3.1) Sentinel-1A土壤水分反演          3.1.1. 水云模型去除植被影响

当研究区内NDVI>0.4时，植被对于雷达后向散射系数的贡献就不可忽略(Neusch和Sties，1999)。利用2016年8月4号的Landsat 8 地表反射率影像计算NDVI值，平均结果达到了0.406，因此植被对于雷达后向散射系数的影响不可忽略。Attema和Ulaby(1978)根据分离植被和地表对于后向散射系数贡献的思想，建立了水云模型，后来经过不断修正和改进，在给定入射角 $\theta $ 下，模型可表达为(Bindlish和Barros，2001)

$\left\{ {\begin{array}{*{20}{l}}\!\!\! {{\sigma ^{\rm{0}}} = \sigma _{{\rm{veg}}}^0 + {\tau ^{\rm{2}}}\sigma _{{\rm{soil}}}^{\rm{0}}}\\[7pt]\!\!\! {{\tau ^{\rm{2}}} = \exp \left( { - 2B{m_{\rm{v}}}\sec \theta } \right)}\\[7pt]\!\!\! {\sigma _{{\rm{veg}}}^0 = A{m_{\rm{v}}}\cos \theta \left( {1 - {\tau ^2}} \right)}\end{array}} \right.$

(1)

式中， ${m_{\rm{v}}}$ 是植被含水量； ${\sigma ^{\text{0}}}$ 、 $\sigma _{{\rm{veg}}}^0$ 和 $\sigma _{{\rm{soil}}}^0$ 分别代表总后向散射系数、植被后向散射系数和土壤后向散射系数； ${\tau ^{\text{2}}}$ 是雷达波穿透植被层的衰减因子； $A\text{、}B$ 是依赖于植被类型的参数，可通过已有研究得到，如表3所示。

本文选取所有植被类型参数，即A=0.0012，B=0.091。

张友静等人(2010)利用NDWI、EVI和NDVI分别反演了小麦含水量，发现NDWI反演精度要优于后两者。因此，本文利用Landsat 8地表反射率影像求取NDWI(Normalized Difference Water Index)，然后反演植被含水量。NDWI(Jackson 等，2004)可通过ENVI软件中的波段运算工具求取

${\rm{NDWI}} = \left({{R_{{\rm{NIR}}}} - {R_{{\rm{SWIR}}}}} \right)/\left({{R_{{\rm{NIR}}}} + {R_{{\rm{SWIR}}}}} \right)$

(2)

式中， ${R_{{\rm{NIR}}}}$ 和 ${R_{{\rm{SWIR}}}}$ 分别代表近红外和短波红外地表反射率，对应于Landsat 8影像的第5和第6波段。考虑到研究区内多为低矮植被，根据Jackson等人(2004)的研究，可建立研究区植被含水量反演模型

${m_{\rm{v}}} = {\rm{1}}{\rm{.44NDW}}{{\rm{I}}^2} + 1.36{\rm{NDWI}} + 0.34$

(3)

根据式(1)、(2)、(3)以及A和B值可去除植被对于雷达后向散射系数的影响，求出裸土对于总雷达后向散射系数的贡献值(图4)。

         3.1.2. Sentinel-1A土壤水分反演模型

在雷达系统参数确定的条件下，后向散射系数值主要受到地表粗糙度、土壤含水量以及入射角的影响，可以表示为

$\sigma _{{\rm{pq}}}^0 = f\left({s, l, {M_{\rm{v}}}, \theta } \right)$

(4)

式中， $\sigma _{{\rm{pq}}}^0$ 是雷达后向散射系数； ${\rm{pq}}$ 是极化方式； $s$ ， $l$ 是地表均方根高度和相关长度； ${M_{\rm{v}}}$ 是土壤含水量； $\theta $ 是入射角。

对于同景多极化雷达数据，任鑫(2004)和李震等人(2006)结合以往的研究再根据大量的模拟分析，建立了地表参数反演模型。利用AIEM模型模拟的后向散射系数值代入该模型反演的水分值与AIEM模型输入水分值相关系数以及该模型反演值与实测土壤水分值相关系数都达到了0.9以上，从而验证了模型的可靠性，该模型可表示为

$\left\{ {\begin{array}{*{20}{l}}\!\!\! {\Delta \sigma _{{\rm{VV}}}^0 = {A_{\rm{V}}}\left(\theta \right)\ln \left({{Z_{\rm{S}}}} \right) + {B_{\rm{V}}}\left(\theta \right)}\\\!\!\! {\sigma _{{\rm{VV}}}^0 = {A_{{\rm{VV}}}}\left(\theta \right)\ln \left({{M_{\rm{V}}}} \right) + {B_{{\rm{VV}}}}\left(\theta \right)\ln \left({{Z_{\rm{S}}}} \right) + {C_{{\rm{VV}}}}\left(\theta \right)}\end{array}} \right.$

(5)

式中， ${Z_{\rm{S}}}$ 是组合粗糙度，可以达到减少未知数，简化模型的目的，可根据Zribi和Dechambre(2003)的研究表达为 ${Z_{\rm{S}}} = {s^2}/l$ ； $\Delta \sigma _{{\rm{VV}}}^0 = \sigma _{{\rm{VH}}}^0 - \sigma _{{\rm{VV}}}^0$ ，是两种极化的后向散射系数差； ${A_{\rm{V}}}\left(\theta \right)$ 、 ${B_{\rm{V}}}\left(\theta \right)$ 和 ${A_{{\rm{VV}}}}\left(\theta \right)$ 、 ${B_{{\rm{VV}}}}\left(\theta \right)$ 、 ${C_{{\rm{VV}}}}\left(\theta \right)$ 是根据模拟结果拟合出的值，只与入射角有关，表达式可参考文献(李震 等，2006)。

模型(5)只能反演VV同极化数据的土壤水分含量，韩桂红(2013)则结合AIEM和Oh模型，给出了VH交叉极化数据反演土壤含水量的模型。综上，可得到研究区地表参数反演方法

$\left\{ {\begin{array}{*{20}{l}}\!\!\!\!\! {{Z_{\rm{S}}} = \exp \left({\left({\Delta \sigma _{{\rm{VV}}}^0 - {B_{\rm{V}}}\left(\theta \right)} \right)/{A_{\rm{V}}}\left(\theta \right)} \right)}\\\!\!\!\!\! {{M_{\rm{V}}}\left({{\rm{VV}}} \right) = \exp \left({\left({\sigma _{{\rm{VV}}}^0 - {B_{{\rm{VV}}}}\left(\theta \right)\ln \left({{Z_{\rm{S}}}} \right) - {C_{{\rm{VV}}}}\left(\theta \right)} \right)/{A_{{\rm{VV}}}}\left(\theta \right)} \right)}\\\!\!\!\!\! {{M_{\rm{V}}}\left({{\rm{VH}}} \right) = \exp \left({\left({\sigma _{{\rm{VH}}}^0 - {B_{{\rm{VH}}}}\left(\theta \right)\ln \left({{Z_{\rm{S}}}} \right) - {C_{{\rm{VH}}}}\left(\theta \right)} \right)/{A_{{\rm{VH}}}}\left(\theta \right)} \right)}\end{array}} \right.$

(6)

    (3.2) Landsat 8土壤水分反演

刘英等人(2013)利用TM/ETM+影像不同波段的地表反射率数据，建立了不同的2维光谱特征空间，并提出了土壤水分监测指数SMMI(Soil Moisture Monitoring index)，且基于SWIR(band 7)–NIR(band 4)特征空间的SMMI(7, 4)在监测地表0—5 cm深度的土壤水分是最优的(R²=0.5410)。

为了验证模型(6)的可靠性，本文借助Landsat 8 OLI地表反射率数据利用光谱特征空间法，同时对研究区进行土壤了水分反演。为了保证与Sentinel-1A数据的时间一致性，本文选取2016-02-10、2016-05-16、2016-09-21、2016-12-10共4期OLI地表反射率影像进行处理。通过查询天气资料可知，2016-02-10、2016-05-16两天均无降水，2016-12-06和2016-12-10前一周也无降水，2016-09-21和2016-09-25两天虽然都有小雨，但可以视为降水一致，因此对比数据的时间差对于两种数据反演结果的对比造成的影响并不大。

因为Landsat 8 OLI影像的第5、7波段分别对应Landsat 5 TM影像、Landsat 7 ETM+影像的第4、7波段，因此，可构建土壤湿度指数SMMI(7, 5)。SMMI计算公式为(刘英 等，2013)

${\rm{SMMI}}\left({i, j} \right) = \left| {OE} \right|{\rm{/}}\left| {OD} \right| = \sqrt {r_i^2 + r_j^2} /\sqrt 2 $

(7)

式中， ${r_i}$ 和 ${r_j}$ 分别为影像的第 $i$ 和第 $j$ 波段地表反射率。在2维光谱特征空间中， $\left| {OE} \right|$ 长度的变化可一定程度上反映土壤水分的变化(图5)。

根据式(7)计算的4期SMMI(7, 5)影像(图6)。刘英等人(2013)利用SMMI影像点和实地对应点的同步实测0—5 cm深土壤水分(SM)数据构建了SMMI-SM散点图，并拟合出了土壤水分反演模型。本文利用得到的SMMI影像基于此模型对研究区土壤水分做反演。

4、结果与分析     (4.1) Sentinel-1A与Landsat 8水分反演结果分析

利用Sentinel-1A雷达数据与Landsat 8地表反射率数据反演的研究区土壤水分空间分布如图7、图8。图7中两种极化数据反演结果相差较大，尤其是在研究区东部的黄土丘陵沟壑区，认为主要是由于两种极化方式对于地形的敏感程度不一样导致。VV极化对于地形更加敏感，从而导致反演结果波动性更大。

研究区内主要包括沙漠裸地、荒漠化草原和黄土丘陵3种地表覆盖类型。分别在每种地表类型内随机采取样点，剔除水分体积含量大于0.5的异常值点以及距离近的相关性点，最终每种地表类型保留30个样点。样点用于分析微波与光学反演结果的相关性，验证模型(6)的可靠性。然后，结合DInSAR结果，在沉陷区及相应对比区采样。沉陷区与对比区采取等间距网格采样，在每个沉陷区及对比区内分别采取9个样本点。对比区的选取主要根据以下几个准则：(1)不能与沉陷区相距太远；(2)地表类型与沉陷区相似；(3)不能在采区工作面推进方向上；(4)云量条件与沉陷区相似。采样点用于分析沉陷区内土壤水分是否受开采沉陷影响。最终的样点分布如图9所示。

对采样点水分值分析可知Landsat 8数据反演的沙漠裸地区的土壤水分大部分为负值，说明依据刘英等人给出的模型反演的土壤水分含量值在沙漠裸地区是偏低的，需要进一步改进和优化。因此，本文只采用荒漠化草原区和黄土丘陵沟壑区的采样数据，对两种方法的结果做对比(图10)。

统计分析60个样点值可知，4期VH极化与OLI数据水分反演结果绝对误差在2%以内的点所占比例分别为48%、80%、80%和50%，相关系数最大为0.20。4期VV极化与OLI数据水分反演结果绝对误差在2%以内的点所占比例分别为20%、32%、22%和17%，相关系数最大为–0.03，表明两种结果基本不相关。此结果说明相对于VV极化，VH极化更适用于土壤水分反演，这与赵昕等人(2016)的研究结果一致。另外，从图10中可以看出，荒漠化草原区2016-02-10与2016-12-06两期VH极化的反演值偏低，认为主要原因是两种模型反映的均是地表浅层土壤水分含量，研究区冬季气温较低时会形成冻土层，对于微波的穿透性影响较大，因此造成土壤水分含量反演值偏低。

进一步对VH、VV极化和OLI数据的荒漠化草原区和黄土丘陵沟壑区反演结果进行分析(图11)。对VH极化与OLI数据的4期反演结果进行比较，荒漠化草原区土壤水分绝对误差在2%以内的点分别占80%、97%、100%和77%，相关系数最大为0.52，黄土丘陵沟壑区的土壤水分绝对误差在2%以内的点分别占17%、63%、60%和23%，相关系数最大为0.19。对VV极化与OLI数据的4期反演结果进行比较，荒漠化草原区土壤水分绝对误差在2%以内的点分别占20%、47%、30%和17%，相关系数最大为0.06，黄土丘陵沟壑区的土壤水分绝对误差在2%以内的点分别占20%、17%、13%和17%，相关系数最大为–0.05。此结果说明模型(6)更适用于平坦地区的土壤水分反演，地形对于后向散射的影响不可忽略(Pasolli 等，2011)。

    (4.2) 沉陷区与对比区土壤水分分析

对VH、VV极化反演的土壤水分结果进行沉陷区与对比区的采样分析(一个样区取9个样点值)，分别去除沉陷区与对比区内9个采样点中水分含量大于0.5的异常值点以及剩余点中的最大值和最小值点然后取平均，得到每个开采沉陷区以及相应对比区2016年12期土壤水分面域平均值，然后做对比(图12)。

统计每个开采沉陷区与相应对比区2016年12期的土壤水分含量大小(表4)。

由表4可知布尔台矿区地表土壤水分受开采沉陷的负面影响最大，呼和乌素矿受开采沉陷影响不大，而上湾矿区内沉陷区土壤水分含量则偏大，这可能是由于不同矿区的地下水埋深造成，需进一步验证。从总体而言，VH极化反演结果对比区土壤水分含量大于沉陷区的有41期，所占比例为57%；VV极化反演结果对比区土壤水分含量大于沉陷区的有36期，所占比例为50%，说明开采沉陷会导致地表土壤水分含量降低，造成沉陷区内土壤干化，但不同矿区之间的这种现象也会有差异。

5、结　论

土壤水分是矿区生态环境评价的一个重要参考标准，也是制约矿区可持续发展的一个主要问题。本文联合Sentinel-1A双极化雷达数据和Landsat 8光学数据反演了内蒙—陕西交界处的干旱—半干旱区的土壤表层水分，并对微波和光学两种反演结果做了对比分析；然后基于Sentinel-1A的DInSAR结果重点分析了开采沉陷对于地表土壤浅层水分含量的影响，得到的主要结论有：

(1)将VH和VV极化反演结果分别与基于OLI的土壤水分监测指数的土壤水分含量反演结果比较，发现VH极化更适用于土壤水分反演。4组对比数据中，反演结果绝对误差在2%以内的点所占比例最高为80%。另外，地形对后向散射的影响不可忽略，荒漠化草原区反演精度更好。

(2)对神东矿区内6个开采沉陷区以及对比区进行采样对比，布尔台矿采样区地表土壤水分受开采沉陷的影响最大，呼和乌素矿采样区受开采沉陷影响不大，而上湾矿采样区土壤水分含量则偏大。在2016年内72期数据中，VH极化反演结果对比区土壤水分含量大于沉陷区的有41期，所占比例为57%；VV极化反演结果对比区土壤水分含量大于沉陷区的有36期，所占比例为50%。说明对于整体而言，开采沉陷通过增加地表粗糙度影响地表土壤水分含量，但也有部分下沉盆地土壤水分含量偏大，可能与当地地下水埋深等因素有关，需进一步做验证。