一种土壤含水量产品降尺度方法与流程

本发明涉及地理测绘技术领域，特别是涉及一种土壤含水量产品降尺度方法。

背景技术：

土壤含水量作为气候与环境干旱化的重要指示因子之一，影响着土壤的理化性质与植被的生长，进而影响我国的粮食产量。同时，土壤含水量是地表能量平衡和水循环的重要组成部分，是全球变化研究中的重要监测因子。目前，各种基于全球观测站点资料建立的土壤含水量数据集，由于观测点的密度和空间代表性不足，模拟和预报的精度难以满足应用需求。

先进微波扫描辐射计AMSR-2(Advanced Microwave Scanning Radiometer for the Earth Observing System)，土壤水分和海洋盐分传感器SMOS(Soil Moisture and Ocean Salinity)，SMAP(Soil Moisture Active Passive)以及风云三号气象卫星均具有全天时、全天候、观测尺度大、重访周期短等突出优势，能够提供高覆盖度高时效性的全球土壤含水量数据。但这些数据的空间分辨率(9-40km)较低，无法满足流域尺度土壤含水量时空动态监测的需求。光学遥感数据空间分辨率可以达到1km以下，高空间分辨率、低时间分辨率、易受天气影响的特点，与被动微波遥感数据恰恰相反。

为了解决被动微波土壤含水量产品的空间分辨率低的问题，国内外学者也提出了各种不同的降尺度方法：一类降尺度模型是以遥感技术获得的土壤物理参数为基础，如利用光学遥感数据反演土壤的蒸散量建立降尺度算法，但是这种算法的局限性在于没有考虑微波遥感数据反演的土壤含水量和光学遥感数据反演的土壤有效蒸发量之间存在很强的非线性关系。还有一类模型是四维变分同化方法，就是在微波遥感数据的第四维尺度上分析土壤含水量。但是这种方法对地表资料的获得要求较高。

由此可见，目前迫切需要发展一种新的土壤含水量产品降尺度方法，建立一种可行的空间降尺度模型，逐步提高土壤含水量数据的空间分辨率，推动流域尺度土壤含水量的时空动态监测。

技术实现要素：

本发明要解决的技术问题是提供一种土壤含水量产品降尺度方法，使其提高被动微波土壤含水量数据产品的空间分辨率，满足流域尺度水资源和农业管理的应用需求，克服目前土壤含水量产品数据集空间代表性不足、空间分辨率低等的不足。

为解决上述技术问题，本发明提供一种土壤含水量产品降尺度方法，所述方法包括以下步骤：

A.获取待研究区的被动微波土壤含水量产品和同一时间的光学遥感影像数据；

B.基于多端元混合像元分解方法对所述光学遥感影像数据进行土壤光谱的提取；

C.利用GA-PLS建立所述土壤光谱反射特征与从所述被动微波土壤含水量产品中获取的土壤含水量之间的定量关系模型；

D.基于步骤C建立的所述定量关系模型，利用泰勒级数展开形式构建土壤含水量降尺度模型，获得高空间分辨率的土壤含水量数据。

作为本发明的一种改进，所述步骤A中被动微波土壤含水量产品采用SMAP土壤含水量数据；所述光学遥感影像数据采用MODIS影像数据。

进一步改进，所述步骤B中提取所述土壤光谱的方法为：将高空间分辨率的MODIS图像重采样到与SMAP数据同样的低空间分辨率，分别对该高空间分辨率的MODIS影像和重采样后低空间分辨率的MODIS应用MESMA方法进行土壤光谱的提取。

进一步改进，所述MESMA方法包括光谱库创建、最优光谱库选取和多端元混合像元分析步骤，

所述光谱库创建包括基于ROI创建光谱库、光谱库元数据制作和光谱库管理；

所述最优光谱库的选取包括创建方形阵列和光谱库优选；

所述多端元混合像元分析采用植被-不透水面-土壤模型，将优选的植被、不透水面、土壤光谱集组合构成2EM、3EM、4EM混合像元分析模型，基于所述最优光谱库的优选结果对MESMA结果进行阴影归一化处理，得到各端元丰度值和表示结果精度的均方根误差，再利用所述各端元丰度值和下列公式得到研究区土壤端元光谱，

其中，R_s(λ)为土壤光谱在波段λ的反射率，R(λ)为像元在波段λ上的反射率，R(i，λ)为第i个端元在波段λ上的反射率，f_i为第i个端元丰度值，N为端元个数，ε_λ是残差，所有端元组分的丰度值之和定义为1。

进一步改进，所述步骤C建立所述土壤光谱反射特征与所述土壤含水量定量关系模型是基于所述土壤光谱计算得到的每个像元的土壤光谱中各波段的反射率、波段比值、曲率与所述土壤含水量的定量关系。

进一步改进，所述步骤D中泰勒级数展开形式的表达式为：

其中，θ_n-1和θ_n分别代表低空间分辨率和高空间分辨率土壤含水量，R_n-1(λ_i)_s和R_n(λ_i)_s分别表示低空间分辨率和高空间分辨率的土壤光谱在第i波段的反射率，Ratio_n-1(j)_s)和Ratio_n(j)_s分别表示低空间分辨率和高空间分辨率的土壤光谱波段比值，Curv_n-1(k)_s和Curv_n(k)_s分别表示高空间分辨率和低空间分辨率的土壤光谱曲率，M、N和L分别代表反射率i、波段比值j和曲率k的变量总数。

进一步改进，所述步骤D中降尺度模型的建立方法采用逐步递减方式进行降尺度运算，逐步达到高空间分辨率的要求。

进一步改进，所述步骤D中高空间分辨率的土壤含水量数据是指500m空间分辨率的土壤含水量数据。

进一步改进，所述步骤A中还同时获取所述待研究区的土壤辅助数据，所述土壤辅助数据包括：土地利用或土地覆盖分类数据，以及DEM数据。

进一步改进，所述步骤C建立的GA-PLS模型还可以增加每个像元的DEM数据或坡度数据作为GA-PLS模型的输入参数。

采用上述技术方案，本发明至少具有以下优点：

(1)本发明基于光学遥感和被动微波遥感数据，获取高空间分辨率的土壤含水量方法，可满足大范围流域尺度区域研究，准确度高，易于建立，省时省力。

(2)本发明可扩展性高，在应用的过程中，可根据实际情况，进行辅助数据或土壤含水量定量反演模型关系项的增减，还可以采用逐级回归式的降尺度方法，以逐步推进的方式提高被动微波土壤含水量产品的空间分辨率，不断提高模型的计算精度。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1是本发明土壤含水量产品降尺度方法的原理流程图。

图2是本发明中土壤最优光谱选取结果图。

图3是本发明中植被最优光谱选取结果图。

图4是本发明中9km空间分辨率MODIS图像MESMA结果图。

图5是本发明中样本点分布图。

图6是本发明中9km尺度GA-PLS模型精度评价结果。

图7是本发明中5km空间分辨率MODIS图像MESMA结果图。

图8是本发明中5km土壤含水量降尺度精度评价结果。

图9是本发明中500m空间分辨率MODIS图像MESMA结果图。

图10是本发明中5km尺度GA-PLS模型精度评价结果。

图11是本发明中500m土壤含水量降尺度精度评价结果。

具体实施方式

本发明在现有技术的基础上，利用土壤的反射光谱特性，即土壤含水量是影响土壤光谱反射率的一个重要因素，又由于土壤含水量具有非常大的时空变异性，高空间分辨率的土壤光谱信息能够更好的反映土壤含水量在时空范围上的变化特征，则本发明利用土壤的光谱特征与土壤含水量之间的定量关系来完成土壤含水量的降尺度研究，为土壤含水量的高时空分辨率监测提供了一个全新的思路。

又鉴于被动微波遥感在全球土壤含水量数据获取方面的优势，创设一种可行的空间降尺度模型，逐步提高土壤含水量数据的空间分辨率，通过光学遥感和被动微波遥感的综合应用，提高被动微波土壤含水量数据产品的实用性，最终推动流域尺度土壤含水量的时空动态监测。其具体的土壤含水量产品降尺度方法如下。

参照附图1所示，本发明土壤含水量产品降尺度方法，主要包括以下步骤：

A.获取被动微波土壤含水量产品和同一时间的光学遥感影像数据；

其具体的数据获取方法为：

获得待研究区SMAP土壤含水量产品。SMAP卫星是于2015年1月由美国NASA发射，其上搭在一个L波段的雷达和一个L波段的辐射计。下述实施例选取的是SMAP L4空间分辨率为9km的土壤水日平均数据。下载地址为：https://ns idc.org/。

目前高空间分辨率的光学图像时间分辨率往往较低，受气象条件的影响，影像获取成功率低，相比之下，中低空间分辨率的光学影像因时间分辨率高更易获取。下述实施例光学遥感影像数据选取的是与SMAP土壤含水量数据同步的空间分辨率为500m的MODIS陆地产品中的MOD09地面反射率数据。其中采用该产品的1-7波段(620～670，841～876，459～479，545～565，1230～1250，1628～1652，2105～2155nm)数据建立土壤含水量反演模型。

该步骤中同时还获取待研究区土壤辅助数据，其包括：(1)土地利用/土地覆盖分类数据，可以获得MODIS土地覆盖类型产品，用于混合像元分解过程中端元组分的选取和验证，最大程度去除植被等其它地物对农业土壤光谱的影响；(2)DEM数据，在受到地形影响的地区，该数据可用于本发明降尺度模型的修正和验证。

B.基于多端元混合像元分解方法进行土壤光谱的提取；

多端元混合像元分解法(Multiple Endmember Spectral Mixture Analysis，MESMA)，首先以端元的物理意义为理论基础，为每类地物选取多条光谱，并以此生成多个端元组合(每个端元组合由不同地物中的某一条光谱组成)，接着对每个像元寻找最小二乘法误差最小的端元组合，进而求出每个像元的端元比例。

本发明中，将500m空间分辨率的MODIS图像重采样到与SMAP数据同样的空间分辨率9km，分别对500m的MODIS影像和重采样后9km的MODIS应用MESMA方法进行土壤光谱的提取，获得9km空间分辨率和500m空间分辨率的土壤光谱。

其中，MESMA处理分析过程主要包括光谱库创建、最优光谱库选取、多端元混合像元分析三部分。其光谱库创建包括基于ROI创建光谱库、光谱库元数据制作和光谱库管理；其最优光谱库的选取是多端元混合像元分解成功的关键，该最优光谱库选取过程包括创建方形阵列和光谱库优选，其优选的方法是结合COB(Count-based Endmember Selection)、EAR(Endmember Average RMES)、MASA(Minimum Average Spectral Angle)计算，得到最能代表各类地物的光谱库；其混合像元分解过程主要采用V-I-S(植被-不透水面-土壤)模型，将优选的植被、不透水面、土壤光谱集组合构成2EM、3EM、4EM混合像元分析模型，基于优选结果对MESMA结果进行阴影归一化处理，最终得到不同端元组合的混合像元分解结果，即各端元丰度值和表示结果精度的均方根误差(RMSE)。

再利用上述MESMA得到的端元丰度值和下列公式(1)得到研究区土壤端元光谱。

其中，R_s(λ)为土壤光谱在波段λ的反射率，R(λ)为像元在波段λ上的反射率，R(i，λ)为第i个端元在波段λ上的反射率，f_i为第i个端元丰度值。N为端元个数，ε_λ是残差。所有端元组分的丰度值之和定义为1。

本实施例还结合ENVI和Matlab设计了纯净像元光谱的自动化提取方法。利用ENVI软件提取各个采样点端元组合模型的丰度值、RMSE和残差，确定采样点选取的端元组合模型并输出采样点组合丰度值f_i。

首先在随机选取的采样点中根据土地利用分类类型剔除掉包含不透水面的采样点，对影像只建立S-V模型，对S-V光谱的组合模型赋属性值，提取各组合模型中土壤端元和植被端元的坐标值，提取光谱库各条光谱坐标值，并对两者进行匹配，得到S-V光谱组合模型的组合光谱；提取各个采样组合模型属性值，与S-V光谱组合模型库中属性值匹配得到各个采样点的S-V端元组合光谱R(i,λ)。根据上述混合像元分解公式(1)计算得到每一个采样点土壤光谱。

C.利用GA-PLS建立土壤光谱和被动微波获取的土壤含水量之间的定量关系模型；

在上述9km空间分辨率的土壤光谱获取的基础上，基于土壤光谱随土壤含水量变化的特征规律，引入偏最小二乘－遗传算法(GA-PLS)，构建土壤含水量与光谱的定量关系模型，作为降尺度的基础。

其中PLS是多因变量对多自变量的回归建模方法，可以应用在全谱数据或部分谱数据的分析。它将数据矩阵的分解和回归相结合，得出与预测组分相关的特征值向量，这使其可以应用于复杂的分析体系，得到更稳健的结果。遗传算法(GA)的基本思想是在演化过程中进行自然选择，该演化过程由基因重组和变异来实现。该算法最显著的优点是避免了初始值选择的问题。本发明为了减小计算量，提高PLS模型的精度，使用遗传算法(GA)对光谱参数进行筛选，剔除对含水量变化反应不明显的变量。

本发明使用GA-PLS模型建立土壤含水量反射光谱特征模型，如将9km的MODIS反射率、波段比值、曲率及SMAP L4土壤水产品作为GA-PLS模型的输入变量，得到9km尺度上的MODIS反射率、波段比值、曲率与土壤含水量之间的GA-PLS模型，并进行精度评价。

D.基于步骤C建立的定量关系模型，利用泰勒级数展开的形式构建土壤含水量降尺度模型，获得空间分辨率为500m的土壤含水量数据。

本发明降尺度的依据是土壤含水量与土壤光谱反射率之间的定量关系，即在步骤C中阐述的基于被动微波土壤水产品与MODIS土壤反射特征如各波段的反射率、波段比值(如MODIS3/MODIS1)、曲率(如：MODIS3×MODIS1/(MODIS2)²)的GA-PLS模型。在空间尺度的变化过程中，土壤含水量与反射率、波段比值、曲率等的关系变化是非线性的，所以本发明利用泰勒级数展开的形式，如公式(2)，来构建不同尺度下的土壤含水量之间的关系。

其中，θ_n-1和θ_n分别代表低空间分辨率(如：9km)和高空间分辨率(如：500m)土壤含水量，R_n-1(λ_i)_s和R_n(λ_i)_s分别表示低空间分辨率和高空间分辨率的土壤光谱在第i波段的反射率，Ratio_n-1(j)_s)和Ratio_n(j)_s分别表示低空间分辨率和高空间分辨率的土壤光谱波段比值，Curv_n-1(k)_s和Curv_n(k)_s分别表示高空间分辨率和低空间分辨率的土壤光谱曲率。M、N和L分别代表反射率(i)、波段比值(j)和曲率(k)的变量总数。关系式f表示由GA-PLS模型分别得到的关于低空间分辨率土壤光谱反射率、波段比值和曲率与土壤含水量的关系式。

该步骤D可根据情况，采用逐步递减的方式进行降尺度运算，如先选择一个尺度间隔t(如4km)将9km土壤含水量数据降到(9-t)尺度，通过反复使用公式(2)，将前一步计算出的高空间分辨率数据作为下一次公式低空间分辨率数据再次计算得到更高空间分辨率的土壤含水量，最终逐步逼近达到最后500m空间分辨率的要求。其中最优尺度间隔t，可以根据实际数据获取情况和结果精度来确定。

该分步降尺度方法的另一个目的是在降尺度计算过程中，不断控制土壤含水量结果的反演精度。每降一级尺度反演的土壤含水量结果，将与同一空间分辨率土壤含水量实测数据建立相关关系，若验证数据相关性较高，降尺度计算继续；若验证数据相关性较差，降尺度运算适时停止，并进行影响因素分析。如果降尺度的精度受到地形、植被等因素的影响，可为降尺度模型增加新的关系项，及时进行调整和优化。最后，模型的降尺度结果将通过同步野外土壤含水量试验数据得到验证。

下述以对美国中部平原地区的SMAP土壤含水量数据进行降尺度的具体实例来说明本发明基于土壤光学特性的土壤含水量产品降尺度方法的过程。

在该实例中，被动微波土壤含水量产品选择9km空间分辨率SMAP L4数据，光学图像选择500m空间分辨率的MODIS 09A1数据。

获取数据后，首先对MODIS数据利用ENVI软件进行投影转换、镶嵌、裁剪、重采样等预处理。

然后进行最优光谱库的选取，通过目视识别的方式分别建立土壤、植被类型光谱库，利用COB(Count-based Endmember Selection)、EAR(Endmember Average RMES)、MASA(Minimum Average Spectral Angle)计算优选最能代表土壤和植被的光谱曲线，建立最优光谱库，最终选择出土壤最优光谱13条，如附图2所示，植被最优光谱6条，如附图3所示。

再将优选的植被、土壤光谱集组合构成3EM混合像元分析模型，对结果进行阴影归一化处理，最终得到不同端元组合的混合像元分解结果。将500m空间分辨率的MODIS图像进行重采样，获得9km、5km空间分辨的图像。

对9km分辨率的MODIS图像进行MESMA分析，如附图4所示。计算得到9km分辨率土壤光谱的反射率、波段比值和曲率。随机选取51个采样点，如附图5所示，其中，36个点作为校正集，15个点作为验证集，分别提取51个采样点的土壤光谱和SMAP土壤含水量，并将其作为参数输入建立土壤光谱和含水量之间的GA-PLS模型。结果显示，校正集决定系数为0.78，验证集决定系数为0.72，如附图6所示。

GA-PLS模型作为降尺度模型的基础，建立逐级空间降尺度模型。逐级降尺度选择空间间隔t＝4km，先将土壤含水量降至5km，需要对重采样至5km分辨率MODIS数据进行MESMA处理分析，如附图7所示，进一步计算得到5km尺度的土壤光谱的反射率、波段比值和曲率，利用上述公式(2)得到5km空间分辨率的土壤含水量，经验证5km分辨率含水量数据与实测数据相关性较9km并未降低，与站点得到的实测数据进行对比分析的精度R²为0.54，如附图8所示。

重复上一步操作，在5km尺度下建立土壤光谱和含水量的GA-PLS模型，共选取了89个采样点，如附图5所示，其中62个校正集样本，27个验证集样本，结果显示，5km尺度下校正集决定系数为0.71，验证集决定系数为0.67，如附图9所示。将500m分辨率MODIS数据进行MESMA处理分析，如附图10所示，并计算得到500m尺度的土壤光谱的反射率、波段比值和曲率，利用5km尺度建立的GA-PLS模型，再次运用公式(2)得到500m空间分辨率土壤含水量数据。利用实测土壤含水量与500m土壤含水量进行相关性分析，与站点得到的实测数据进行对比分析的精度R²为0.49，如附图11所示。

本发明利用被动微波遥感数据和光学遥感数据二者在时空分辨率上的优势，将二者有效整合进行土壤含水量产品的降尺度研究，实现了流域尺度土壤含水量实时或准实时的动态监测。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。