基于双抛物线型NDVI‑Ts特征空间的土壤湿度遥感监测方法与流程

本发明涉及土壤检测技术领域，尤其涉及一种基于双抛物线型NDVI-T_s特征空间的土壤湿度遥感监测方法。

背景技术：

干旱是人类面临的主要自然灾害之一，而土壤湿度作为陆面水资源形成、转化、消耗过程的基本参数，是反映陆地干旱的重要指标。多时相、多光谱、高光谱遥感数据反映了大面积的地表信息，使得快速、及时、动态监测旱情成为可能。

国内外学者基于多源、多光谱、多时相遥感数据，提出了多种土壤湿度遥感监测的方法和模型。总体上，这些方法主要可分为可见光——红外法和微波遥感法。利用可见光和红外波段遥感监测土壤湿度的具体方法有距平植被指数法、条件植被指数法、条件温度指数法、归一化温度指数法、条件植被温度指数法、植被供水指数法、作物缺水指数法、热惯量和表观热惯量法、光谱特征空间法等。温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)是一种应用较为广泛的方法。Price和Carlson等发现，当研究区域的植被覆盖和土壤湿度变化范围较大时，根据遥感资料得到的以归一化植被指数(Normalized Difference Vegetation Index，NDVI)为横轴和以地表辐射温度为纵轴的散点图呈三角形特征；Moran等认为NDVI-T_s(地表温度)之间为梯形特征空间；Sandholt等提出了TVDI的计算方法；Patel等基于MODIS产品数据的NDVI-T_s特征空间计算TVDI，并估算了印度亚湿润地区的土壤湿度，结果表明TVDI与实测土壤水分之间存在强负相关关系；Schnur等采用250m分辨率MODIS产品NDVI和增强型植被指数(Enhanced Vegetation Index，EVI)数据评估美国西南区域的土壤水分，发现NDVI与土壤湿度的相关性略高于EVI和土壤湿度的相关性。

国内外相关研究大多采用1km或5km分辨率MODIS产品数据来构建NDVI-T_s特征空间，并认为此特征空间呈三角形或梯形，但是，采用NDVI-T_s三角形或梯形特征空间计算TVDI的原理为T_smax与NDVI呈线性关系，随着NDVI的增加，T_smax线性减小，且一般认为，当NDVI＜0.15时陆地表面是裸地，线性拟合时不予考虑。

技术实现要素：

针对上述缺陷或不足，本发明的目的在于提供一种基于双抛物线型NDVI-T_s特征空间的土壤湿度遥感监测方法。

本发明的技术方案为：

一种基于双抛物线型NDVI-T_s特征空间的土壤湿度遥感监测方法，包括：

1)获取研究地区遥感影像数据；

2)对遥感影像数据进行预处理，所述预处理包括对遥感影像数据进行辐射校正、大气校正、以及几何校正，得到预处理后的遥感影像图；

3)根据遥感影像图，获取研究区归一化植被指数(NDVI)与地表温度(T_s)数据，利用温度植被干旱指数(TVDI)计算公式获取每个遥感影像图中像元的TVDI值，TVDI值表达式为：

$TVDI=\frac{T_s-T_{smin}}{T_{smax}-T_{smin}}---\left(1\right)$

其中，T_s表示地表温度；T_smin表示相同NDVI值对应的最小地表温度，是NDVI-T_s特征空间中的湿边；T_smax为研究区相同NDVI值对应的最大地表温度，代表特征空间中的干边；

4)构建双抛物线型NDVI-T_s特征空间散点图，并根据双抛物线型NDVI-T_s特征空间获取TVDI中T_smax、T_smin算法表达式：

T_smax＝a₁×NDVI²+b₁×NDVI+c₁

T_smin＝a₂×NDVI²+b₂×NDVI+c₂ (2)

其中，a₁、b₁、c₁、a₂、b₂、c₂为方程拟合系数；

5)根据公式(1)和公式(2)得到研究区每个像元的TVDI值；

6)根据TVDI值的大小，依据标准将土壤湿度进行划分，得到待监测地区的土壤湿度分布状况。

所述步骤6)中根据获得的TVDI，结合前人研究成果，将土壤湿度划分为5类，具体为：

当0＜TVDI≤0.2时，则为极湿润；当0.2＜TVDI≤0.4时，则为湿润；当0.4＜TVDI≤0.6时，则为正常；当0.6＜TVDI≤0.8时，则为干旱；当0.8＜TVDI≤I时，则为极干旱。

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

本发明提供了一种基于双抛物线型NDVI-T_s特征空间的土壤湿度遥感监测方法，分析了基于1km、500m、250m分辨率MODIS数据的NDVI-T_s特征空间，发现NDVI-T_s特征空间均具有双抛物线型特征。与基于1km、500m、250m分辨率MODIS数据的NDVI-T_s三角形特征空间进行对比分析，并将得到的TVDI数据与实测土壤湿度进行相关性分析，结果表明：基于双抛物线型NDVI-T_s特征空间的TVDI_c在反映地表浅层土壤湿度方面要优于基于三角形NDVI-T_s特征空间的TVDI_t；且250m双抛物线型NDVI-T_s特征空间在监测地表0-5cm土壤湿度方面是最好的，而500m双抛物线型NDVI-T_s特征空间在监测10cm深土壤湿度时更有优势。

附图说明

图1是本发明基于双抛物线型NDVI-T_s特征空间的土壤湿度遥感监测方法流程图；

图2是TVDI在三角形NDVI-T_s特征空间中的定义图；

图3是TVDI在双抛物线型NDVI-T_s特征空间中的定义图；

图4是本发明实施案例中10cm深度土壤湿度与TVDI的关系图。

具体实施方式

下面结合附图对本发明做详细描述。

如图1所示，本发明提供了一种基于双抛物线型NDVI-T_s特征空间的土壤湿度遥感监测方法，包括：

1)获取研究地区遥感影像数据；

2)对遥感影像数据进行预处理，所述预处理包括对遥感影像数据进行辐射校正、大气校正、以及几何校正，得到预处理后的遥感影像图；

3)根据遥感影像图，获取研究区归一化植被指数(NDVI)与地表温度(T_s)数据，利用温度植被干旱指数(TVDI)计算公式获取每个遥感影像图中像元的TVDI值，TVDI值表达式为：

$TVDI=\frac{T_s-T_{smin}}{T_{smax}-T_{smin}}---\left(1\right)$

其中，T_s表示地表温度；T_smin表示相同NDVI值对应的最小地表温度，是NDVI-T_s特征空间中的湿边；T_smax为研究区相同NDVI值对应的最大地表温度，代表特征空间中的干边；TVDI值的取值范围是0到1。TVDI值越大，T_s越接近干边，越干旱；反之TVDI值越小，T_s越接近湿边，越湿润。

在现有技术中，三角形或梯形NDVI-T_s特征空间计算TVDI原理是：随着NDVI增加，T_smax呈线性减小趋势；当NDVI＜0.15时，认为陆地表面是裸地，无植被覆盖，在线性拟合T_smax时不予考虑。本发明发现，NDVI-T_s特征空间呈双抛物线型，随着NDVI增加，T_smax呈非线性减小趋势，如图2所示，与NDVI-T_s三角形或梯形特征空间不一致。但若将双抛物线型特征空间中的干湿边在拟合时不考虑NDVI＜0.15部分，则NDVI-T_s散点图呈三角形。

利用上述两种特征空间分别计算TVDI，经对比发现两种特征空间得到的结果相似，但与实测土壤湿度数据(Soil Moisture，SM)进行相关分析时发现，由双抛物线型特征空间计算的TVDI与0-5cm和10cm深SM相关性要优于由三角形特征空间计算的TVDI与SM的相关性。因而，认为NDVI＜0.15部分应被包含在双抛物线型NDVI-T_s特征空间中，不应舍弃。

本发明将NDVI扩展到0.15以内，构建双抛物线型NDVI-T_s特征空间如图3所示。

4)构建双抛物线型NDVI-T_s特征空间散点图，并根据双抛物线型NDVI-T_s特征空间获取TVDI中T_smax、T_smin算法表达式：

T_smax＝a₁×NDVI²+b₁×NDVI+c₁

T_smin＝a₂×NDVI²+b₂×NDVI+c₂ (2)

其中，a₁、b₁、c₁、a₂、b₂、c₂为方程拟合系数；

5)根据公式(1)和公式(2)得到研究区每个像元的TVDI值；

6)根据TVDI值的大小，依据标准将土壤湿度进行划分，得到待监测地区的土壤湿度分布状况，为土壤湿度的监测提供依据。

根据获得的TVDI，将土壤湿度进行划分，具体为：

将土壤湿度划分为5类：

当0＜TVDI≤0.2时，则为极湿润；当0.2＜TVDI≤0.4时，则为湿润；当0.4＜TVDI≤0.6时，则为正常；当0.6＜TVDI≤0.8时，则为干旱；当0.8＜TVDI≤1时，则为极干旱。

示例性的，本发明实施案例中，以神东矿区环境与灾害遥感监测与分析为对象，根据MODIS/AQUA卫星产品的归一化植被指数NDVI和地表温度T_s数据，分析了基于1km、500m、250m分辨率数据的NDVI-T_s特征空间，发现NDVI-T_s特征空间均具有双抛物线型特征。与基于1km、500m、250m分辨率数据的NDVI-T_s三角形特征空间进行对比分析，并将得到的温度植被干旱指数TVDI数据与实测土壤湿度进行相关性分析，结果表明：基于双抛物线型NDVI-T_s特征空间的TVDIc在反映地表浅层土壤湿度方面要优于基于三角形NDVI-T_s特征空间的TVDIt；且250m双抛物线型NDVI-T_s特征空间在监测地表0-5cm土壤湿度方面是最好的，而500m双抛物线型NDVI-T_s特征空间在监测10cm深土壤湿度时更有优势。

具体验证过为：

筛选出影像单元均一、面积约1km²的代表性样方30个，并在神东矿区进行与卫星同步的土壤样点采样，并根据编号分类保存，现场测定土壤湿重，精确记录。在各样地按0-5cm、10cm深度分别取样，每个样点分别取样23次。将采集的土样带回实验室采用烘干称重法进行含水量测试，烘干温度为105℃，烘干时间约为12h，取平均值得到土壤湿度数据。

利用神东矿区实测0-5cm、10cm深度土壤湿度数据对TVDI进行验证。以10cm土壤湿度(Soil Moisture，SM)为横坐标，以双抛物线型NDVI-T_s特征空间反演得到神东矿区TVDI_c为横坐标，构建SM-TVDI散点图，如图4所示，并与三角形NDVI-T_s特征空间得到的TVDIt进行对比，计算其相关系数，如表1所示。

表1神东矿区TVDI与土壤湿度线性相关系数(R²)

注：***，**，*分别表示通过99％、95％、90％显著性检验。

利用MODIS数据产品NDVI和地表温度T_s数据，对神东矿区土壤湿度状况进行了遥感监测与分析，表明NDVI-T_s三角形特征空间并不能完整描述植被指数与地表辐射温度的关系。本发明将NDVI扩展到0.15以内，提出NDVI-T_s双抛物线型特征空间，其监测土壤湿度状况优于NDVI-T_s三角形特征空间。250m双抛物线型NDVI-T_s特征空间在监测地表0-5cm土壤湿度方面是最好的，而500m双抛物线型NDVI-T_s特征空间在监测10cm深土壤湿度时更有优势。