一种基于遥感信息的地表水、热和碳通量耦合估算方法与流程

本发明涉及一种地表碳水通量耦合估算方法，尤其涉及一种基于遥感信息的地表水、热和碳通量耦合估算方法，属于遥感技术领域。

背景技术：

陆地和大气之间的水热交换是水文学、气象学和全球变化研究领域关心的重要物理过程。植被和大气之间的碳交换是生态系统碳循环的核心内容，也对全球气候变化起着关键作用。因此，地表显热通量、潜热通量和碳通量的研究对深入理解水循环过程机理，理解陆地生态系统功能、二氧化碳的源汇过程及它们与气候变化的关系极其重要。由于近三十年来遥感观测能力的显著提高，并且由于遥感观测具有传统方法不可比拟的大面积动态监测能力，利用遥感技术进行地表水、热和碳通量的研究目前已经成为重要的方法和手段。目前，基于遥感的机理性的碳水耦合估算方法仍然很少，相关研究仍处于起步阶段。现有的耦合方法没有区分光照和阴影下光合速率的差异，没有考虑土壤水分对冠层导度的影响，土壤蒸发计算没有直接引入遥感信息较繁琐，没有同时估算地表显热通量。

现有方法存在以下缺陷：第一，光照和阴影下植被的气孔导度差异较大，导致光合作用差异较大，现有的方法没有考虑这种差异，因此对光合速率的估算造成误差；第二，土壤水分胁迫对植被气孔导度有显著影响，现有的耦合方法没有考虑土壤水分限制的影响，给气孔导度的计算带来误差，进而给植被蒸腾的估算带来误差；第三，现有的方法对土壤蒸发的计算需要降雨数据或者额外的土壤水分估算模型，较繁琐；第四，现有方法没有地表显热通量的计算过程，而显热通量是陆气之间能量交换的重要过程。

技术实现要素：

为了解决上述技术所存在的不足之处，本发明提供了一种基于遥感信息的地表水、热和碳通量耦合估算方法。

为了解决以上技术问题，本发明采用的技术方案是：一种基于遥感信息的地表水、热和碳通量耦合估算方法，方法的步骤包括：

步骤s1、地面观测数据获取；获取研究区内所有地表通量观测台站的地面观测数据；

步骤s2、地面数据预处理；按各台站的植被类型将台站归类，相同类型划分为一类；按年月日的顺序排列数据，数据格式统一化、数据异常值剔除、数据缺失值插补，以保证数据的连续性，有效性；

步骤s3、地面输入数据制备；计算以下地面输入数据：土壤净辐射、植被净辐射、土壤热通量、空气动力学阻抗、摩擦风速、土壤边界层阻抗、整体叶片边界层阻抗、植被吸收的光合有效辐射、光照和阴影叶片吸收的光合有效辐射、温度限制因子、台站土壤水分限制因子、植被覆盖率、空气饱和水汽压、叶片表面co2浓度、碳通量以及冠层导度；

步骤s4、模型参数率定；率定设定的多目标函数为：

f＝2-(nseet+nsegpp)公式28

其中，ac是碳通量，le为潜热通量，nse为纳什效率系数，下标et和gpp分别表示蒸散发和生态系统总初级生产力；n是观测值的个数，下标sim和obs代表模拟和观测；f代表目标函数；最小化的f值所对应的参数即是最后率定得到的参数结果值；

步骤s5、遥感数据获取；下载时间和空间分辨率相同的遥感数据产品、与地表反照率产品时间最接近的植被无明显变化时期的叶面积指数产品lai、茎面积指数sai、植被聚集指数产品ω以及下垫面地表类型没有明显变化时期的土地覆盖土地利用产品；

步骤s6、气象数据获取；根据下载的地表反照率产品的时间和空间分辨率，下载时间和空间分辨率最接近的气象数据；

步骤s7、其它参数获取；根据研究区下垫面的情况给出占主导的植被类型的冠层高度hc和叶片直径s的经验值；

步骤s8、区域数据预处理；根据下载的地表反照率产品的空间分辨率，空间范围和投影坐标，将下载的遥感数据做几何裁剪、几何纠正和重投影，使得这些数据能够完全配准到一起；将下载的气象数据做空间重采样或空间插值，使得数据与地表反照率产品具有相同的空间分辨率，之后再做几何裁剪、几何纠正和重投影，使得气象数据与其它数据能够完全配准到一起；

步骤s9、区域输入数据制备；利用下载的数据计算所需要的输入参数，包括地表净辐射rn，土壤净辐射rns，植被净辐射rnc，土壤热通量g，空气动力学阻抗ra，土壤边界层阻抗rs，整体叶片边界层阻抗rx，植被吸收的光合有效辐射apar，光照下植被最大光能利用率βmsl，阴影下植被最大光能利用率βmsh，地表水分指数lswi，温度限制因子ft，土壤水分限制因子fw，叶片饱和水汽压ec^*，冠层内空气水汽压eac；

步骤s10、碳通量计算；

ac＝(βmsl×aparsl+βmsh×aparsh)×ft×fw公式33

其中，ac是碳通量，βmsl为光照下植被最大光能利用率，βmsh为阴影下植被最大光能利用率，aparsl和aparsh是光照和阴影叶片吸收的光合有效辐射，ft和fw是温度限制因子和土壤水分限制因子；

步骤s11、气孔导度计算；

其中，m和d0是经验系数，通过率定得到；rc是冠层阻力；是空气饱和水汽压，ea是空气实际水汽压，cb是叶片表面co2浓度；

步骤s12、地表潜热通量和显热通量的计算；

其中，lec是植被蒸腾，δ是饱和水汽压曲线斜率，ρ是空气密度，cp是定压比热，γ是干湿表常数；

les＝fw·lep公式36

其中，les是土壤蒸发；lep是土壤潜在蒸散发；

le＝les+lec公式38

其中，le是地表总蒸散发。lec是植被蒸腾；les是土壤蒸发；

hc＝rnc-lec公式39

其中，hc是植被显热通量；rnc是植被净辐射；lec是植被蒸腾；

hs＝rns-les-g公式40

其中，hs是土壤显热通量；rns是土壤净辐射；les是土壤蒸腾；g是土壤热通量；

h＝hs+hc公式41

其中，h是地表总蒸散发；hc是植被显热通量；hs是土壤显热通量。

进一步地，步骤s3中地面输入数据的计算方法为：

a、土壤净辐射：

rns＝rn×exp[-ωlai]公式1

其中，ω是消光系数，等于0.45，rn是净辐射，rns是土壤净辐射；

b、植被净辐射：

rnc＝rn{1-exp[-ωlai]}公式2

其中，rnc是植被净辐射，lai是叶面积指数；

c、土壤热通量：

g＝0.35rns公式3

其中，rns是土壤净辐射，g是土壤热；

d、空气动力学阻抗：

其中ra是空气动力学阻抗，zu和zt是风速和气温的观测高度；d是零平面位移；zom和zov分别是动量传输和能量传输粗糙度长度，um是风速；k是冯卡曼常数；

e、摩擦风速：

其中，u*为摩擦风速，其他各参数含义与上式相同。

进一步地，步骤s3中地面输入数据的计算方法为：

f、土壤边界层阻抗：

cscs，barews+cs，dense(1-ws)公式7

ws＝e^-(lai+sai)公式9

其中，rs是土壤边界层阻抗；cs是湍流传输系数；u*是摩擦风速；cs，bare是裸土表面的湍流传输系数；cs，dense是浓密植被表面的湍流传输系数；zomg是地表动量粗糙度长度；v是空气运动粘滞长度；sai是茎面积指数；a是系数；ws代表权重；

g、整体叶片边界层阻抗：

其中，rb是整体叶片边界层阻抗；u*是摩擦风速；

h、植被吸收的光合有效辐射：

apar＝fpar×par公式11

fpar＝1-exp(-kl×lai)公式12

其中fpar是光合有效辐射吸收比率；apar是植被吸收的光合有效辐射；kl是消光系数。

进一步地，步骤s3中地面输入数据的计算方法为：

i、光照和阴影叶片吸收的光合有效辐射：

pardif＝par(0.7527+3.8453r-16.316r²+18962r³-7.0802r⁴)公式15

pardif，u＝pardifexp[-0.5ωlai/(0.537+0.025lai)]公式16

c＝0.07ωpardir(1.1-0.1lai)exp(-cosθ)公式17

pardir＝par-pardif公式18

laish＝lai-laisl公式20

其中σ是地表反照率，由遥感产品得到，pardif和pardir是散射和直射光合有效辐射，由公式15-18计算；r是天空清晰度指数；s0是太阳常数；pardif，u是冠层下的散射光合有效辐射；c是单位面积多次散射对冠层内散射辐射的贡献率；θ是太阳天顶角；laisl和laish是光照和阴影叶面积指数；ω是聚集指数。

进一步地，步骤s3中地面输入数据的计算方法为：

j、温度限制因子：

其中ta，min和ta，max是能够进行光合作用的最小和最大植被温度；ta，opt是能够实现最大光能利用率的最佳植被温度；它们通过参数率定获得；

k、台站土壤水分限制因子：

其中，sw是土壤水分；swmax和swmin是土壤田持水量和凋萎系数；fw是土壤水分限制因子；

l、植被覆盖率：

fc＝1-exp(-0.5lai)公式23

其中，fc是植被覆盖率，lai是叶面积指数；

m、空气饱和水汽压：

其中，是空气饱和水汽压；ta是近地表空气温度；

n、叶片表面co2浓度：

cb＝ca-ac(1.3rb+ra)公式25

其中，cb是叶片表面co2浓度；ca是大气co2浓度；ac是碳通量；rb是纠正叶片气孔分布不均匀后的叶片整体边界层阻抗；ra是空气动力学阻抗；

o、碳通量：

ac＝(βmsl×aparsl+βmsh×aparsh)×ft×fw公式26其中，ac是碳通量；βmsl和βmsh分别是光照和阴影下的最大光能利用率；aparsl和aparsh是光照和阴影叶片吸收的光合有效辐射；ft和fw是温度限制因子和土壤水分限制因子；

p、冠层导度：

其中，m和d0是经验系数，通过率定得到；rc是冠层阻力；是空气饱和水汽压，ea是空气实际水汽压，cb是叶片表面co2浓度。

进一步地，地面观测数据包括日平均的净辐射rn、地表反照率σ，气温ta、风速um、空气水汽压ea、大气co2浓度ca、土壤水分sw；日总量的光合有效辐射par、潜热通量le和gpp；当日的冠层高度hc、叶片直径s、叶面积指数lai、茎面积指数sai；台站的下垫面类型，饱和土壤水分sws和凋萎系数swmin。

进一步地，步骤s5中遥感数据产品包括地表反照率产品σ，比辐射率产品ε，par产品；步骤s6中气象数据包括近地表空气温度ta、近地表空气湿度ea、近地表风速μ、近地表太阳辐射s0和近地表大气co2浓度ca。

进一步地，步骤s9中地表水分指数的计算方法为：

其中，lswi是土壤水分指数；vnir和ρswir分别是地表红光波段和短波红外波段的反射率；

土壤水分限制因子的计算方法为：

其中，fw是土壤水分限制因子；lswimax是地表水分指数的年最大值。

本发明具有的有益效果为：本发明提出了一种基于遥感信息的地表水、热和碳通量耦合估算方法，实现了地表潜热通量(le)，显热通量(h)和碳通量(ac)的联合估算。本发明主要解决的技术问题为四点。第一，融合双叶遥感光能利用率模型，区分了光照和阴影下光合速率的差异，降低了只把冠层作为大叶处理造成的光合速率估算误差；第二，考虑了土壤水分对冠层导度的影响，通过引入土壤水分限制因子定量化土壤水分对冠层导度的胁迫作用，降低了由于没有考虑土壤水分影响造成的冠层导度估算误差；第三，考虑了土壤水分对土壤蒸发的影响，并且利用遥感信息实现了土壤水分对土壤蒸发影响的估算，提高了可操作性和应用性；第四，能够同时估算地表显热通量，通过使用能量平衡方程实现地表潜热通量和显热通量的同时估算，弥补了现有方法中缺少地表显热通量估算的不足。

附图说明

图1为本发明的方法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1所示的一种基于遥感信息的地表水、热和碳通量耦合估算方法，方法的步骤包括：

步骤s1、地面观测数据获取；

获取研究区内所有地表通量观测台站的日平均的地表净辐射rn、地表反照率σ，气温ta、风速um、空气水汽压ea、大气co2浓度ca、土壤水分sw；获取日总量的光合有效辐射par、潜热通量le和gpp；获取当日的冠层高度hc、叶片直径s、叶面积指数lai、茎面积指数sai；按igbp的土地覆盖/土地利用类别确定台站的下垫面类型，饱和土壤水分sws和凋萎系数swmin；

步骤s2、地面数据预处理；

根据igbp划分的土地覆盖/土地利用类型，按各台站的植被类型将台站归类，相同类型划分为一类；按年月日的顺序排列数据，数据格式统一化、数据异常值剔除、数据缺失值插补，以保证数据的连续性，有效性；

步骤s3、地面输入数据制备；

利用获取的观测数据，计算以下地面输入数据，包括：

a、土壤净辐射：

rns＝rn×exp[-ωlai]公式1

其中，ω是消光系数，等于0.45，rn是净辐射，rns是土壤净辐射。

b、植被净辐射：

rnc＝rn{1-exp[-ωlai]}公式2

其中，rnc是植被净辐射，lai是叶面积指数。

c、土壤热通量：

g＝0.35rns公式3

其中，rns是土壤净辐射，g是土壤热。

d、空气动力学阻抗：

其中ra是空气动力学阻抗。zu和zt是风速和气温的观测高度；d是零平面位移，d＝2hc/3，hc是冠层高度；zom和zov分别是动量传输和能量传输粗糙度长度，zom＝0.123hc，zov＝0.1zom；um是风速；k是冯卡曼常数＝0.41。

e、摩擦风速：

u*为摩擦风速，其它各参数含义与上式相同。

f、土壤边界层阻抗：

cs＝cs，barews+cs，dense(1-ws)公式7

ws＝e^-(lai+sai)公式9

其中，rs是土壤边界层阻抗；cs是湍流传输系数；u*是摩擦风速；cs，bare是裸土表面的湍流传输系数；cs，dense是浓密植被表面的湍流传输系数，cs，dense′＝0.004；zomg是地表动量粗糙度长度，zomg＝0.01m；v是空气运动粘滞长度，v＝1.5×10^-5m²s^-1；sai是茎面积指数；a是系数，a＝0.13；ws代表权重。

g、整体叶片边界层阻抗：

其中，rb是整体叶片边界层阻抗；u*是摩擦风速，与前文相同。

h、植被吸收的光合有效辐射：

apar＝fpar×par公式11

fpar＝1-exp(-kl×iai)公式12

其中fpar是光合有效辐射吸收比率；apar是植被吸收的光合有效辐射；kl是消光系数，通常设定为0.5。

i、光照和阴影叶片吸收的光合有效辐射：

pardif＝par(0.7527+3.8453r-16.316r²+18.962r³-7.0802r⁴)公式15

pardif，u＝pardifexp[-0.5ωlai/(0.537+0.025lai)]公式16

c＝0.07ωpardir(1.1-0.1lai)exp(-cosθ)公式17

pardir＝par-pardif公式18

laish＝lai-laisl公式20

其中σ是地表反照率，由遥感产品得到，pardif和pardir是散射和直射光合有效辐射，由公式15-18计算；r是天空清晰度指数＝par/(0.5s0cosθ)；s0是太阳常数＝1367w.m^-2；pardif，u是冠层下的散射光合有效辐射；c是单位面积多次散射对冠层内散射辐射的贡献率；θ：是太阳天顶角；laisl和laish是光照和阴影叶面积指数；ω是聚集指数。

j、温度限制因子：

其中ta，min和ta，max是能够进行光合作用的最小和最大植被温度；ta，opt是能够实现最大光能利用率的最佳植被温度；它们通过参数率定获得。

k、台站土壤水分限制因子：

其中，sw是土壤水分；swmax和swmin是土壤田持水量和凋萎系数；fw是土壤水分限制因子。

l、植被覆盖率：

fc＝1-exp(-0.5lai)公式23

其中，fc是植被覆盖率，lai是叶面积指数。

m、空气饱和水汽压：

其中，是空气饱和水汽压；ta是近地表空气温度。

n、叶片表面co2浓度：

cb＝ca-ac(1.3rb+ra)公式25

其中，cb是叶片表面co2浓度；ca是大气co2浓度；ac是碳通量；rb是纠正叶片气孔分布不均匀后的叶片整体边界层阻抗；ra是空气动力学阻抗。

o、碳通量：

ac＝(βmsl×aparsl+βmsh×aparsh)×ft×fw公式26

其中，ac是碳通量；βmsl和βmsh分别是光照和阴影下的最大光能利用率；aparsl和aparsh是光照和阴影叶片吸收的光合有效辐射；ft和fw是温度限制因子和土壤水分限制因子。

p、冠层导度：

其中，m和d0是经验系数，通过率定得到；rc是冠层阻力；是空气饱和水汽压，ea是空气实际水汽压，cb是叶片表面co2浓度。

步骤s4、模型参数率定；

采用matlab中全局优化的遗传算法进行参数率定，所有相同植被功能类型的观测站用于率定这个类型下的参数值。率定设定的多目标函数为：

f＝2-(nseet+nsegpp)公式28

其中，ac是碳通量，le为潜热通量，nse为纳什效率系数，下标et和gpp分别表示蒸散发和生态系统总初级生产力；n是观测值的个数，下标sim和obs代表模拟和观测；f代表目标函数；最小化的f值所对应的参数即是最后率定得到的参数结果值。

由此，可以率定出不同植被类型下的模型参数包括βmsl，βmsh，m，d0，ta，min，ta，max，ta，opt。每种植被类型具有相同的参数值。

步骤s5、遥感数据获取；

下载时间和空间分辨率相同的遥感数据产品，包括地表反照率产品σ，比辐射率产品ε，par产品；下载与地表反照率产品时间最接近的植被无明显变化时期的叶面积指数产品lai，茎面积指数sai，植被聚集指数产品ω，以及下垫面地表类型没有发生明显变化时期的土地覆盖土地利用产品；

步骤s6、气象数据获取；

根据下载的地表反照率产品的时间和空间分辨率，下载时间和空间分辨率最接近的气象数据，包括近地表空气温度ta、近地表空气湿度ea、近地表风速μ、近地表太阳辐射s0和近地表大气co2浓度ca；

步骤s7、其它参数获取；

根据研究区下垫面的情况给出占主导的植被类型的冠层高度hc和叶片直径s的经验值；

步骤s8、区域数据预处理；

根据下载的地表反照率产品的空间分辨率，空间范围和投影坐标，将下载的地表反照率产品，比辐射率产品，par产品、叶面积指数产品，植被聚集指数产品和叶面积指数产品，土地覆盖/土地利用产品分别做几何裁剪、几何纠正和重投影，使得这些数据能够完全配准到一起；将下载的气象数据(近地表空气温度、近地表空气湿度、近地表风速、近地表太阳辐射和近地表大气co2浓度)做空间重采样或空间插值，使得数据与地表反照率产品具有相同的空间分辨率，之后再做几何裁剪、几何纠正和重投影，使得气象数据与其它数据能够完全配准到一起；

步骤s9、区域输入数据制备；

利用下载的数据计算所需要的输入参数，包括地表净辐射rn，土壤净辐射rns，植被净辐射rnc，土壤热通量g，空气动力学阻抗ra，土壤边界层阻抗rs，整体叶片边界层阻抗rx，植被吸收的光合有效辐射apar，光照下植被最大光能利用率βmsl，阴影下植被最大光能利用率βmsh，地表水分指数lswi，温度限制因子ft，土壤水分限制因子fw，叶片饱和水汽压ec^*，冠层内空气水汽压eac；

对于区域尺度的计算，输入数据的计算方法与上文相同，不同的包括以下：

a、地表水分指数：

其中，lswi是土壤水分指数；ρnir和ρswir分别是地表红光波段和短波红外波段的反射率。

b、土壤水分限制因子：

其中，fw是土壤水分限制因子；lswimax是地表水分指数的年最大值。

步骤s10、碳通量计算；

ac＝(βmsl×aparsl+βmsh×aparsh)×ft×fw公式33

其中，ac是碳通量，βmsl为光照下植被最大光能利用率，βmsh为阴影下植被最大光能利用率，aparsl和aparsh是光照和阴影叶片吸收的光合有效辐射，ft和fw是温度限制因子和土壤水分限制因子；

步骤s11、气孔导度计算；

其中，m和d0是经验系数，通过率定得到；rc是冠层阻力；是空气饱和水汽压，ea是空气实际水汽压，cb是叶片表面co2浓度；

步骤s12、地表潜热通量和显热通量的计算；

a、penman-moteith公式计算植被蒸腾：

其中，lec是植被蒸腾，δ是饱和水汽压曲线斜率，ρ是空气密度，cp是定压比热，γ是干湿表常数；其它参数与前文相同。

b、计算地表潜热通量：

les＝fw·lep公式36

其中，les是土壤蒸发；lep是土壤潜在蒸散发；其它参数与前文相同。

c、计算地表潜热通量：

le＝les+lec公式38

其中，le是地表总蒸散发。lec是植被蒸腾；les是土壤蒸发；

d、计算植被显热通量：

hc＝rnc-lec公式39

其中，hc是植被显热通量；rnc是植被净辐射；lec是植被蒸腾；

e、计算土壤显热通量：

hs＝rns-les-g公式40

其中，hs是土壤显热通量；rns是土壤净辐射；les是土壤蒸腾；g是土壤热通量；

f、计算地表显热通量：

h＝hs+hc公式41

其中，h是地表总蒸散发；hc是植被显热通量；hs是土壤显热通量。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。