一种基于特征空间的土壤水分微波产品降尺度方法

技术领域

本发明涉及遥感技术领域，具体涉及一种基于特征空间的土壤水分微波产品降尺度方法。

背景技术

土壤水分是地表水热平衡的重要决定参数，其时空分布信息对于农作物估产、水文预报以及旱情监测等具有重要意义。传统的地面观测技术（譬如土钻烘干法、时域反射仪、中子仪）虽能准确测量站点尺度的土壤水分，但存在着仪器昂贵、费时费力等明显缺点，难以有效反映区域尺度土壤水分时空变化细节信息。在此背景下，卫星遥感成为大范围土壤水分监测的主流方法，特别是微波遥感，由于具有穿透云雾和地表的能力，可以实现土壤水分的全天候重复监测，现有的土壤水分遥感产品基本都是基于微波遥感研发的。然而受限于传感器设计技术，当前土壤水分微波产品的空间分辨率普遍偏低（多为几十千米），这在很大程度上影响了其在气象、水文和农业领域的实际应用。在此背景下，如何研发土壤水分微波产品降尺度方法，获取高空间分辨率土壤水分信息，成为业界关注的热点。

研究发现，基于光学遥感反演的地表温度（LST）和植被指数(VI)，不但具有较高的空间分辨率，而且与土壤水分本身具有很好的相关性，基于特定空间范围两者的散点图，可以进行高空间分辨率土壤水分指数的构建，进而实现土壤水分微波产品的降尺度，该方法称为特征空间法。特征空间法的基本假设是：若选定的空间范围内存在着足够多可以充分反映地表土壤水分和植被覆盖度变化情况的像元，则由这些像元构成的LST-VI特征空间将形成具有物理意义的三角形或者梯形；三角形或梯形的上下边界分别代表土壤水分的两种极端条件，每个像元的土壤湿度指数可以根据其在上下边界中的相对位置，通过插值的方法求得，如图1所示。可以看出，特征空间法成功应用的关键在于上下边界的正确构建。在早期的特征空间方法中，上下边界是完全基于LST-VI的散点分布进行拟合的，很明显，这种方法具有很大的随机性和经验性。因此，后续很多研究基于地表能量平衡原理，提出了上下边界的理论计算方法以克服上述缺陷。虽然这些理论方法能够在很大程度上提高土壤水分的监测精度，但由于地表能量平衡原理过于复杂，需要结合较多的地面实测参数才能开展，限制了其在地面实测资料匮乏地区的应用，亟需研发一种脱离地面实测数据依赖、完全基于遥感数据开展的特征空间构建方法。

发明内容

针对现有技术中的上述不足，本发明提供了一种基于特征空间的土壤水分微波产品降尺度方法。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于特征空间的土壤水分微波产品降尺度方法，包括以下步骤：

S1、将土壤水分微波产品与地表温度和植被指数光学产品进行像元匹配，获取低空间分辨率土壤水分与高空间分辨率地表温度和植被指数的一对多特征参数；

S2、依次遍历各个土壤水分微波像元，采用特征空间法获取对应土壤水分微波像元内地表温度最大值和地表温度最小值，并在年尺度上拟合地表温度最大值和地表温度最小值分别与太阳天顶角的余弦函数关系；

S3、根据地表温度最小值与太阳天顶角的余弦函数关系确定特征空间下边界，并基于土壤水分微波像元内光学像元的土壤湿度指数与土壤水分值对地表温度最大值进行修正，根据修正后的地表温度最大值确定特征空间上边界；

S4、根据特征空间下边界和特征空间上边界计算对应土壤水分微波像元与光学产品空间分辨率相匹配的土壤水分，得到所有土壤水分微波像元的降尺度结果。

可选择地，步骤S1具体包括以下分步骤：

S11、提取土壤水分微波产品的粗分辨率像元边界；

S12、利用土壤水分微波产品的粗分辨率像元边界分别提取地表温度和植被指数光学产品的高空间分辨率光学像元；

S13、将植被指数光学像元中的植被指数转化为植被覆盖度。

可选择地，步骤S13中将植被指数光学像元中的植被指数转化为植被覆盖度的方法为：

其中，Fc为植被覆盖度，NDVI为植被指数，NDVI

可选择地，步骤S2中在年尺度上拟合地表温度最大值和地表温度最小值分别与太阳天顶角的余弦函数关系的方法为：

其中，LSTmax为地表温度最大值，Amax为地表温度最大值与太阳天顶角的余弦函数的振幅，

可选择地，步骤S3中基于土壤水分微波像元内光学像元的土壤湿度指数与土壤水分值对地表温度最大值进行修正的方法为：

S31、设定地表温度最大值与太阳天顶角的余弦函数的振幅的修正迭代区间和迭代步长；

S32、根据每天的太阳天顶角和地表温度最大值与太阳天顶角的余弦函数计算日尺度下的地表温度最大值；

S33、根据地表温度光学像元中的地表温度、植被指数光学像元转化的植被覆盖度、空气温度、日尺度下的地表温度最大值和地表温度最小值计算每个光学像元的土壤湿度指数；

S34、根据每个光学像元的土壤湿度指数和土壤水分微波像元的土壤水分计算每个光学像元的土壤水分；

S35、计算同一土壤水分微波像元在年尺度下获取的有效光学像元总量，并分别计算各个光学像元的土壤湿度指数与土壤水分的相关系数；

S36、判断当前迭代条件下计算的所有相关系数的最大值是否大于前一迭代条件下计算的所有相关系数的最大值；

若是，则进入下一迭代条件，跳转至步骤S32；

否则，以当前迭代条件下相关系数最大值对应的振幅修正值作为上边界最优振幅，迭代结束。

可选择地，步骤S33中计算每个光学像元的土壤湿度指数的方法为：

其中，SMI为光学像元的土壤湿度指数，LST为地表温度光学像元中的地表温度，Fc为植被指数光学像元转化的植被覆盖度，Ta为空气温度，LSTmax为地表温度最大值，LSTmin为地表温度最小值。

可选择地，步骤S34中计算每个光学像元的土壤水分的方法为：

其中，SM

可选择地，步骤S35中计算各个光学像元的土壤湿度指数与土壤水分的相关系数的方法为：

其中，Corr为光学像元的土壤湿度指数与土壤水分的相关系数，n为一年内的有效天数，

可选择地，步骤S4中根据特征空间下边界和特征空间上边界计算对应土壤水分微波像元与光学产品空间分辨率相匹配的土壤水分的方法为：

S41、根据每天的太阳天顶角、特征空间下边界和特征空间上边界计算日尺度下的地表温度最小值和地表温度最大值；

S42、根据地表温度光学像元中的地表温度、植被指数光学像元转化的植被覆盖度、空气温度、日尺度下的地表温度最大值和地表温度最小值计算每个光学像元的土壤湿度指数；

S43、根据每个光学像元的土壤湿度指数和土壤水分微波像元的土壤水分计算每个光学像元的土壤水分。

本发明具有以下有益效果：

本发明通过土壤水分微波产品与地表温度、植被指数光学产品的有机结合，提出了一种迭代统计寻优的特征空间构建方法，该方法既能消除原始特征空间法的经验性与随机性，又能降低地表能量平衡方法对地面观测数据的依赖性，而且直接对标土壤水分微波遥感产品，在保证了土壤水分微波产品降尺度精度的同时，实现了特征空间法在区域尺度的纯遥感高效运行。

附图说明

图1为特征空间法原理示意图；

图2为本发明实施例中一种基于特征空间的土壤水分微波产品降尺度方法的流程示意图；

图3为本发明实施例中上下边界余弦函数拟合示意图；

图4为本发明实施例中程序运行结果示意图；

图5为本发明实施例中土壤水分降尺度精度示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

本发明的土壤水分微波产品降尺度方法，首先通过微波像元与光学像元的空间匹配，提取特征空间构建的关键参数，而后通过迭代统计寻优的方法构建特征空间上下边界，最后在最优边界的约束下完成土壤水分微波产品的空间降尺度，实施方式简单概括为像元空间匹配——关键参数提取——原始边界拟合——边界迭代修正——最优边界确立——统计降尺度。

如图2所示，本实施例提供的一种基于特征空间的土壤水分微波产品降尺度方法，包括以下步骤S1至S4：

S1、将土壤水分微波产品与地表温度和植被指数光学产品进行像元匹配，获取低空间分辨率土壤水分与高空间分辨率地表温度和植被指数的一对多特征参数；

在本发明的一个可选实施例中，步骤S1具体包括以下分步骤：

S11、提取土壤水分微波产品的粗分辨率像元边界；

S12、利用土壤水分微波产品的粗分辨率像元边界分别提取地表温度和植被指数光学产品的高空间分辨率光学像元；

S13、将植被指数光学像元中的植被指数转化为植被覆盖度。

具体而言，以SMAP土壤水分微波产品和MODIS光学产品为例，SMAP土壤水分微波产品的空间分辨率为36*36 km，MODIS光学产品的空间分辨率为1*1 km，因此利用土壤水分微波产品的36*36 km空间分辨率像元边界分别提取地表温度和植被指数光学产品的1*1 km空间分辨率的光学像元，每个土壤水分（SM）微波像元对应着1296个地表温度（LST）和植被指数（NDVI）光学像元；进而将植被指数光学像元中的植被指数转化为植被覆盖度，表示为：

其中，Fc为植被覆盖度，NDVI为植被指数，NDVI

S2、依次遍历各个土壤水分微波像元，采用特征空间法获取对应土壤水分微波像元内地表温度最大值和地表温度最小值，并在年尺度上拟合地表温度最大值和地表温度最小值分别与太阳天顶角的余弦函数关系；

在本发明的一个可选实施例中，本实施例所针对的研究区域内可能存在多个土壤水分微波像元，因此依次遍历各个土壤水分微波像元分别进行后续处理。

本实施例采用特征空间法获取对应土壤水分微波像元内地表温度最大值和地表温度最小值；并在年尺度上拟合地表温度最大值和地表温度最小值分别与太阳天顶角的余弦函数关系，具体方法为：

其中，LSTmax为地表温度最大值，Amax为地表温度最大值与太阳天顶角的余弦函数的振幅，

具体而言，本实施例将微波像元范围内1296个地表温度LST数值取最大值（LSTmax）和最小值（LSTmin），由于在均一下垫面条件下，地表温度LST年内变化主要受控于太阳辐射，因此地表温度最大值LSTmax和地表温度最小值LSTmin的年内变化可用太阳天顶角θ的余弦函数非线性拟合反映，太阳天顶角θ可以通过MOD03产品获取，拟合效果示例如图3所示。

S3、根据地表温度最小值与太阳天顶角的余弦函数关系确定特征空间下边界，并基于土壤水分微波像元内光学像元的土壤湿度指数与土壤水分值对地表温度最大值进行修正，根据修正后的地表温度最大值确定特征空间上边界；

在本发明的一个可选实施例中，特征空间下边界代表土壤水分饱和条件，由于降水等原因，特征空间的下边界是真实存在的，因此采用地表温度最小值拟合的函数关系可以直接反映下边界。特征空间上边界代表极端干旱条件，在现实中往往是不存在的，这意味着通过上述公式拟合的Amax值偏小，因此需要通过迭代寻优的方法修正采用地表温度最大值拟合余弦函数的振幅。

本实施例基于土壤水分微波像元内光学像元的土壤湿度指数与土壤水分值对地表温度最大值进行修正的方法为：

S31、设定地表温度最大值与太阳天顶角的余弦函数的振幅的修正迭代区间和迭代步长；

具体而言，本实施例设置振幅Amax的修正迭代区间为[Amax,

S32、根据每天的太阳天顶角和地表温度最大值与太阳天顶角的余弦函数计算日尺度下的地表温度最大值；

具体而言，本实施例可以根据每天的太阳天顶角和地表温度最大值与太阳天顶角的余弦函数确定每天的地表温度最大值。

S33、根据地表温度光学像元中的地表温度、植被指数光学像元转化的植被覆盖度、空气温度、日尺度下的地表温度最大值和地表温度最小值计算每个光学像元的土壤湿度指数；

具体而言，本实施例在地表温度最大值LSTmax和地表温度最小值LSTmin已知的条件下，根据特征空间模型计算每个光学像元的土壤湿度指数SMI，计算方法为：

其中，SMI为光学像元的土壤湿度指数，LST为地表温度光学像元中的地表温度，Fc为植被指数光学像元转化的植被覆盖度，Ta为空气温度，LSTmax为地表温度最大值，LSTmin为地表温度最小值。其中空气温度Ta可以通过MOD07_L2产品获取。

S34、根据每个光学像元的土壤湿度指数和土壤水分微波像元的土壤水分计算每个光学像元的土壤水分；

具体而言，本实施例在每个光学像元的土壤湿度指数SMI和土壤水分微波像元的土壤水分SM已知的条件下，计算每个光学像元的土壤水分SM

其中，SM

S35、计算同一土壤水分微波像元在年尺度下获取的有效光学像元总量，并分别计算各个光学像元的土壤湿度指数与土壤水分的相关系数；

具体而言，由于SMAP和MODIS卫星的过境日期不统一，并且光学遥感不同于微波遥感，其数据质量受云影响较大，仅在大面积晴天条件下才有充分保障。因此，为了解决在不同天的样本匹配上存在的质量控制问题，本实施例判断当SMAP像元内的MODIS有效样本量大于1296*80%时，通过步骤S32至S34计算得到年尺度下不同天数的每个光学像元的土壤水分SM

本实施例通过上述判断，计算同一土壤水分微波像元在年尺度下获取的有效光学像元总量为

其中，Corr为光学像元的土壤湿度指数与土壤水分的相关系数，n为一年内的有效天数，

S36、判断当前迭代条件下计算的所有相关系数的最大值是否大于前一迭代条件下计算的所有相关系数的最大值；

若是，则进入下一迭代条件，跳转至步骤S32；

否则，以当前迭代条件下相关系数最大值对应的振幅修正值作为上边界最优振幅，迭代结束。

具体而言，随着幅值Amax在修正区间[Amax,∞)内的变化，每个幅值Amax将求得一个相关系数Corr，而且由于样本量是一致的，相关系数Corr之间是可比较的，最优的幅值Amax应对应最大的相关系数Corr。据此，为了提高幅值Amax迭代修正的运行效率，本实施例判断当Amax的继续增加不足以进一步提高Corr时，即

S4、根据特征空间下边界和特征空间上边界计算对应土壤水分微波像元与光学产品空间分辨率相匹配的土壤水分，得到所有土壤水分微波像元的降尺度结果。

在本发明的一个可选实施例中，本实施例根据特征空间下边界和特征空间上边界计算对应土壤水分微波像元与光学产品空间分辨率相匹配的土壤水分的方法为：

S41、根据每天的太阳天顶角、特征空间下边界和特征空间上边界计算日尺度下的地表温度最小值和地表温度最大值；

S42、根据地表温度光学像元中的地表温度、植被指数光学像元转化的植被覆盖度、空气温度、日尺度下的地表温度最大值和地表温度最小值计算每个光学像元的土壤湿度指数；

S43、根据每个光学像元的土壤湿度指数和土壤水分微波像元的土壤水分计算每个光学像元的土壤水分。

具体而言，本实施例将得到的最优Amax带入余弦函数，即可求得任意一天的特征空间上边界，下边界依然采用步骤S2中LSTmin的拟合函数，据此求得SMAP土壤水分微波像元逐日所有MODIS光学像元的土壤湿度指数SMI，进而将其转化为土壤水分SM。

通过上述步骤遍历完研究区域内的所有土壤水分微波像元即可获得研究区所有土壤水分微波像元的降尺度结果，模型输出的空间分辨率为1 km。

上述以SMAP土壤水分微波产品和MODIS光学产品为例的程序运行结果如图4所示，左侧为原始SMAP土壤水分数据，右侧为土壤水分降尺度结果；基于站点实测值的精度验证结果如图5所示。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。