技术领域
[0001] 本
发明属于遥感技术领域，具体地说，涉及一种联合遥感和气象数据获取全天候蒸散发的方法。
背景技术
[0002] 蒸散发是陆表
能量平衡的重要组成部分，对于
气候变化、
水文循环和生态环境等诸多领域具有十分重要的意义。传统获取蒸散发的方法主要是基于点的观测，但是点观测数据的空间代表性较差，难以反映区域蒸散发的
时空变化。目前在遥感蒸散发监测中，主要是利用光学与热红外遥感数据。然而，光学与热红外遥感要求晴空无
云的天气条件，一旦遇到云的情况则无法获取正常的地表参数，从而导致无法反演全天候的蒸散发信息。另一方面，目前国内外已有部分格网气象数据产品发布，这些数据产品能够弥补光学与热红外遥感有云像元无法获取蒸散发信息的
缺陷，为反演全天候蒸散发提供了新的信息来源。
[0003] 现有的蒸散发获取方法主要分为传统的地表观测方法和遥感反演方法。其中，传统的地表观测方法通常是基于单点、小尺度的测量，其缺点是难以推广到区域尺度，无法满足当前应用各领域对区域尺度蒸散发信息的要求；遥感反演的方法主要是利用光学与热红外信息，它最大的缺点是容易受到天气条件的影响，无法获取全天候的蒸散发信息。
发明内容
[0004] 有鉴于此，本发明针对
现有技术中存在的不能够获取全天候蒸散发的问题，提供了一种联合遥感和气象数据获取全天候蒸散发的方法，本发明联合现有遥感数据和格网气象数据产品直接获取全天候蒸散发，能够有效地为区域尺度气候、水文和生态环境等领域研究和应用提供实时、可靠的参数，具有十分重要的应用价值。
[0005] 为了解决上述技术问题，本发明公开了一种联合遥感和气象数据获取全天候蒸散发的方法，包括以下步骤：
[0006] 步骤1、采集遥感数据：遥感数据包括1km空间
分辨率的地表
温度产品、地表反照率产品、8天合成的叶面积指数产品、16天合成归一化植被指数产品；遥感数据来自MODIS卫星数据产品；
[0007] 步骤2、采集气象数据：气象数据包括气温、
风速、短波
辐射和比湿这四个气象要素；气象数据来自CLDAS大气驱动场产品；
[0008] 步骤3、经纬度匹配：计算每个MODIS遥感像元归属的CLDAS数据格网
位置；
[0009] 步骤4、估算MODIS遥感像元蒸散发：遍历MODIS遥感像元，判断MODIS的地表温度产品像元是否为有云，根据判断结果，估算MODIS有云像元和MODIS无云像元的蒸散发。
[0010] 进一步地，步骤3中经纬度匹配中的计算每个MODIS遥感像元归属的CLDAS数据格网位置具体为：
[0011] 对每个MODIS遥感像元，读取其经纬度(Lat_M，Lon_M)，对每个CLDAS格网气象数据，读取其经纬度(Lat_C，Lon_C)，遍历所有的CLDAS格网，计算两个经纬度的距离d，公式为：
[0012]
[0013] 其中，d最小值出现的格网即为MODIS遥感像元归属CLDAS的格网。
[0014] 进一步地，步骤4中估算MODIS有云像元的蒸散发具体为：当MODIS的地表温度产品像元有云，其值为空值，这种情况下，根据步骤3所述方法找到该像元所在的CLDAS归属格网，读取该格网的风速，计算摩擦风速，公式为：
[0015]
[0016] 其中，K是von Karman常数，取值为0.41；z是CLDAS风速的观测高度，其值为10m。d是零平面位移，z0是表面粗糙长度；植被高度为1m，则d＝0.67，z0＝0.1；
[0017] 计算2m处风速u2，公式为：
[0018]
[0019] 计算空气阻抗ra，公式为：
[0020]
[0021] 其中，Z0M和Z0H分别是动量传输和能量传输粗糙度长度；植被高度为1m，则Z0M＝0.1，Z0H＝0.01；
[0022] 地表阻抗rs的计算公式为:
[0023]
[0024] 其中，us是地表粗糙度影响最小的高度处的风速，利用叶面积指数，叶子直径和植被高度进行估算；当植被高度为1m，叶子直径为0.1m，则us可以近似为：
[0025]
[0026] 其中，LAI是叶面积指数，从MODIS数据直接读取；MODIS有云像元蒸散发ET的计算公式为：
[0027]
[0028] 其中，Δ是饱和水汽压随温度变化的斜率，计算公式为：
[0029]
[0030] 其中，Ta是气温，从CLDAS数据中读取；Rn是净辐射，利用CLDAS数据读取下行短波辐射Sd、气温Ta以及MODIS的地表反照率albedo来估算，G是
土壤热通量，表示为Rn的函数；ρ是空气
密度，取值1.29kg/m3；Cp是空气定压
比热，取值1004J/(kg·K)；es(Ta)是空气温度为Ta时的实际水汽压；ea为饱和水汽压；γ是
干湿表常数，取值为0.066；净辐射Rn用以下公式估算：
[0031] Rn＝(1-albedo)Sd+εεaσTa4-εσTa4 (9)
[0032] 其中，σ是Stefan-Boltzmann常数，取值为5.67×10-8；εa是天空发射率，表示为气-6 2温的函数εa＝9.2×10 ×Ta；ε是地表发射率，从MODIS温度产品中读取；
[0033] 土壤热通量G的计算公式为：
[0034] G＝[0.05+0.265×(1-FVC)]Rn (10)
[0035] 其中，FVC是植被
覆盖度，根据MODIS的归一化植被指数NDVI进行计算：
[0036]
[0037] 其中，NDVI直接从MODIS数据中读取。
[0038] 进一步地，步骤4中估算MODIS无云像元的蒸散发具体为：
[0039] 首先计算干燥裸土的温度Ts,max和受胁迫植被的温度Tc,max，计算公式为：
[0040]
[0041]
[0042] 其中，αs和αc分别是干燥裸土和受胁迫植被反照率，分别取值0.3和0.2；Sd为短波辐射，从CLDAS数据中读取；εs和εc分别为干燥裸土和受胁迫植被的发射率，分别取值0.95和0.98；σ是Stefan-Boltzmann常数，取值为5.67×10-8；Ta是气温，从CLDAS数据中读取；ρ是空气密度，取值1.29kg/m3；Cp是空气定压比热，取值1004J/(kg·K)；εa是天空发射率，表示为-6 2
气温的函数(εa＝9.2×10 Ta )；c是裸土条件下联系土壤热通量与净辐射的系数取值
0.315；ra,s和ra,c分别是干燥裸土和受
水分胁迫的全植被覆盖对应的阻抗，它们均表达为风速的函数，风速直接由CLDAS数据提供；
[0043]
[0044] 其中， 和 分别是动量传输和能量传输的大气稳定度函数；z0m和z0h分别是动量传输和能量传输粗糙度长度，对于ra,s计算，d取值0，z0m取值0.005，z0h取值0.0005；对于ra,c计算，d取值0.67，z0m取值0.1，z0h取值0.01；
[0045] 对于MODIS无云像元，蒸散发ET估算公式为：
[0046]
[0047] 与现有技术相比，本发明可以获得包括以下技术效果：
[0048] 1)本发明充分利用现有光学与热红外遥感数据，并充分利用格网气象数据信息，能够弥补现有技术方法中单纯使用遥感数据在有云条件下无法反演蒸散发的缺陷。
[0049] 2)本发明融合遥感和气象数据信息，能够获取空间上连续的全天候蒸散发，为区域研究提供完整的蒸散发数据。
[0050] 3)本发明满足了定量遥感领域研究对地面数据信息完整性的需求。
[0051] 当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。
附图说明
[0052] 此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性
实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0053] 图1是本发明联合遥感和气象数据获取全天候蒸散发的方法的
流程图；
[0054] 图2是本发明2015年4月29日10:00河南省的风速(单位：m/s)和气温(单位：K)，其中，(a)为风速，(b)为气温；
[0055] 图3是本发明2015年4月29日河南省MODIS地表温度(单位：K)；
[0056] 图4是本发明2015年4月29日蒸散发(单位：W/m2)；
[0057] 图5是本发明2015年5月12日10:00河南省的风速(单位：m/s)和气温(单位：K)，其中，(a)为风速，(b)为气温；
[0058] 图6是本发明2015年5月12日河南省MODIS地表温度(单位：K)；
[0059] 图7是本发明2015年5月12日蒸散发(单位：W/m2)，其中，(a)传统遥感反演方法，(b)本发明的方法。
具体实施方式
[0060] 以下将配合实施例来详细说明本发明的实施方式，藉此对本发明如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。
[0061] 本发明提供一种联合遥感和气象数据获取全天候蒸散发的方法，如图1所示，包括以下步骤：
[0062] 步骤1、遥感数据下载和处理：
[0063] 本发明的遥感数据来自MODIS(Moderate Resolution Imaging Spectroradiometer)卫星数据产品，主要包括1km空间分辨率的地表温度(LST)产品(MOD11A1)、地表反照率(Albedo)产品(MCD43A3)、8天合成的叶面积指数(LAI)产品(MOD15A2)、16天合成归一化植被指数(NDVI)产品(MOD13A2)。数据可以从
网站免费下载(https://ladsweb.nascom.nasa.gov/search/)。
[0064]
数据处理过程包括：投影、拼接、裁剪等，最后输出为tiff格式，作为本发明的输入数据之一。
[0065] 步骤2、气象数据下载和处理：
[0066] 本发明的气象数据来自中国气象局陆面数据同化系统(CLDAS，China Meteorological Administration Land Data Assimilation System)大气驱动场产品，主要包括气温、风速、短波辐射和比湿这四个气象要素。这些数据为0.0625°×0.0625°等经纬度格网产品，可以从网站免费下载(http://data.cma.cn/data/detail/dataCode/NAFP_CLDAS_RT.html)。数据下载后裁剪出研究区所在区域，作为本发明的另一输入数据，无需做其他处理。
[0067] 步骤3、经纬度匹配：
[0068] 计算每个MODIS遥感像元归属的CLDAS数据格网位置，方法是：
[0069] 对每个MODIS遥感像元，读取其经纬度(Lat_M，Lon_M)，对每个CLDAS格网气象数据，读取其经纬度(Lat_C，Lon_C)，遍历所有的CLDAS格网，计算两个经纬度的距离d，公式为：
[0070]
[0071] 其中，d最小值出现的格网即为MODIS遥感像元归属CLDAS的格网。
[0072] 步骤4、估算MODIS遥感像元蒸散发：
[0073] 遍历MODIS遥感像元，判断MODIS的地表温度产品像元是否为有云，如果像元有云，其值为空值，这种情况下，根据步骤3所示方法找到该像元所在的CLDAS归属格网，读取该格网的风速，计算摩擦风速，公式为：
[0074]
[0075] 其中，K是von Karman常数，取值为0.41。z是CLDAS风速的观测高度，其值为10m。d是零平面位移，z0是表面粗糙长度。假设植被高度为1m，则d＝0.67，z0＝0.1。
[0076] 计算2m处风速u2，公式为：
[0077]
[0078] Z0M和Z0H分别是动量传输和能量传输粗糙度长度。假设植被高度为1m，则Z0M＝0.1，Z0H＝0.01。计算空气阻抗ra，公式为：
[0079]
[0080] 地表阻抗rs的计算公式为:
[0081]
[0082] 其中，us是地表粗糙度影响最小的高度处的风速，可以利用叶面积指数，叶子直径和植被高度进行估算。假设植被高度为1m，叶子直径为0.1m，则us可以近似为：
[0083]
[0084] 其中，LAI是叶面积指数，从MODIS数据直接读取。MODIS有云像元蒸散发ET的计算公式为：
[0085]
[0086] 其中，Δ是饱和水汽压随温度变化的斜率，计算公式为：
[0087]
[0088] 其中，Ta是气温，从CLDAS数据中读取。Rn是净辐射，利用CLDAS数据读取下行短波辐射(Sd)、气温(Ta)以及MODIS的地表反照率(albedo)来估算，G是土壤热通量，可以表示为Rn的函数。ρ是空气密度(取值1.29kg/m3)；Cp是空气定压比热(取值1004J/(kg·K))；es(Ta)是空气温度为Ta(从CLDAS数据读取)时的实际水汽压，ea为饱和水汽压；γ是干湿表常数(取值0.066)。净辐射Rn用以下公式估算：
[0089] Rn＝(1-albedo)Sd+εεaσTa4-εσTa4 (9)
[0090] 其中，σ是Stefan-Boltzmann常数(5.67×10-8)；εa是天空发射率，表示为气温的函数εa＝9.2×10-6×Ta2；ε是地表发射率，可以从MODIS温度产品(LST)中读取。
[0091] 土壤热通量G的计算公式为：
[0092] G＝[0.05+0.265×(1-FVC)]Rn (10)
[0093] 其中，FVC是植被覆盖度，可以根据MODIS的归一化植被指数(NDVI)进行计算：
[0094]
[0095] 其中，NDVI可直接从MODIS数据中读取。
[0096] 对于MODIS有云像元，根据公式(7)可以计算得到该像元对应的蒸散发。
[0097] 对于MODIS无云像元，首先计算干燥裸土的温度Ts,max和受胁迫植被的温度Tc,max，计算公式为：
[0098]
[0099]
[0100] 其中，αs和αc分别是干燥裸土和受胁迫植被反照率(分别取值0.3和0.2)；Sd为短波辐射，从CLDAS数据中读取；εs和εc分别为干燥裸土和受胁迫植被的发射率(分别取值0.95和0.98)；σ是Stefan-Boltzmann常数(5.67×10-8)；Ta是气温，从CLDAS数据中读取；ρ是空气密度，(取值1.29kg/m3)；Cp是空气定压比热(取值1004J/(kg·K))；εa是天空发射率，可以表示为气温的函数(εa＝9.2×10-6Ta2)；c是裸土条件下联系土壤热通量与净辐射的系数(取值
0.315)。ra,s和ra,c分别是干燥裸土和受水分胁迫的全植被覆盖对应的阻抗，它们均可以表达为风速的函数，风速可以直接由CLDAS数据提供。
[0101]
[0102] 其中， 和 分别是动量传输和能量传输的大气稳定度函数。z0m和z0h分别是动量传输和能量传输粗糙度长度，对于ra,s计算，d取值0，z0m取值0.005，z0h取值0.0005；对于ra,c计算，d取值0.67，z0m取值0.1，z0h取值0.01。
[0103] 对于MODIS无云像元，蒸散发ET估算公式为：
[0104]
[0105] 实施例1
[0106] 以河南省为研究区，图2是CLDAS格网气象数据下载后裁剪出的研究区2015年4月29日10:00气象数据，图3是同一天MODIS的地表温度(LST)产品，从图3可以看出，这一天很多地方，尤其是东部地区由于云等不利因素的影响，表现出大片白色空白区域，没有正常的温度值，呈现出温度不连续的状况，这些地方无法正常利用遥感方法反演蒸散发，从而无法利用传统遥感反演方法获取全天候的蒸散发。图4是蒸散发反演结果。从图中可以看出，传统方法在没有像元值的地方无能为
力，无法获取蒸散发；本发明方法能够有效弥补传统光学与热红外遥感的这种缺陷，获取到传统光学与热红外因云等因素导致无值的地区的蒸散发，实现全天候的蒸散发反演。从结果来看，本发明的方法获得的蒸散发结果连续性良好。
[0107] 实施例2
[0108] 以河南省为研究区，图5是CLDAS格网气象数据下载后裁剪出的研究区2015年5月12日10:00气象数据。图6是当天MODIS数据，从图中可以看出，在河南中部地区有部分区域为空白无值地区。图7是蒸散发反演结果。图7(a)是传统遥感方法反演的结果，从图中可以看出，传统方法在没有像元值的地方无法获取蒸散发；图7(b)是本发明的方法的结果。
[0109] 从实施例2结果可以看出，本发明方法能够有效弥补传统光学与热红外遥感因为云等不利气象因素导致的无值区域不能反演蒸散发的问题，能够实现全天候的蒸散发反演，有效获取全天候蒸散发。
[0110] 上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、
修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附
权利要求的保护范围内。
