首页 / 专利库 / 物理 / 饱和蒸汽 / 基于遥感地表温度与植被盖度两阶段空间探测地表蒸散的方法

基于遥感地表温度与植被盖度两阶段空间探测地表蒸散的方法

阅读:506发布:2023-01-24

专利汇可以提供基于遥感地表温度与植被盖度两阶段空间探测地表蒸散的方法专利检索,专利查询,专利分析的服务。并且一种基于遥感地表 温度 与植被盖度两阶段空间探测地表蒸散的方法,其包括:步骤1,利用遥感数据产品获取地表温度和植被盖度;步骤2,构建地表温度与植被盖度两阶段特征空间;步骤3,利用地表温度与植被盖度两阶段特征空间提取植被和裸土的组分温度;步骤4,结合地表温度与植被盖度两阶段特征空间和Priestley‑Taylor公式计算植被和裸土组分的 蒸发 比;步骤5,根据 辐射 收支平衡方程计算植被组分的可利用 能量 和裸土组分的可利用能量;步骤6,通过地表蒸散的二源模式计算输出植被蒸腾、 土壤 蒸发、及地表蒸散。,下面是基于遥感地表温度与植被盖度两阶段空间探测地表蒸散的方法专利的具体信息内容。

1.一种基于遥感地表温度与植被盖度两阶空间探测地表蒸散的方法,其包括:
步骤1,利用遥感数据产品获取地表温度和植被盖度;
步骤2,构建地表温度与植被盖度两阶段特征空间;
步骤3,利用地表温度与植被盖度两阶段特征空间提取植被和裸土的组分温度;
步骤4,结合地表温度与植被盖度两阶段特征空间和Priestley-Taylor公式计算植被和裸土组分的蒸发比;
步骤5,根据辐射收支平衡方程计算植被组分的可利用能量和裸土组分的可利用能量;
步骤6,通过地表蒸散的二源模式计算输出植被蒸腾、土壤蒸发、及地表蒸散;
步骤1中,利用MODIS遥感数据产品,根据MODIS产品的质量文件,通过十进制向二进制的转换,实现低质量和无效的MODIS地表温度(LST)和植被指数(NDVI)数据的自动过滤;
然后,利用数据产品转换公式,获取地表温度(LST)和植被指数(NDVI);
最后,利用植被指数(NDVI)通过式1计算植被覆盖度(fv):
式中,NDVImin对应于裸土;NDVImax对应于全覆盖植被,两者的取值分别设为0.2和0.86;
其中NDVI由MODIS产品获得;
步骤2中,所述地表温度与植被盖度两阶段特征空间由干边和湿边确定;其中,所述湿边为一平直线,其对应的地表温度等于近地表气温(Ta),近地表气温(Ta)由气象站实测得到;
所述干边为一倾斜直线,其由裸土表面温度(Tsmax)和植被表面温度(Tvmax)确定;其中,对于位于干边上的裸土,裸土表面温度(Tsmax)通过式2计算:
max 4 max 3 max 2 max
as×(Ts ) +bs×(Ts ) +cs×(Ts ) +ds×(Ts )+es=0   (式2);
式2中的系数设置如下:
式3中,as~es为针对裸土的系数;
εs为裸土发射率,为0.95;σ为斯蒂芬-玻尔兹曼常数,为5.67×10-8W·m-1·K-4;
ρ为空气密度,为1.293kg·m-3;cp为空气的质量定压热容,为1005.0J·kg-1·K-1;ras为裸土表面的空气动学阻抗,由式8确定;n为土壤热通量与地表净辐射通量的比值,为
0.35;
εa为空气的发射率,由式6确定;αs为裸土表面反照率,为0.3;Sd代表下行太阳短波辐射,由气象站实测得到;Ta为近地表气温,实测获得;
对于位于干边上的植被,植被表面温度(Tvmax),通过式4计算:
av×(Tvmax)4+bv×(Tvmax)3+cv×(Tvmax)2+dv×(Tvmax)+ev=0   (式4);
式4中的系数设置如下:
式5中,av~ev为针对植被的系数;
εv为植被发射率,为0.98;
rav为植被表面的空气动力学阻抗,由式12确定;αv为植被表面反照率,为0.2;εa为空气的发射率,通过式6确定:
式6中,Lv=2.5×106J·kg-1;Rv=461J·kg-1·K-1;T0=273K;Ta为近地表气温,实测获得;Ta与Td的单位均为K;Td为露点温度,Td通过式7计算:
式7中,a=17.27;b=237.7℃;RH为相对湿度,单位为百分比,实测获得;式7中,Ta与Td的单位均为摄氏度;
式3中,ras通过式8计算:
式8中,u1m为裸土表面1米高处的速,单位m/s;u*为摩擦速度,单位m/s,实测获得;k为von Karman常数,为0.41;zm风速的参考高度,为1m;d为零平面位移,对于裸土为0m;zom为动量传输的粗糙度长度,对于裸土为0.005m;Ψm(1)和 分别是稳定度校正系数,由式10或式11确定;
稳定度校正系数Ψm(1)和 的计算与裸土表面Monin-Obukhov长度(Ls)的取值有关,其中Ls可通过式9确定:
式9中,g为重力加速度,取值9.8m/s2;Hs为裸土表面的感热通量;
当Ls<0时,通过下式计算稳定度校正系数:
当Ls>0时,通过下式计算稳定度校正系数:
式5中,rav由式12确定:
式12中,zT为气温观测的参考高度,为2m;zoh与植被高度(hc)有关,zoh=hc/70m,植被高度(hc)实测获得;  和  为稳定度校正系数,通过式14或15确定,与植被表面Monin-Obukhov长度(Lv)的取值有关;
针对植被表面,Lv通过式13确定:
当Lv<0时,通过式14计算 和
当Lv>0时,通过式15计算 和
以Tsmax=Tvmax=Ta为初始条件,利用迭代方式,通过上述式2-式15计算得到新的Tsmax和Tvmax值,判断新值和旧值之间的差异,设置判断阈值,当差异大于阈值时,用新值代替旧值,继续迭代计算;当小于阈值时,输出最终结果的Tsmax和Tvmax值;
计算得到最终结果的Tsmax和Tvmax值后,得到干边的数学式,如式16所示:
LST=(Tvmax-Tsmax)×fv+Tsmax  (式16)
max max max
式16中,Tv -Ts 表示干边的斜率;Ts 表示干边的截距;fv为植被覆盖度。
2.如权利要求1所述的基于遥感地表温度与植被盖度两阶空间探测地表蒸散的方法,其特征在于:
在步骤3中,将所示地表温度与植被盖度的两阶段特征空间分为上三形ABC、下三角形ACD;地表温度与植被盖度的散点如位于下三角形ABC之中,则为P点,地表温度与植被盖度的散点如位于上三角形ACD之中,则为P’点;
对于P点,植被组分温度(Tv)和裸土组分温度(Ts)通过式17确定:
式17中,OP和ON分别代表O点到P点以及O点到N点的垂直距离,OP/ON通过式18计算:
式18中,LST代表P点的地表温度;
对于P’点,Tv和Ts通过式19确定:
式19中,MP′和MO分别代表M点到P’点以及M点到O点的垂直距离,MP′/MO可通过式20计算:
式20中,LST′代表P’点的地表温度;
其中O点为地表温度与植被盖度两阶空间的原点。
3.如权利要求2所述的基于遥感地表温度与植被盖度两阶空间探测地表蒸散的方法,其特征在于:
在步骤4中,
如果地表温度与植被盖度的散点位于下三角形ABC空间内,即为P点,那么其所对应的裸土组分的蒸发比(EFs)和植被组分的蒸发比(EFv)可通过式23计算:
式23中,OP/ON通过式18计算;
如果地表温度与植被盖度的散点位于上三角形ACD空间内,即为P’点,那么其所对应的裸土组分的蒸发比(EFs)和植被组分的蒸发比(EFv)可通过式24计算:
式24中,MP′/MO通过式20计算;
其中,Δ为饱和蒸汽压曲线在某空气温度处的斜率;γ为湿度计常数;
Δ/(Δ+γ)与近地表气温密切相关,采取式22确定:
4.如权利要求3所述的基于遥感地表温度与植被盖度两阶空间探测地表蒸散的方法,其特征在于:
在步骤5中,植被组分的可利用能量Qv和裸土组分的可利用能量Qs,可通过式25计算:
5.如权利要求4所述的基于遥感地表温度与植被盖度两阶空间探测地表蒸散的方法,其特征在于:
在步骤6中,植被蒸腾通过植被蒸腾潜热通量(LEv)表现,植被蒸腾潜热通量LEv通过式
26计算:
LEv=EFv×Qv  (式26)
土壤蒸发通过土壤蒸发潜热通量(LEs)表示,土壤蒸发潜热通量LEs通过式27计算:
LEs=EFs×Qs  (式27)
地表蒸散通过地表蒸散潜热通量(LE)表现,地表蒸散潜热通量(LE)通过式28计算:
LE=fv×LEv+(1-fv)×LEs  (式28)。

说明书全文

基于遥感地表温度与植被盖度两阶段空间探测地表蒸散的

方法

技术领域

[0001] 本发明涉及地表蒸散的遥感探测领域,具体来说涉及基于地表温度与植被盖度特征空间的地表蒸散遥感探测方法。

背景技术

[0002] 地表蒸散是地表循环过程中的关键参量,直接或间接地影响着地表水分和热量的收支状况。地表蒸散主要由植被蒸腾和土壤蒸发两个部分组成,其中土壤蒸发反映了土壤中的水分通过非生物作用进入大气的过程,与表层土壤湿度密切相关;而植被蒸腾则反映土壤中的水分通过生物作用进入大气的过程,与植被根层的土壤湿度密切相关。一般在农业应用中,土壤蒸发被认为是水分浪费,而植被蒸腾则被认为是真正的水分利用,并且植被蒸腾及根层土壤湿度与植物的生长、发育、以及产量形成等直接相关。因此,对于节水农业、农业干旱监测、灌溉需水量估算、农作物长势监测、以及农作物产量预测等,准确的土壤蒸发和植被蒸腾数据以及表层和根层土壤湿度信息,具有非常重要的实用价值。
[0003] 传统的土壤蒸发和植被蒸腾观测方式,包括涡动相关系统、波文比系统、蒸渗仪、树干液流观测、稳定同位素等通常具有代价高昂、操作复杂、或空间代表性不足等缺陷,不适宜大范围、长时间观测。遥感技术可在一定程度上,弥补传统观测方式的缺陷,其一个像元代表地面的一区域,可以很方便地实现由点到面的转换。因此,遥感成为监测地表蒸散和土壤湿度的重要手段。自上世纪70-80年代至今,基于遥感技术的地表蒸散和土壤湿度监测取得了长足进展。在这些方法中,一类利用地表温度与植被盖度(Land Surface Temperature-Fractional vegetation cover,LST-Fv)特征空间的方法,吸引了国内外众多学者的关注,获得了广泛应用,因为它具有许多显著的优点,比如:(1)LST-Fv特征空间易于反演植被或裸土组分温度、实现地表蒸散的分离;(2)LST-Fv特征空间与深层土壤湿度的变异密切相关;(3)该方法简单易用,对蒸散和土壤湿度的估算精度与某些复杂模型的估算精度相当;(4)LST-Fv特征空间可用于计算蒸发比,易于由遥感瞬时测算的蒸散向日尺度转换;(5)LST-Fv特征空间主要使用遥感数据,较少依赖其他辅助数据,并且主要利用地表温度的相对变化信息,并不要求使用其绝对值,等等。
[0004] 目前,学者们基于LST-Fv特征空间发展了多个估算地表蒸散的一源模型和二源模型,一源模型包括Jiang and Islam(1999,2001)的三算法和Wang(2006)的矩形算法;二源模型包括Nishida’s two-source model、TTME模型、以及HTEM模型等。一源模型仅能计算地表蒸散,而二源模型不仅可以计算地表蒸散,还可以实现地表蒸散的分离,输出土壤蒸发和植被蒸腾。然而现有的二源模型均基于图1(a)所示的地表温度与植被盖度的传统特征空间模式。
[0005] 需要注意的是,在土壤由湿润到干燥的驱动过程下,由于植被根系作用,植被与裸土的辐射温度存在着响应速度和程度的差异。图1(a)所示的传统特征空间模式考虑了响应程度的差异,裸土响应程度大,植被响应程度小,所以才形成梯形形状;但是没有考虑响应速度的差异。不考虑响应速度的差异,会导致植被和裸土组分温度的反演存在较大误差,进而导致植被蒸腾、土壤蒸发、和地表蒸散的估算存在较大误差,直接影响了农业干旱监测、灌溉需水量估算、以及农作物产量预测等的应用效果。

发明内容

[0006] 本发明的目的在于建立一种基于地表温度与植被盖度两阶段特征空间,探测植被蒸腾、土壤蒸发、以及地表蒸散的方法。
[0007] 本发明的基于遥感地表温度与植被盖度两阶空间探测地表蒸散的方法,其包括:
[0008] 步骤1,利用遥感数据产品获取地表温度和植被盖度;
[0009] 步骤2,构建地表温度与植被盖度两阶段特征空间;
[0010] 步骤3,利用地表温度与植被盖度两阶段特征空间提取植被和裸土的组分温度;
[0011] 步骤4,结合地表温度与植被盖度两阶段特征空间和Priestley-Taylor公式计算植被和裸土组分的蒸发比;
[0012] 步骤5,根据辐射收支平衡方程计算植被组分的可利用能量和裸土组分的可利用能量;
[0013] 步骤6,通过地表蒸散的二源模式计算输出植被蒸腾、土壤蒸发、及地表蒸散。
[0014] 优选地,步骤1中,利用MODIS遥感数据产品,根据MODIS产品的质量文件,通过十进制向二进制的转换,实现低质量和无效的MODIS地表温度(LST)和植被指数(NDVI)数据的自动过滤;
[0015] 然后,利用数据产品转换公式,获取地表温度(LST)和植被指数(NDVI);
[0016] 最后,利用植被指数(NDVI)通过式1计算植被覆盖度(fv):
[0017] (式1);
[0018] 式中,NDVImin对应于裸土;NDVImax对应于全覆盖植被,两者的取值分别设为0.2和0.86;其中NDVI由MODIS产品获得。
[0019] 优选地,步骤2中,所述地表温度与植被盖度两阶段特征空间由干边和湿边确定;其中,所述湿边为一水平直线,其对应的地表温度等于近地表气温(Ta),近地表气温(Ta)由气象站实测得到;
[0020] 所述干边为一倾斜直线,其由裸土表面温度(Tsmax)和植被表面温度(Tvmax)确定;其max中,对于位于干边上的裸土,裸土表面温度(Ts )通过式2计算:
[0021] as×(Tsmax)4+bs×(Tsmax)3+cs×(Tsmax)2+ds×(Tsmax)+es=0  (式2);
[0022] 式2中的系数设置如下:
[0023]   (式3);
[0024] 式3中,as~es为针对裸土的系数;
[0025] εs为裸土发射率,为0.95;σ为斯蒂芬-玻尔兹曼常数,为5.67×10-8W·m-1·K-4;
[0026] ρ为空气密度,为1.293kg·m-3;cp为空气的质量定压热容,为1005.0J·kg-1·K-1;ras为裸土表面的空气动学阻抗,由式8确定;n为土壤热通量与地表净辐射通量的比值,为
0.35;
[0027] εa为空气的发射率,由式6确定;αs为裸土表面反照率,为0.3;Sd代表下行太阳短波辐射,由气象站实测得到;Ta为近地表气温,实测获得;
[0028] 对于位于干边上的植被,植被表面温度(Tvmax),通过式4计算:
[0029] av×(Tvmax)4+bv×(Tvmax)3+cv×(Tvmax)2+dv×(Tvmax)+ev=0  (式4);
[0030] 式4中的系数设置如下:
[0031]   (式5);
[0032] 式5中,av~ev为针对植被的系数;
[0033] εv为植被发射率,为0.98;
[0034] rav为植被表面的空气动力学阻抗,由式12确定;αv为植被表面反照率,为0.2;εa为空气的发射率,通过式6确定:
[0035]   (式6);
[0036] 式6中,Lv=2.5×106J·kg-1;Rv=461J·kg-1·K-1;T0=273K;Ta为近地表气温,实测获得;Ta与Td的单位均为K;Td为露点温度,Td通过式7计算:
[0037]   (式7);
[0038] 式7中,a=17.27;b=237.7℃;RH为相对湿度,单位为百分比,实测获得;式7中,Ta与Td的单位均为摄氏度;
[0039] 式3中,ras通过式8计算:
[0040]   (式8);
[0041] 式8中,u1m为裸土表面1米高处的速,单位m/s;u*为摩擦速度,单位m/s,实测获得;k为von Karman常数,为0.41;zm风速的参考高度,为1m;d为零平面位移,对于裸土为0m;zom为动量传输的粗糙度长度,对于裸土为0.005m;Ψm(1)和 分别是稳定度校正系数,由式10或式11确定;
[0042] 稳定度校正系数Ψm(1)和 的计算与裸土表面Monin-Obukhov长度(Ls)的取值有关,其中Ls可通过式9确定:
[0043]   (式9);
[0044] 式9中,g为重力加速度,取值9.8m/s2;Hs为裸土表面的感热通量;
[0045] 当Ls<0时,通过下式计算稳定度校正系数:
[0046]   (式10);
[0047] 当Ls>0时,通过下式计算稳定度校正系数:
[0048]   (式11);
[0049] 式5中,rav由式12确定:
[0050]   (式12);
[0051] 式12中,zT为气温观测的参考高度,为2m;zoh与植被高度(hc)有关,zoh=hc/70m,植被高度(hc)实测获得; 和 为稳定度校正系数,通过式14或15确定,与植被表面Monin-Obukhov长度(Lv)的取值有关;
[0052] 针对植被表面,Lv通过式13确定:
[0053]   (式13);
[0054] 当Lv<0时,通过式14计算 和
[0055]   (式14);
[0056] 当Lv>0时,通过式15计算 和
[0057]   (式15);
[0058] 以Tsmax=Tvmax=Ta为初始条件,利用迭代方式,通过上述式2-式15计算得到新的Tsmax和Tvmax值,判断新值和旧值之间的差异,设置判断阈值(0.001K),当差异大于阈值时,用新值代替旧值,继续迭代计算;当小于阈值时,输出最终结果的Tsmax和Tvmax值。
[0059] 计算得到最终结果的Tsmax和Tvmax值后,得到干边的数学式,如式16所示:
[0060] LST=(Tvmax-Tsmax)×fv+Tsmax  (式16)
[0061] 式16中,Tvmax-Tsmax表示干边的斜率;Tsmax表示干边的截距;fv为植被覆盖度。
[0062] 优选地,在步骤3中,将所示地表温度与植被盖度的两阶段特征空间分为上三角形ABC、下三角形ACD;地表温度与植被盖度的散点如位于下三角形ABC之中,则为P点,地表温度与植被盖度的散点如位于上三角形ACD之中,则为P’点;
[0063] 对于P点,植被组分温度(Tv)和裸土组分温度(Ts)通过式17确定:
[0064]   (式17)
[0065] 式17中,OP和ON分别代表O点到P点以及O点到N点的垂直距离,OP/ON通过式18计算:
[0066]   (式18)
[0067] 式18中,LST代表P点的地表温度;
[0068] 对于P’点,Tv和Ts通过式19确定:
[0069]   (式19)
[0070] 式19中,MP′和MO分别代表M点到P’点以及M点到O点的垂直距离,MP′/MO可通过式20计算:
[0071]   (式20)
[0072] 式20中,LST′代表P’点的地表温度;
[0073] 其中O点为地表温度与植被盖度两阶空间的原点。
[0074] 优选地,如果地表温度与植被盖度的散点位于下三角形ABC空间内,即为P点,那么其所对应的裸土组分的蒸发比(EFs)和植被组分的蒸发比(EFv)可通过式23计算:
[0075]   (式23)
[0076] 式23中,OP/ON通过式18计算;
[0077] 如果地表温度与植被盖度的散点位于上三角形ACD空间内,即为P’点,那么其所对应的裸土组分的蒸发比(EFs)和植被组分的蒸发比(EFv)可通过式24计算:
[0078]   (式24)
[0079] 式24中,MP′/MO通过式20计算;
[0080] 其中,Δ为饱和蒸汽压曲线在某空气温度处的斜率;γ为湿度计常数;
[0081] Δ/(Δ+γ)与近地表气温密切相关,采取式22确定:
[0082]   (式22)。
[0083] 优选地,在步骤5中,植被组分的可利用能量Qv和裸土组分的可利用能量Qs,可通过式25计算:
[0084]   (式25)。
[0085] 优选地,在步骤6中,植被蒸腾通过植被蒸腾潜热通量(LEv)表现,植被蒸腾潜热通量LEv通过式26计算:
[0086] LEv=EFv×Qv  (式26)
[0087] 土壤蒸发通过土壤蒸发潜热通量(LEs)表示,土壤蒸发潜热通量LEs通过式27计算:
[0088] LEs=EFs×Qs  (式27)
[0089] 地表蒸散通过地表蒸散潜热通量(LE)表现,地表蒸散潜热通量(LE)通过式28计算:
[0090] LE=fv×LEv+(1-fv)×LEs  (式28)。
[0091] 本发明构建的地表蒸散探测方法,基于地表温度与植被盖度的两阶段特征空间模式(如图1(b)所示),兼顾了植被与裸土辐射温度存在的响应速度和响应程度差异,可以提高植被和裸土组分温度的反演精度,并进一步改善地表蒸散的监测效果,有望为农作物灌溉需水量估算、农业干旱监测、农作物产量预测等应用需求,提供更为精确的地表蒸散估算及其分离的方法模型,也有助于提高我国定量遥感的水平、支撑相关领域的科学研究进展。附图说明
[0092] 图1.(a)地表温度与植被盖度的传统特征空间模式;(b)地表温度与植被盖度的两阶段特征空间模式
[0093] 图2.利用地表温度与植被盖度两阶段特征空间进行地表参数反演
[0094] 图3.涡动相关仪观测站空间分布
[0095] 图4.遥感估算EF与实测EF的散点图比较,图中(a)~(f)分别对应于7、10~14等六个涡动相关观测站
[0096] 图5.遥感估算EF与实测EF的时序比较,图中(a)~(f)分别对应于7、10~14等六个涡动相关观测站,DOY为天数,年份为2012
[0097] 图6.不同初始场景下MTME模型计算的LST-SVI理论特征空间,(a)场景1,(b)场景2[0098] 图7.当P点在三角形ABC以内时,MTME模型对输入参数的敏感性
[0099] 图8.当P点在三角形ACD以内时,MTME模型对输入参数的敏感性
[0100] 图9.本发明方法的流程图

具体实施方式

[0101] 一种基于地表温度与植被盖度两阶段特征空间探测地表蒸散的方法,其特征在于包含以下功能或步骤:
[0102] A、利用遥感数据产品获取地表温度和植被盖度
[0103] 以MODIS遥感数据产品为例,首先根据MODIS产品的质量文件,通过十进制向二进制的转换功能,实现低质量和无效的MODIS地表温度(LST)和植被指数(NDVI)数据的自动过滤。然后利用数据产品转换公式,获取地表温度和植被指数,最后利用植被指数通过下式计算植被覆盖度(fv):
[0104]   (式1)
[0105] 式中,NDVImin对应于裸土;NDVImax对应于全覆盖植被,两者的取值分别设为0.2和0.86。
[0106] B、构建地表温度与植被盖度两阶段特征空间
[0107] 地表温度与植被盖度两阶段特征空间是由代表地表温度变化的纵轴和代表植被盖度变化的横轴所组成的二维散点图,该散点图被抽象为由上下两个共边三角形组成的梯形,每个三角形均由一组从底部到顶部逐渐干燥的土壤湿度等值线组成,下三角形反映表层土壤湿度的变异,上三角形则反映植被根层土壤湿度的变异;两阶段特征空间边界由四个特征顶点组成,即裸土覆盖地表潜热最小/感热最大、裸土覆盖地表潜热最大/感热最小、植被全覆盖地表潜热最大/感热最小、以及植被全覆盖地表潜热最小/感热最大。通过在地表能量平衡方程式中,设置一些极值状态,求解两阶段特征空间的理论边界。
[0108] 地表温度与植被盖度两阶段特征空间是由代表地表温度变化的纵轴和代表植被盖度变化的横轴所组成的二维散点图,该散点图的边界轮廓形似为一个倾斜的梯形,其内部是由上下两个共边三角形组成,每一个三角形内部均由一组渐变的土壤湿度等值线组成,从底部到顶部代表土壤湿度越来越干燥,下三角形反映表层土壤湿度的变异,上三角形则反映植被根层土壤湿度的变异。两阶段特征空间边界由4个特征顶点组成,分别代表地表能量平衡中的一些极值状态,即裸土覆盖地表潜热最小/感热最大(A点)、裸土覆盖地表潜热最大/感热最小(B点)、植被全覆盖地表潜热最大/感热最小(C点)、以及植被全覆盖地表潜热最小/感热最大(D点)等4种极值状态。BC的连线定义为湿边,代表地表蒸散达到潜在最大状态。AD的连线定义为干边,代表地表蒸散达到理论最小。AC的连线为两阶段特征空间的分界线,代表土壤蒸发为理论最小而植被蒸腾为潜在最大。从CB到CA,一组围绕C点逐渐抬升的土壤湿度等值线构成两阶段特征空间的下三角形ABC。从AC到AD,一组围绕A点逐渐抬升的土壤湿度等值线构成两阶段特征空间的上三角形ACD。
[0109] 由于湿边代表地表蒸散达到潜在最大状态,根据能量平衡方程,湿边对应的地表温度接近于近地表气温Ta。因此,地表温度与植被盖度两阶段特征空间的湿边BC为一水平直线,其对应的地表温度等于Ta,Ta由气象站实测得到。
[0110] 地表温度与植被盖度两阶段特征空间的干边,为一倾斜直线。对于位于干边上的裸土,裸土表面温度Tsmax通过下式计算:
[0111] as×(Tsmax)4+bs×(Tsmax)3+cs×(Tsmax)2+ds×(Tsmax)+es=0  (式2)[0112] 式中的系数设置如下:
[0113]   (式3)
[0114] 式中,as~es为针对裸土的a~e等系数;εs为裸土发射率(取0.95);σ为斯蒂芬-玻尔兹曼常数(5.67×10-8 W·m-1·K-4);ρ为空气密度(1.293kg·m-3);cp为空气的质量定压-1 -1 s热容(1005.0J·kg ·K );ra为裸土表面的空气动力学阻抗;n为土壤热通量与地表净辐射通量的比值(取0.35);εa为大气的比辐射率;αs为裸土表面反照率(取0.3);Sd代表下行太阳短波辐射,由气象站实测得到。
[0115] 对于位于干边上的植被,植被表面温度Tvmax,通过下式计算:
[0116] av×(Tvmax)4+bv×(Tvmax)3+cv×(Tvmax)2+dv×(Tvmax)+ev=0  (式4)[0117] 公式中的系数设置如下:
[0118]   (式5)
[0119] 式中,av~ev为针对植被的a~e等系数;εv为植被发射率(取0.98);rav为植被表面的空气动力学阻抗;αv为植被表面反照率(取0.2)。
[0120] 空气的发射率εa则通过下式确定:
[0121]   (式6)
[0122] 式中,Lv=2.5×106J·kg-1;Rv=461J·kg-1·K-1;T0=273K;Ta为实测气温;Td为露点温度,可利用百分比湿度RH的观测量进行估算,估算式如下:
[0123]   (式7)
[0124] 式中,a=17.27;b=237.7℃;RH为相对湿度,单位为百分比;Ta与Td的单位均为摄氏度。值得注意的是,Ta与Td的单位均为K。
[0125] 在确定了上述参数后,剩余两个关键的参数是植被和裸土表面的空气动力学阻抗rav和rsv。利用如下系列表达式,计算这两个空气动力学阻抗。
[0126] (a)裸土表面空气动力学阻抗
[0127]   (式8)
[0128] 式中,u1m为裸土表面1米高处的风速,单位m/s;u*为摩擦速度,单位m/s,可由通量观测站提供;k为von Karman常数,取值0.41;zm风速的参考高度,这里取1m;d为零平面位移,对于裸土取0m;zom为动量传输的粗糙度长度,对于裸土取0.005m;Ψm(1)和 分别是稳定度校正系数。稳定度校正系数的计算与裸土表面Monin-Obukhov长度(Ls)的取值有关,其中Ls可通过下式确定:
[0129]   (式9)
[0130] 式中,g为重力加速度,取值9.8m/s2;Hs为裸土表面的感热通量。当Ls<0时,通过下式计算稳定度校正系数:
[0131]   (式10)
[0132] 当Ls>0时,通过下式计算稳定度校正系数:
[0133]   (式11)
[0134] (b)植被表面空气动力学阻抗
[0135] 对于植被表面空气动力学阻抗rac,则通过下式确定:
[0136]   (式12)
[0137] 式中,zT为气温观测的参考高度,取2m;zoh与植被高度(hc)有关,假设植被高度为hc=1m,zoh=hc/70m; 和 为稳定度校正系数,与植被表面Monin-Obukhov长度(Lv)的取值有关。针对植被表面,Lv通过下式确定:
[0138]   (式13)
[0139] 当Lv<0时,通过下式计算 和
[0140]   (式14)
[0141] 当Lv>0时,通过下式计算 和
[0142]   (式15)
[0143] 在利用公式2-公式15计算图1(b)中A点的地表温度Tsmax和D的地表温度Tvmax的过程中,Monin-Obukhov长度(L)的计算很重要,然而L的计算又需要Tsmax和Tvmax作为输入参数之一。因此需要迭代求解算法,初始条件下设Tsmax=Tvmax=Ta,然后通过以上系列公式计算得到新的Tsmax和Tvmax值,判断新值和旧值之间的差异,设置判断阈值(0.001K),当大于阈值时,用新值代替旧值,迭代计算;当小于阈值时,输出最终结果。计算得到Tsmax和Tvmax之后,干边AD可用如下方程表示:
[0144] LST=(Tvmax-Tsmax)×fv+Tsmax  (式16)
[0145] 式中,Tvmax-Tsmax表示为干边的斜率;Tsmax为干边的截距;fv为植被覆盖度。
[0146] C、利用地表温度与植被盖度两阶段特征空间提取植被和裸土的组分温度[0147] 需要首先判断地表温度与植被盖度散点在两阶段特征空间中的位置,当散点位于下三角形,植被组分温度与近地表气温相同,裸土组分温度则通过散点相对于下三角形的边界位置,内插得到;当地表温度与植被盖度散点位于上三角形,裸土组分温度取裸土表面潜热最小/感热最大的极值状态,植被组分温度则利用散点相对于上三角形的边界位置,内插得到。
[0148] 由于地表温度与植被盖度的两阶段特征空间由上下两个三角形组成,这两个三角形具有不同的物理意义,因此需要首先判断地表温度与植被盖度散点在两阶段特征空间中的位置。如图2所示,地表温度与植被盖度散点可能位于三角形ABC之中(即P点),也可能位于三角形ACD之中(即P’点)。对于P点,植被组分温度(Tv)和裸土组分温度(Ts)通过下式确定:
[0149]   (式17)
[0150] 式中,OP和ON分别代表O点到P点以及O点到N点的垂直距离,OP/ON通过下式计算:
[0151]   (式18)
[0152] 式中,LST代表P点的地表温度。
[0153] 对于P’点,Tv和Ts通过下式确定:
[0154]   (式19)
[0155] 式中,MP′和MO分别代表M点到P’点以及M点到O点的垂直距离,MP′/MO可通过下式计算:
[0156]   (式20)
[0157] 式中,LST′代表P’点的地表温度。
[0158] D、结合地表温度与植被盖度两阶段特征空间和Priestley-Taylor公式计算植被和裸土组分的蒸发比
[0159] 需要首先判断地表温度与植被盖度散点在两阶段特征空间中的位置,当地表温度与植被盖度散点位于下三角形,则利用下三角形内插Priestley-Taylor参数,计算土壤表面蒸发,此时植被蒸腾处于潜在最大状态,取Priestley-Taylor参数为1.26计算得到;当地表温度与植被盖度散点位于上三角形,则利用上三角形内插Priestley-Taylor参数,计算植被蒸腾,此时土壤表面蒸发处于理论最小状态,取Priestley-Taylor参数为0计算得到。
[0160] 根据Priestley-Taylor公式,蒸发比(EF)的计算表达式如下:
[0161]   (式21)
[0162] 式中φ代表Priestley-Taylor参数,与地表湿润状况密切相关;Δ为饱和蒸汽压曲线在某空气温度处的斜率;γ为湿度计常数。Δ/(Δ+γ)与近地表气温密切相关,采取如下的经验方程确定:
[0163]   (式22)
[0164] 以图2为例,如果地表温度与植被盖度散点位于三角形ABC空间内,即P点,那么其所对应的裸土组分的蒸发比(EFs)和植被组分的蒸发比(EFv)可通过下式计算:
[0165]   (式23)
[0166] 式中,OP/ON式18计算。
[0167] 如果地表温度与植被盖度散点位于三角形ACD空间内,即P’点,那么其所对应的EFs和EFv可通过下式计算:
[0168]   (式24)
[0169] 式中,MP′/MO式20计算。
[0170] E、根据辐射收支平衡方程计算植被和裸土组分的可利用能量
[0171] 植被和裸土组分的可利用能量Qv和Qs,可通过下式计算:
[0172]   (式25)
[0173] 式中,Rnv为植被组分的净辐射能量,αv为植被组分的反照率,Sd为入射太阳辐射,εv为植被组分的发射率,εa为空气的发射率,σ为斯特潘-玻尔兹曼常数,Ta为空气温度,Tv为植被组分温度;n为裸土土壤热通量占裸土净辐射能量Rns的比值,αs为裸土组分的反照率,Ts为裸土组分的温度。
[0174] F、通过地表蒸散的二源模式计算输出植被蒸腾、土壤蒸发、及地表蒸散[0175] 二源模式的基本原理是分别估算裸土组分和植被组分的蒸发比,再通过加权平均得到混合像元的蒸发比。假设混合地表由植被和裸土两种简单组分组成,在不考虑两种组分之间能量交互的基础上,混合地表的潜热通量LE等于植被和裸土组分潜热通量的加权平均,权重值为植被覆盖度。
[0176] 植被蒸腾潜热通量LEv通过下式计算得到:
[0177] LEv=EFv×Qv  (式26)
[0178] 土壤蒸发潜热通量LEs通过下式计算得到:
[0179] LEs=EFs×Qs  (式27)
[0180] 地表蒸散潜热通量则通过下式计算得到:
[0181] LE=fv×LEv+(1-fv)×LEs  (式28)
[0182] 式中,LEv和LEs分别为植被和裸土组分的潜热通量(W/m2);fv为植被覆盖度。
[0183] 实例:
[0184] 本发明实施区域选择为甘肃省张掖市黑河流域中游的盈科绿洲,遥感数据为2012年4月份~2012年10月份的MODIS地表温度和植被指数产品,具体包括(1)Terra MODIS地表温度产品MOD11A1,时间分辨率为daily,空间分辨率1km;(2)Aqua MODIS地表温度产品MYD11A1,时间分辨率为daily,空间分辨率1km;(3)Terra MODIS植被指数产品MOD13A2,时间分辨率为16-day,空间分辨率1km;(4)Aqua MODIS植被指数产品MYD13A2,时间分辨率为16-day,空间分辨率1km。气象观测数据为2012年4月份~2012年10月份左右的近地表气温(℃)、空气湿度(%)、太阳下行短波辐射(W/m2)、摩擦速度(m/s)。选用的测试站点均位于张掖市盈科灌区,下垫面均为玉米,各站点所在空间位置如图3所示。这些站点其空间位置和海拔,以及涡动相关仪的架高等信息如表1所示。
[0185] 表1.地表蒸散观测站点位置及海拔信息
[0186]
[0187]
[0188] 1.系统流程
[0189] A、数据预处理
[0190] 根据MODIS产品的质量文件,通过十进制向二进制的转换功能,实现低质量和无效数据的自动过滤,以及符合质量要求数据的自动提取。
[0191] B、获取地表温度和植被盖度
[0192] MODIS数据产品为HDF格式,首先对其进行格式转换、投影转换、裁切等预处理,然后根据地表温度转换系数,计算得到以K为单位的地表温度,最后通过植被盖度的计算公式,将MODIS植被指数转换为植被盖度。
[0193] C、计算地表温度与植被盖度两阶段特征空间的理论边界
[0194] 根据气象站数据,通过迭代运算的方式,计算地表温度与植被盖度两阶段特征空间的理论干边和湿边。
[0195] D、提取植被和裸土的组分温度
[0196] 基于地表温度与植被盖度两阶段特征空间的理论边界,在每一个通量观测站对应的像素上,对遥感观测的地表温度进行分解,得到植被和裸土的组分温度。
[0197] E、计算植被和裸土组分的蒸发比
[0198] 根据植被和裸土组分温度以及地表温度与植被盖度两阶段特征空间的理论边界,结合Priestley-Taylor公式,计算植被和裸土组分的蒸发比。
[0199] F、计算植被和裸土组分的可利用能量
[0200] 根据地表辐射收支平衡方程,分别计算植被和裸土组分的可利用能量。
[0201] G、计算输出植被蒸腾、土壤蒸发、及地表蒸散
[0202] 通过地表蒸散的二源模式,计算输出植被蒸腾、土壤蒸发、及地表蒸散。
[0203] 2.结果分析及本方法的优越性
[0204] A.本发明模型与其他模型及实测数据的比较
[0205] 图4以散点图的形式,在每个涡动相关仪观测站上,比较了站点实测EF和基于多种方法估算的EF,这些模型包括本发明方法(简称MTME模型)、Long and Singh(2012)发展的TTME二源模型、Long et al.(2012)发展的梯形算法,简称Trapezoid、Wang et al.(2006)发展的矩形算法,简称Rectangle、以及Jiang and Islam发展的三角形算法,简称Triangle。图(a)~图(f)分别对应6个EF观测站点:7号站、10号站、11号站、12号站、13号站、以及14号站,每个站点上多个时相数据的比较,构成了图中的散点样式。图中横轴为实测值,纵轴为遥感估算值。图中黑色斜线代表实测值与估算值的1:1线,各散点越接近该1:1线,越表明估算值与实测值的一致性较高。图4表明:在多数站点上,MTME模型的估算结果比其他算法更接近于实测值,表现在图4中,红色圆点更接近1:1等值线。
[0206] 表2.遥感估算EF与实测EF的定量比较
[0207]
[0208] 选择若干统计指标,进一步比较了遥感估算EF与实测EF。这些统计指标包括,均方根误差(Root Mean Square Error,RMSE)、平均相对误差(Mean Relative Error,MRE)、以及相关系数R,结果如表2所示。表2的统计结果表明,在大多数站点上,MTME模型的RMSE要比其他模型的RMSE小,MRE也要比其他模型的MRE小。从平均值来看,MTME模型的RMSE为0.14,R为0.71,MRE为17%。其次是Rectangle和TTME模型的RMSE及MRE较小,Trapezoid模型的较大。总之,表2的统计结果从站点数据的比较上,证明了本发明方法的优越性,其平均估算误差为17%,精度达到83%。
[0209] 图5比较了遥感估算EF与实测EF在时间序列上的变化,图中EC表示涡动相关仪(Eddy Covariance)观测的EF,图中横轴表示天数(Day of year,DOY),年份是2012年。其中第200天对应于2012年7月18号;第265天对应于2012年9月21号。该图表明:从DOY 200(7月中旬)到DOY 250(9月上旬),植被生长比较旺盛,地表的EF值较高,自进入9月份之后,EF逐渐下降。MTME模型估算的EF与实测EF具有较为一致的时间变化信息,在大多数站点上,MTME模型的估算结果比其他算法的估算结果更接近EF的实测值。总之,图5从遥感估算值与实测值的时序变好比较上,也说明了本发明方法的优越性。
[0210] B.本发明模型的敏感性分析
[0211] 模型的敏感性分析有助于了解模型的误差来源及其不确定性,选定两个初始场景(两组初始的MTME输入参数),初始场景选定的依据是像素点在两阶段特征空间中的位置。表3展示了这两个初始场景的参数值,包括气象站实测的空气温度Ta、空气湿度RH、太阳下行短波辐射Sd、和摩擦速度U_star,以及MODIS提供的植被指数NDVI和地表温度LST,其中场景1对应于10号站、第250天、Terra MODIS过境时刻(上午11:30)的环境;场景2对应于11号站、第207天、Terra MODIS过境时刻(上午10:12)的环境。
[0212] 表3.MTME初始场景的参数设置
[0213]
[0214]
[0215] 将表3所示的初始场景参数代入MTME模型,可计算得到这两个场景所对应的两阶段特征空间。结果如图6所示,图6(a)和(b)分别为场景1和2对应的两阶段特征空间及遥感观测像素在特征空间的位置,图中ABCD为计算得到的LST-SVI特征空间理论边界上的特征点,P点为MODIS观测像素在两阶段特征空间中的位置,O点为P点对应的AC边上的特征点,M点为P点对应在AD边上的特征点。在场景1中P点位于三角形ABC以内,在场景2中P点位于三角形ACD以内。
[0216] 在图6所示的场景中,逐个调节每个输入参数,然后计算新场景下的EF与初始EF的相对大小,以此可评价模型对该调节参数的敏感性。其中将LST和Ta的调节范围设为[LST-2K,LST+2K]和[Ta-2K,Ta+2K],调节步长设为0.5K;将其余参数的调节范围设为其本身大小的±20%,调节步长设为其本身大小的5%,比如Sd的调节范围为[Sd-Sd×20%,Sd+Sd×
20%],调节步长为Sd×5%。EF的相对变化(Relative Change,RC)定义为:
[0217]   (式29)
[0218] 式中,EF±和EF0分别为调节参数下和初始参数下的MTME模型的输出值EF。
[0219] 图7为场景1环境下,即P点在三角形ABC以内时,MTME模型对输入参数的敏感性。图8为场景2环境下,即P点在三角形ACD以内时,MTME模型对输入参数的敏感性。对比这两个图可知:
[0220] (1)MTME模型对NDVI的敏感性,在不同的场景下有所不同。当P点在三角形ABC以内时,EF随着NDVI的增大而增大,当NDVI从-20%到+20%发生相对变化时,EF的相对变化为-5%到+5%。但是当P点在三角形ACD以内时,EF随着NDVI的增大而减小,而且EF针对NDVI的相对变化比较敏感,NDVI增加50%时,EF减小40%左右。
[0221] (2)无论是在场景1还是在场景2的环境下,MTME模型对RH都不太敏感。在两个场景中,EF均随着RH的增大而增大。在场景1中,当RH从-20%到+20%发生相对变化时,EF的相对变化在±1%以内。在场景2中,当RH从-20%到+20%发生相对变化时,EF的相对变化在±3%以内。
[0222] (3)在场景1和场景2中,EF均随着Sd的增大而增大。当Sd从-20%到+20%发生相对变化时,EF在场景1中的相对变化在±5%以内,但在场景2中的相对变化在±30%以内。说明在场景2中,MTME模型对Sd更敏感。
[0223] (4)在两个场景中,EF均随着U_star的增大而减小,这可能是因为U_star的增大使得空气动力学阻抗减小,从而增大了感热通量,使得蒸发比减小。MTME模型对U_star的敏感性也在场景1中较小,在场景2中较大。
[0224] (5)在两个场景中,EF随着LST的增大而减小,随着Ta的增大而增大,这可能是因为LST的增大将P点拉向干边,从而EF减小;而Ta的增大使得P点向湿边相对移动,从而EF增大。MTME对LST和Ta的敏感性,在场景1中较小,在场景2中较大。
[0225] 总之,MTME模型对输入参数Ta、RH、Sd、U_star、NDVI、和LST的敏感性与地表像素P在两阶段特征空间中的位置有关,当P在三角形ABC以内时,MTME对这些输入参数较不敏感,当P点在三角形ACD以内时,MTME模型对这些参数较敏感。P点在三角形ABC以内,对应着植被蒸散接近潜在最大的状态,此时植被不受干旱的胁迫。而P点在三角形ACD以内的情形,对应着植被受到干旱胁迫的情景。这意味着,本发明方法在植被不发生干旱时,稳定性较高,但在植被发生干旱的时期,可能具有一定的不稳定性。
高效检索全球专利

专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台,是提供专利分析,专利查询,专利检索等数据服务功能的知识产权数据服务商。

我们的产品包含105个国家的1.26亿组数据,免费查、免费专利分析。

申请试用

分析报告

专利汇分析报告产品可以对行业情报数据进行梳理分析,涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度,系统通过AI智能系统对图表进行解读,只需1分钟,一键生成行业专利分析报告。

申请试用

QQ群二维码
意见反馈