CN113343408A - 基于多源遥感与能量平衡模型造林响应地表温度的方法 - Google Patents

Abstract

本发明公开了一种基于多源遥感与能量平衡模型造林响应地表温度的方法，包括以下步骤：1)收集可用的遥感数据、气候数据，以及土地覆盖数据；2)利用空间分析的方法生成高效分辨率的未变化人工林、天然林以及开阔地的时空分布图；3)将收集到的数据进行尺度转换与匹配，统一分辨率，并基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析；4)构建地表能量平衡方程与内在生物物理理论集成的造林响应地表温度辨别的模型。本申请集成多源遥感数据、气候数据和能量平衡模型评估造林的区域温度响应机制，使得响应温度变化的模型更为全面与可靠，为合理造林的经营举措的实施提供重要科学评价依据。

Description

基于多源遥感与能量平衡模型造林响应地表温度的方法

技术领域

本发明属于森林变化响应温度效应模型技术领域，具体涉及一种基于多源遥感与能量平衡模型造林响应地表温度的方法。

背景技术

如何精确地量化造林时空分布的生物物理机制响应区域温度变化有助于了解森林变化的气候调控作用。传统的基于土地利用变化和野外观测的统计方法通常被用于探测森林响应温度变化。然而由于气象站或者野外布设站点分布较为稀疏，不能准确地定位森林覆盖动态位置的温度时空变化规律。其作为一种最具说服力的方法，特别适用于小区域，但其缺少机理分析。气候模式已被广泛应用于造林的生物物理反馈机制对于气候影响的研究，但全球气候模式更适合于大尺度的研究，因为气候模型的输入数据存在空间分辨率较粗，以及气候模型结构、物理过程和参数设定的差异，不能以能量变化解释的温度响应，会给模型的结果带来很多不确定性。区域气候模式对局地尺度有更好的适应性，但区域模式未能充分认识森林变化的内在地表生物物理效应和大气反馈对地表温度变化的贡献。以遥感观测为基础的研究已经成功的将森林覆盖数据与陆地表面产品结合，实现监测不同尺度的温度变化。将原有的点尺度监测扩展到了更为广阔的时空尺度。另外，基于空间分析与比对的遥感观测与气候模型模拟的结果存在较大差异。无论采用哪种方法，人们逐渐关注其物理机制如何以能量收支和能量平衡解释，以及这些关键生物物理参数和气候参数如何通过遥感手段获取并应用。

基于MODIS陆表数据结合简单建模或者经验数据分析方法是分析省级及以上尺度森林与气候关联采取的主要手段。MODIS数据及其衍生的地面产品种类丰富且全面，有较一致与连续的时空分辨率。但是MODIS土地覆盖产品在区域尺度因精度较低不能很好地反映真实的地表覆盖类型，特别是在具有异质性土地覆被的地区。基于Landsat或者多源数据结合的产品已用于量化森林变化的地表温度变化。但是方法中关于造林时的地表状况改变及其区域与局地反照率辐射反馈、能量再分配等生物物理能量平衡机制以及气候反馈是如何变化的还有待评估。总之，一方面中高分辨率造林数据的有限可用性一直阻碍着对温度变化影响的精准量化，另一方面还不清楚究竟具体哪些生物物理因素导致了区域温度的变化。

常规的观测方法无法实现造林及其区域温度变化响应的时空尺度扩展，一是缺乏精确的地表变化时空信息，二是响应温度变化的观测手段和方法的单一性，比如仅使用遥感数据结合一些空间统计分析量化手段，忽略了很多地表能量平衡的关键参数以及内在的生物物理机制过程对于地面温度的作用，即使使用了一些站点数据结合地表能量平衡模型进行森林响应温度变化的分析，但是站点数据的有限覆盖度以及不能与森林变化的分布很好的契合，这都阻碍了省级以及更高尺度范围内的动态监测与响应。依据现存的很多气候模型与模式用于森林响应温度变化，其中全球气候模式在局地应用中局限性较多，其分辨率太粗不适于局地的估测；区域气候模式需要调整很多参数才能本地化，单独使用这些方法，都凸显出了各自的不足之处，因而很难在局地尺度充分发挥各种观测手段与方法的特色之处。粗分辨率的MODIS土地覆盖数据因数据分辨率低，通常被用于大尺度的研究，对于小区域尺度发生的森林变化无法准确定位。

发明内容

针对现有技术存在的上述缺陷，本发明所要解决的技术问题在于提供一种基于多源遥感与能量平衡模型造林响应地表温度的方法。

为了解决上述技术问题，本发明所采用的技术方案如下：

基于多源遥感与能量平衡模型造林响应地表温度的方法，包括以下步骤：

1)收集可用的遥感数据、气候数据，以及土地覆盖数据；

2)利用空间分析的方法生成高效分辨率的未变化人工林、天然林以及开阔地的时空分布图；

3)将收集到的数据进行尺度转换与匹配，统一分辨率，并基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析；

4)构建地表能量平衡方程与内在生物物理理论集成的造林响应地表温度辨别的模型。

进一步的，步骤1)中，可用的遥感数据、气候数据，以及土地覆盖数据包括：MODIS地表温度MYD11A2数据、MODIS蒸散发MOD16A2数据、MODIS反照率MCD43B3数据、气温数据、MREEA-2再分析数据和高程数据。

进一步的，MODIS反照率MCD43B3数据包括黑空反照率与白空反照率。

进一步的，步骤3)中，在将收集到的数据进行尺度转换与匹配，统一分辨率时，利用双线性内插法将收集到的数据的分辨率升尺度或者降尺度为1km，实现尺度统一。

进一步的，步骤3)中，将收集到的数据进行尺度转化与匹配具体包括：将黑空反照率与白空反照率进行均值处理；对按月份、季节、年份的所有MODIS数据进行求均值处理，生成按月份、年际的地表温度白天值、地表温度夜晚值、平均最高气温、平均最低气温、蒸散发与反照率；将MERRA-2再分析数据用于提取晴空的地表下行长波辐射、下行短波辐射；利用随机森林模型反演站点观测的温度数据，重采样到1km，基于温度数据与高程值的线性回归方程，修正由海拔差异造成的温度差异。

进一步的，步骤3)中，在基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析时，遵循以下准则：人工林覆盖至少在5％以上，天然林或者开阔地覆盖在80％以上。

进一步的，步骤3)中，在基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析时，参照造林的反照率差异表达方程(1)：

Δα＝α_PF-α_OP (1)：

其中Δα表示造林后的反照率空间变化，α_PF是造林后的人工林分布区域的反照率均值，α_OP是造林前的开阔地分布区域的反照率均值。

进一步的，步骤3)中，在基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析时，对造林前后每个土地覆盖类型所在区域计算月与年能量重分配因子f的分布按照如下公式(2)计算：

其中T_s为卫星观测的地表温度；T_a为空气温度；

为净辐射；

SW_net为净地表短波辐射，SW_net＝(1-α)SW↓，单位为W m^-2；LW↓为入射长波辐射，单位为W m^-2；SW↓为入射短波辐射，单位为W m^-2；G为土壤热通量，G＝0.14(T_a，n-T_a，n-1)，n表示月份数1-12；

表示长波辐射反馈的温度敏感性，单位为K(Wm^-2)^-1；其中λ₀、f、

G分别表示造林前农田、草地和天然林对应变量的值，若不考虑大气反馈，则

Δα为反照率变化值；ΔG为土壤热通量扩散变化；Δf是湍流能量重分配的平均变化。

进一步的，在对造林前后每个土地覆盖类型所在区域计算月与年能量重分配因子f的分布按照如下公式(2)计算后，将计算所得的上述值输入地表能量平衡方程构建的地表温度响应模型，模型如公式(3-7)：

该公式依据公式(1)被改进为如下公式(4-7)：

ΔT_s＝ΔT_{s_α}+ΔT_{s_G}+ΔT_{s_f} (4)；

其中ΔT_s为建模地表温度；ΔT_{s_α}为辐射强迫应对反照率变化的地表温度响应；ΔT_{s_G}为土壤热通量扩散的地表温度响应；ΔT_{s_f}为湍流能量再分配的地表温度响应；当ΔT_s为正值表明造林带来增温影响，否者则为降温影响，零值则表明影响不明显。

有益效果：相比于现有技术，本发明的优点为：

(1)本申请将高分辨率表征造林时空模式的土地覆盖数据作为地表能量平衡模型的重要输入数据，不同于以往采用粗分辨率数据导致地表变化监测的较大误差，从而影响造林响应气候变化的判定。

(2)本申请集成多源遥感数据、气候数据和能量平衡模型评估造林的区域温度响应机制，充分发挥了多源数据的特性，使得响应温度变化的模型更为全面与可靠，为合理造林的经营举措的实施提供重要科学评价依据。

附图说明

图1为本申请的技术流程图；

图2为造林前后(未变化CR/GR/NF-＞PF)能量重分配因子f、净辐射

反照率Δα或者ΔAlbedo、下行短波辐射SW_net、蒸散发ET和土壤热通量G的月份差异值；图中，每条竖线代表学生t检验估计的95％的置信区间；图中，PF-CR、PF-GR、PF-NF分别表示农田、草地与天然林转为人工林后对应值的前后差异；

图3为造林前后ΔAlbedo与

的关联关系图；图中，直线为线性回归线；

图4为造林的建模地表温度数据(ΔT_{s_m})、f驱动地表温度数据(ΔT_{s_f})、反照率驱动地表温度数据(ΔT_{s_α})、土壤热通量驱动地表温度数据(ΔT_{s_G})；图中，每条竖线代表学生t检验估计的95％的置信区间；

图5为造林的建模地表温度数据(ΔT_{s_m})与Δf的关联关系图；图中，直线为线性回归线；

图6为造林的建模地表温度数据(ΔT_{s_m})与ΔT_{s_f}的关联关系图；图中，直线为线性回归线。

具体实施方式

下面结合具体实施例对本发明进一步进行描述。

实施例1

如图1所示的基于多源遥感与能量平衡模型的人工林响应地表温度变化的方法，具体包括如下步骤：首先收集可用的遥感与气候数据，以及土地覆盖数据，利用空间分析的方法生成高分辨率的未变化造林时空分布图，包括人工林、天然林与开阔地(草地与农田)，然后将获取到的所有数据进行尺度转换与匹配，统一分辨率，并通过移动窗口空间采样方法，实现基于辐射、土壤热通量、卫星地表温度与空气温度的能量重分配因子估测，最后，基于造林引起的能量重分配因子、反照率、土壤热通量的变化，计算出各因子对地表温度的影响，将三者求和即为造林对地表温度的影响。

1)数据获取

以广东省为例，获取30m分辨率的2000年和2010年的Shen等(2019)中的“SGB-NDVI森林非森林产品”以及国家测绘局的GlobeLand30数据用于提取土地覆盖分布用于识别造林的时空分布。此外，还获取地表温度数据，具体参数：每8天一期的MODIS MYD11A2数据产品，时间跨度为2002-2010年；反照率数据来源于MODLAND团队，具体参数：每8天一期的MCD43B3，时间跨度为2001-2010年，包括黑空反照率和白空反照率；ET数据来源蒙大拿大学数字地形动力学模拟组，具体参数：每8天一期的MOD16A2数据产品，所有数据的分辨率为1km。需要强调的是这些被选取的地表数据的年份都晚于土地覆盖数据1年，是因为地表这些特性的改变通常都滞后于土地覆盖的变化。然后从MERRA-2再分析数据集(第2版本，https：//disc.gsfc.nasa.gov/daac-bin/FTPSubset2.p1)获取月地表净下行长波辐射通量和下行短波辐射通量(2001-2010年)，分辨率为0.5°×0.625°。还获取了地上2m处的月空气温度数据(最高温度、最低温度)，来源于中国气象数据服务中心(http：//data.cma.cn/en)。

2)数据分析

基于上述两种土地覆盖数据生成的广东省人工林、天然林与开阔地(草地、农田)，2000-2010年间未变化的人工林、天然林与开阔地像元被作为未变化的人工林、天然林与开阔地，主要用于表征假设的造林时空模式用于后续的分析其潜在的影响。因为通常只选取当年的土地类型数据会不可避免导致一些分类误差。另外，这些30m分辨率的数据都利用最邻近采样方法聚合成1km分辨率并与其他气候变化进行相互匹配。

对于获取的所有地表数据，都先进行了数据质量的筛查，确保数据的有效性，其中获取的黑空和白空反照率数据值进行了均值处理用于表示某一时间的反照率。进而所有的MODIS数据首先进行月均值处理，然后整合为季节均值，其中白天和夜晚的地表温度数据、平均最高与最低气温、ET、反照率进行了年平均处理。此外，平均白天与夜晚地表与空气温度被作为昼夜平均地表与空气温度。还对MERRA-2下行长波与短波辐射数据进行双线性内插处理保证所有变量降尺度后与其他变量进行统一。

3)建模假设造林的生物物理机制对温度的影响

a-基于移动窗口空间采样法的造林及其相关参数的变化分析

人工林假设变化即一种假设的森林变化概念，并没有在现实中发生，带来的是潜在的影响。所有的变化发生在单一年份的采样格网内。采样格网采用5×5km尺寸覆盖整个研究区域，在每个格网内量化2000-2010年天然林、草地与农田转为人工林的时空模式。对于每个格网及其土地覆盖变化的量化需遵照准则如下：人工林覆盖至少在5％以上，天然林或者开阔地覆盖在80％以上，基于此准则将一些不符合的格网进一步去除。除了5×5km尺寸以外，其他的采样格网尺寸也被证明能够得到相似的结论，因而不做其他的格网选择。

对所有格网进行量化分析，计算未变化人工林与开阔地对应的多年平均LST、反照率、空气温度、下行长波辐射与下行短波辐射(2001-2010年)并得到两者的差值。即基于格网采样尺度的移动窗口空间采样方法实现假设变化。人工林对于区域观测地表温度、空气温度、生物物理参数的影响，可参照造林的反照率差异表达方程(1)：

Δα＝α_PF-α_OP (1)；

其中Δα表示造林后的反照率空间变化，α_PF是造林后的人工林分布区域的反照率均值，α_OP是造林前的开阔地分布区域的反照率均值。

b-利用地表能量平衡模型和内在生物物理机制理论建模评估造林对地表的影响

准备内在生物物理机制理论(IBM)以及地表能量平衡方程中所需要的参数。首先，利用造林相关的土地变化数据、净辐射、土壤热通量、卫星观测地表温度、空气温度，基于3a)中的移动窗口空间采样方法对造林前后每个土地覆盖类型所在区域计算月与年能量重分配因子f的分布按照如下公式(2)计算：

其中T_s为卫星观测的地表温度；T_a为空气温度；

为净辐射

G＝0.14(T_a，n-T_a，n-1)，G为土壤热通量(n表示月份数1-12)；

表示长波辐射反馈的温度敏感性(ε_s为表面发射率，默认为1)，单位为(K(Wm^-2)^-1)；SW_net为净地表短波辐射，SW_net＝(1-α)SW↓，单位为W m^-2；LW↓为入射长波辐射，单位为W m^-2；SW↓为入射短波辐射，单位为W m^-2。其中λ₀、f、

G分别表示造林前农田、草地和天然林对应变量的值，是控制建模T_s变化的重要参量。若不考虑大气反馈，

Δα为反照率变化值；ΔG为土壤热通量扩散变化；Δf是湍流能量重分配的平均变化。另外，需要注意因为采用的数据较多数据单位不同。

再利用3a)移动窗口空采样方法，将计算所得的上述值输入地表能量平衡方程构建的地表温度响应模型，模型如公式(3-7)：

该公式依据公式(1)被改进为如下公式(4-7)：

ΔT_s＝ΔT_{s_α}+ΔT_{s_G}+ΔT_{s_f} (4)；

其中上述ΔT_s为建模地表温度；ΔT_{s_α}为辐射强迫应对反照率变化的地表温度响应；ΔT_{s_G}为土壤热通量扩散的地表温度响应；ΔT_{s_f}为湍流能量再分配的地表温度响应。ΔT_s为正值表明造林带来增温影响，否者则为降温影响，零值则表明影响不明显。

4)结果与分析

a--造林对于地表生物物理参数与地表通量的评价

图2显示了造林前后不同地表生物物理参数与通量按月份变化的差异值。从造林后的能量重分配因子f值的季节差异发现绝大多数造林的f值基本高于零值，这也表明造林的非辐射过程起到了主导作用，值越高表明作用越大，特别是暖季f值是最高的。然而天然林替换为人工林的f值并不是很明显。造林前后ET的变化显示农田造林后的ET值最低，天然林转为人工林后的ET值最高，但总体都低于0，而且ΔET和Δf没有呈现明显的关联关系。造林前后

与反照率以及下行短波辐射有相关关系，主要受反照率变化的影响，呈现负相关关系(图3)，特别是草地造林。造林前后G值与反照率有一定关系，但与短波辐射关系并不明显。

b-造林对于地表温度影响的分析与评价

图4指出了建模后的地表温度以及主导地表温度变化的其他三个主要要素的分布情况。农田造林和天然林位置处造林带来的是年净增温，具体为0.083±0.034K与0.058±0.004K。相反，草地造林带来的是年净降温-0.058±0.019K，降温的趋势并不是很明显。其中由f变化驱动的建模地表温度变化与总体的趋势保持一致，并且比其他的驱动因素对总体的贡献更大。总的来说，由反照率驱动的辐射过程对于建模地表温度的贡献不足，而以f变化为代表的非辐射过程作为一个主要分离变量，从农田造林与天然林造林主导增温影响，从草地造林主导降温影响(图5-6)。

本发明基于多源遥感数据、气候数据和能量平衡模型集成，估测了造林的生物物理机制对地表温度的影响，该方法不仅尝试了以高分辨率的土地覆盖数据为基础提高对地表观测的精确性，还构建了详实的森林响应温度变化模型，为评价区域造林响应气候提供了科学的方法论。

上述实施例仅为本发明专利较佳实施方式，但本发明专利的实施方式不仅限于此，其他任何未背离本发明专利原理下作的改变修饰、替代、简化均为等效的置换方式，都包含在本发明专利技术方案的包含范围之内中。

Claims (9)

1.基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，包括以下步骤：

1)收集可用的遥感数据、气候数据，以及土地覆盖数据；

2)利用空间分析的方法生成高效分辨率的未变化人工林、天然林以及开阔地的时空分布图；

3)将收集到的数据进行尺度转换与匹配，统一分辨率，并基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析；

4)构建地表能量平衡方程与内在生物物理理论集成的造林响应地表温度辨别的模型。

2.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，步骤1)中，可用的遥感数据、气候数据，以及土地覆盖数据包括：MODIS地表温度MYD11A2数据、MODIS蒸散发MOD16A2数据、MODIS反照率MCD43B3数据、气温数据、MREEA-2再分析数据和高程数据。

3.根据权利要求2所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，所述MODIS反照率MCD43B3数据包括黑空反照率与白空反照率。

4.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，步骤3)中，在将收集到的数据进行尺度转换与匹配，统一分辨率时，利用双线性内插法将收集到的数据的分辨率升尺度或者降尺度为1km，实现尺度统一。

5.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，步骤3)中，将收集到的数据进行尺度转化与匹配具体包括：将黑空反照率与白空反照率进行均值处理；对按月份、季节、年份的所有MODIS数据进行求均值处理，生成按月份、年际的地表温度白天值、地表温度夜晚值、平均最高气温、平均最低气温、蒸散发与反照率；将MERRA-2再分析数据用于提取晴空的地表下行长波辐射、下行短波辐射；利用随机森林模型反演站点观测的温度数据，重采样到1km，基于温度数据与高程值的线性回归方程，修正由海拔差异造成的温度差异。

6.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，步骤3)中，在基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析时，遵循以下准则：人工林覆盖至少在5％以上，天然林或者开阔地覆盖在80％以上。

7.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，步骤3)中，在基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析时，参照造林的反照率差异表达方程(1)：

Δα＝α_PF-α_OP (1)；

其中Δα表示造林后的反照率空间变化，α_PF是造林后的人工林分布区域的反照率均值，α_OP是造林前的开阔地分布区域的反照率均值。

8.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，步骤3)中，在基于格网采样尺度的移动窗口空间采样法进行造林的时空变化模拟及相关参数的变化分析时，对造林前后每个土地覆盖类型所在区域计算月与年能量重分配因子f的分布按照如下公式(2)计算：

其中T_s为卫星观测的地表温度；T_a为空气温度；

为净辐射；

SW_net为净地表短波辐射，SW_net＝(1-α)SW_↓，单位为W m^-2；LW_↓为入射长波辐射，单位为W m^-2；SW_↓为入射短波辐射，单位为W m^-2；G为土壤热通量，G＝0.14(T_a，n-T_a，n-1)，n表示月份数1-12；

表示长波辐射反馈的温度敏感性，单位为K(Wm^-2)^-1；其中λ₀、f、

G分别表示造林前农田、草地和天然林对应变量的值，若不考虑大气反馈，则

Δα为反照率变化值；ΔG为土壤热通量扩散变化；Δf是湍流能量重分配的平均变化。

9.根据权利要求1所述的基于多源遥感与能量平衡模型造林响应地表温度的方法，其特征在于，在对造林前后每个土地覆盖类型所在区域计算月与年能量重分配因子f的分布按照如下公式(2)计算后，将计算所得的上述值输入地表能量平衡方程构建的地表温度响应模型，模型如公式(3-7)：

该公式依据公式(1)被改进为如下公式(4-7)：

ΔT_s＝ΔT_{s_α}+ΔT_{s_G}+ΔT_{s_f} (4)；

其中ΔT_s为建模地表温度；ΔT_{s_α}为辐射强迫应对反照率变化的地表温度响应；ΔT_{s_G}为土壤热通量扩散的地表温度响应；ΔT_{s_f}为湍流能量再分配的地表温度响应；当ΔT_s为正值表明造林带来增温影响，否者则为降温影响，零值则表明影响不明显。

Images