本發(fā)明屬于遙感技術(shù)領(lǐng)域，具體地說(shuō)，涉及一種聯(lián)合遙感和氣象數(shù)據(jù)獲取全天候蒸散發(fā)的方法。

背景技術(shù)：

蒸散發(fā)是陸表能量平衡的重要組成部分，對(duì)于氣候變化、水文循環(huán)和生態(tài)環(huán)境等諸多領(lǐng)域具有十分重要的意義。傳統(tǒng)獲取蒸散發(fā)的方法主要是基于點(diǎn)的觀測(cè)，但是點(diǎn)觀測(cè)數(shù)據(jù)的空間代表性較差，難以反映區(qū)域蒸散發(fā)的時(shí)空變化。目前在遙感蒸散發(fā)監(jiān)測(cè)中，主要是利用光學(xué)與熱紅外遙感數(shù)據(jù)。然而，光學(xué)與熱紅外遙感要求晴空無(wú)云的天氣條件，一旦遇到云的情況則無(wú)法獲取正常的地表參數(shù)，從而導(dǎo)致無(wú)法反演全天候的蒸散發(fā)信息。另一方面，目前國(guó)內(nèi)外已有部分格網(wǎng)氣象數(shù)據(jù)產(chǎn)品發(fā)布，這些數(shù)據(jù)產(chǎn)品能夠彌補(bǔ)光學(xué)與熱紅外遙感有云像元無(wú)法獲取蒸散發(fā)信息的缺陷，為反演全天候蒸散發(fā)提供了新的信息來(lái)源。

現(xiàn)有的蒸散發(fā)獲取方法主要分為傳統(tǒng)的地表觀測(cè)方法和遙感反演方法。其中，傳統(tǒng)的地表觀測(cè)方法通常是基于單點(diǎn)、小尺度的測(cè)量，其缺點(diǎn)是難以推廣到區(qū)域尺度，無(wú)法滿足當(dāng)前應(yīng)用各領(lǐng)域?qū)^(qū)域尺度蒸散發(fā)信息的要求；遙感反演的方法主要是利用光學(xué)與熱紅外信息，它最大的缺點(diǎn)是容易受到天氣條件的影響，無(wú)法獲取全天候的蒸散發(fā)信息。

技術(shù)實(shí)現(xiàn)要素：

有鑒于此，本發(fā)明針對(duì)現(xiàn)有技術(shù)中存在的不能夠獲取全天候蒸散發(fā)的問(wèn)題，提供了一種聯(lián)合遙感和氣象數(shù)據(jù)獲取全天候蒸散發(fā)的方法，本發(fā)明聯(lián)合現(xiàn)有遙感數(shù)據(jù)和格網(wǎng)氣象數(shù)據(jù)產(chǎn)品直接獲取全天候蒸散發(fā)，能夠有效地為區(qū)域尺度氣候、水文和生態(tài)環(huán)境等領(lǐng)域研究和應(yīng)用提供實(shí)時(shí)、可靠的參數(shù)，具有十分重要的應(yīng)用價(jià)值。

為了解決上述技術(shù)問(wèn)題，本發(fā)明公開(kāi)了一種聯(lián)合遙感和氣象數(shù)據(jù)獲取全天候蒸散發(fā)的方法，包括以下步驟：

步驟1、采集遙感數(shù)據(jù)：遙感數(shù)據(jù)包括1km空間分辨率的地表溫度產(chǎn)品、地表反照率產(chǎn)品、8天合成的葉面積指數(shù)產(chǎn)品、16天合成歸一化植被指數(shù)產(chǎn)品；遙感數(shù)據(jù)來(lái)自modis衛(wèi)星數(shù)據(jù)產(chǎn)品；

步驟2、采集氣象數(shù)據(jù)：氣象數(shù)據(jù)包括氣溫、風(fēng)速、短波輻射和比濕這四個(gè)氣象要素；氣象數(shù)據(jù)來(lái)自cldas大氣驅(qū)動(dòng)場(chǎng)產(chǎn)品；

步驟3、經(jīng)緯度匹配：計(jì)算每個(gè)modis遙感像元?dú)w屬的cldas數(shù)據(jù)格網(wǎng)位置；

步驟4、估算modis遙感像元蒸散發(fā)：遍歷modis遙感像元，判斷modis的地表溫度產(chǎn)品像元是否為有云，根據(jù)判斷結(jié)果，估算modis有云像元和modis無(wú)云像元的蒸散發(fā)。

進(jìn)一步地，步驟3中經(jīng)緯度匹配中的計(jì)算每個(gè)modis遙感像元?dú)w屬的cldas數(shù)據(jù)格網(wǎng)位置具體為：

對(duì)每個(gè)modis遙感像元，讀取其經(jīng)緯度(lat_m，lon_m)，對(duì)每個(gè)cldas格網(wǎng)氣象數(shù)據(jù)，讀取其經(jīng)緯度(lat_c，lon_c)，遍歷所有的cldas格網(wǎng)，計(jì)算兩個(gè)經(jīng)緯度的距離d，公式為：

其中，d最小值出現(xiàn)的格網(wǎng)即為modis遙感像元?dú)w屬cldas的格網(wǎng)。

進(jìn)一步地，步驟4中估算modis有云像元的蒸散發(fā)具體為：當(dāng)modis的地表溫度產(chǎn)品像元有云，其值為空值，這種情況下，根據(jù)步驟3所述方法找到該像元所在的cldas歸屬格網(wǎng)，讀取該格網(wǎng)的風(fēng)速，計(jì)算摩擦風(fēng)速，公式為：

其中，k是vonkarman常數(shù)，取值為0.41；z是cldas風(fēng)速的觀測(cè)高度，其值為10m。d是零平面位移，z0是表面粗糙長(zhǎng)度；植被高度為1m，則d＝0.67，z0＝0.1；

計(jì)算2m處風(fēng)速u2，公式為：

計(jì)算空氣阻抗ra，公式為：

其中，z0m和z0h分別是動(dòng)量傳輸和能量傳輸粗糙度長(zhǎng)度；植被高度為1m，則z0m＝0.1，z0h＝0.01；

地表阻抗rs的計(jì)算公式為:

其中，us是地表粗糙度影響最小的高度處的風(fēng)速，利用葉面積指數(shù)，葉子直徑和植被高度進(jìn)行估算；當(dāng)植被高度為1m，葉子直徑為0.1m，則us可以近似為：

其中，lai是葉面積指數(shù)，從modis數(shù)據(jù)直接讀?。籱odis有云像元蒸散發(fā)et的計(jì)算公式為：

其中，δ是飽和水汽壓隨溫度變化的斜率，計(jì)算公式為：

其中，ta是氣溫，從cldas數(shù)據(jù)中讀??；rn是凈輻射，利用cldas數(shù)據(jù)讀取下行短波輻射sd、氣溫ta以及modis的地表反照率albedo來(lái)估算，g是土壤熱通量，表示為rn的函數(shù)；ρ是空氣密度，取值1.29kg/m³；cp是空氣定壓比熱，取值1004j/(kg·k)；es(ta)是空氣溫度為ta時(shí)的實(shí)際水汽壓；ea為飽和水汽壓；γ是干濕表常數(shù)，取值為0.066；凈輻射rn用以下公式估算：

rn＝(1-albedo)sd+εεaσta⁴-εσta⁴(9)

其中，σ是stefan-boltzmann常數(shù)，取值為5.67×10^-8；εa是天空發(fā)射率，表示為氣溫的函數(shù)εa＝9.2×10^-6×ta²；ε是地表發(fā)射率，從modis溫度產(chǎn)品中讀??；

土壤熱通量g的計(jì)算公式為：

g＝[0.05+0.265×(1-fvc)]rn(10)

其中，fvc是植被覆蓋度，根據(jù)modis的歸一化植被指數(shù)ndvi進(jìn)行計(jì)算：

其中，ndvi直接從modis數(shù)據(jù)中讀取。

進(jìn)一步地，步驟4中估算modis無(wú)云像元的蒸散發(fā)具體為：

首先計(jì)算干燥裸土的溫度ts,max和受脅迫植被的溫度tc,max，計(jì)算公式為：

其中，αs和αc分別是干燥裸土和受脅迫植被反照率，分別取值0.3和0.2；sd為短波輻射，從cldas數(shù)據(jù)中讀取；εs和εc分別為干燥裸土和受脅迫植被的發(fā)射率，分別取值0.95和0.98；σ是stefan-boltzmann常數(shù)，取值為5.67×10^-8；ta是氣溫，從cldas數(shù)據(jù)中讀??；ρ是空氣密度，取值1.29kg/m³；cp是空氣定壓比熱，取值1004j/(kg·k)；εa是天空發(fā)射率，表示為氣溫的函數(shù)(εa＝9.2×10^-6ta²)；c是裸土條件下聯(lián)系土壤熱通量與凈輻射的系數(shù)取值0.315；ra,s和ra,c分別是干燥裸土和受水分脅迫的全植被覆蓋對(duì)應(yīng)的阻抗，它們均表達(dá)為風(fēng)速的函數(shù)，風(fēng)速直接由cldas數(shù)據(jù)提供；

其中，和分別是動(dòng)量傳輸和能量傳輸?shù)拇髿夥€(wěn)定度函數(shù)；z0m和z0h分別是動(dòng)量傳輸和能量傳輸粗糙度長(zhǎng)度，對(duì)于ra,s計(jì)算，d取值0，z0m取值0.005，z0h取值0.0005；對(duì)于ra,c計(jì)算，d取值0.67，z0m取值0.1，z0h取值0.01；

對(duì)于modis無(wú)云像元，蒸散發(fā)et估算公式為：

與現(xiàn)有技術(shù)相比，本發(fā)明可以獲得包括以下技術(shù)效果：

1)本發(fā)明充分利用現(xiàn)有光學(xué)與熱紅外遙感數(shù)據(jù)，并充分利用格網(wǎng)氣象數(shù)據(jù)信息，能夠彌補(bǔ)現(xiàn)有技術(shù)方法中單純使用遙感數(shù)據(jù)在有云條件下無(wú)法反演蒸散發(fā)的缺陷。

2)本發(fā)明融合遙感和氣象數(shù)據(jù)信息，能夠獲取空間上連續(xù)的全天候蒸散發(fā)，為區(qū)域研究提供完整的蒸散發(fā)數(shù)據(jù)。

3)本發(fā)明滿足了定量遙感領(lǐng)域研究對(duì)地面數(shù)據(jù)信息完整性的需求。

當(dāng)然，實(shí)施本發(fā)明的任一產(chǎn)品并不一定需要同時(shí)達(dá)到以上所述的所有技術(shù)效果。

附圖說(shuō)明

此處所說(shuō)明的附圖用來(lái)提供對(duì)本發(fā)明的進(jìn)一步理解，構(gòu)成本發(fā)明的一部分，本發(fā)明的示意性實(shí)施例及其說(shuō)明用于解釋本發(fā)明，并不構(gòu)成對(duì)本發(fā)明的不當(dāng)限定。在附圖中：

圖1是本發(fā)明聯(lián)合遙感和氣象數(shù)據(jù)獲取全天候蒸散發(fā)的方法的流程圖；

圖2是本發(fā)明2015年4月29日10:00河南省的風(fēng)速(單位：m/s)和氣溫(單位：k)，其中，(a)為風(fēng)速，(b)為氣溫；

圖3是本發(fā)明2015年4月29日河南省modis地表溫度(單位：k)；

圖4是本發(fā)明2015年4月29日蒸散發(fā)(單位：w/m²)；

圖5是本發(fā)明2015年5月12日10:00河南省的風(fēng)速(單位：m/s)和氣溫(單位：k)，其中，(a)為風(fēng)速，(b)為氣溫；

圖6是本發(fā)明2015年5月12日河南省modis地表溫度(單位：k)；

圖7是本發(fā)明2015年5月12日蒸散發(fā)(單位：w/m²)，其中，(a)傳統(tǒng)遙感反演方法，(b)本發(fā)明的方法。

具體實(shí)施方式

以下將配合實(shí)施例來(lái)詳細(xì)說(shuō)明本發(fā)明的實(shí)施方式，藉此對(duì)本發(fā)明如何應(yīng)用技術(shù)手段來(lái)解決技術(shù)問(wèn)題并達(dá)成技術(shù)功效的實(shí)現(xiàn)過(guò)程能充分理解并據(jù)以實(shí)施。

本發(fā)明提供一種聯(lián)合遙感和氣象數(shù)據(jù)獲取全天候蒸散發(fā)的方法，如圖1所示，包括以下步驟：

步驟1、遙感數(shù)據(jù)下載和處理：

本發(fā)明的遙感數(shù)據(jù)來(lái)自modis(moderateresolutionimagingspectroradiometer)衛(wèi)星數(shù)據(jù)產(chǎn)品，主要包括1km空間分辨率的地表溫度(lst)產(chǎn)品(mod11a1)、地表反照率(albedo)產(chǎn)品(mcd43a3)、8天合成的葉面積指數(shù)(lai)產(chǎn)品(mod15a2)、16天合成歸一化植被指數(shù)(ndvi)產(chǎn)品(mod13a2)。數(shù)據(jù)可以從網(wǎng)站免費(fèi)下載(https://ladsweb.nascom.nasa.gov/search/)。

數(shù)據(jù)處理過(guò)程包括：投影、拼接、裁剪等，最后輸出為tiff格式，作為本發(fā)明的輸入數(shù)據(jù)之一。

步驟2、氣象數(shù)據(jù)下載和處理：

本發(fā)明的氣象數(shù)據(jù)來(lái)自中國(guó)氣象局陸面數(shù)據(jù)同化系統(tǒng)(cldas，chinameteorologicaladministrationlanddataassimilationsystem)大氣驅(qū)動(dòng)場(chǎng)產(chǎn)品，主要包括氣溫、風(fēng)速、短波輻射和比濕這四個(gè)氣象要素。這些數(shù)據(jù)為0.0625°×0.0625°等經(jīng)緯度格網(wǎng)產(chǎn)品，可以從網(wǎng)站免費(fèi)下載(http://data.cma.cn/data/detail/datacode/nafp_cldas_rt.html)。數(shù)據(jù)下載后裁剪出研究區(qū)所在區(qū)域，作為本發(fā)明的另一輸入數(shù)據(jù)，無(wú)需做其他處理。

步驟3、經(jīng)緯度匹配：

計(jì)算每個(gè)modis遙感像元?dú)w屬的cldas數(shù)據(jù)格網(wǎng)位置，方法是：

對(duì)每個(gè)modis遙感像元，讀取其經(jīng)緯度(lat_m，lon_m)，對(duì)每個(gè)cldas格網(wǎng)氣象數(shù)據(jù)，讀取其經(jīng)緯度(lat_c，lon_c)，遍歷所有的cldas格網(wǎng)，計(jì)算兩個(gè)經(jīng)緯度的距離d，公式為：

其中，d最小值出現(xiàn)的格網(wǎng)即為modis遙感像元?dú)w屬cldas的格網(wǎng)。

步驟4、估算modis遙感像元蒸散發(fā)：

遍歷modis遙感像元，判斷modis的地表溫度產(chǎn)品像元是否為有云，如果像元有云，其值為空值，這種情況下，根據(jù)步驟3所示方法找到該像元所在的cldas歸屬格網(wǎng)，讀取該格網(wǎng)的風(fēng)速，計(jì)算摩擦風(fēng)速，公式為：

其中，k是vonkarman常數(shù)，取值為0.41。z是cldas風(fēng)速的觀測(cè)高度，其值為10m。d是零平面位移，z0是表面粗糙長(zhǎng)度。假設(shè)植被高度為1m，則d＝0.67，z0＝0.1。

計(jì)算2m處風(fēng)速u2，公式為：

z0m和z0h分別是動(dòng)量傳輸和能量傳輸粗糙度長(zhǎng)度。假設(shè)植被高度為1m，則z0m＝0.1，z0h＝0.01。計(jì)算空氣阻抗ra，公式為：

地表阻抗rs的計(jì)算公式為:

其中，us是地表粗糙度影響最小的高度處的風(fēng)速，可以利用葉面積指數(shù)，葉子直徑和植被高度進(jìn)行估算。假設(shè)植被高度為1m，葉子直徑為0.1m，則us可以近似為：

其中，lai是葉面積指數(shù)，從modis數(shù)據(jù)直接讀取。modis有云像元蒸散發(fā)et的計(jì)算公式為：

其中，δ是飽和水汽壓隨溫度變化的斜率，計(jì)算公式為：

其中，ta是氣溫，從cldas數(shù)據(jù)中讀取。rn是凈輻射，利用cldas數(shù)據(jù)讀取下行短波輻射(sd)、氣溫(ta)以及modis的地表反照率(albedo)來(lái)估算，g是土壤熱通量，可以表示為rn的函數(shù)。ρ是空氣密度(取值1.29kg/m³)；cp是空氣定壓比熱(取值1004j/(kg·k))；es(ta)是空氣溫度為ta(從cldas數(shù)據(jù)讀取)時(shí)的實(shí)際水汽壓，ea為飽和水汽壓；γ是干濕表常數(shù)(取值0.066)。凈輻射rn用以下公式估算：

rn＝(1-albedo)sd+εεaσta⁴-εσta⁴(9)

其中，σ是stefan-boltzmann常數(shù)(5.67×10^-8)；εa是天空發(fā)射率，表示為氣溫的函數(shù)εa＝9.2×10^-6×ta²；ε是地表發(fā)射率，可以從modis溫度產(chǎn)品(lst)中讀取。

土壤熱通量g的計(jì)算公式為：

g＝[0.05+0.265×(1-fvc)]rn(10)

其中，fvc是植被覆蓋度，可以根據(jù)modis的歸一化植被指數(shù)(ndvi)進(jìn)行計(jì)算：

其中，ndvi可直接從modis數(shù)據(jù)中讀取。

對(duì)于modis有云像元，根據(jù)公式(7)可以計(jì)算得到該像元對(duì)應(yīng)的蒸散發(fā)。

對(duì)于modis無(wú)云像元，首先計(jì)算干燥裸土的溫度ts,max和受脅迫植被的溫度tc,max，計(jì)算公式為：

其中，αs和αc分別是干燥裸土和受脅迫植被反照率(分別取值0.3和0.2)；sd為短波輻射，從cldas數(shù)據(jù)中讀??；εs和εc分別為干燥裸土和受脅迫植被的發(fā)射率(分別取值0.95和0.98)；σ是stefan-boltzmann常數(shù)(5.67×10^-8)；ta是氣溫，從cldas數(shù)據(jù)中讀取；ρ是空氣密度，(取值1.29kg/m³)；cp是空氣定壓比熱(取值1004j/(kg·k))；εa是天空發(fā)射率，可以表示為氣溫的函數(shù)(εa＝9.2×10^-6ta²)；c是裸土條件下聯(lián)系土壤熱通量與凈輻射的系數(shù)(取值0.315)。ra,s和ra,c分別是干燥裸土和受水分脅迫的全植被覆蓋對(duì)應(yīng)的阻抗，它們均可以表達(dá)為風(fēng)速的函數(shù)，風(fēng)速可以直接由cldas數(shù)據(jù)提供。

其中，和分別是動(dòng)量傳輸和能量傳輸?shù)拇髿夥€(wěn)定度函數(shù)。z0m和z0h分別是動(dòng)量傳輸和能量傳輸粗糙度長(zhǎng)度，對(duì)于ra,s計(jì)算，d取值0，z0m取值0.005，z0h取值0.0005；對(duì)于ra,c計(jì)算，d取值0.67，z0m取值0.1，z0h取值0.01。

對(duì)于modis無(wú)云像元，蒸散發(fā)et估算公式為：

實(shí)施例1

以河南省為研究區(qū)，圖2是cldas格網(wǎng)氣象數(shù)據(jù)下載后裁剪出的研究區(qū)2015年4月29日10:00氣象數(shù)據(jù)，圖3是同一天modis的地表溫度(lst)產(chǎn)品，從圖3可以看出，這一天很多地方，尤其是東部地區(qū)由于云等不利因素的影響，表現(xiàn)出大片白色空白區(qū)域，沒(méi)有正常的溫度值，呈現(xiàn)出溫度不連續(xù)的狀況，這些地方無(wú)法正常利用遙感方法反演蒸散發(fā)，從而無(wú)法利用傳統(tǒng)遙感反演方法獲取全天候的蒸散發(fā)。圖4是蒸散發(fā)反演結(jié)果。從圖中可以看出，傳統(tǒng)方法在沒(méi)有像元值的地方無(wú)能為力，無(wú)法獲取蒸散發(fā)；本發(fā)明方法能夠有效彌補(bǔ)傳統(tǒng)光學(xué)與熱紅外遙感的這種缺陷，獲取到傳統(tǒng)光學(xué)與熱紅外因云等因素導(dǎo)致無(wú)值的地區(qū)的蒸散發(fā)，實(shí)現(xiàn)全天候的蒸散發(fā)反演。從結(jié)果來(lái)看，本發(fā)明的方法獲得的蒸散發(fā)結(jié)果連續(xù)性良好。

實(shí)施例2

以河南省為研究區(qū)，圖5是cldas格網(wǎng)氣象數(shù)據(jù)下載后裁剪出的研究區(qū)2015年5月12日10:00氣象數(shù)據(jù)。圖6是當(dāng)天modis數(shù)據(jù)，從圖中可以看出，在河南中部地區(qū)有部分區(qū)域?yàn)榭瞻谉o(wú)值地區(qū)。圖7是蒸散發(fā)反演結(jié)果。圖7(a)是傳統(tǒng)遙感方法反演的結(jié)果，從圖中可以看出，傳統(tǒng)方法在沒(méi)有像元值的地方無(wú)法獲取蒸散發(fā)；圖7(b)是本發(fā)明的方法的結(jié)果。

從實(shí)施例2結(jié)果可以看出，本發(fā)明方法能夠有效彌補(bǔ)傳統(tǒng)光學(xué)與熱紅外遙感因?yàn)樵频炔焕麣庀笠蛩貙?dǎo)致的無(wú)值區(qū)域不能反演蒸散發(fā)的問(wèn)題，能夠?qū)崿F(xiàn)全天候的蒸散發(fā)反演，有效獲取全天候蒸散發(fā)。

上述說(shuō)明示出并描述了發(fā)明的若干優(yōu)選實(shí)施例，但如前所述，應(yīng)當(dāng)理解發(fā)明并非局限于本文所披露的形式，不應(yīng)看作是對(duì)其他實(shí)施例的排除，而可用于各種其他組合、修改和環(huán)境，并能夠在本文所述發(fā)明構(gòu)想范圍內(nèi)，通過(guò)上述教導(dǎo)或相關(guān)領(lǐng)域的技術(shù)或知識(shí)進(jìn)行改動(dòng)。而本領(lǐng)域人員所進(jìn)行的改動(dòng)和變化不脫離發(fā)明的精神和范圍，則都應(yīng)在發(fā)明所附權(quán)利要求的保護(hù)范圍內(nèi)。