基于区域特征相似度的微波土壤水分反演结果可信度评价

2018-03-16陈鲁皖秦小宝王文娟

水资源与水工程学报 2018年1期

陈鲁皖，韩玲，秦小宝，王文娟

(长安大学地质工程与测绘学院，陕西西安 710064)

1 研究背景

土壤水分是农业、林业、气象学、水文学和生态学等领域研究的重要参数[1]，利用微波遥感反演土壤水分是一种很有效的手段。利用SAR反演土壤水分的原理是土壤表面的后向散射系数与土壤含水量、地表粗糙度、地表物质介电特性、土壤物理特性、植被特性以及雷达系统参数密切相关[2]。目前通常利用野外地表实测数据结合同期SAR数据来拟合或训练土壤水分反演经验方程(关系)，以此得到更高精度的土壤水分反演结果[3-5]。反演结果的精度评价一般是基于实测数据的，即选取部分采样点，将采样点的土壤水分实测数据和反演数据进行比较，通常选择均方根误差(RMSE)和决定系数(R2)作为精度指标[6-9]。

野外实测布设样区时，受多种因素的影响(卫星过顶时间、测量人员数目和测量经费等)，往往样方的设置范围和采样点的数量都比较有限。根据主动微波探测土壤水分的物理机制，影响地表后向散射系数的各因子在不同地区数值不同，即反演区域内存在各影响因子的异质性。所以有限范围的样区中实测数据拟合或训练得到的土壤水分反演经验方程(关系)，在较大范围的区域进行反演，不同地区的反演结果精度必然不同。目前常用的基于实测数据的土壤水分反演精度评价方法实际上只能反映实测数据所在样区的反演精度，不能反映所有反演区域的反演精度。

本文提出了一种基于区域特征相似度的微波土壤水分反演结果可信度评价方法，该方法可以有效反映所有反演区域的土壤水分反演精度。这里的可信度，是反演区域内土壤水分反演结果的可信程度，表示了反演结果与土壤水分真实值之间的接近程度。本文方法的原理是，反演区域中不同地区与样区在土壤湿度分布、地表粗糙度、地表物质介电特性、土壤物理特性、植被特性和雷达系统参数等方面具有不同的相似性，而土壤水分反演经验方程(关系)又是基于样区的实测数据得到的，因此可以通过不同反演地区与样区的相似程度来反映土壤水分反演结果的可信度。通过分析TM影像、DEM数据、HWSD土壤质地数据集和EnvisatASAR数据，结合实际地表特征，选择土壤湿度、地表温度、NDVI、STI(SoilTextureIndex，土壤质地指数)、地形指数、雷达入射角和LandsatTM影像的4个波段(b3、b4、b5、b7)共10个影响因子，通过主成分分析得到前3个主成分，将前3个主成分合成为RGB影像，再基于该RGB影像使用分水岭算法进行分割，构建各分割区域的特征向量，计算分割后各区域与样本区域的马氏距离，基于特征相似度数据集计算反映土壤水分反演结果的可信度。

2 试验数据

本文所使用的试验数据来源于中国科学院寒区旱区科学数据中心所提供的“黑河综合遥感联合试验”。

2.1 SAR数据来源

本文选用欧空局的Envisat-1卫星上ASAR传感器获取的SAR影像作为土壤含水量反演的数据源。采用了2008年7月11日的ASAR数据，该影像的入射波段为C波段(f=5.331GHZ)，入射角范围31.0°～36.3°，经度范围99°28′E～100°43′E，纬度范围38°42′N～39°48′N，地面分辨率为12.5m×12.5m，工作模式为AlternatingPolarization，极化方式为VV和VH两种。

2.2 地面数据来源

为验证本文所提出的方法的有效性，特地选取了具有多个采样区的甘肃省黑河中游临泽地区作为研究区。该地区有4个样区，分别为360m×360m的样方B、C、D、E，样点间距为60m，B和C样方地表类型是带稀疏杂草的盐碱地，D样方地表类型是苜蓿，E样方地表类型是大麦地。与ASAR影像同期的共187组地面样方观测数据，观测数据包括地表土壤含水量、均方根高度和相关长度。为了对比验证区域特征相似度是否能够有效反映反演结果的可信度，分别基于B样区和E样区中的部分实测数据构建土壤水分反演经验方程，其余的B、C、D、E样区中的实测数据为验证数据。

2.3 其他数据来源

本文选用了Landsat5卫星的TM影像，空间分辨率为30m×30m；还选用了中国科学院寒区旱区科学数据中心提供的“黑河流域HWSD土壤质地数据集”和“黑河流域ASTERGDEM数据集”。

3 研究方法

土壤水分反演结果可信度评价方法的流程如图1所示。

图1 研究方法流程图

具体步骤：(1)利用多元数据获取10个与地表土壤后向散射系数密切相关的影响因子；(2)将各影响因子的影像进行波段合成；(3)利用主成分分析法[10](PCA)对步骤(2)的结果图像处理提取出主成分，将前3个主成分合成为RGB影像；(4)利用分水岭算法[11]对步骤(3)的结果影像做分割，获得过分割影像；(5)构建各分割区域和样区的特征向量，分别计算不同分割区域与各样本区域间的特征相似度；(6)对步骤(5)的结果进行处理，得到反演结果的可信度。

3.1 获取影响因子

基于SAR反演土壤水分的物理机制，本文选择了土壤湿度、地表温度、NDVI、STI、地形指数、雷达入射角这6种影响因子。

TM影像4个波段(b3、b4、b5、b7)中，通过构建各波段之间二维光谱特征空间[12]，发现TM红外、近红外、短波红外波段的光谱特征与土壤水分存在复杂而紧密的联系，因而可以用来监测地表土壤湿度。

NDVI、地表温度和土壤湿度可由TM反演计算得到；STI为本文提出的一种反映土壤质地的新指数，由“黑河流域HWSD土壤质地数据集” 中的土壤砂土、黏土、淤泥含量计算得到；地形指数由“黑河流域ASTERGDEM数据集”通过提取地形的坡度计算得到；雷达入射角由辐射定标后的SAR影像中得到。3.1.1 地表温度和土壤湿度利用TM影像和单窗算法可以反演地表温度[13-14]。统计分析后发现土壤水分含量与其温度之间存在相关性较好的线性关系，利用线性方程(1)可直接反演土壤湿度的分布。

Mv=cT+d

(1)

式中：Mv为土壤水分，cm3/cm3； T为地表土壤温度，℃； c、d为方程的经验系数。

用通过反演得到的地表土壤温度来估算土壤湿度，可以反映出土壤湿度的分布规律。

3.1.2 土壤质地指数(SoilTextureIndex,STI) 土壤质地是根据土壤的颗粒组成划分的土壤类型，通常使用砂土含量、黏土含量和淤泥含量3种指标来量化土壤质地。为综合考虑这3种指标，减少参与后续主成分分析的因子数量，本文提出了土壤质地指数STI。

本文使用AIEM模型建立基于不同土壤砂土含量sv和黏土含量cv的后向散射系数模拟数据集，其中sv和cv的取值参照“黑河流域HWSD土壤质地数据集”中提供的黑河地区的sv和cv的数值。AIEM模型的地表参数和雷达系统参数为固定值，做如下设定：雷达入射角给定为30.0°，入射频率5.331GHz，土壤水分Mv=30%，均方根高度S=1.2cm，相关长度l=12cm(粗糙度参数通过统计分析“黑河综合遥感联合试验”、美国NASA与农业部的合作试验“Washita92”可得)。

(2)

相关系数R2=0.999819，整理式(2)可得：

(3)

若设：

STI=3sv+cv

(4)

则可得后向散射系数与STI的线性方程：

(5)

式中：a、b为方程的经验系数。

3.1.3 地形指数地形指数是指反演区域地球表面积与投影面积之比，可通过坡度的余弦得到[15]。

Z=1/cos(S)

(6)

式中： Z为地形指数； S为坡度，rad。

3.1.4 雷达入射角对分布目标而言，雷达回波在近距离较强，随着向远距离移动将逐渐减弱。雷达入射角可对微波遥感的反演结果产生较大影响[16]。而本文选取的EnvisatASAR影像，沿单幅图像距离向入射角大约有5°的变化。因此，雷达入射角也是影响土壤后向散射系数的一个因素[2]。

本文通过使用NEST(NestESASARToolbox)软件处理EnvisatASAR数据，可得到雷达入射角影像。

3.2 基于主成分分析的分水岭算法影像分割

3.2.1 主成分分析主成分分析法又称K-L(Karhunen-Loeve)变换，可实现在尽可能不丢失信息的同时，用几个综合性分量代表多波段的原图像[10]。一般情况下，前3个主成分(PC1、PC2、PC3)包含所有波段中90%以上的方差信息。

在对10个后向散射系数影响因子的图层进行波段合成后，利用ENVI5.1软件中的PrincipalComponents模块进行主成分分析。

3.2.2 基于分水岭算法的影像分割分水岭算法是一种建立在数学形态学基础之上算法。传统分水岭算法只能用于灰度图像的分割，无法分割彩色图像。本文把主成分分析后得到的前3个主成分合成为RGB彩色影像，选择基于形态学梯度的分水岭彩色图像分割方法，首先计算彩色图像的梯度图像，然后结合多尺度形态学和巴特沃斯低通滤波对彩色梯度进行修正得到新的彩色梯度图，最后对修正后的梯度图像进行分水岭变换得到分割结果[11-17]。

3.3 计算区域特征相似度

在分割后得到的区域和样区之间进行特征相似度的判断，使用数值来有效地表示各分割区域和样区间地表情况的相似程度，这便是区域特征相似度的计算问题。特征相似度大多可以表示成向量的形式，本文构建特征向量时选用了参与主成分分析的6个因子(土壤湿度、地表温度、NDVI、STI、地形指数、雷达入射角)作为特征向量的分量。

通常采用距离法来表示特征的相似度量，常用的有欧氏距离、马氏距离[18-19]。对后向散射系数有影响的因素中，土壤湿度、地形指数和NDVI值这几个因子相对于其他因子对后向散射系数的影响更大[20]。马氏距离主要用于特征向量的各分量具有相关性或各分量的权重不等的情况。

因此本文采用马氏距离来计算区域特征相似度，马氏距离[19]表示为：

(7)

式中：x和y为两个准备计算的区域特征向量;E为影像的特征向量协方差矩阵;d(x,y)为向量x和y之间的马氏距离。

3.4 计算反演结果可信度

反演结果可信度基于区域特征相似度进行计算，反映了土壤水分反演结果的可信程度，可通过下式计算：

(8)

式中：Re为反演结果可信度;d为马氏距离，min为1/d中的最小值，max为1/d中的最大值。

3.5 土壤水分反演

根据文献[3-5]，给定入射角时，裸土地表的后向散射系数与地表粗糙度、土壤含水量之间的关系可以表示为[4]：

(9)

Zs=S3/l2

(10)

4 试验对比

4.1 获取影响因子的试验

本文使用Landsat5的TM影像，反演得到NDVI、地表温度和土壤湿度影像。图2(a)和2(b)分别是反演区域的地表温度和土壤湿度影像。图2(a)中红色区域地表温度最高，其土地覆盖类型为裸土；其次是地表温度较高的黄色区域，其土地覆盖类型为稀疏草地；再次是温度较低的绿色区域，土地覆盖类型为小麦地、大麦地、玉米地和湿地。土壤湿度图也有类似的分布情况。

根据“黑河流域HWSD土壤质地数据集”得到的研究区土壤砂土含量、黏土含量，采用公式(4)计算土壤质地指数STI。地形指数由“黑河流域ASTERGDEM数据集”得到；入射角由辐射定标后的SAR影像中得到。

4.2 基于主成分分析的分水岭分割试验

对10个影响因子进行主成分分析，得到表1。

从表1可以看出，前3个主成分的样本方差累计贡献率已经达到96.65%，把前3个主成分合成为RGB彩色图。使用分水岭算法对RGB影像进行分割，图3为分割后的结果图。从图3中可以看出，裸地、草地、农田和水体都得到了较好的同质分割。

表1 主成分因子贡献率

4.3 反演结果可信度计算

选择临泽地区的B样区，计算每个分割区域特征向量与B样区特征向量之间的马氏距离，然后根据公式(9)计算得到反演结果可信度。为验证反演结果可信度的有效性，对比基于不同样区的实测数据得到的反演结果，又选择E样区，同样计算每个分割区域特征向量与E样区特征向量之间的马氏距离和反演结果可信度。表2为部分分割区域与B、E样区之间的马氏距离和反演结果可信度，表2中各分割区域的地物类型参考“黑河生态水文遥感试验：黑河流域土地利用覆被数据集(2011年7月)”。

表2 部分分割区域与样本区域之间的特征相似度

B样方地表类型是盐碱草地，表2中与B样方相似程度较高的区域为F4、F43、F5与F11。E样方地表类型是大麦地，表2中与E样方较为相似的区域为F24、F40、F22、F48。

各分割区域中，相似度较好分割区域的土地覆盖类型基本与B样方一致，这从一定程度上证明了本文方法的有效性。

为评估本文提出方法的有效性，采用包括临泽地区4个样方共187个样点的地面观测数据参与验证试验。选择B样方中的部分地面观测数据作为训练数据结合同期EnvisatASAR观测影像，得到土壤水分反演经验方程。

将B样区训练数据包括土壤含水量、均方根高度和相关长度等数据，代入到公式(9)、(10)中，可得到公式(11)、(12)，两式联立消参求解可得反演区域的土壤含水量。

0.96lnZslnMv+5.31

(11)

0.22lnZslnMv-12.86

(12)

图4(a)为式(11)、(12)的土壤水分反演结果，图4(b)为反演结果可信度。

同理，选择E样方中的部分地面观测数据作为训练数据，可反演得到土壤水分。图5(a)为土壤水分反演结果，图5(b)为反演结果可信度。

4.4 对比分析

首先将B样方的其他土壤水分实测数据以及C、D、E样方的土壤水分实测数据作为验证数据，与基于样区B的土壤水分反演结果中各采样点的土壤水分反演数据进行对比；同样地对基于样区E的土壤水分反演结果进行对比，表3为对比的结果。

图2 部分试验影像及TM影像

图3 RGB合成图分割结果

训练样区验证样区R2分割区域特征相似度ReB0.7956F40.01321.0000BC0.7024F170.46130.3924D0.6121F221.06070.2939E0.4142F242.10760.1881B0.4211F42.10760.3694EC0.3875F172.55890.2812D0.6201F221.04830.4728E0.7442F240.02371.0000