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

一种考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法，包括：1)生成30m空间分辨率的精细尺度植被指数；2)获取与土壤水分数据具有滞后性的NDVI数据，设计植被记忆数据；3)获取目标区域的MODIS数据、SRTM数据、SoilGrids数据、ERA5

【技术实现步骤摘要】
一种考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法


[0001]本专利技术涉及土壤水分降尺度领域，具体是一种考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法。

技术介绍

[0002]土壤水分在水资源管理、水文循环和干旱发展等过程中扮演着重要的角色，作为地表与大气之间水热交换的重要变量，对指导农业生产和保证粮食安全至关重要。对于具有复杂的地形、多变的土地覆盖类型、异质的土壤类型区域，迫切需要高时空分辨率的土壤水分数据提升干旱监测的精度和能力。但是当前全球的土壤水分产品通常空间分辨率较低，无法反映多种地表覆盖、不同气候条件和复杂地形共同作用下所导致的土壤水分较大的空间异质性。因此，如何充分挖掘土壤水分数据与现有辅助数据之前复杂的关系，获得高时空分辨率的土壤水分数据，是一个非常重要的问题。
[0003]土壤水分的产品众多，其中SMAP(Soil Moisture Active and Passive,土壤水分主被动探测计划)同时使用L波段雷达和辐射计开展并行、同步测量，结合了主动(雷达)和被动(辐射计)微波遥感的优势，提供一种能够同时满足高时间分辨率、宽空间覆盖、最佳传感深度以及在中等植被条件下能够高精度的反演土壤水分的要求，被广泛应用于全球土壤水分监测研究中。SMAP项目生产三种土壤水分产品，产品名称分别为SM_P、SM_A和SM_AP，空间分辨率分别为36km、3km和9km，其中由于L波段雷达在2015年7月失效，目前只能获取36km和9km空间分辨率的土壤水分数据，较低空间分辨率的遥感像元中包含了异质的土地覆盖类型和复杂的地形，无法提供土壤水分的详细变化，难以准确的反映局部区域尺度的干旱信息。
[0004]土壤水分的空间分布与气候因子、地形地貌、土壤性质和植被等众多因素相互作用。植被通过影响太阳入射辐射和土壤蒸散发进而影响土壤湿度；地表温度通过影响土壤水蒸散发进而影响土壤水分，即地表温度越高，土壤水蒸散发越大，地表则表现出越干旱；高程通过影响地表温度从而影响土壤湿度，且高程将会影响径流量大小以及降雨的再分配。土壤水分时空分布的差异性是由多种因素共同作用的结果，上述这些因子应同时应用于土壤水分降尺度模型中。同时，光学传感器经常受到云污染影响，导致研究区可用的遥感影像数据量往往较少，但是土壤水分降尺度需要多光谱数据的参与，如何弥补研究区可用数据的空间和时间缺失数据是需要重点关注的问题。此外，现有的土壤水分空间分辨率较低，如何提供更多有关土壤水分的精细尺度信息，提升土壤水分的空间分辨率，也是降尺度研究需要考虑的问题。另一方面，植被与土壤水分之间存在一定的时间滞后性，当前土壤水分的状态需要在一定时间段后才能响应在植被上，通常比土壤水分滞后一个月的时间，现有的土壤水分降尺度模型较少考虑植被后滞性对土壤水分的影响。

技术实现思路

[0005]本专利技术针对现有技术较少考虑土壤水分与植被之间的滞后性、土壤水分与精细尺度植被信息之间的联系、土壤水分数据与多特征辅助数据之间的关系以及弥补光学数据因云雾影响导致无法有效应用于土壤水分降尺度模型中的问题，提出一种考虑植被记忆和精细尺度植被数据的联合多特征的土壤水分降尺度方法，为土壤水分降尺度方法提供更多有关土壤水分的精细尺度信息和植被记忆信息，同时弥补了研究区可用多光谱数据空间和时间的缺失，充分挖掘土壤水分与多特征数据之间的非线性关系，获取了1km空间分辨率的土壤水分数据。
[0006]为实现上述目的，按照本专利技术的一个方面，提供了一种联合多特征的土壤水分降尺度方法，包括以下步骤：
[0007](1)采用增强时空自适应融合算法生成30m空间分辨率的精细尺度植被指数；
[0008](2)利用MOD13A3月尺度NDVI获取与土壤水分数据具有滞后性的NDVI数据，设计土壤水分降尺度辅助因子之一的植被记忆数据，确保植被记忆数据与土壤水分数据时间滞后1个月左右；
[0009](3)获取目标区域的SMAP数据、MODIS数据、SRTM数据、SoilGrids数据、ERA5
‑
Land数据和Landsat数据，对MODIS数据、SRTM数据、SoilGrids数据、ERA5
‑
Land数据和Landsat数据进行预处理使具有相同的投影方式、相同的影像覆盖范围和与SMAP土壤水分具有一致的空间分辨率，然后与SMAP数据进行时空匹配以作为土壤水分降尺度模型的输入数据；
[0010](4)联合植被记忆数据、精细尺度植被指数、SMAP数据、MODIS数据、SRTM数据、SoilGrids数据、ERA5
‑
Land数据和Landsat数据，构建土壤水分降尺度模型以获取1km的高空间分辨率土壤水分产品。
[0011]进一步的，步骤(1)包括：
[0012](1.1)准备数据：为预测日期提供两对与预测日期相近，且确保两对数据具有相同的时间和覆盖范围的MODIS和Landsat反射率影像，以及一组用于被预测的不同时间序列的MODIS影像；
[0013](1.2)数据预处理：将所有的MODIS数据利用ENVI软件进行重采样，重采样方法选择双线性内插以减少地理参考误差的影响，对重采样后的MODIS进行地理校正以获得和Landsat相同的地理配准，之后对所有的数据进行裁剪操作以确保影像具有完全相同的覆盖区域，所有数据都经过预处理以获得表面反射率，预处理后的MODIS和Landsat影像具有相同的投影、空间分辨率和一致的覆盖范围；
[0014](1.3)搜索邻近像元：使用两个高分辨率图像搜索预测期影像局部窗口中心像元的相似像元，相似像元的获取采用滑动窗口的方法获取，相邻像元值需要满足式(1)的条件，即两个像元的标准差较小，以确保在搜索窗口范围内获得与中心像元具有较高光谱相似特性的像元：
[0015]|L(x
i
,y
i
,t
k
,B)
‑
L(x
w/2
,y
w/2
,t
k
,B)|≤σ(B)
·
2/m
ꢀꢀ
(10)
[0016]式中：L为Landsat影像，(x
i
,y
i
)是第i个相似像元的位置；(x
w/2
,y
w/2
)为预测时间的中心像元位置，B表示影像波段，t
k
表示影像时间，L(x
i
,y
i
,t
k
,B)为t
k
时刻Landsat影像B波段在(x
i
,y
i
)像元处的反射率；L(x
w/2
,y
w/2
,t
p
,B)为t
k
时刻Landsat影像B波段搜索窗口中心像元的反射率；σ(B)为波段B反射率值的标准差；m表示估计类别的数量；
[0017](1本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法，其特征在于：包括如下步骤：(1)采用增强时空自适应融合算法生成30m空间分辨率的精细尺度植被指数；(2)利用MOD13A3月尺度NDVI获取与土壤水分数据具有滞后性的NDVI数据，设计土壤水分降尺度辅助因子之一的植被记忆数据，确保植被记忆数据与土壤水分数据时间滞后1个月左右；(3)获取目标区域的SMAP数据、MODIS数据、SRTM数据、SoilGrids数据、ERA5
‑
Land数据和Landsat数据，对MODIS数据、SRTM数据、SoilGrids数据、ERA5
‑
Land数据和Landsat数据进行预处理使具有相同的投影方式、相同的影像覆盖范围和与SMAP土壤水分具有一致的空间分辨率，然后与SMAP数据进行时空匹配以作为土壤水分降尺度模型的输入数据；(4)联合植被记忆数据、精细尺度植被指数、SMAP数据、MODIS数据、SRTM数据、SoilGrids数据、ERA5
‑
Land数据和Landsat数据，构建土壤水分降尺度模型以获取1km的高空间分辨率土壤水分产品。2.如权利要求1所述的考虑植被记忆和精细尺度植被指数的土壤水分降尺度方法，其特征在于：步骤(1)包括：(1.1)准备数据：为预测日期提供两对与预测日期相近，且确保两对数据具有相同的时间和覆盖范围的MODiS和Landsat反射率影像，以及一组用于被预测的不同时间序列的MODIS影像；(1.2)数据预处理：将所有的MODIS数据利用ENVI软件进行重采样，重采样方法选择双线性内插以减少地理参考误差的影响，对重采样后的MODIS进行地理校正以获得和Landsat相同的地理配准，之后对所有的数据进行裁剪操作以确保影像具有完全相同的覆盖区域，所有数据都经过预处理以获得表面反射率，预处理后的MODIS和Landsat影像具有相同的投影、空间分辨率和一致的覆盖范围；(1.3)搜索邻近像元：使用两个高分辨率图像搜索预测期影像局部窗口中心像元的相似像元，相似像元的获取采用滑动窗口的方法获取，相邻像元值需要满足式(1)的条件，即两个像元的标准差较小，以确保在搜索窗口范围内获得与中心像元具有较高光谱相似特性的像元：|L(x
i
，y
i
，t
k
，B)
‑
L(x
w/2
，y
w/2
，t
k
，B)|≤σ(B)
·
2/m
ꢀꢀꢀꢀ
(1)式中：L为Landsat影像，(x
i
，y
i
)是第i个相似像元的位置；(x
w/2
，y
w/2
)为预测时间的中心像元位置，B表示影像波段，t
k
表示影像时间，L(x
i
，y
i
，t
k
，S)为t
k
时刻Landsat影像B波段在(x
i
，y
i
)像元处的反射率；L(x
w/2
，y
w/2
，t
p
，S)为t
k
时刻Landsat影像B波段搜索窗口中心像元的反射率；σ(B)为波段B反射率值的标准差；m表示估计类别的数量；(1.4)计算相似像元权重：每个相似像元与之对应的高分辨率和低分辨率像元之间光谱相似度计算公式如(2)所示：其中：L
i
＝{L(x
i
，y
i
，t
m
，B1)，...，L(x
i
，y
i
，t
m
，B
n
)，L(x
i
，y
i
，t
n
，B1)，...，L(x
i
，y
i
，t
n
，B
n
)}
ꢀꢀꢀꢀ
(3)
M
i
＝{M(x
i
，y
i
，t
m
，B1)，...，M(x
i
，y
i
，t
n
，B
n
)，M(x
i
，y
i
，t
n
，B1)，...，M(x
i
，y
i
，t
n
，B
n
)}
ꢀꢀꢀꢀ
(4)R
i
是描述相似像元i的高分辨率像元和低分辨率像元之间的光谱相关系数；L
i
和M
i
分别表示高空间分辨率和低空间分辨率数据在t
m
和t
n
时间段内每个波段相似像元的反射率集合；E()表示期望值；D(L
i
)和D(M
i
)分别是L
i
和M
i
的方差，di是相似像元i的地理距离；第i个相似像元与中心像元之间的地理距离d
i
如公式(5)所示：结合像元i的光谱相似度和地理距离，计算出指数D
i
：D
i
＝(1
‑
R
i
)
×
d
i
ꢀꢀꢀꢀꢀ
(6)根据D
i
值较大的相似像元对中心像元的算计结果贡献较小，因此相似像元权重W
i
的计算公式如(7)所示，即计算D
i
的归一化倒数，W
i
的范围为0
‑
1，并且所有相似像元的总权重为1：时间权重T
k
依据MODIS影像反射率在时间t
k
(k＝m，n)和预测时间t
p
之间的差异计算，如公式(10)所示：(1.5)计算转换系数：转换系数的计算公式如(9)所示：v(x，y)为高分辨率影像与低分辨率影像的反射率变化的比值，t
m
和t
n
为两个不同时间，(x，y)表示像元位置，B表示波段，L(x，y，t
m
，B)和L(x，y，t
n
，B)分别表示Landsat数据在t
m
和t
n
时间B波段(x，y)像元位置处的像元值，M(x，y，t
m
，B)和M(x，y，t
n
，B)分别表示MODIS数据在t
m
和t
n
时间B波段(x，y)像元位置处的像元值；(1.6)计算时间权重：计算出相似像元权重和转换系数后，利用t
m
和t
n
两个时间点的MODIS数据以及预测时间t
p
时MODIS数据分别代入公式(10)，得到时间加权，为精细分辨率反射率设置更大的时间权重，以确保离预测日期更近的精细分辨率数据具有相近的反射率值：式中：m和n表示不同的时间；W表示搜索窗口大小，(x
i
，y
i
)表示像元位置，B表示影像的波段，M(x
i
，y
i
，t
k
，B)表示t
k
时间B波段(x
i
，y
i
)像元位置的MODIS数据，M(x
i
，y
i
，t
p
，B)表示预测时间t
p
B波段(x
i
，y
i
)像元位置的MODIS数据，T
k
表示时间权重；
(1.7)计算预测期中心像元值：根据时间权重T
m
和T
n
，高空间分辨率预测期t
p
中心像元值计算公式如下：L(x
w/2
，y
w/2
，t
p
，B)＝T
m
×
L
m
(x
w/2
，y<...

【专利技术属性】
技术研发人员：崔长露，向大享，李喆，陈喆，姜莹，陈希炽，赵静，吴仪邦，文雄飞，李经纬，
申请(专利权)人：长江水利委员会长江科学院，
类型：发明
国别省市：