CN110781602B

CN110781602B - 基于特征空间法来获取时空连续的土壤水方法

Info

Publication number: CN110781602B
Application number: CN201911066438.9A
Authority: CN
Inventors: 朱文彬
Original assignee: Institute of Geographic Sciences and Natural Resources of CAS
Current assignee: Institute of Geographic Sciences and Natural Resources of CAS
Priority date: 2019-11-04
Filing date: 2019-11-04
Publication date: 2021-06-15
Anticipated expiration: 2039-11-04
Also published as: CN110781602A

Abstract

本发明提供一种基于特征空间法来获取时空连续的土壤水方法，获取研究区卫星遥感数据和站点实测数据，对数据进行预处理，保证两种数据时间一致；根据太阳辐射与地表温度之间的余弦函数关系，分别构建干、湿边界函数，并以站点实测数据为驱动对干边进行参数优化。为减少优化方法对实测站点数量的依赖，仅采用一个站点进行率定，由多个站点参与验证，得到的精度与传统方法得到的精度相当。本发明实现了土壤水在时间上的连续估算。此外本发明通过使用卫星遥感数据进行土壤水估算，快速地获取陆地表面土壤水的分布信息，实现了土壤水在空间上的连续估算。

Description

基于特征空间法来获取时空连续的土壤水方法

技术领域

本发明属于土壤水测量方法，具体涉及一种基于特征空间法来获取时空连续的土壤水方法。

背景技术

土壤水是陆地和大气能量交换过程中的重要变量，并对陆表蒸散、水的运移和碳循环有控制作用，也是大气、水文、生态、农业等领域衡量土壤干旱程度的重要指标。土壤水的获取方法通常分为两大类：传统观测方法和卫星遥感估算方法。土壤水的传统观测方法是指从地面仪器直接测量，精度较高，但传统技术由于站点密度过于稀疏，无法有效反映其时空分布趋势。卫星遥感技术的发展为快速获取陆地表面土壤水的分布式信息提供了丰富的数据源。具体来说，基于微波遥感的土壤水分观测技术日益成熟，其能够穿透云层，具有“全天候”监测能力。然而，土壤水的微波遥感产品空间分辨率相对较低，其反演精度仍有待提高。相比之下，常用的光学/热红外遥感数据具有较高的空间分辨率，但是其容易受到云的影响，对于获取时空连续的土壤水信息具有一定困难，是当前研究面临的普遍性难题。在此背景下如何实时准确地获取时空连续的土壤水信息是非常紧迫而重要的科学问题。本专利基于中分辨率成像光谱仪获取的卫星遥感数据和大气辐射测量设备获取的实测土壤水数据，提出了一种获取时空连续的土壤水方法，并应用到大型异质区域进行了评价。

地表温度-植被指数特征空间法(简称特征空间法)是土壤湿度遥感监测的主流方法之一。其基本假设是，在均匀的大气条件下，如果研究区有足够多的像元，则由地表温度(T_s)和植被指数(VI)组成的二维空间会形成三角形或梯形形状。三角形/梯形形状的干湿边界表示像元在极大水分胁迫条件下可以达到的最大和最小T_s。因此，使用线性插值法根据每个像元在要素空间中的位置可确定其位置。显然，特征空间方法的准确性很大程度上取决于其干湿边界的正确选取。Sandholt等提出的温度植被干旱指数(TVDI)方法，其干边的确定主要分为经验法和理论法。经验法的检索过程存在极大的主观性和不确定性；而基于地表能量平衡的理论方法足以克服经验法的不足之处。但是理论法需要输入大量复杂参数，这与特征空间法的初衷相矛盾。

发明内容

针对上述技术问题，本发明提出一种改进的温度植被干旱指数(MTVDI)方法，基于裸土干点与湿点的物理属性，根据太阳直接辐射与太阳高度角之间的关系，以晴天条件下单个站点的土壤湿度观测数据为驱动，在年尺度上通过目标函数和约束条件的设定，以参数优化的方法构建干湿边界与太阳高度角的经验函数，率定干点参数方程中的振幅参数，最终将站点获得的最优振幅参数应用到整个研究区土壤水的求解中，以实现对土壤水的时空连续监测。

具体的技术方案为：

基于特征空间法来获取时空连续的土壤水方法，包括以下步骤：

(1)获取研究区卫星遥感数据(太阳天顶角、大气剖面数据、地表温度、植被指数)和站点实测数据(实测的体表土壤含水量)；

(2)对数据进行预处理，所述的预处理主要包括对步骤(1)获取的遥感数据进行数据子集的提取、投影变换、拼接、裁剪以及有效值转换，对实测数据进行数据子集提取，保证两种数据时间一致；

(3)利用研究区的植被归一化指数(NDVI)，进行植被覆盖度(f_c,i)计算；

(4)地表辐射温度(T_s,i)随着植被覆盖度的变化可以用土壤水分等值线上的裸土温度(T_smax,i)与植被温度(T_cmax,i)的加权和表示为公式(1)；此外，随着植被覆盖度的增加，整个像元(i)表面的辐射温度与植被冠层内的空气温度(T_a,i)保持平衡，因此T_smax,i可以推导为公式(2)：

T_s,i＝f_c,i·T_cmax,i+(1-f_c,i)·T_smax,i (1)

(5)假设地表面的情况稳定而均匀，T_smax的季节性变化以接收到的太阳辐射为主，T_smax的变化通过太阳天顶角(θ)来捕获，T_smax算法表达式：

T_smax,j＝A₁·cos(a₁·θ_j+b₁) (3)

式中，j表示第j天观察到所有变量；T_smax,j由晴天条件下，从地表温度数据中检索的日最高地表温度(T_max,j)和定位数据中检索的日平均θ_j来表示，并率定出表达式中的参数a₁和b₁，而振幅参数A₁利用优化算法确定；

(6)振幅参数A₁的优化过程主要由目标函数和约束函数组成:

目标函数：

式中，n代表晴天数量，SM_i,j和MTVDI_i,j分别代表站点i在第j天的实测体表土壤含水量和反演的MTVDI。

约束函数1：针对研究区的所有像元：

0≤MTVDI≤1 (5)

约束函数2：T_smax,j取值范围的限定：

T_max,j≤T_smax,j≤U (6)

其中，U是一个经验值，表示T_smax,j的上限；

(7)相对于干边的确定过程而言，湿边的确定方式更为简单，T_w,j被定义为研究区最低地表温度(T_min,j)，其计算表达式表示为：

T_w,j＝min(T_smin,j,T_min,j) (7)

T_smin,j＝A₂·cos(a₂·θ_j+b₂) (8)

式中，T_smin,j是T_min,j的余弦函数估计值，T_smin,j表达式中的三个参数(A₂，a₂，b₂)通过在晴天条件下，由地表温度数据中检索的日最低地表温度(T_min,j)和定位数据中检索的日平均θ_j来率定得到；

(8)结合步骤(4)、步骤(5)、步骤(7)计算的变量，反演MTVDI，计算公式表示为：

式中，θ_v,i表示像元i的体表土壤含水量，θ_max,i和θ_min,i分别为θ_v,i的上下限；MTVDI与θ_v呈负相关关系；

(9)通过步骤(8)获得MTVDI值，结合前人的研究成果，输入实测的体表土壤含水量与MTVDI，采用最小二乘回归方法拟合出参数θ_max,i和θ_min,i，最后进行土壤水含量的时空连续估算。

与现有技术比较，本发明的有益效果为：

本发明提出了基于特征空间法来获取时空连续的土壤水方法，并将其应用到大型异质区域土壤水的估算，克服了传统特征空间方法在土壤水分反演过程中存在的复杂性和不连续性。本发明根据太阳辐射与地表温度之间的余弦函数关系，分别构建干、湿边界函数，并以站点实测数据为驱动对干边进行参数优化。为减少优化方法对实测站点数量的依赖，仅采用一个站点进行率定，由多个站点参与验证，得到的精度与传统方法得到的精度相当。当研究区大部分区域处于晴天条件下的率定结果对部分有云情况下也同样适用，尽管精度略低于晴天条件下的验证结果，但是相对于传统特征空间法来说是一种突破，由此可见，本发明实现了土壤水在时间上的连续估算。此外，传统技术获取土壤水由于站点空间代表性差，无法获取空间连续的土壤水分布情况，本发明通过使用卫星遥感数据进行土壤水估算，快速地获取陆地表面土壤水的分布信息，实现了土壤水在空间上的连续估算。

附图说明

图1本发明的流程图；

图2地表温度-植被覆盖度特征空间图；

图3(a)E1站点在晴天条件下的率定结果；

图3(b)E2站点在晴天条件下的率定结果；

图3(c)E6站点在晴天条件下的率定结果；

图3(d)E7站点在晴天条件下的率定结果；

图3(e)E9站点在晴天条件下的率定结果；

图3(f)E10站点在晴天条件下的率定结果；

图3(g)E11站点在晴天条件下的率定结果；

图3(h)E12站点在晴天条件下的率定结果；

图3(i)E13站点在晴天条件下的率定结果；

图3(j)E20站点在晴天条件下的率定结果；

图3(k)E22站点在晴天条件下的率定结果；

图4(a)E1站点在晴天条件下的验证结果；

图4(b)E2站点在晴天条件下的验证结果；

图4(c)E6站点在晴天条件下的验证结果；

图4(d)E7站点在晴天条件下的验证结果；

图4(e)E9站点在晴天条件下的验证结果；

图4(f)E10站点在晴天条件下的验证结果；

图4(g)E11站点在晴天条件下的验证结果；

图4(h)E12站点在晴天条件下的验证结果；

图4(i)E13站点在晴天条件下的验证结果；

图4(j)E20站点在晴天条件下的验证结果；

图4(k)E22站点在晴天条件下的验证结果；

图5(a)E1站点在部分有云条件下的验证结果；

图5(b)E2站点在部分有云条件下的验证结果；

图5(c)E6站点在部分有云条件下的验证结果；

图5(d)E7站点在部分有云条件下的验证结果；

图5(e)E9站点在部分有云条件下的验证结果；

图5(f)E10站点在部分有云条件下的验证结果；

图5(g)E11站点在部分有云条件下的验证结果；

图5(h)E12站点在部分有云条件下的验证结果；

图5(i)E13站点在部分有云条件下的验证结果；

图5(j)E20站点在部分有云条件下的验证结果；

图5(k)E22站点在部分有云条件下的验证结果；

图6(a)2004年E1站点土壤湿度估计值与实测值之间的趋势图；

图6(b)2004年E2站点土壤湿度估计值与实测值之间的趋势图；

图6(d)2004年E7站点土壤湿度估计值与实测值之间的趋势图；

图6(e)2004年E9站点土壤湿度估计值与实测值之间的趋势图；

图6(f)2004年E10站点土壤湿度估计值与实测值之间的趋势图；

图6(g)2004年E11站点土壤湿度估计值与实测值之间的趋势图；

图6(h)2004年E12站点土壤湿度估计值与实测值之间的趋势图；

图6(i)2004年E13站点土壤湿度估计值与实测值之间的趋势图；

图6(j)2004年E20站点土壤湿度估计值与实测值之间的趋势图；

图6(k)2004年E22站点土壤湿度估计值与实测值之间的趋势图。

具体实施方式

结合实施例说明本发明的具体技术方案。

以2004年美国南部大平原地区为研究对象，对土壤水观测的时空拓展技术进行了验证，按照图1的流程，具体验证过程如下：

分别下载美国南部大平原2004年的MODIS数据和站点实测土壤水数据，对数据进行预处理，得到用于构建特征空间的植被盖度与地表温度变量和构建MTVDI模型的变量，特征空间的框架可见图2。

选取研究区11个站点进行实验。首先分别输入每个站点晴天条件下的数据进行干边的参数优化，每个站点得到一个最优干边，代入公式(9)计算得到MTVDI。图3(a)～图3(k)分别表示11个站点的率定结果精度，由MTVDI与实测土壤水之间的散点图即图3(a)～图3(k)可知，二者在大部分站点都呈现出了明显的负相关趋势，相关系数的平均值为-0.59。然后将率定结果应用到其他站点进行验证，如图4(a)～图4(k)所示。同样可以观察到，大部分站点呈现出显著的负相关关系，并且相关系数的平均值为-0.58，与率定站点结果精度的差异很小，几乎为0。由此可见，一个站点的率定结果对其他站点同样适用，并且精度相当，实现了土壤水的空间连续估算。

将各个站点在晴天条件下的率定结果应用到其对应站点在部分有云的天气条件下，结果如图5(a)～图5(k)所示。MTVDI与实测土壤水之间的相关性处于-0.35至-0.84之间，平均值为-0.56，同样与晴天的率定结果之间差异很小。因此，可以看出晴天条件下的率定结果对部分有云情况下也同样适用，实现了土壤水的时间连续估算。最后根据最小二乘回归方法，通过输入实测的体积土壤含水量与MTVDI，拟合出公式(9)中的参数θ_max,i和θ_min,i，进行土壤含水量的连续估算。

如图6(a)～图6(k)所示，展示了研究区内11个站点土壤水的估计值与实测值的趋势图，可以看出二者在大部分站点的吻合度较高。下表1列出了每个站点土壤水的估计精度。为比较起见，引入了其他研究方法的结果进行对比，以验证本发明所提出方法的准确性。总体而言，尽管所有三个统计数据都显示本发明验证结果产生的精度略高一些，但是这两种方法的精度没有显著差异。

表1每个站点土壤水的估计精度对比

在所有站点中，本发明产生的r、MAE和RMSE分别为0.78、0.017m³m^-3和0.022m³m^-3，而其他研究产生的r、MAE和RMSE分别为0.075、0.013m³m^-3和0.025m³m^-3。因此，本发明提出的土壤水时空连续估算精度足以达到或超过前人的研究精度。

Claims

1.基于特征空间法来获取时空连续的土壤水方法，其特征在于，包括以下步骤：

(1)获取研究区卫星遥感数据和站点实测数据，站点实测数据即实测的体表土壤含水量；

(2)对数据进行预处理，保证卫星遥感数据和站点实测数据两种数据时间一致；

(3)利用研究区的植被归一化指数NDVI，进行植被覆盖度f_c，i计算；

(4)地表辐射温度T_s，i随着植被覆盖度的变化可以用土壤水分等值线上的裸土温度T_smax，i与植被温度T_cmax，i的加权和表示为公式(1)；此外，随着植被覆盖度的增加，整个像元i表面的辐射温度与植被冠层内的空气温度T_a，i保持平衡，T_smax，i推导为公式(2)：

T_s，i＝f_c，i·T_cmax，i+(1-f_c，i)·T_smax，i (1)

(5)假设地表面的情况稳定而均匀，T_smax，j的季节性变化以接收到的太阳辐射为主，T_smax，j的变化通过太阳天顶角θ_j来捕获，T_smax，j算法表达式：

T_smax，j＝A₁·cos(a₁·θ_j+b₁) (3)

式中，j表示第j天观察到所有变量；T_smax，j由晴天条件下，从地表温度数据中检索的日最高地表温度T_max，j和定位数据中检索的日平均θ_j来表示，并率定出表达式中的参数a₁和b₁，而振幅参数A₁利用优化算法确定；

(6)振幅参数A₁的优化过程由目标函数和约束函数组成:

目标函数：

式中，n代表晴天数量，SM_i，j和MTVDI_i，j分别代表站点i在第j天的实测体表土壤含水量和反演的MTVDI；

约束函数1：针对研究区的所有像元：

0≤MTVDI≤1 (5)

约束函数2：T_smax，j取值范围的限定：

T_max，j≤T_smax，j≤U (6)

其中，U是一个经验值，表示T_smax，j的上限；

(7)T_w，j被定义为研究区最低地表温度T_min，j，其计算表达式表示为：

T_w，j＝min(T_smin，j，T_min，j ) (7)

T_smin，j＝A₂·cos(a₂·θ_j+b₂) (8)

式中，T_smin，j是T_min，j的余弦函数估计值，T_smin，j表达式中的三个参数即A₂，a₂，b₂，通过在晴天条件下，由地表温度数据中检索的日最低地表温度T_min，j和定位数据中检索的日平均θ_j来率定得到；

(8)结合步骤(4)、步骤(5)、步骤(7)计算的变量，反演MTVDI；计算公式表示为：

式中，θ_v，i表示像元i的体表土壤含水量，θ_max，i和θ_min，i分别为θ_v，i的上下限；MTVDI与θ_v，i呈负相关关系；

(9)通过步骤(8)获得MTVDI值，输入实测的体表土壤含水量与MTVDI，采用最小二乘回归方法拟合出参数θ_max，i和θ_min，i，最后进行土壤水含量的时空连续估算。

2.根据权利要求1所述的基于特征空间法来获取时空连续的土壤水方法，其特征在于，步骤(1)中所述的研究区卫星遥感数据，包括太阳天顶角、大气剖面数据、地表温度、植被指数。

3.根据权利要求1所述的基于特征空间法来获取时空连续的土壤水方法，其特征在于，步骤(2)中所述的预处理包括，对步骤(1)获取的遥感数据进行数据子集的提取、投影变换、拼接、裁剪以及有效值转换，对实测数据进行数据子集提取，保证两种数据时间一致。