1.融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,包括如下步骤:步骤1,获取目标山区除水体区域之外的其他区域的遥感数据,包括旋风全球导航卫星系统接收数据、数字高程模型数据、中分辨率成像光谱仪观测数据以及SoilGrids世界土壤信息数据;同时获取目标山区各观测站的土壤湿度数据;对获取的遥感数据和土壤湿度数据进行预处理;
步骤2,选取海拔高度、坡度、坡向、经纬度、植被指数、粘土含量以及信噪比作为山区GNSS‑R土壤湿度反演的影响因子,从预处理后的数字高程模型数据中获取海拔高度、坡度、坡向及经纬度数据,利用预处理后的中分辨率成像光谱仪观测数据计算植被指数,从预处理后的SoilGrids世界土壤信息数据中获取粘土含量数据;
步骤3,对预处理后的旋风全球导航卫星系统接收数据中的信噪比进行校正,利用预处理后的数字高程模型数据计算地形湿度指数;将海拔高度、坡度、坡向、经纬度、植被指数、粘土含量、校正后的信噪比、地形湿度指数以及预处理后的观测站土壤湿度数据时空匹配为1km分辨率网格下的日数据;
步骤4,获取信噪比不为空值的1km分辨率网格所对应的信噪比和地形湿度指数日数据,将其中的信噪比作为因变量,地形湿度指数作为自变量,建立每日对应的空间拟合函数模型,利用每日对应的空间拟合函数模型对当日信噪比为空值的1km分辨率网格的信噪比进行空间化,得到每日所有1km分辨率网格的信噪比;
步骤5,将所有1km分辨率网格对应的海拔高度、坡度、坡向、经纬度、植被指数、粘土含量、校正后的信噪比以及预处理后的观测站土壤湿度日数据作为数据集,将数据集按时间顺序从前往后排序,从数据集最前面开始,选取一段时间中信噪比为校正值且土壤湿度不为空值的1km分辨率网格对应的日数据作为训练集和验证集,将剩余时间对应的所有日数据作为测试集,测试集包括信噪比为空值和信噪比不为空值的所有1km分辨率网格对应的日数据,其中信噪比为空值的1km分辨率网格的信噪比通过步骤4空间化获得;
步骤6,构建山区GNSS‑R土壤湿度反演模型,包括随机森林模型、XGBoost模型、支持向量机模型以及Stacking集成模型,海拔高度、坡度、坡向、经纬度、植被指数、粘土含量和校正后的信噪比均作为随机森林模型、XGBoost模型和支持向量机模型的输入,将随机森林模型、XGBoost模型和支持向量机模型输出的土壤湿度预测结果作为Stacking集成模型的输入,由Stacking集成模型对随机森林模型、XGBoost模型和支持向量机模型输出的土壤湿度预测结果分配权重,并输出最终的土壤湿度预测结果;利用训练集和验证集对山区GNSS‑R土壤湿度反演模型进行训练和验证,得到训练好的山区GNSS‑R土壤湿度反演模型;
步骤7,将测试集输入训练好的山区GNSS‑R土壤湿度反演模型,得到1km分辨率网格下土壤湿度预测日数据,将每个月的土壤湿度预测日数据取平均,得到1km分辨率网格下土壤湿度预测月数据。
2.根据权利要求1所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,所述步骤1中,对获取的遥感数据和土壤湿度数据进行预处理,具体如下:对旋风全球导航卫星系统接收数据进行预处理,包括:依次剔除入射角大于65°、信噪比小于2dB或大于14dB、接收机增益小于0dB以及反射信号的峰值功率大于‑147dB的接收数据;
对数字高程模型数据进行预处理,包括:剔除数字高程模型数据中的噪声数据;
对中分辨率成像光谱仪观测数据进行预处理,包括:格式转换、投影转换和数据重采样;
对土壤湿度数据进行预处理,包括:剔除各观测站观测到的土壤湿度数据中的缺失值和异常值。
3.根据权利要求1所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,所述步骤2中,从预处理后的数字高程模型数据中获取海拔高度和经纬度;
使用ArcGIS软件中的坡度计算工具计算坡度;其中,坡度的计算公式如下:,
其中, 为坡度,为地形表面上的高程值, 和 表示水平方向上空间位置的变化;
使用ArcGIS软件中的坡向计算工具计算坡向,按顺时针方向计算,正北方向为0度;
植被指数的计算公式如下:
,
其中, 为植被指数, 为红光波段的反射率, 为近红外波段的反射率;
从预处理后的SoilGrids世界土壤信息数据中,选取L波段对应的表层粘土含量数据。
4.根据权利要求1所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,所述步骤3中,校正后的信噪比的计算公式如下:,
其中, 为校正后的信噪比, 为全球导航卫星系统信号的发射功率, 为发射天线的增益, 为采样点与全球导航卫星系统发射机之间的距离, 为采样点与旋风全球导航卫星系统接收机之间的距离,为全球导航卫星系统信号的载波波长, 为预处理后的旋风全球导航卫星系统接收数据中的信噪比;
地形湿度指数的计算公式如下:
,
其中, 为地形湿度指数,为单位汇水面积, 为栅格单元的坡度,且
,
,
其中, 为进入栅格单元的上游总面积,栅格单元即数字高程模型数据的栅格单元;
为有效等高线长度,即单元栅格的长度; 为相邻栅格单元的高度差; 为相邻栅格单元中心的长度。
5.根据权利要求1所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,所述步骤3中,时空匹配具体如下:将海拔高度、坡度、坡向、经纬度、植被指数、粘土含量以及地形湿度指数以最邻近法重采样到1km分辨率网格下的日数据;
将信噪比数据格式转换到1km分辨率网格下,即对于每个网格,将24小时内落在网格内的所有采样点计算得到的信噪比取平均,得到1km分辨率网格下的日信噪比;将24小时内无采样点的网格的日信噪比,记为空值;
将预处理后的土壤湿度数据转换到1km分辨率网格下,即将各观测站预处理后的土壤湿度作为各观测站所在网格的土壤湿度;将无观测站的网格的土壤湿度记为空值。
6.根据权利要求1所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,所述步骤4中,建立每日对应的空间拟合函数模型后,对于空间拟合函数模型计算决定系数,当决定系数大于预设的第三阈值时,表示建立的空间拟合函数模型有效,保留当日建立的空间拟合函数模型拟合出的信噪比;否则将当日建立的空间拟合函数模型拟合出的信噪比删除。
7.根据权利要求1所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法,其特征在于,所述步骤6中,对山区GNSS‑R土壤湿度反演模型进行训练和验证时,采用均方根误差和相关系数表征反演模型的精度,在均方根误差小于预设的第一阈值且相关系数大于预设的第二阈值时,山区GNSS‑R土壤湿度反演模型训练完成;
其中,均方根误差和相关系数的计算公式分别如下:
,
,
其中, 为均方根误差, 为相关系数, 为样本数,为反演模型对第 个样本预测的土壤湿度值, 为第 个样本的对应站点土壤湿度观测值,为所有样本观测值的平均值,为所有样本预测值的平均值。
8.一种计算机设备,包括存储器、处理器,以及存储在所述存储器中并能够在所述处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现如权利要求1至7任一项所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法的步骤。
9.一种计算机可读存储介质,所述计算机可读存储介质存储有计算机程序,其特征在于,所述计算机程序被处理器执行时实现如权利要求1至7任一项所述的融合地形湿度指数的山区GNSS‑R土壤湿度反演方法的步骤。