[发明专利]一种联合遥感和气象数据获取全天候蒸散发的方法

权利要求书

1.一种联合遥感和气象数据获取全天候蒸散发的方法，其特征在于，包括以下步骤：

步骤1、采集遥感数据：遥感数据包括1km空间分辨率的地表温度产品、地表反照率产品、8天合成的叶面积指数产品和16天合成归一化植被指数产品；遥感数据来自MODIS卫星数据产品；

步骤2、采集气象数据：气象数据包括气温、风速、短波辐射和比湿这四个气象要素；气象数据来自CLDAS大气驱动场产品；

步骤3、经纬度匹配：计算每个MODIS遥感像元归属的CLDAS数据格网位置；

步骤4、估算MODIS遥感像元蒸散发：遍历MODIS遥感像元，判断MODIS的地表温度产品像元是否为有云，根据判断结果，估算MODIS有云像元和MODIS无云像元的蒸散发；

所述步骤3中经纬度匹配中的计算每个MODIS遥感像元归属的CLDAS数据格网位置具体为：

对每个MODIS遥感像元，读取其经纬度(Lat_M，Lon_M)，对每个CLDAS格网气象数据，读取其经纬度(Lat_C，Lon_C)，遍历所有的CLDAS格网，计算两个经纬度的距离d，公式为：

其中，d最小值出现的格网即为MODIS遥感像元归属CLDAS的格网；

所述步骤4中估算MODIS有云像元的蒸散发具体为：当MODIS的地表温度产品像元有云，其值为空值，这种情况下，根据步骤3所述方法找到该像元所在的CLDAS归属格网，读取该格网的风速，计算摩擦风速，公式为：

其中，K是von Karman常数，取值为0.41；z是CLDAS风速的观测高度，其值为10m，d是零平面位移，z₀是表面粗糙长度；植被高度为1m，d＝0.67，z₀＝0.1；

计算2m处风速u₂，公式为：

计算空气阻抗r_a，公式为：

其中，Z_0m和Z_0h分别是动量传输和能量传输粗糙度长度；植被高度为1m，Z_0m＝0.1，Z_0h＝0.01；

地表阻抗r_s的计算公式为:

其中，u_s是地表粗糙度影响最小的高度处的风速，利用叶面积指数，叶子直径和植被高度进行估算；植被高度为1m，叶子直径为0.1m，u_s近似为：

其中，LAI是叶面积指数，从MODIS数据直接读取；MODIS有云像元蒸散发ET的计算公式为：

其中，Δ是饱和水汽压随温度变化的斜率，计算公式为：

其中，T_a是气温，从CLDAS数据中读取；R_n是净辐射，利用CLDAS数据读取下行短波辐射S_d、气温T_a以及MODIS的地表反照率albedo来估算，G是土壤热通量，表示为R_n的函数；ρ是空气密度，取值1.29kg/m³；C_p是空气定压比热，取值1004J/(kg·K)；e_s(T_a)是空气温度为T_a时的实际水汽压；e_a为饱和水汽压；γ是干湿表常数，取值为0.066；净辐射R_n用以下公式估算：

R_n＝(1-albedo)S_d+εε_aσT_a⁴-εσT_a⁴ (9)

其中，σ是Stefan-Boltzmann常数，取值为5.67×10^-8；ε_a是天空发射率，表示为气温的函数ε_a＝9.2×10^-6×T_a²；ε是地表发射率，从MODIS温度产品中读取；

土壤热通量G的计算公式为：

G＝[0.05+0.265×(1-FVC)]R_n (10)

其中，FVC是植被覆盖度，根据MODIS的归一化植被指数NDVI进行计算：

其中，NDVI直接从MODIS数据中读取。

2.根据权利要求1所述的联合遥感和气象数据获取全天候蒸散发的方法，其特征在于，所述步骤4中估算MODIS无云像元的蒸散发具体为：

首先计算干燥裸土的温度T_s,max和受胁迫植被的温度T_c,max，计算公式为：

其中，α_s和α_c分别是干燥裸土和受胁迫植被反照率，分别取值0.3和0.2；S_d为短波辐射，从CLDAS数据中读取；ε_s和ε_c分别为干燥裸土和受胁迫植被的发射率，分别取值0.95和0.98；σ是Stefan-Boltzmann常数，取值为5.67×10^-8；T_a是气温，从CLDAS数据中读取；ρ是空气密度，取值1.29kg/m³；C_p是空气定压比热，取值1004J/(kg·K)；ε_a是天空发射率，表示为气温的函数ε_a＝9.2×10^-6T_a²；c是裸土条件下联系土壤热通量与净辐射的系数取值0.315；r_a,s和r_a,c分别是干燥裸土和受水分胁迫的全植被覆盖对应的阻抗，它们均表达为风速的函数，风速直接由CLDAS数据提供；

其中，和分别是动量传输和能量传输的大气稳定度函数；z_0m和z_0h分别是动量传输和能量传输粗糙度长度，对于r_a,s计算，d取值0，z_0m取值0.005，z_0h取值0.0005；对于r_a,c计算，d取值0.67，z_0m取值0.1，z_0h取值0.01；

对于MODIS无云像元，蒸散发ET估算公式为：