本发明公开基于遥感数据的区域地表感热/潜热通量反演方法，包括以下步骤：确定研究区域和遥感数据；根据所述遥感数据制备遥感地表和区域气象参数，其中，所述遥感地表参数包括：归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts、叶面积指数LAI；所述区域气象参数包括：空气温度Ta和相对湿度RH；根据所述遥感地表和区域气象参数进行净辐射通量Rn反演；根据所述叶面积指数LAI、空气温度Ta以及净辐射通量Rn进行土壤热通量G反演；根据所述净辐射通量Rn、土壤热通量G、值被覆盖度f和地表温度Ts的理论二维空间以及感热通量估算的温度廓线方程进行辐射‑对流阻抗rae反演；根据所述辐射‑对流阻抗rae实现区域地表感热H/潜热LE通量反演。
基于遥感数据的区域地表感热/潜热通量反演方法及系统
技术领域
[0001] 本发明涉及区域地表感热/潜热通量的反演方法及系统。更具体地，涉及基于遥感数据的区域地表感热/潜热通量反演方法及系统。
背景技术
[0002] 区域地表感热/潜热通量的反演是流域水资源管理、农作物灌溉制度制定、农业用水效率评价、区域旱情监测、农业节水等应用领域的基础，同时也是构建土壤-植被-大气连续体能量和水量迁移转换，进行陆-气耦合模拟研究的前提。区域地表感热/潜热通量的反演结果能够为水利、农业、林业专业技术人员提供重要的下垫面状况信息，为辅助决策提供重要支持信息。
[0003] 遥感作为地面观测的重要手段，为区域地表感热/潜热通量的反演提供了重要的技术手段。现有的基于遥感的地表感热/潜热通量的方法主要分为如下几类：
[0004] 统计经验模型、地表能量平衡模型、地表温度-植被指数特征空间模型、地表阻抗类模型。统计经验模型将站点观测数据同遥感数据结合，建立地表感热/潜热通量的经验回归关系，统计经验模型依赖于地表观测数据，可移植性差，难以应用于区域地表感热/潜热通量的估算。地表能量平衡模型包括单层模型和双层模型，单层模型将下垫面植被和土壤统一考虑，代表性模型如SEBAL模型和SEBS模型。
[0005] 由于遥感反演的地表温度和地表空气动力学温度的差异，导致单层模型的应用和推广存在较大的不确定性。双层模型则将下垫面土壤和植被分开考虑，分别计算土壤和植被的感热/潜热通量，代表性模型如TSEB模型。双层模型参数众多，实际操作和应用受到限制。地表温度-植被指数特征空间模型基于遥感反演的地表温度和植被指数构成的二维空间来反演下垫面的感热/潜热通量信息。但特征空间的边界确定存在较大的主观性和不确定性。地表阻抗类模型是将遥感反演的下垫面参数同传统的感热/潜热模型相结合，实现区域尺度感热/潜热通量的估算。水汽阻抗的参数化是地表阻抗类模型的关键，区域尺度水汽阻抗的刻画仍然存在较多不足。
[0006] 目前，基于遥感的地表感热/潜热通量模型方法仍然具有较大的不确定性，现有的模型方法中，仍然没有一种能够有效的解决无资料区、复杂非均匀下垫面的地表感热/潜热通量估算问题。尤其是无资料区、复杂非均匀下垫面的地表感热/潜热通量缺乏有效的观测手段，是目前研究和实际应用的难点。
[0007] 因此，需要提供一种基于遥感数据的区域地表感热/潜热通量反演方法和系统，以解决上述问题。
发明内容
[0008] 本发明的目的在于提供基于遥感数据的区域地表感热/潜热通量反演方法及系统，在不依赖风速和不需要地表空气动力学阻抗相关参数的模型情况下实现区域地表感热/潜热通量的反演。
[0009] 为达到上述目的，本发明采用下述技术方案：
[0010] 基于遥感数据的区域地表感热/潜热通量反演方法，包括以下步骤：
[0011] S1：确定研究区域和遥感数据，所述遥感数据为Landsat、ASTER或MODIS遥感数据；
[0012] S2：根据所述遥感数据或/和研究区域站点观测数据制备遥感地表参数和区域气象参数，其中，
[0013] 所述遥感地表参数包括：归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts、叶面积指数LAI；
[0014] 所述区域气象参数包括：空气温度Ta和相对湿度RH；
[0015] S3：根据所述遥感地表和区域气象参数进行净辐射通量Rn反演；
[0016] S4：根据所述叶面积指数LAI、空气温度Ta以及净辐射通量Rn进行土壤热通量G反演；
[0017] S5：根据所述植被覆盖度f和地表温度Ts的理论二维空间以及感热通量H估算的温度廓线方程进行辐射-对流阻抗rae反演；
[0018] 所述感热通量H估算的温度廓线方程为：
[0019]
[0020] 其中，ρ为空气密度；Cp为空气定压比热；Taero为动力学温度；Ta为空气温度；ra为空气动力学阻抗；Ts为地表温度；rex为调和地表温度和空气动力学温度的附加阻抗；rae为辐射-对流阻抗；
[0021] 基于所述净辐射通量Rn、土壤热通量G、值被覆盖度f和地表温度Ts的理论二维空间，感热通量H估算公式为：
[0022] H＝(Rn-G)-LE＝(Rn-G)-(1-WDI)Ep (2)
[0023] 其中，LE为潜热通量；WDI为水分亏缺指数；Ep为潜在蒸散；
[0024] 所述潜在蒸散Ep采用下述P-T公式：
[0025]
[0026] 其中，Δ为饱和水汽压曲线斜率，γ为干湿表常数；
[0027] 所述水分亏缺指数WDI计算公式为：
[0028]
[0029] 其中，a为植被覆盖度f和地表温度Ts的理论二维空间干边同估算点的温差；b为植被覆盖度f和地表温度Ts的理论二维空间湿边同估算点的温差；计算公式分别为：
[0030]
[0031] 其中，
[0032]
[0033] 其中，Ta为空气温度；VPD为水汽亏缺，其由空气相对湿度RH计算得到；rcp为植被最小冠层阻抗；rcx为冠层最大阻抗；
[0034] 联立上述公式(1)(2)(3)(4)(5)(6)可得辐射-对流阻抗rae方程：
[0035]
[0036] 其中，c4、c3、c2、c1以及c0为联立公式(1)(2)(3)(4)(5)(6)所得系数，所述辐射-对流阻抗rae为辐射-对流阻抗rae方程(7)中实数解。
[0037] S6：根据所述辐射-对流阻抗rae以及感热通量H估算的温度廓线方程(1)和感热通量H估算公式(2)进行区域地表感热H/潜热LE通量反演，其中，在极端干燥情况下，所述辐射-对流阻抗rae为所述实数解中最小值。
[0038] 优选的，S3中净辐射通量Rn计算公式为：
[0039] Rn＝[Rsd×(1.0-albedo)+Emiss×Rld-Rlu]
[0040] 其中，Rsd为短波入射辐射，albedo为反照率，Emiss为地表比辐射率，Rld为大气下行辐射，Rlu为地表上行辐射；
[0041] 所述短波入射辐射Rsd、大气下行辐射Rld和地表上行辐射Rlu计算公式分别为：
[0042] Rsd＝Gsc×dr×Tsw×cosθ
[0043]
[0044] Rlu＝δ×Emiss×Ts4
[0045] 其中，Gsc为太阳常数，大小为1367瓦/平方米；dr为日地距离修订系数；Tsw为大气透过率；θ为太阳天顶角；ea为水汽压；Ta为空气温度；δ为斯特凡波尔兹曼常数；Ts为地表温度，Emiss为地表比辐射率；
[0046] 所述日地距离修订系数dr和所述大气透过率Tsw计算公式为：
[0047]
[0048]
[0049] 其中，nDays为一年中第n天；PI等于3.1415；W为大气含水量，单位为mm；Kt为浊度系数，范围在0到1之间；θ为太阳天顶角；Pa为大气压强。
[0050] 优选的，S4中，所述土壤热通量G反演过程如下：当所述叶面积指数LAI≥0.5时，所述土壤热通量G估算公式为：
[0051] G＝[0.05+0.18exp(-0.521×LAI)]×Rn
[0052] 当所述叶面积指数LAI＜0.5时，所述土壤热通量G估算公式为：
[0053] G＝1.80(Ts-273.15)+0.084×Rn
[0054] 其中，Ts为地表温度。一种实现上述方法的反演系统，包括：初始化模块，遥感地表参数和区域气象参数制备模块、净辐射通量Rn反演模块、土壤热通量G反演模块、辐射-对流阻抗rae反演模块以及区域地表感热H/潜热LE通量反演模块；其中，
[0055] 所述初始化模块被配置为确定研究区域和遥感数据，所述遥感数据为Landsat、ASTER或MODIS遥感数据；
[0056] 所述遥感地表参数和区域气象参数制备模块被配置为根据所述遥感数据制备遥感地表参数和区域气象参数，所述遥感地表参数包括：归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts、叶面积指数LAI，所述区域气象参数包括：空气温度Ta和相对湿度RH；
[0057] 所述净辐射通量Rn反演模块被配置为根据所述遥感地表和区域气象参数确定净辐射通量Rn；
[0058] 所述土壤热通量G反演模块被配置为根据所述叶面积指数LAI、空气温度Ta以及净辐射通量Rn确定土壤热通量G；
[0059] 所述辐射-对流阻抗rae反演模块被配置为根据方程(7)确定所述辐射-对流阻抗rae；
[0060] 所述区域地表感热H/潜热LE通量反演模块被配置为根据所述辐射-对流阻抗rae以及感热通量H估算的温度廓线方程(1)和感热通量H估算公式(2)确定区域地表感热H/潜热LE通量，其中，在极端干燥情况下，所述辐射-对流阻抗rae为所述辐射-对流阻抗rae方程(7)实数解中最小值。
[0061] 优选的，所述净辐射通量Rn反演模块根据短波入射辐射Rsd、反照率albedo、地表比辐射率Emiss、大气下行辐射Rld以及地表上行辐射Rlu计算净辐射通量Rn，计算公式为：
[0062] Rn＝[Rsd×(1.0-albedo)+Emiss×Rld-Rlu]
[0063] 所述短波入射辐射Rsd、大气下行辐射Rld和地表上行辐射Rlu计算公式分别为：
[0064] Rsd＝Gsc×dr×Tsw×cosθ
[0065]
[0066] Rlu＝δ×Emiss×Ts4
[0067] 其中，Gsc为太阳常数，大小为1367瓦/平方米；dr为日地距离修订系数；Tsw为大气透过率；θ为太阳天顶角；ea为水汽压；Ta为空气温度；δ为斯特凡波尔兹曼常数；Ts为地表温度，Emiss为地表比辐射率；
[0068] 所述日地距离修订系数dr和所述大气透过率Tsw计算公式为：
[0069]
[0070]
[0071] 其中，nDays为一年中第n天；PI等于3.1415；W为大气含水量，单位为mm；Kt为浊度系数，范围在0到1之间；θ为太阳天顶角；Pa为大气压强。
[0072] 优选的，当所述叶面积指数LAI≥0.5时，所述净辐射通量Rn反演模块中土壤热通量G估算公式为：
[0073] G＝[0.05+0.18exp(-0.521×LAI)]×Rn
[0074] 当所述叶面积指数LAI＜0.5时，所述净辐射通量Rn反演模块中土壤热通量G估算公式为：
[0075] G＝1.80(Ts-273.15)+0.084×Rn。
[0076] 本发明的有益效果如下：
[0077] 提供基于遥感数据的区域地表感热/潜热通量反演遥感方法及系统，在不依赖风速的情况下，利用地表感热通量估算的廓线方程和地表温度和植被覆盖度所构成的理论特征空间四个角点方程，可同时实现在无资料区和复杂非均匀下垫面的地表情况下，求解地表感热通量和辐射-对流阻抗具有区域应用的有实用性和可操作性，满足区域水资源管理、农业用水效率评价、区域旱情监测、农业节水等方面的应用领域的需求。
附图说明
[0078] 下面结合附图对本发明的具体实施方式作进一步详细的说明。
[0079] 图1示出本发明所述方法流程图。
[0080] 图2示出植被覆盖度f和地表温度Ts的理论二维空间。
[0081] 图3示出本发明所述系统结构示意图。
[0082] 图4示出黑河地区地表感热H和地表潜热LE通量估算结果同观测结果的对比。
具体实施方式
[0083] 为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。
[0084] 如图1所示，基于遥感数据的区域地表感热/潜热通量反演遥感方法包括以下步骤：
[0085] S1：根据研究区域所需，下载Landsat、ASTER或MODIS遥感卫星产品数据；其中，Landsat和ASTER满足分辨率为100米的通量反演计算，MODIS遥感产品数据满足分辨率为1km的通量反演计算。
[0086] S2：根据所述遥感数据制备遥感地表参数和区域气象参数，所述遥感地表参数包括：归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts、叶面积指数LAI；所述区域气象参数包括：空气温度Ta和相对湿度RH。
[0087] 基于Landsat卫星遥感数据的归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts以及叶面积指数LAI可参照Allen,R.G.,M.Tasumi,and R.Trezza(2007),Satellite-based energy balance for mapping evapotranspiration with internalized calibration(METRIC)-Model,Journal of Irrigation and Drainage Engineering-Asce,133(4),380-394.的文献。
[0088] 基于ASTER卫星遥感数据的归一化植被指数、植被覆盖度、反照率、比辐射率、地表温度、叶面积指数可参照French,A.N.,et al.(2005),Surface energy fluxes with the Advanced Spaceborne Thermal Emission and Reflection radiometer(ASTER)at the Iowa 2002SMACEX site(USA)(vol 99,pg 55,2005),Remote Sensing of Environment,99(4),471-471.的文献。
[0089] 基于MODIS卫星遥感数据，归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts以及叶面积指数LAI等遥感参数产品可以直接计算获得。
[0090] 具体的，归一化植被指数NDVI的提取，运用MODIS遥感产品数据第一波段和第二波段的反射率进行计算，计算方法为：
[0091]
[0092] 其中，b1和b2分别为MODIS遥感产品数据第一和第二波段的波段反射率。
[0093] 植被覆盖度f为特定地域植物垂直投影面积与该地域面积之比，使用归一化植被指数NDVI计算如下：
[0094]
[0095] 其中，NDVImax和NDVImin分别为研究区内NDVI的最大值和最小值。
[0096] 反照率albedo的遥感反演，采用Liang,S.L.,2001.Narrowband to broadband conversions of land surface albedo I Algorithms.Remote Sensing of Environment,76(2):213-238.提出的加权方法：
[0097] albedo＝0.160×b1+0.291×b2+0.243×b3+0.116×b4+0.112×b5+0.018×b7-0.0015
[0098] 其中，b1，b2，b3，b4，b5，b7分别为MODIS波段1，2，3，4，5，7的波段反射率。
[0099] 比辐射率Emiss的反演如下：
[0100] 对于裸土区，归一化植被指数NDVI小于0.2，
[0101] Emiss＝0.9825-0.051×b1
[0102] 对于植被和裸土混合区，归一化植被指数NDVI介于0.2和0.5之间，[0103] Emiss＝0.971-0.018×f
[0104] 对于植被区，归一化植被指数NDVI大于0.5，
[0105] Emiss＝0.990
[0106] 对于水体，
[0107] Emiss＝0.995
[0108] 地表温度Ts：
[0109] Ts＝DNTS×0.02
[0110] 其中，DNTS为MODIS遥感产品数据中MOD11A1产品中的白天地表温度数据层。
[0111] 叶面积指数LAI:
[0112] LAI＝Lai1km×0.1
[0113] 其中，Lai1km为MODIS遥感产品数据中MOD15A2产品中的叶面积指数数据层。
[0114] 所述区域气象参数中空气温度Ta和相对湿度RH的确定，如果选择基于Landsat遥感产品数据，则需要获得研究区域内站点观测数据的支持，利用研究区域内的站点观测数据进行空间插值得到区域的空气温度Ta和相对湿度RH；如果选择基于MODIS遥感数据产品，则不需要区域内的气象观测数据支持，MODIS遥感卫星的MOD07L2产品提供卫星过境时刻的大气廓线产品，可以从MODIS大气廓线产品中获得近地层的压强Pa、空气温度Ta和相对湿度RH。
[0115] 具体的，当选择基于MODIS遥感数据产品时，首先要对MOD07L2数据进行处理，采用MOD07L2产品中20个层的空气温度Ta和相对湿度RH中最近地表层的空气温度Ta和相对湿度RH参与计算，对MOD07L2的20个层从地表层进行循环计算，直到取得有效值为止。
[0116] 大气压强Pa的计算如下(单位为kPa)：
[0117] Pa＝0.01×SurfacePressure
[0118] 其中，SurfacePressure为MOD07L2中的压强数据层；
[0119] 空气温度Ta的计算如下(单位为K)：
[0120] Ta＝(AirTemperature+15000)×0.01
[0121] 其中，AirTemperature为MOD07L2中的空气温度数据层；
[0122] 露点温度Tdew的计算如下(单位为℃)：
[0123] Tdew＝(Mositure+15000)×0.01-273.15
[0124] 其中，Mositure为MOD07L2中的露点温度数据层。
[0125] 水汽压ea的计算如下(单位为hpa)：
[0126]
[0127] 大气含水量W的计算如下(单位为mm)：
[0128] W＝0.01×WaterVapor
[0129] 其中，WaterVapor为MOD07L2中的大气含水量数据层。
[0130] S3：根据所述遥感地表和区域气象参数进行净辐射通量Rn反演；
[0131] 具体的，净辐射通量Rn计算公式为：
[0132] Rn＝[Rsd×(1.0-albedo)+Emiss×Rld-Rlu]
[0133] 其中，Rsd为短波入射辐射，albedo为反照率，Emiss为地表比辐射率，Rld为大气下行辐射，Rlu为地表上行辐射；
[0134] 进一步，所述短波入射辐射Rsd、大气下行辐射Rld和地表上行辐射Rlu计算公式分别为：
[0135] Rsd＝Gsc×dr×Tsw×cosθ
[0136]
[0137] Rlu＝δ×Emiss×Ts4
[0138] 其中，Gsc为太阳常数，大小为1367瓦/平方米；dr为日地距离修订系数；Tsw为大气透过率；θ为太阳天顶角；ea为水汽压；Ta为空气温度；δ为斯特凡波尔兹曼常数；Ts为地表温度，Emiss为地表比辐射率。
[0139] 进一步，所述日地距离修订系数dr和所述大气透过率Tsw计算公式为：
[0140]
[0141]
[0142] 其中，nDays为一年中第n天；PI等于3.1415；W为大气含水量，单位为mm；Kt为浊度系数，范围在0到1之间；θ为太阳天顶角；Pa为大气压强。
[0143] S4：根据所述叶面积指数LAI、空气温度Ta以及净辐射通量Rn进行土壤热通量G反演；
[0144] 具体的，当所述叶面积指数LAI≥0.5时，所述土壤热通量G估算公式为：
[0145] G＝[0.05+0.18exp(-0.521×LAI)]×Rn
[0146] 当所述叶面积指数LAI＜0.5时，所述土壤热通量G估算公式为：
[0147] G＝1.80(Ts-273.15)+0.084×Rn
[0148] 其中，在极端干燥情况下，所述辐射-对流阻抗rae为所述实数解中最小值。
[0149] S5：根据所述值被覆盖度f和地表温度Ts的二维空间以及感热通量H估算的温度廓线方程进行辐射-对流阻抗rae反演；
[0150] 具体的，所述感热通量H估算的温度廓线方程为：
[0151]
[0152] 其中，ρ为空气密度；Cp为空气定压比热；Taero为动力学温度；Ta为空气温度；ra为空气动力学阻抗；Ts为地表温度；rex为调和地表温度和空气动力学温度的附加阻抗；rae为辐射-对流阻抗。
[0153] 基于所述净辐射通量Rn、土壤热通量G、值被覆盖度f和地表温度Ts的二维空间，感热通量H估算公式为：
[0154] H＝(Rn-G)-LE＝(Rn-G)-(1-WDI)Ep (2)
[0155] 其中，LE为潜热通量；WDI为水分亏缺指数；Ep为潜在蒸散。
[0156] 所述潜在蒸散Ep采用下述P-T公式：
[0157]
[0158] 其中，Δ为饱和水汽压曲线斜率，γ为干湿表常数。
[0159] 所述水分亏缺指数WDI计算公式为：
[0160]
[0161] 其中，如图2所示，a为植被覆盖度f和地表温度Ts二维空间干边同估算点的温差；b为植被覆盖度f和地表温度Ts二维空间湿边同估算点的温差；计算公式分别为：
[0162]
[0163] 其中，
[0164]
[0165] 联立上述公式(1)(2)(3)(4)(5)(6)可得：
[0166]
[0167] 其中，c4、c3、c2、c1以及c0为联立公式(1)(2)(3)(4)(5)(6)所得系数，所述辐射-对流阻抗rae为公式(7)中实数解。
[0168] S6：根据所述辐射-对流阻抗rae进行区域地表感热H/潜热LE通量反演；
[0169] 具体的，根据所述公式(1)(2)估算所述区域地表感热通量H和区域地表感热潜热通量LE；其中，在极端干燥情况下，所述辐射-对流阻抗rae为所述实数解中最小值。
[0170] 如图3所示，基于遥感数据的区域地表感热/潜热通量的反演系统，包括：初始化模块，遥感地表参数和区域气象参数制备模块、净辐射通量Rn反演模块、土壤热通量G反演模块、辐射-对流阻抗rae反演模块以及区域地表感热H/潜热LE通量反演模块。
[0171] 其中，所述初始化模块被配置为确定研究区域和遥感数据，所述遥感数据为Landsat、ASTER或MODIS遥感数据；
[0172] 所述遥感地表参数和区域气象参数制备模块被配置为根据所述遥感数据或/和研究区域站点观测数据制备遥感地表参数和区域气象参数；所述遥感地表参数包括：归一化植被指数NDVI、植被覆盖度f、反照率albedo、地表比辐射率Emiss、地表温度Ts、叶面积指数LAI，所述区域气象参数包括：空气温度Ta和相对湿度RH；
[0173] 具体的，选择使用Landsat遥感数据时，需要获得研究区的站点观测数据，选择使用MODIS遥感数据时，不需要区域内的气象观测数据支持。
[0174] 所述净辐射通量Rn反演模块被配置为根据所述遥感地表和区域气象参数确定净辐射通量Rn；
[0175] 所述土壤热通量G反演模块被配置为根据所述叶面积指数LAI、空气温度Ta以及净辐射通量Rn确定土壤热通量G；
[0176] 所述辐射-对流阻抗rae反演模块被配置为根据公式(7)确定所述辐射-对流阻抗rae；
[0177] 所述区域地表感热H/潜热LE通量反演模块被配置为根据所述辐射-对流阻抗rae确定区域地表感热H/潜热LE通量。
[0178] 实施例
[0179] 以黑河流域中游为研究区，使用ASTER卫星遥感数据基于本发明方法反演了地表感热H/潜热LE通量，并将反演结果同研究区内的涡度相关仪器测量的感热H/潜热LE通量值进行对比。
[0180] 具体的，如图4所示，地表感热H/潜热LE通量估算的均方根误差为50W/m2和60W/m2，这一结果和其他模型方法在这一区域的验证结果相近，说明本发明提出的模型方法精度较高，结果合理可信。
[0181] 显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。
