一种考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法

本发明涉及土壤水分降尺度领域，具体是一种考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法。

背景技术：

1、土壤水分在水资源管理、水文循环和干旱发展等过程中扮演着重要的角色，作为地表与大气之间水热交换的重要变量，对指导农业生产和保证粮食安全至关重要。对于具有复杂的地形、多变的土地覆盖类型、异质的土壤类型区域，迫切需要高时空分辨率的土壤水分数据提升干旱监测的精度和能力。但是当前全球的土壤水分产品通常空间分辨率较低，无法反映多种地表覆盖、不同气候条件和复杂地形共同作用下所导致的土壤水分较大的空间异质性。因此，如何充分挖掘土壤水分数据与现有辅助数据之前复杂的关系，获得高时空分辨率的土壤水分数据，是一个非常重要的问题。

2、土壤水分的产品众多，其中smap(soil moisture active and passive,土壤水分主被动探测计划)同时使用l波段雷达和辐射计开展并行、同步测量，结合了主动(雷达)和被动(辐射计)微波遥感的优势，提供一种能够同时满足高时间分辨率、宽空间覆盖、最佳传感深度以及在中等植被条件下能够高精度的反演土壤水分的要求，被广泛应用于全球土壤水分监测研究中。smap项目生产三种土壤水分产品，产品名称分别为sm_p、sm_a和sm_ap，空间分辨率分别为36km、3km和9km，其中由于l波段雷达在2015年7月失效，目前只能获取36km和9km空间分辨率的土壤水分数据，较低空间分辨率的遥感像元中包含了异质的土地覆盖类型和复杂的地形，无法提供土壤水分的详细变化，难以准确的反映局部区域尺度的干旱信息。

3、土壤水分的空间分布与气候因子、地形地貌、土壤性质和植被等众多因素相互作用。植被通过影响太阳入射辐射和土壤蒸散发进而影响土壤湿度；地表温度通过影响土壤水蒸散发进而影响土壤水分，即地表温度越高，土壤水蒸散发越大，地表则表现出越干旱；高程通过影响地表温度从而影响土壤湿度，且高程将会影响径流量大小以及降雨的再分配。土壤水分时空分布的差异性是由多种因素共同作用的结果，上述这些因子应同时应用于土壤水分降尺度模型中。同时，光学传感器经常受到云污染影响，导致研究区可用的遥感影像数据量往往较少，但是土壤水分降尺度需要多光谱数据的参与，如何弥补研究区可用数据的空间和时间缺失数据是需要重点关注的问题。此外，现有的土壤水分空间分辨率较低，如何提供更多有关土壤水分的精细尺度信息，提升土壤水分的空间分辨率，也是降尺度研究需要考虑的问题。另一方面，植被与土壤水分之间存在一定的时间滞后性，当前土壤水分的状态需要在一定时间段后才能响应在植被上，通常比土壤水分滞后一个月的时间，现有的土壤水分降尺度模型较少考虑植被后滞性对土壤水分的影响。

技术实现思路

1、本发明针对现有技术较少考虑土壤水分与植被之间的滞后性、土壤水分与精细尺度植被信息之间的联系、土壤水分数据与多特征辅助数据之间的关系以及弥补光学数据因云雾影响导致无法有效应用于土壤水分降尺度模型中的问题，提出一种考虑植被记忆和精细尺度植被数据的联合多特征的土壤水分降尺度方法，为土壤水分降尺度方法提供更多有关土壤水分的精细尺度信息和植被记忆信息，同时弥补了研究区可用多光谱数据空间和时间的缺失，充分挖掘土壤水分与多特征数据之间的非线性关系，获取了1km空间分辨率的土壤水分数据。

2、为实现上述目的，按照本发明的一个方面，提供了一种联合多特征的土壤水分降尺度方法，包括以下步骤：

3、(1)采用增强时空自适应融合算法生成30m空间分辨率的精细尺度植被指数；

4、(2)利用mod13a3月尺度ndvi获取与土壤水分数据具有滞后性的ndvi数据，设计土壤水分降尺度辅助因子之一的植被记忆数据，确保植被记忆数据与土壤水分数据时间滞后1个月左右；

5、(3)获取目标区域的smap数据、modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据，对modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据进行预处理使具有相同的投影方式、相同的影像覆盖范围和与smap土壤水分具有一致的空间分辨率，然后与smap数据进行时空匹配以作为土壤水分降尺度模型的输入数据；

6、(4)联合植被记忆数据、精细尺度植被指数、smap数据、modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据，构建土壤水分降尺度模型以获取1km的高空间分辨率土壤水分产品。

7、进一步的，步骤(1)包括：

8、(1.1)准备数据：为预测日期提供两对与预测日期相近，且确保两对数据具有相同的时间和覆盖范围的modis和landsat反射率影像，以及一组用于被预测的不同时间序列的modis影像；

9、(1.2)数据预处理：将所有的modis数据利用envi软件进行重采样，重采样方法选择双线性内插以减少地理参考误差的影响，对重采样后的modis进行地理校正以获得和landsat相同的地理配准，之后对所有的数据进行裁剪操作以确保影像具有完全相同的覆盖区域，所有数据都经过预处理以获得表面反射率，预处理后的modis和landsat影像具有相同的投影、空间分辨率和一致的覆盖范围；

10、(1.3)搜索邻近像元：使用两个高分辨率图像搜索预测期影像局部窗口中心像元的相似像元，相似像元的获取采用滑动窗口的方法获取，相邻像元值需要满足式(1)的条件，即两个像元的标准差较小，以确保在搜索窗口范围内获得与中心像元具有较高光谱相似特性的像元：

11、|l(xi,yi,tk,b)-l(xw/2,yw/2,tk,b)|≤σ(b)·2/m  (10)

12、式中：l为landsat影像，(xi,yi)是第i个相似像元的位置；(xw/2,yw/2)为预测时间的中心像元位置，b表示影像波段，tk表示影像时间，l(xi,yi,tk,b)为tk时刻landsat影像b波段在(xi,yi)像元处的反射率；l(xw/2,yw/2,tp,b)为tk时刻landsat影像b波段搜索窗口中心像元的反射率；σ(b)为波段b反射率值的标准差；m表示估计类别的数量；

13、(1.4)计算相似像元权重：每个相似像元与之对应的高分辨率和低分辨率像元之间光谱相似度计算公式如(2)所示：

14、

15、其中：

16、li＝{l(xi，yi，tm，b1)，...，l(xi，yi，tm，bn)，l(xi，yi，tn，b1)，...，l(xi，yi，tn，bn)}  (12)

17、mi＝{m(xi，yi，tm，b1)，...，m(xi，yi，tm，bn)，m(xi，yi，tn，b1)，...，m(xi，yi，tn，bn)}  (13)

18、ri是描述相似像元i的高分辨率像元和低分辨率像元之间的光谱相关系数；li和mi分别表示高空间分辨率和低空间分辨率数据在tm和tn时间段内每个波段相似像元的反射率集合；e()表示期望值；d(li)和d(mi)分别是li和mi的方差，di是相似像元i的地理距离；

19、第i个相似像元与中心像元之间的地理距离di如公式(5)所示：

20、

21、结合像元i的光谱相似度和地理距离，计算出指数di：

22、di＝(1-ri)×di  (15)

23、根据di值较大的相似像元对中心像元的算计结果贡献较小，因此相似像元权重wi的计算公式如(7)所示，即计算di的归一化倒数，wi的范围为0-1，并且所有相似像元的总权重为1：

24、

25、时间权重tk依据modis影像反射率在时间tk(k＝m，n)和预测时间tp之间的差异计算，如公式(10)所示。

26、

27、(1.5)计算转换系数：转换系数的计算公式如(9)所示：

28、

29、v(x，y)为高分辨率影像与低分辨率影像的反射率变化的比值，tm和tn为两个不同时间，(x，y)表示像元位置，b表示波段，l(x，y，tm，b)和l(x，y，tn，b)分别表示landsat数据在tm和tn时间b波段(x，y)像元位置处的像元值，m(x，y，tm，b)和m(x,y，tn，b)分别表示modis数据在tm和tn时间b波段(x，y)像元位置处的像元值；

30、(1.6)计算时间权重：计算出相似像元权重和转换系数后，利用tm和tn两个时间点的modis数据以及预测时间tp时modis数据分别代入公式(10)，得到时间加权，为精细分辨率反射率设置更大的时间权重，以确保离预测日期更近的精细分辨率数据具有相近的反射率值：

31、

32、式中：m和n表示不同的时间；w表示搜索窗口大小，(xi，yi)表示像元位置，b表示影像的波段，m(xi，yi，tk，b)表示tk时间b波段(xi，yi)像元位置的modis数据，m(xi，yi，tp，b)表示预测时间tpb波段(xi，yi)像元位置的modis数据，tk表示时间权重；

33、(1.7)计算预测期中心像元值：根据时间权重tm和tn，高空间分辨率预测期tp中心像元值计算公式如下：

34、l(xw/2，yw/2，tp，b)＝tm×lm(xw/2，yw/2，tp，b)+tn×ln(xw/2，yw/2，tp，b)  (18)

35、式中：tm和tn为时间权重；w表示搜索窗口大小，(xw/2，yw/2)表示中心像素，tp表示预测时间，b表示预测波段，lm(xw/2，yw/2，tp，b)和ln(xw/2，yw/2，tp，b)分别为以高分辨率影像b波段的中心像元tm和tn时刻为基准日计算得到的预测日期的结果；

36、(1.8)根据获取到的预测日期的结果，即estarfm算法融合得到的30m空间分辨率的多光谱数据，获取精细尺度植被指数ndvifused，具体计算公式如下：

37、ndvifused＝(nir-r)/(nir+r)

38、式中：nir表示融合得到的30m空间分辨率多光谱数据的近红外波段，r表示融合得到的30m空间分辨率多光谱数据的红波段。

39、进一步的，步骤(2)中与土壤水分数据时间滞后1个月左右的植被记忆数据选择modis13a2，该数据为月尺度ndvi数据，预测的smap土壤水分数据b月，选择的modis13a2数据为相同年份的b+1月数据，即ndvilagged。

40、进一步的，步骤(3)中的smap数据选择为9km每日spl3smp_e数据，具体选择土壤水分和亮温数据，亮温数据包括水平极化亮温和垂直极化亮温；

41、modis数据包括植被、温度、蒸散、初级生产力、地表反照率、土地覆盖类型，植被数据包括归一化植被指数(ndvi,mod13a3)、增强植被指数(evi,mod13a2)和叶面积指数(lai,mcd15a3)，温度数据为地表温度(lst,mod11a1)，蒸散发数据包括蒸散发(et,mod16a2)和潜在蒸散(pet,mod16a2)，初级生产力数据为总初级生产力(gpp,mod17a2)，地表反照率数据为(albedo,mcd43a3)，土地覆盖类型为(landcover,mcd12q1)；

42、srtm数据包括高程，通过高程数据在google earth engine平台计算坡度、坡向和山体阴影；

43、soilgrids数据选择0-5厘米深度250米空间分辨率的平均土壤砂粒含量、平均粉粒含量、平均粘粒含量和土壤平均ph值；

44、era5-land数据在google earth engine(gee)的公共数据档案中下载空间分辨率为11132米的era5-land小时数据集，具体获取数据集中0-7厘米表层土壤水分；

45、landsat数据利用gee云平台获取，数据集为landsat8 oli/tirs传感器，包含经过正射校正处理为表面反射率的6个可见和近红外波段和2个短波红外波段，处理为正射校正亮温数据的2个热红外波段和1个全色波段，所有的数据已在gee平台进行辐射校正、地形校正、几何校正预处理工作，并利用lasrc方法进行了大气校正，包括使用cfmask算法生成的云、阴影、水和雪的掩膜以及每个像素饱和度掩模，投影转化为epsg:32649，输出格式为geotiff。

46、进一步的，步骤(3)中对modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据进行预处理具体包括：对smap数据、modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据在gee平台进行投影转换、裁剪和重采样，确保上述数据具有相同的投影方式、相同的影像覆盖范围和与smap土壤水分具有一致的空间分辨率。

47、进一步的，其中投影转换使用的投影坐标系为wgs84地理坐标，裁剪用到研究区的矢量数据，重采样方式为最邻近法。

48、进一步的，步骤(3)中与smap数据进行时空匹配以作为土壤水分降尺度模型的输入数据，具体包括：

49、smap数据的水平极化亮温(tbh)和垂直极化亮温(tbv)时间分辨率与smap土壤水分一致，改数据只进行空间降尺度；

50、mod13a3时间分辨率为月，考虑土壤水分与植被指数之间滞后约1个月时间，选择土壤水分数据月份的下一个月植被指数，作为植被滞后数据；

51、mcd15a3时间分辨率为4天，选择土壤水分当日前4天或后4天中最近的数据，优先选择后4天数据；

52、mod13a2时间分辨率为16天，选择土壤水分当日前16天或后16天中最近的数据，优先选择后16天数据；

53、mcd15a3时间分辨率为4天，选择土壤水分当日前4天或后4天中最近的数据，优先选择后4天数据；

54、mod11a1时间分辨率为每日，与smap土壤水分一致，改数据只进行空间降尺度；

55、mod16a2时间分辨率为8天，选择土壤水分当日前8天或后8天中最近的数据，优先选择后8天数据；

56、mod17a2时间分辨率为8天，选择土壤水分当日前8天或后8天中最近的数据，优先选择后8天数据；

57、mcd43a3时间分辨率为每日，与smap土壤水分一致，改数据只进行空间降尺度；

58、srtm数据为地形数据，考虑到地形的稳定性，所有土壤水分降尺度数据均选择相同的srtm数据；

59、soilgrids数据为土壤质地数据，考虑到地形的稳定性，所有土壤水分降尺度数据均选择相同的soilgrids数据；

60、era5-land数据粗尺度土壤水分数据，时间分辨率为小时，选择与smap土壤水分过境时间一致的上午9点数据；

61、modis和landsat数据经过estarfm时空融合算法，获取30米反射率数据，与landsat数据一起计算ndvi，构成精细尺度植被指数，时间分辨率小于16天，选择土壤水分当日前16天或后16天中最近的数据，优先选择后16天数据。

62、进一步的，所述步骤(4)具体包括：

63、(4.1)利用时空匹配后的与smap土壤水分具有相同空间分辨率的smap数据、modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据作为输入数据，采用随机森林作为降尺度的训练模型，输出数据为低空间分辨率的土壤水分，获得土壤水分降尺度模型；

64、(4.2)将时空匹配后的原始smap数据、modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据进行重采样至1km空间分辨率，获取高空间分辨率的对应数据；

65、(4.3)将重采样后的高空间分辨率smap数据、modis数据、srtm数据、soilgrids数据、era5-land数据和landsat数据输入土壤水分降尺度模型中，以获取高空间分辨率的土壤水分数据。

66、进一步的，所述土壤水分降尺度模型为：

67、

68、其中，sm为土壤水分，rf为基于随机森林的土壤水分降尺度模型，ndvifused为estarfm算法融合的ndvi数据，ndvilagged为与土壤水分降尺度时间滞后1个月时间尺度的ndvi数据，evi为增强植被指数，lai为叶面积指数，lst为地表温度数据，et为蒸散发数据，pet为潜在蒸散发数据，gpp为总初级生产力，albedo为地表反照率，landcover为土地覆盖类型，elevation为高程数据，slope为坡度数据，aspect为坡向数据，hillshade为山体阴影数据，sand为土壤砂粒含量百分比、silt为土壤粉粒含量百分比、clay为土壤黏粒含量百分比，ph为土壤水分中的ph值，era5选用表层土壤水分0-7厘米数据，tbh为smap水平极化亮温，tbv为smap垂直极化亮温；其中ndviiused为estarfm算法融合得到的30m空间分辨率的多光谱数据，光谱波段包括计算ndvi所用的近红外波段nir和红光波段r，ndvifused具体计算公式如下：

69、ndvifused＝(nir-r)/(nir+r)。

70、总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

71、1、构建了一种利用多源遥感数据联合多特征的土壤水分降尺度模型，实现了低空间分辨率土壤水分数据的空间降尺度：该模型考虑到光学传感器受到云污染影响而导致可用的数据量减少，为了提供更多有关土壤水分的精细尺度信息，本发明借助时空融合算法estarfm生成了精细尺度植被指数ndvifused数据，弥补研究区可用数据的空间和时间缺失数据；考虑到植被对土壤水分存在一定的时间滞后性，当前土壤水分的状态需要在一定时间段后才能响应在植被上，通常比土壤水分滞后一个月，本发明设计了与土壤水分降尺度数据对应时间滞后月尺度的ndvilagged数据，解决了土壤水分与植被存在滞后关系从而影响土壤水分降尺度模型的精度；本发明同时采用了土壤水分降尺度常用的植被数据和温度数据，地表植被通过影响太阳入射辐射和土壤蒸散发进而影响土壤湿度；地表温度通过影响土壤水蒸散发进而影响土壤水分，即地表温度越高，土壤水蒸散发越大，地表则表现出越干旱；同时本发明也考虑到了反照率、蒸散、地形、亮温、土壤质地、土地利用、初级生产力因素对降尺度的影响，另外模型加入了粗尺度的土壤水分数据，低空间分辨率的土壤水分可以作为粗尺度下的平均土壤水分含量，促进生成高空间分辨率且可靠的土壤水分；

72、2、本发明利用降尺度后的高空间分辨率土壤水分数据，实现了区域尺度的土壤水分干旱监测：土壤水分通常从p波段、l波段、c波段和x波段数据中反演，p波段土壤水分的研究较少，仅有个别基于遥感建模与机载飞行观测的研究，l波段比c波段和x波段对土壤水分更敏感性，l波段的土壤水分数据应用更为广泛。常用的l波段的土壤水分数据包括smap和smos数据，两种数据的空间分辨率分别在9km和25km以上。低空间分辨率的遥感像元中包含了异质的土地覆盖类型和复杂的地形，无法提供土壤水分的详细变化，难以准确的反映局部区域尺度的干旱信息。土壤水分对农作物和自然植被的生长有着决定性的作用，原因在于土壤水分的动态决定了农业生态系统中可用的水资源。本发明获取到的高空间分辨率土壤水分数据较为详细地展示了区域尺度的土壤水分空间分布情况，可以提供每日连续时间尺度的土壤水分状况，为区域尺度尤其是复杂下垫面干旱监测提供了精细化的土壤水分数据。