CN105841847A - 一种估算地表潜热通量的方法 - Google Patents

Abstract

本发明公开一种估算地表潜热通量的方法，该方法基于Penman‑Monteith模型，运用遥感数据以及地表通量塔观测数据，对模型参数进行获取，并在参数获取时加入地表辐射通量以及大气二氧化碳浓度，最终根据模型以及参数获取地表潜热通量。本发明能够实现通过在获取模型参数的同时加入地表辐射通量以及大气二氧化碳浓度，从而获得提高模型计算精度的效果。

Description

一种估算地表潜热通量的方法

技术领域

本发明涉及一种估算地表潜热通量的方法。

背景技术

地表的潜热通量是陆表能量平衡和水量平衡的重要组成部分，估算地表潜热通量对研究气候变化、水循环过程以及这些过程对全球生态系统的影响有着重要的科学意义(Wang&Dickinson,2012)。同时在实际应用中，农业灌溉、水资源管理等也需要对地表潜热通量进行较高精度的估计(Anderson et al.,2012)。传统的估算地表潜热通量的方法，通常只能利用有限的站点观测数据来估计田间尺度的通量值，所需要的输入参数较多，且估算的空间范围有限。随着遥感观测技术的发展，利用遥感数据估计地表潜热通量成为可能。其观测的范围广，时效性高，观测精度较好，正成为一项先进的估算地表潜热通量的方法(Li et al.,2009)。利用遥感估算地表潜热通量的模型可以分为四类：

1)经验统计的方法，在地面观测的潜热通量与遥感植被指数之间建立回归关系(即统计模型)，通过数据拟合获取模型中的经验参数，再将这一模型应用至其他区域上，估算地表潜热通量(Wang et al.,2010)。这类方法比较依赖训练样本的选择，物理机制不够明确。

2)基于地表能量平衡的方法，在地表能量平衡中，净辐射通量通常(R_n)表达为潜热通量(LE)、感热通量(H)以及土壤热通量(G)的和，因此把潜热通量看作是净辐射通量的一个剩余项，利用观测或模型计算的地表净辐射通量、土壤热通量以及感热通量来估算潜热通量(Kustas&Norman,1997)，计算公式如下：

LE＝R_n-H-G

这类模型中，估算感热通量需要利用遥感地表温度，同时其估算结果对地表温度数据十分敏感。因此当遥感地表温度数据误差较大时，模型估算结果的可靠性不高(Timmermans et al.,2007)。

3)基于数据同化的方法，将遥感反演的各类地表参数，如地表温度、叶面积指数等，与地表过程模型相结合，利用数据同化技术估算地表潜热通量(Xu etal.,2011)。这些过程模型可以包括：辐射传输模型、地表能量平衡模型、水量平衡模型等。基于数据同化的方法尽管使用了物理机制较强的过程模型进行估算，但其模型复杂度较高，计算耗时长，结果容易受到各类误差的影响。

4)基于传统方法的半经验半物理机制模型，利用遥感、气象等各类观测数据反演或者直接获取相应的地表参数(如植被指数、空气饱和水气压差、净辐射通量等)，将这些数据应用到传统模型(如Penman-Monteith模型)中进行地表潜热通量的估算(Zhang et al.,2010；Mu et al.,2011)。主要模型公式如下：

LE_t＝LE_c+LE_s

${\mathrm{LE}}_c=\frac{{\mathrm{\ΔA}}_c+{\mathrm{\ρC}}_p\·VPD\·G_a}{\mathrm{\Δ}+\mathrm{\γ}\·(1+G_a/G_c)}$

${\mathrm{LE}}_s={\mathrm{RH}}^{VPD}\frac{\mathrm{\Δ}}{\mathrm{\Δ}+\mathrm{\γ}}{A}_s$

其中，LE_t为总潜热通量；LE_C为植被冠层的潜热通量；LE_S为土壤的潜热通量；A_C为植被冠层可利用能量；A_S为土壤的可利用能量；VPD为饱和水气压差；Δ为温度－饱和水气压差曲线斜率；γ为温湿计常数；RH为相对湿度；ρ为空气密度；C_p为空气比热容；G_a为空气动力学导度；G_C为植被冠层导度。模型中，植被冠层导度是一个重要的生物物理参数，可以利用一些地表观测数据对其进行估计，如相对湿度、空气温度、空气饱和水气压差、遥感观测的植被指数和叶面积指数等。

发明内容

本发明的目的是提供一种估算地表潜热通量的方法，以解决计算植被冠层导度时误差较大从而导致估算结果误差较大的技术问题。

为实现以上发明目的，本发明提供一种估算地表潜热通量的方法，其特征在于，包括如下步骤：

步骤1：获取遥感数据以及地表通量塔观测数据；

步骤2：对所述遥感数据进行预处理；

步骤3：通过所述预处理的遥感数据以及地表通量塔观测数据获取计算模型所需参数；

步骤4：根据模型计算地表潜热通量。

进一步地，所述步骤1中获取遥感数据的方法包括：

通过中分辨率成像光谱仪获取叶面积指数以及植被指数。

进一步地，所述步骤1中获取植被指数的具体公式如下：

$NDVI=\frac{{\mathrm{\ρ}}_n-{\mathrm{\ρ}}_r}{{\mathrm{\ρ}}_n+{\mathrm{\ρ}}_r}$

其中，“NDVI”是植被指数数据，“ρ_n”是近红外波段反射率，“ρ_r”是红波段反射率。

进一步地，所述步骤2中对所述遥感数据进行预处理的方法包括如下步骤：

步骤201：通过光谱平滑滤波方法消除所述遥感数据中信号噪声；

步骤202：对所述消除信号噪声的遥感数据进行线性内插。

进一步地，所述步骤3中计算模型所需参数包括：

植被与土壤可利用能、空气动力学导度、土壤的热潜热通量以及植被冠层导度。

进一步地，所述步骤3中获取植被与土壤可利用能的方法包括：

根据叶面积指数以及地表净辐射数据获取植被可利用能以及土壤可利用能；具体公式如下：

A_c＝R_n(1-e^-0.6·LAI)

A_s＝R_n·e^-0.6·LAI-G

其中，“Ac”是植被可利用能，“As”是土壤可利用能，“LAI”是叶面积指数数据，“Rn”是日平均地表净辐射数据，“G”是地表土壤热通量；日平均土壤热通量近似取0。

进一步地，所述步骤3中获取空气动力学导度的方法包括：

根据日平均风速以及植被高度数据计算空气动力学导度；具体公式如下：

$G_a=\frac{k^2{u}_m}{ln\[(z_m-d)/z_{om}\]ln\[(z_m-d)/z_{ov}\]}$

其中，“k”是卡曼系数，“Zm”是风速观测参考高度，“Zom”是动量传输粗糙度长度，“Zov”是水汽传输粗糙度长度，“d”是零平面位移；

其中，k＝0.41，d＝2h/3，Zom＝0.123h，Zov＝0.1Zom；“h”是植被高度数据。

进一步地，所述步骤3中获取土壤的热潜热通量的方法包括：

根据相对湿度以及水汽压差计算土壤的热潜热通量；具体公式如下：

${\mathrm{LE}}_s={\mathrm{RH}}^{VPD}\frac{\mathrm{\Δ}}{\mathrm{\Δ}+\mathrm{\γ}}{A}_s$

其中，“Δ”是温度-饱和水汽压斜率，“γ”是温度计常数。

进一步地，所述步骤3中获取植被冠层导度的计算公式如下：

G_c＝(c₁NDVI+c₂)·(VPD^-0.5)·R_n/C_a

其中，“c₁”和“c₂”分别是待定系数，“C_a”是大气二氧化碳浓度；

其中，对于温带阔叶林，所述待定系数C₁＝0.016，C₂＝0。

进一步地，所述步骤4中根据模型计算地表潜热通量的模型公式如下：

${\mathrm{LE}}_t=\frac{{\mathrm{\ΔA}}_c+{\mathrm{\ρC}}_p\·VPD\·G_a}{\mathrm{\Δ}+\mathrm{\γ}\·(1+G_a/\[(c_{1}NDVI+c_2)\·\left({\mathrm{VPD}}^{-0.5}\right)\·R_n/C_a\])}+{\mathrm{LE}}_s$

其中，“LE_t”是地表总的潜热通量，“ρ”是空气密度，“C_p”是空气比热容。

与现有技术相比，本发明的有益效果是：

1.运用将地表辐射通量以及大气二氧化碳浓度代入模型计算的技术方案，获得了提高原有模型适用性的技术效果。

2.运用对遥感数据进行预处理的技术方案，获得了提高遥感数据精确度的技术效果。

附图说明

图1是本发明的估算地表潜热通量的方法的流程框图；

图2是本发明的估算地表潜热通量的方法的一个具体实施例流程图

图3是通量塔站点空间分布图；

图4是模型估计值与潜热通量观测值间的散点图；

图5是使用本发明方法模拟潜热通量结果在站点US-WCr上的时间序列图；

图6是使用LE_Leuning方法模拟潜热通量结果在站点US-WCr上的时间序列图；。

图7是使用LE_Zhang方法模拟潜热通量结果在站点US-WCr上的时间序列图；

图8是使用本发明方法模拟潜热通量结果在站点IT-Ro1上的时间序列图；

图9是使用LE_Leuning方法模拟潜热通量结果在站点IT-Ro1上的时间序列图；

图10是使用LE_Zhang方法模拟潜热通量结果在站点IT-Ro1上的时间序列图；

图11是使用本发明方法模拟潜热通量结果在站点IT-Ro2上的时间序列图

图12是使用LE_Leuning方法模拟潜热通量结果在站点IT-Ro2上的时间序列图；

图13是使用LE_Zhang方法模拟潜热通量结果在站点IT-Ro2上的时间序列图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

实施例1：

如图1所示，本发明的估算地表潜热通量的方法，包括如下步骤：

步骤1：获取遥感数据以及地表通量塔观测数据；

步骤2：对所述遥感数据进行预处理；

步骤3：通过所述预处理的遥感数据以及地表通量塔观测数据获取计算模型所需参数；

步骤4：根据模型计算地表潜热通量。

具体来说，对地表潜热通量的估算主要是在于精度上，背景技术中的四种方式，精度相对较高的是第四种方式，即运用Penman-Monteith模型进行估算，运用该模型进行估算所需的参数可以通过遥感数据以及地表通量塔观测数据获取，由于遥感数据以及地表通量塔观测数据获取的数据存在噪声，从计算精度上考虑，必须对这些存在噪声的数据进行预处理。预处理之后才可用于计算估算模型所需的参数，最终将参数代入模型进行估算。

实施例2：

如图2所示，本发明的估算地表潜热通量的方法，包括如下步骤：

步骤A1：通过中分辨率成像光谱仪获取叶面积指数以及植被指数；获取地表通量塔观测数据；

步骤A201：消除所述遥感数据中信号噪声；

步骤A202：对所述消除信号噪声的遥感数据进行线性内插；

步骤A3：通过所述预处理的遥感数据以及地表通量塔观测数据获取植被与土壤可利用能、空气动力学导度、土壤的热潜热通量以及植被冠层导度；

步骤A4：根据模型计算地表潜热通量。

具体来说，本发明采用基于传统模型:Penman-Monteith模型；对地表潜热通量进行计算。该模型计算地表潜热通量需要植被与土壤可利用能、空气动力学导度、土壤的热潜热通量以及植被冠层导度这些参数。遥感数据以及地表通量塔观测数据包含计算上述参数的一些基础量。所以我们需要先从遥感数据以及地表通量塔观测数据获取这些基础量，一般我们采用中分辨率成像光谱仪(简称：“MODIS”)获取遥感数据。

我们通过MODIS地表反射率数据获取8天的植被指数(简称：“NDVI”)，计算公式如下：

$NDVI=\frac{{\mathrm{\ρ}}_n-{\mathrm{\ρ}}_r}{{\mathrm{\ρ}}_n+{\mathrm{\ρ}}_r}$

其中，ρ_n为近红外波段反射率，ρ_r为红波段反射率

叶面积指数(简称：“LAI”)以及植被指数(简称：“NDVI”)。

由于MODIS遥感数据会受到云层等数据影响，而使得数据不精确，所以我们对初步获取的8天植被指数以及叶面积指数进行预处理。先通过光谱平滑(英文名称：Savitzky-Golay)滤波方法消除数据中云等信号噪声的影响，再对处理过的数据进行线性内插，最终获得更加准确的每天的遥感数据。基础量获取后，就可以根据基础量对模型参数进行计算获取。

植被与土壤可利用能的获取方法如下：

根据叶面积指数以及地表净辐射数据获取植被可利用能以及土壤可利用能；具体公式如下：

A_c＝R_n(1-e^-0.6·LAI)

A_s＝R_n·e-^0.6·LAI-G

其中，“Ac”是植被可利用能，“As”是土壤可利用能，“LAI”是叶面积指数数据，“Rn”是日平均地表净辐射数据，“G”是地表土壤热通量；日平均土壤热通量近似取0。

空气动力学导度的获取方法如下：

根据日平均风速以及植被高度数据计算空气动力学导度；具体公式如下：

$G_a=\frac{k^2{u}_m}{ln\[(z_m-d)/z_{om}\]ln\[(z_m-d)/z_{ov}\]}$

其中，“k”是卡曼系数，“Zm”是风速观测参考高度，“Zom”是动量传输粗糙度长度，“Zov”是水汽传输粗糙度长度，“d”是零平面位移；

我们优选：k＝0.41；d＝2h/3；Zom＝0.123h；Zov＝0.1Zom；其中“h”是植被高度数据。

土壤的热潜热通量的获取方法如下：

根据相对湿度以及水汽压差计算土壤的热潜热通量；具体公式如下：

${\mathrm{LE}}_s={\mathrm{RH}}^{VPD}\frac{\mathrm{\Δ}}{\mathrm{\Δ}+\mathrm{\γ}}{A}_s$

其中，“Δ”是温度-饱和水汽压斜率，“γ”是温度计常数。

植被冠层导度的获取方法如下：

G_c＝(c₁NDVI+c₂)·(VPD^-0.5)·R_n/C_a

其中，“c₁”和“c₂”分别是待定系数，“C_a”是大气二氧化碳浓度。

我们在获取植被可利用能以及土壤可利用能时将地表净辐射数据加入计算，提高了植被可利用能以及土壤可利用能计算结果的精确性。

我们在获取植被冠层导度时将大气二氧化碳浓度加入计算，提高了植被冠层导度计算结果的精确性。由于植被冠层导度对计算地表潜热通量的精度起到较为关键的作用，所以提高植被冠层导度计算进度是相当必要的。

在所有参数都计算获取完后，便可生成计算模型：

${\mathrm{LE}}_t=\frac{{\mathrm{\ΔA}}_c+{\mathrm{\ρC}}_p\·VPD\·G_a}{\mathrm{\Δ}+\mathrm{\γ}\·(1+G_a/\[(c_{1}NDVI+c_2)\·\left({\mathrm{VPD}}^{-0.5}\right)\·R_n/C_a\])}+{\mathrm{LE}}_s$

其中，“LE_t”是地表总的潜热通量，“ρ”是空气密度，“C_p”是空气比热容。

以下提供一个案例：

如图3所示，图中，实心圆形数据仅用于模型验证。实心三角形数据部分用于模型率定，其余部分用于模型验证。实验采用统计方法对模型的效果进行评估。主要采用的指标有：相关系数、决定系数和均方根误差。相关系数数值在-1到1之间，表示两个变量之间的相关程度，负值代表负相关，正值代表正相关，其绝对值越接近于1两者的相关程度越大。决定系数数值也在0到1之间，表示模型中自变量对因变量的解释程度，其数值越接近于1表示这种解释程度越高。均方根误差(RMSE)则反映了模型预测值与观测值的偏差程度，其数值越小结果的精度越高。

实验采用的通量数据来自La Thuile数据集，它是由全球250多个通量站点的数据构成，该观测数据集包括了二氧化碳通量、水热通量、空气温度、空气湿度等多类地表生物物理属性。

实验首先对模型参数进行率定，率定过程采用的通量塔数据如表1所示：

表1

其中，通量塔名称代表了不同森林地区的观测站点(如表1所示)。一个“站点年”代表某个站点一年的观测数据。通量塔US-WCr、IT-Ro1和IT-Ro2中的部分数据应用于模型的率定，剩余部分用于模型的验证(如表1所示)。模型率定共包含14个站点年，经过率定后，模型的待定系数为C₁＝0.015，C₂约为0，因此在应用时可以忽略常数项。

潜热通量模型的表达式为：

${\mathrm{LE}}_t=\frac{{\mathrm{\ΔA}}_c+{\mathrm{\ρC}}_p\·VPD\·G_a}{\mathrm{\Δ}+\mathrm{\γ}\·(1+G_a/\[(c_{1}NDVI+c_2)\·\left({\mathrm{VPD}}^{-0.5}\right)\·R_n/C_a\])}+{\mathrm{LE}}_s$

用于模型验证的通量塔数据如表所示：

表2

模型验证共包含9个站点年。利用模型估算每日的地表潜热通量，将模型预测值与通量塔观测值做回归分析，二者的相关系数为0.84，属于显著相关；决定系数(R²)为0.7，均方根误差为26.8W/m²，这与已有估算结果精度相当(Zhang et al.,2010)。但当地表观测的潜热通量较低时，模型估算结果可能存在一定的高估。

图如4所示，该图表示了模型估计值与潜热通量观测值间的散点图。图中散点基本分布于1:1线附近，表明模型估计值与观测值相当。

最后，在不同站点上对预测值与观测值的季节变化进行比较，同时也对比了不同参数化方法对模型模拟结果的影响，进行比较的参数化方法包括：Leuning等人提出的潜热通量模型(简称：LE_Leuning)(Leuning et al.,2008)和Zhang等人提出的潜热通量模型(简称：LE_Zhang)(Zhang et al.,2010)。

站点US-WCr结果如图5-7所示：

图中蓝色点代表模型预测值，红色点代表站点观测值。模型模拟了站点2004－2006年3年的每日地表潜热通量变化。图1为本发明的结果，图2为LE_Leuning的结果，图3为LE_Zhang的结果。从图中可以看到，在这个站点上三种方法都能比较好地描述每日地表潜热通量季节变化。从统计指标来看，本发明模型预测值与观测值的相关系数为0.77，RMSE为26.0W/m²；LE_Leuning的相关系数为0.83，RMSE为23.4W/m²；LE_Zhang的相关系数为0.54，RMSE为27.3W/m²。因此在此站点上，本发明结果精度较好，与已有结果相近。

站点IT-Ro1结果如图8-10所示：

图中标记含义与上述内容一致。从图中可以看到，三种方法结果也能较好地描述地表潜热通量的季节变化，不过在2005年三种方法对观测的潜热通量都有一定程度的高估。从统计指标来看，本发明模型预测值与观测值的相关系数为0.85，RMSE为29.7W/m²；LE_Leuning的相关系数为0.83，RMSE为39.3W/m²；LE_Zhang的相关系数为0.84，RMSE为46.7W/m²。因此在此站点上，本发明的各个统计指标均优于其它方法。

站点IT-Ro2结果如图11-13所示：

图中标记含义与上述内容一致。从图中可以看到，三种方法结果都能较好地描述地表潜热通量的季节变化。从统计指标来看，本发明模型预测值与观测值的相关系数为0.93，RMSE为23.8W/m²；LE_Leuning的相关系数为0.91，RMSE为31.1W/m²；LE_Zhang的相关系数为0.90，RMSE为40.3W/m²。因此在此站点上，本发明的各个统计指标均优于其它方法，获得了较好的精度。

除上述实施例外，本发明还可以有其他实施方式，凡采用等同替换或等效变换形成的技术方案，均落在本发明的保护范围内。

Claims (10)

1.一种估算地表潜热通量的方法，其特征在于，包括如下步骤：

步骤1：获取遥感数据以及地表通量塔观测数据；

步骤2：对所述遥感数据进行预处理；

步骤3：通过所述预处理的遥感数据以及地表通量塔观测数据获取计算模型所需参数；

步骤4：根据模型计算地表潜热通量。

2.如权利要求1所述的估算地表潜热通量的方法，其特征在于，所述步骤1中获取遥感数据的方法包括：

通过中分辨率成像光谱仪获取叶面积指数以及植被指数。

3.如权利要求2所述的估算地表潜热通量的方法，其特征在于，所述步骤1中获取植被指数的具体公式如下：

$NDVI=\frac{{\mathrm{\ρ}}_n-{\mathrm{\ρ}}_r}{{\mathrm{\ρ}}_n+{\mathrm{\ρ}}_r}$

其中，“NDVI”是植被指数数据，“ρ_n”是近红外波段反射率，“ρ_r”是红波段反射率。

4.如权利要求1所述的估算地表潜热通量的方法，其特征在于，所述步骤2中对所述遥感数据进行预处理的方法包括如下步骤：

步骤201：通过光谱平滑滤波方法消除所述遥感数据中信号噪声；

步骤202：对所述消除信号噪声的遥感数据进行线性内插。

5.如权利要求1所述的估算地表潜热通量的方法，其特征在于，所述步骤3中计算模型所需参数包括：

植被与土壤可利用能、空气动力学导度、土壤的热潜热通量以及植被冠层导度。

6.如权利要求5所述的估算地表潜热通量的方法，其特征在于，所述步骤3中获取植被与土壤可利用能的方法包括：

根据叶面积指数以及地表净辐射数据获取植被可利用能以及土壤可利用能；具体公式如下：

A_c＝R_n(1-e^-0.6·LAI)

A_s＝R_n·^-0.6·LAI-G

其中，“Ac”是植被可利用能，“As”是土壤可利用能，“LAI”是叶面积指数数据，“Rn”是日平均地表净辐射数据，“G”是地表土壤热通量；日平均土壤热通量近似取0。

7.如权利要求5所述的估算地表潜热通量的方法，其特征在于，所述步骤3中获取空气动力学导度的方法包括：

根据日平均风速以及植被高度数据计算空气动力学导度；具体公式如下：

$G_a=\frac{k^2{u}_m}{ln\[(z_m-d)/z_{om}\]ln\[(z_m-d)/z_{ov}\]}$

其中，“k”是卡曼系数，“Zm”是风速观测参考高度，“Zom”是动量传输粗糙度长度，“Zov”是水汽传输粗糙度长度，“d”是零平面位移；

其中，k＝0.41，d＝2h/3，Zom＝0.123h，Zov＝0.1Zom；“h”是植被高度数据。

8.如权利要求5所述的估算地表潜热通量的方法，其特征在于，所述步骤3中获取土壤的热潜热通量的方法包括：

根据相对湿度以及水汽压差计算土壤的热潜热通量；具体公式如下：

${\mathrm{LE}}_s={\mathrm{RH}}^{VPD}\frac{\mathrm{\Δ}}{\mathrm{\Δ}+\mathrm{\γ}}{A}_s$

其中，“Δ”是温度-饱和水汽压斜率，“γ”是温度计常数。

9.如权利要求5-8中任一所述的估算地表潜热通量的方法，其特征在于，所述步骤3中获取植被冠层导度的计算公式如下：

G_c＝(c₁NDVI+c₂)·(VPD^-0.5)·R_n/C_a

其中，“c₁”和“c₂”分别是待定系数，“C_a”是大气二氧化碳浓度；

其中，对于温带阔叶林，所述待定系数C₁＝0.016，C₂＝0。

10.如权利要求9所述的估算地表潜热通量的方法，其特征在于，所述步骤4中根据模型计算地表潜热通量的模型公式如下：

${\mathrm{LE}}_t=\frac{{\mathrm{\ΔA}}_c+{\mathrm{\ρC}}_p\·VPD\·G_a}{\mathrm{\Δ}+\mathrm{\γ}\·(1+G_a/\[(c_{1}NDVI+c_2)\·\left({\mathrm{VPD}}^{-0.5}\right)\·R_n/C_a\])}+{\mathrm{LE}}_s$

其中，“LE_t”是地表总的潜热通量，“ρ”是空气密度，“C_p”是空气比热容。