基于静止气象卫星的植被冠层可燃物含水率反演方法及系统

1.本发明涉及遥感反演技术领域，特别是一种基于静止气象卫星的植被冠层可燃物含水率反演方法及系统。


背景技术：

2.野火是生态系统中一种复杂的自然干扰过程，不仅破坏生态环境，释放温室气体，影响全球气候，还威胁人民生命财产安全。野火风险评估与预警作为森林和草原管理的一项重要举措，对于火灾的有效预防和管控工作具有重要意义。
3.可燃物含水率(fuel moisture content,fmc)定义为植被冠层含水量与其干重的比率，与火灾的发生及蔓延过程密切相关，是评估野火风险的关键指标因子。因此，准确、及时获取fmc对于野火风险评估与预警至关重要。卫星遥感技术能够进行长时间序列、高分辨率的地表观测，为大范围、近实时的fmc监测提供了有效途径。目前，fmc反演通常使用太阳同步轨道卫星获取的光学遥感影像，包括landsat系列卫星、sentinel-2系列卫星、中等分辨率成像光谱仪modis产品等。这些光学遥感影像虽然具有较高空间分辨率，但其时间分辨率较低，导致反演获取的fmc时效性较差，无法满足近实时的野火风险评估与预警需求。相比之下，静止气象卫星himawari-8能够以10分钟的重访周期进行对地观测，可形成丰富的遥感观测数据，在近实时fmc获取及野火风险评估与预警中具有较大的应用潜力。
4.fmc反演研究可采用经验统计和基于物理模型的方法，前者严重依赖于实测数据，难以大范围应用，后者考虑植被辐射传输机制，具有较高的普适性。然而，目前针对基于静止卫星数据的fmc反演研究鲜少，如何基于静止气象卫星的多时相、多角度观测数据反演fmc，实现近实时获取fmc，是进一步提高野火风险评估与预警时效性关键。


技术实现要素：

5.有鉴于此，本发明的目的在于提供基于静止气象卫星的植被冠层可燃物含水率反演方法，该方法用于近实时估算植被冠层可燃物含水率。
6.为达到上述目的，本发明提供如下技术方案：
7.本发明提供的基于静止气象卫星的植被冠层可燃物含水率反演方法，包括以下步骤：
8.步骤1：获取静止气象卫星观测数据，生成用于植被指数计算的地表反射率数据；
9.步骤2：构建植被辐射传输模型；构建关于归一化差分植被指数ndvi、增强植被指数evi、归一化差分红外指数ndii的代价函数t_rmse；
10.步骤3：获取代价函数t_rmse关于自由变量的偏导数，用于计算梯度信息；
11.步骤4：利用既约梯度算法分两阶段迭代求解fmc，在第一阶段反演敏感参数等水分厚度ewt、叶面积指数lai和土壤因子psoil，在第二阶段反演弱敏感参数可燃物含水率fmc和叶绿素a和b含量c
ab
；
12.步骤5：针对不同初值组合下使用两阶段既约梯度算法迭代求解出的所有可燃物
含水率fmc结果。
13.进一步，所述归一化差分植被指数ndvi、增强植被指数evi、归一化差分红外指数ndii及其代价函数t_rmse，按照以下公式进行计算：
[0014][0015][0016][0017][0018]
其中，ρ
blue
，ρ
red
，ρ
nir
，ρ
swir2
分别是蓝光、红光、近红外和短波红外波段的光谱反射率；ndvio，evio，ndiio为从卫星数据提取的观测植被指数；ndvim，evim，ndiim为prosail模型模拟得到的植被指数，也是关于自由变量c
ab
，ewt，fmc、lai和psoil的函数，t_rmse表示当代价函数。
[0019]
进一步，所述步骤3中的自由变量的偏导数是利用自动微分工具tapenade进行计算的。
[0020]
进一步，所述步骤4中的利用既约梯度算法分两阶段迭代求解fmc，具体步骤如下：
[0021]
步骤4.1：设置敏感参数ewt、lai和psoil初值；
[0022]
步骤4.2：利用既约梯度算法迭代求解敏感参数ewt、lai和psoil，当代价函数t_rmse1满足公式(5)时，停止第一阶段迭代，同时输出ewt、lai和psoil，以及fmc和c
ab
，用于第二阶段数值迭代的初值输入；
[0023][0024]
其中，所述代价函数t_rmse1以公式(4)为原型，且固定参数fmc和c
ab
，仅将ewt、lai和psoil作为自由变量；
▽
t_rmse1为t_rmse1的梯度，由步骤3计算得出；ε1为允许的误差限；
[0025]
步骤4.3：引入基于fmc与lai和ewt的生态学规则，通过公式(6)动态更新fmc取值区间；
[0026][0027]
其中，g1(lai)为fmc关于lai的线性拟合方程，rmse
fmc-lai
为拟合方程的均方根误差；g2(ewt)为fmc关于ewt的线性拟合方程，rmse
fmc-ewt
为拟合方程的均方根误差；
[0028]
步骤4.4：利用既约梯度算法迭代求解弱敏感参数fmc和c
ab
；当代价函数t_rmse1不满足公式(5)时，转入步骤4.2；当代价函数t_rmse2满足公式(11)时，停止第二阶段迭代，输出fmc和c
ab
的迭代结果；
[0029][0030]
其中，所述代价函数t_rmse2以公式(4)为原型，且固定参数ewt、lai和psoil，仅将fmc和c
ab
作为自由变量；
▽
t_rmse2为t_rmse2的梯度，由步骤3计算得出；ε2为允许的误差限。
[0031]
进一步，所述步骤5中的针对不同初值组合下使用两阶段既约梯度算法迭代求解
出的所有可燃物含水率fmc结果，具体步骤如下：
[0032]
按照与其对应的代价函数函数值t_rmse及其梯度的模|
▽
t_rmse|进行主排序和次排序，首先选出代价函数最小的m组可燃物含水率fmc反演结果，然后再从其中选出梯度的模最小的n组可燃物含水率fmc反演结果，取这n组可燃物含水率fmc的均值作为最终反演结果。
[0033]
进一步，所述ε1的取值区间为0.1-1；或所述ε2的取值区间为10-6-10-5
；或所述m的取值为15；或所述n的取值为9。
[0034]
进一步，步骤1中的所述卫星观测数据包括全圆盘观测模式下的标准反照率产品和太阳天顶角及方位角、卫星天顶角及方位角的几何参数产品；并完成对数据的云掩模、大气校正。
[0035]
进一步，所述步骤2中的植被辐射传输模型按照以下方式进行：
[0036]
选取叶绿素a和b含量c
ab
、等水分厚度ewt、可燃物含水率fmc、叶面积指数lai和土壤因子psoil作为prosail植被辐射传输模型的自由变量，模型的其它参数为常量。
[0037]
本发明提供的基于静止气象卫星的植被冠层可燃物含水率反演系统，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述方法。
[0038]
本发明的有益效果在于：
[0039]
本发明提供的基于静止气象卫星的可燃物含水率反演方法，涉及遥感反演技术领域。本发明针对静止气象卫星himawari-8数据，基于prosail植被辐射传输模型，首先利用植被指数构建以均方根误差为原型的代价函数，然后利用自动微分工具计算代价函数的梯度信息，再基于两阶段的数值迭代过程快速反演可燃物含水率，最后结合多迭代初值选取策略，缓解弱敏感参数fmc的病态反演问题。
[0040]
本方法操作简单，仅需输入遥感植被指数，即可快速反演可燃物含水率，且能考虑太阳天顶角及方位角等几何参数的动态变化，适用于基于静止气象卫星多时相、多角度观测数据的fmc反演，这对于近实时的森林火灾风险管控具有十分重要的意义。
[0041]
本方法操作简单，仅需输入静止气象卫星观测的植被指数和几何参数信息，即可基于两阶段的迭代反演过程快速估测可燃物含水率。该方法可考虑太阳天顶角及方位角等几何参数的动态变化，较适用于基于静止气象卫星的fmc反演，对于近实时的森林火灾风险管控具有十分重要的意义。
[0042]
本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。
附图说明
[0043]
为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：
[0044]
图1为本发明的整体方法流程示意图。
[0045]
图2为本发明具体实施方案中的研究地块地理位置示意图。
[0046]
图3为本发明具体实施方案中反演中国和澳大利亚可燃物含水率的结果。
具体实施方式
[0047]
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
[0048]
实施例1
[0049]
如图1所示，本实施例提供基于静止气象卫星的植被冠层可燃物含水率反演方法，包括以下步骤：
[0050]
步骤1：获取静止气象卫星观测数据，包括全圆盘观测模式下的标准反照率产品和太阳天顶角及方位角、卫星天顶角及方位角等几何参数产品；并完成对数据的云掩模、大气校正，生成用于植被指数计算的地表反射率数据；
[0051]
步骤2：构建植被辐射传输模型和代价函数t_rmse，选取叶绿素a和b含量c
ab
、等水分厚度ewt、可燃物含水率fmc、叶面积指数lai和土壤因子psoil作为prosail植被辐射传输模型的自由变量，模型的其它参数为常量；构建关于归一化差分植被指数ndvi(公式1)、增强植被指数evi(公式2)、归一化差分红外指数ndii(公式3)的代价函数t_rmse，其计算方法如公式(4)所示。
[0052][0053][0054][0055][0056]
其中，ρ
blue
，ρ
red
，ρ
nir
，ρ
swir2
分别是蓝光、红光、近红外和短波红外波段的光谱反射率；ndvio，evio，ndiio为从卫星数据提取的观测植被指数；ndvim，evim，ndiim为prosail模型模拟得到的植被指数，也是关于自由变量c
ab
，ewt，fmc、lai和psoil的函数，t_rmse表示当代价函数；
[0057]
步骤3：利用自动微分工具tapenade获取代价函数t_rmse关于上述5个自由变量的偏导数，用于计算梯度信息；
[0058]
步骤4：利用既约梯度算法分两阶段迭代求解fmc，在第一阶段反演敏感参数ewt、lai和psoil，在第二阶段反演弱敏感参数fmc和c
ab
；
[0059]
具体地包括4个部分：
[0060]
步骤4.1：敏感参数ewt、lai和psoil初值设定；
[0061]
ewt、lai和psoil初值选取范围通过实测数据确定，均从每个变量的最小值、均值、最大值中选取，故总共拥有3
×3×
3＝27种初值组合；
[0062]
步骤4.2：利用既约梯度算法迭代求解敏感参数ewt、lai和psoil，当代价函数t_rmse1满足公式(5)时，停止第一阶段迭代，同时输出ewt、lai和psoil，以及fmc和c
ab
，用于第二阶段数值迭代的初值输入；
[0063][0064]
其中，所述代价函数t_rmse1以公式(4)为原型，且固定参数fmc和c
ab
，仅将ewt、lai和psoil作为自由变量；为t_rmse1的梯度，由步骤3计算得出；ε1为允许的误差限；
[0065]
所述ε1的取值区间为0.1-1；
[0066]
步骤4.3：引入基于fmc与lai和ewt的生态学规则，通过公式(6)动态更新fmc取值区间；
[0067][0068]
其中，g1(lai)为fmc关于lai的线性拟合方程，rmse
fmc-lai
为拟合方程的均方根误差；g2(ewt)为fmc关于ewt的线性拟合方程，rmse
fmc-ewt
为拟合方程的均方根误差；
[0069]
针对中国实测草地冠层可燃物含水率数据：
[0070][0071]
针对中国实测森林冠层可燃物含水率数据：
[0072][0073]
针对澳大利亚实测草地冠层可燃物含水率数据：
[0074][0075]
针对澳大利亚实测森林冠层可燃物含水率数据：
[0076][0077]
其中，为拟合的可燃物含水率，lai为实测的叶面积指数，ewt为实测的等水分厚度；
[0078]
步骤4.4：利用既约梯度算法迭代求解弱敏感参数fmc和c
ab
；当代价函数t_rmse1不满足公式(5)时，转入步骤4.2；当代价函数t_rmse2满足公式(11)时，停止第二阶段迭代，输出fmc和c
ab
的迭代结果；
[0079][0080]
其中，所述代价函数t_rmse2以公式(4)为原型，且固定参数ewt、lai和psoil，仅将fmc和c
ab
作为自由变量；为t_rmse2的梯度，由步骤3计算得出；ε2为允许的误差限；
[0081]
所述ε2的取值区间为10-6-10-5
。
[0082]
步骤5：针对不同初值组合下使用两阶段既约梯度算法迭代求解出的所有fmc结果，按照与其对应的代价函数函数值(t_rmse)及其梯度的模进行主排序和次排序，首先选出代价函数最小的m组fmc反演结果，然后再从其中选出梯度的模最小的n组fmc反演结果，取这n组fmc的均值作为最终反演结果；
[0083]
所述m的取值为15；所述n的取值为9。
[0084]
实施例2
[0085]
如图1所示，下面结合具体实施例和说明书附图对本发明提供的基于静止气象卫星的植被冠层可燃物含水率反演方法作进一步说明，包括以下步骤：
[0086]
步骤1：数据准备
[0087]
中国实测草地冠层可燃物含水率(fmc)数据，于2015年7月28日至8月2日期间在中国青海湖流域草原进行野外调查所得；中国实测森林冠层可燃物含水率数据，于2020年12月5日至10日和2021年3月16日至20日期间在中国四川省西南地区进行野外调查所得。澳大利亚实测森林和草地冠层可燃物含水率数据，于2015年10月1日至2016年4月20日期间在澳大利亚堪培拉进行野外调查测得。如图2所示，即为本实施例中所用实测数据的地理位置分布情况。本实施例中所用himawari-8卫星全圆盘观测模式下的标准反照率产品和太阳天顶角、太阳方位角、卫星天顶角和卫星方位角等几何参数产品。对该数据进行预处理操作，包括云掩模和大气校正，生成地表反射率产品。
[0088]
步骤2：代价函数构建
[0089]
对prosail植被辐射传输模型进行参数化，如表1所示；
[0090]
表1 prosail模型输入参数(g：草地，f：森林)
[0091][0092]
选取叶绿素a和b含量c
ab
(μg/cm2)、等水分厚度ewt(g/cm2)、干物质重量dmc(g/cm2)、叶面积指数lai和土壤因子psoil作为prosail植被辐射传输模型的自由变量，模型的其它参数为常量；其中，可燃物含水率构建关于归一化差分植被指数ndvi(公式1)、增强植被指数evi(公式2)、归一化差分红外指数ndii(公式3)的代价函数t_rmse，其计算方法如公式(4)所示。
[0093]
步骤3：代价函数的偏导数获取
[0094]
利用自动微分工具tapenade(下载地址：https://team.inria.fr/ecuador/en/
tapenade，版本：3.6)获取代价函数t_rmse关于c
ab
，ewt，fmc、lai和psoil的偏导数，用于计算梯度信息；
[0095]
步骤4：数值迭代求解fmc
[0096]
利用既约梯度算法分两阶段迭代求解fmc，在第一阶段反演敏感参数ewt、lai和psoil，在第二阶段反演弱敏感参数fmc和c
ab
；
[0097]
具体地包括4个部分：
[0098]
步骤4.1：敏感参数ewt、lai和psoil初值设定；
[0099]
单组敏感参数初值组合下，迭代点容易陷入局部最优值，这可能导致fmc反演的不确定性，因此需要设定多组初值重复数值寻优过程，有效避免局部最优问题。ewt、lai和psoil初值选取范围通过实测数据确定，均从每个变量的最小值、均值、最大值中选取，故总共拥有3
×3×
3＝27种初值组合；
[0100]
步骤4.2：数值迭代求解敏感参数ewt、lai和psoil；
[0101]
利用既约梯度算法迭代求解敏感参数ewt、lai和psoil。当代价函数t_rmse1满足公式(5)时，停止第一阶段迭代，同时输出ewt、lai和psoil，以及fmc和c
ab
，用于第二阶段数值迭代的初值输入。
[0102]
步骤4.3：引入生态学规则；
[0103]
针对野外实测可燃物含水率数据，拟合实测fmc与实测lai和ewt间的经验方程，结合公式(6)～(10)，构建生态学规则，具体如(12)～(15)所示，并以此动态更新fmc取值区间；
[0104]
针对中国实测草地冠层可燃物含水率反演：
[0105][0106]
针对中国实测森林冠层可燃物含水率反演：
[0107]
|fmc-56.08*lai-61.34||≤14.09％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0108]
针对澳大利亚实测草地冠层可燃物含水率反演：
[0109]
|fmc-108.5*lai-4.52||≤29.95％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0110]
针对澳大利亚实测森林冠层可燃物含水率反演：
[0111]
|fmc-9.25*lai-98.34||≤11.67％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0112]
步骤4.4：数值迭代求解弱敏感参数fmc和c
ab
；
[0113]
利用既约梯度算法迭代求解弱敏感参数fmc和c
ab
；当代价函数t_rmse1不满足公式(5)时，转入步骤4.2；当代价函数t_rmse2满足公式(11)时，停止第二阶段迭代，输出fmc和c
ab
的迭代结果；
[0114]
利用网格化搜索，统计代价函数关于变量的偏导数的区间分布，如表2所示，误差限ε2取值应小于10-5
，为提高反演精度和计算效率，ε2优选为10-6
；
[0115]
表2代价函数关于输入参数的偏导数的值域
[0116][0117]
步骤5：可燃物含水率反演
[0118]
针对不同初值组合下使用两阶段既约梯度算法迭代求解出的所有fmc结果，按照与其对应的代价函数函数值(t_rmse)及其梯度的模进行主排序和次排序，首先选出代价函数最小的m组fmc反演结果，然后再从其中选出梯度的模最小的n组fmc反演结果，取这n组fmc的均值作为最终反演结果；为了缓解病态反演问题，利用网格搜索法进行调参，即分别设置m＝6，9，
…
，27；n＝3，6，
…
，24(m＞n)，通过重复实验确定产生最佳反演结果的参数m和n的值，可知m取值优选18，n取值优选9；
[0119]
根据上述反演方法，求解出可燃物含水率。如图3所示，即为本发明在本实施例中对可燃物含水率的反演结果。
[0120]
以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。