專利名稱::從modis數(shù)據(jù)反演地表溫度方法
技術領域:
:本發(fā)明涉及一種利用對地觀測衛(wèi)星上MODIS傳感器獲得的地面熱輻射信息反演地表溫度的方法,該方法克服了以往反演方法中參數(shù)獲取不同步和復雜難以實用的缺點。能夠應用在氣象、農(nóng)業(yè)、環(huán)境監(jiān)測和旱情監(jiān)測等遙感部門。
背景技術:
:地表溫度是指地表表層的溫度,英文是skintemperature。地表溫度在區(qū)域資源環(huán)境研究中的重要性已經(jīng)使熱紅外遙感成為遙感研究的一個重要領域,目前已經(jīng)開發(fā)了很多實用的地表溫度遙感反演方法,如熱輻射傳輸方程法、劈窗算法、單窗算法和多通道算法。許反演算法是針對具體的傳感器開發(fā)的,例如Qinetal[2001]針對NOAA-AVHRR開發(fā)了一個劈窗算法[QinZ.H.,OlmoGD.,andKarnieliA.,DerivationofsplitwindowalgorithmanditssensitivityanalysisforretrievinglandsurfacetemperaturefromNOAA-advancedveryhighresolutionradiometerdata,Gec^A".i饑,2001,22,655-22,670.],由于AVHRR傳感器上沒有近紅外波段反演大氣水汽含量,只能通過氣象站點或者其他傳感器獲得大氣參數(shù)。1999、2002年搭載MODIS遙感器的對地觀測衛(wèi)星發(fā)射成功,為全球和區(qū)域資源環(huán)境動態(tài)監(jiān)測開辟了又一新的途徑。MODIS是一個擁有36個波段的中分辨率遙感系統(tǒng)(如圖1),每12天可獲得一次全球觀測數(shù)據(jù),其飛行與太陽同步,每天同一區(qū)域至少可獲得晝夜兩景圖像,并且是免費接收,因此非常適合于中大尺度的區(qū)域資源環(huán)境動態(tài)監(jiān)測。在MODIS的36個波段中有8個是熱紅外波段(如表O,因而非常合適于區(qū)域尺度的地表熱量空間差異分析。但是,目前針對MODIS遙感數(shù)據(jù)的地表溫度反演方法不是很多,萬正明等在1996和1997年在MODIS傳感器還沒有搭載衛(wèi)星上天之前,針對這個傳感器提出了兩種反演算法[WanZ.andDozierJ.,Ageneralizedsplit-windowalgorithmforretrievinglandsurfacetemperaturemeasurementfromspace,rra肌Gmsc/,/emoMSe打s.,1996,34:892-905.;WanZ.M.,LiZ..L.,APhysics-BasedAlgorithmforRetrievingland-surfaceemissivityandtemperaturefromEOS/MODISdata,Gmyc/.及e附We&肌,1997,35:980-996.]。其中一個是劈窗算法,只反演地表溫度;另外一個算法是能同時反演地表溫度和發(fā)射率,該算法同時需要白天和晚上的數(shù)據(jù),由于數(shù)據(jù)匹配原因和該算法需要的一些參數(shù)需要其它同時搭載的傳感器或者其它反演產(chǎn)品獲得,這兩個算法被美國宇航局(NASA)采用作為產(chǎn)品算法。但算法比較復雜,一般科研人員難以實現(xiàn),比較復雜。國內許多研究者針對MODIS傳感器的地表溫度反演做過一些研究,比如毛克彪和覃志豪等[毛克彪,覃志豪,施建成,宮鵬,針對MODIS數(shù)據(jù)的劈窗算法研究,武漢大學學報(信息科學版),2005(8):703-708.;毛克彪,覃志豪,施建成,用MODIS影像和劈窗算法反演山東半島的地表溫度,中國礦業(yè)大學學報(自然科學版),2005(1):46-50.]。雖然MODIS衛(wèi)星傳感器搭載在美國衛(wèi)星上,但我國MODIS數(shù)據(jù)地面接收站比較多,比如農(nóng)業(yè)部資源遙感與數(shù)字農(nóng)業(yè)重點室,中國氣象局衛(wèi)星氣象中心,中國科學院地理所等都有接收站點。雖然美國宇航局向全球發(fā)布MODIS溫度產(chǎn)品,其精度在美國本土比較準確,但在世界其它地區(qū)有些地方精度還是不夠,需要進一步提高[WAN,Z.,ZHANG,Y.,ZHANG,Q.andLI,Z.-L.,2002,Validationoftheland-surfacetemperatureproductsretrievedfromTerraModerateResolutionImagingSpectroradiometerdata.'RemoteSensingofEnvironment,2002,83,163-180.]。需要我們根據(jù)世界各地本地的實際情況,進一步研究新算法或者做進一步的校正,以提高精度。表1MODIS遙感器技術參數(shù)<table>tableseeoriginaldocumentpage5</column></row><table><table>tableseeoriginaldocumentpage6</column></row><table><table>tableseeoriginaldocumentpage7</column></row><table>
發(fā)明內容本發(fā)明的目的在于提供一種從遙感數(shù)據(jù)MODIS反演地表溫度的方法,以克服現(xiàn)有地表溫度反演方法復雜,且難以滿足科研人員個人幵發(fā)使用的需要,以及更加充實國內研究人員針對熱紅外遙感器開發(fā)的產(chǎn)品算法,為我國2020之前計劃發(fā)射100顆衛(wèi)星上熱紅外傳感器提供地表溫度反演方法參考,而且還能進一步提高目前農(nóng)業(yè)部業(yè)務運行中地表溫度反演精度,提高旱情監(jiān)測和農(nóng)作物估產(chǎn)精度[已經(jīng)擬定作為農(nóng)業(yè)部農(nóng)情監(jiān)測中的地表溫度反演方法選擇的方法之一]。為實現(xiàn)上述目的,本發(fā)明提供的從遙感數(shù)據(jù)MODIS反演地表溫度方法步驟為第一步,通過大氣水汽含量,計算MODIS第31和32波段的透過率1-1)利用MODIS第19波段和第2波段計算比值T,然后利用這個比值計算大氣水汽含量W[KaufmanY.J.,GaoBo-Cai.,RemoteSensingofWaterVaporintheNearIRfromEOS/MODIS,壓五五Traw.Geosc/.Wewo化Se朋.,1992,30:871掘.〗<formula>formulaseeoriginaldocumentpage8</formula>1-2)利用大氣輻射傳輸模型MODTRAN4模擬計算得到大氣水汽含量與MODIS第31和32波段透過率的關系,將l-l)中計算得到的大氣水汽含量代入下面公式,計算得到第31波段和32波段的透過率(r31,)。<formula>formulaseeoriginaldocumentpage8</formula>第二步、通過植被指數(shù)NDVI(NormalizedDifferenceVegetationIndex)計算MODIS傳感器第31和32波段的發(fā)射率。2-1)讀入MODIS第l波段(近紅外波段)和第2波段(紅光波段),計算<formula>formulaseeoriginaldocumentpage8</formula>其中A表示第1波段,A表示第二波段。2-2)利用NDVI來判斷1平方公里分辨率像元的地物類型,然后根據(jù)美國噴汽推進實驗室(JPL)測試的ASTER波譜庫計算MODIS傳感器第31波段和32波段的發(fā)射率(f31,£32分別表示第31和32波段發(fā)射率)。計算如下當NDVK0時,是水體和雪,f31=0.992,£32=0.988當0<NDVI<0.05,是裸土,s31=0.986,s32=0.991當0.05<NDVI<0.65,為植被和裸土的混合像元。利用PV指數(shù)混合像元中植被的覆蓋的比率[KERR,Y.H.,LAGOUARDE,J.P.andIMBERNON,J.,1992,AccuratelandsurfacetemperatureretrievalfromAVHRRdatawithuseofanimprovedsplitwindowalgorithm.RemoteSensingofEnvironment,1992,41,197-209.]。計算公式PV=(NDVI-0.05)/0.6(式5)s31=0.986*(l-PV)+0.972*PV,s32=0.991*(1-PV)+0.976*PV;當NDVI〉0.65,為植被f31=0.972,f32=0.976第三步簡化輻射傳輸方程3-1)簡化普朗克(Planck)函數(shù)。在能量平衡方程中,每一項都包含普朗克(Planck)函數(shù),A,r=,(e^-1)-1(式6)5"是分譜輻射通量密度,單位是『.m—2.,、A是波長,單位輝;A是普朗克常數(shù)(6.6256xl0-34J*s);c是光速(3xl08m/s);/fc是玻耳茲曼常數(shù)(1.38x10-23J/K);T是絕對溫度(K)。普朗克(Planck)函數(shù)是一個非線性函數(shù),為了簡化方程組,需要對普朗克(Planck)函數(shù)進行線性簡化。分別對MODIS的第31波段C10.78011.280戶)和32波段(11.77-12.27,))的熱輻射與溫度在273K-322K區(qū)間內的變化關系進行計算,得到如圖所示的計算結果。從圖中可以看出,熱輻射強度隨溫度的變化接近于線性關系。因此,對散點圖建立線性回歸方程,得到對第31波段331(7)=0.13787731-31.65677,i2=0.9971對第32波段532(7>0.11849r32-26.50036,及2=0.99783-2)簡化輻射傳輸方程組,9對于MODIS的第31和32波段,簡化的輻射傳輸方程組可以寫成如下[毛克彪,針對M0DIS數(shù)據(jù)的地表溫度反演方法研究,碩士學位論文,凌二貧乂學,2004.5.]:a31(t31)=r31的&(6)萬31(t;)+[i-r31的][i+(i-f31⑨^⑨]^(t;)(式7A)A2(4)=&(e)s32(0)532(rs)+[i_r32的][i+(i-f32(爭32(-,32(r。)(式7B)上式中r,.(P)是第31和32波段在角度0方向的透過率,s.(0是第31和32波段在角度0方向的發(fā)射率,r。是大氣平均作用溫度,t;是地表溫度,左邊第一項s,(;.)是第31和32波段的星上亮度溫度,等式右邊第一項是地表輻射,等式右邊第二項是大氣影響。把3-1)中第31和32波段的簡化方程S31(7)=0.13787r31-31.65677和532(rH).11849r32-26.50036分別代入輻射傳輸方程組,得到0.13787£"31r317;=0.13787:r31+31.65677£"31r31-(l-r31)[l+(a-f31A31](0.13787ra-31.65677)-31.65677(式8A)0.11849s32T327;=0.11849:r32+26.50036£"32z"32-(l-ir32)[l+(l-s32)r32](0.118497>26.50036)-26.50036(式8B)為了便于計算,將上面方程組中的系數(shù)分別記為A尸O.13787*r3/+31.65677*r31*£"31墨31.65677&尸(1〈31)*(1+(1<31)*r31)*0.13787332=0.11849*r32+26.50036*T"32*f32-26.50036C32=(l-r32)*(l+(l-s32)*r32)*0.11849""=(l-r32)*(l+(l-Q2)*r32)*26.50036輻射傳輸方程組可以寫成r"(式9A)&(式9B)解方程組,陸地表面溫度的計算公式為K=(C"^+i^)-C3//(C^;;-C"》(式10)第四步反演地表溫度,4-1)利用普朗克函數(shù)(Planck)反函數(shù)計算第31和32波段的星上亮度溫度T31和T32,下式中B31和B32分別是MODIS第31和32波段的星上亮度。T31=14380/(11.03*ln(2*59500000/(B31*11.03A5)+1))T32=14380/(12.02*ln(2*59500000/(B32*12.02A5)+l))4-2)將第一步和第二步計算得到的透過率和發(fā)射率,以及4-l)計算得到的星上亮度溫度分別代入系數(shù)爿w,,C^,C,利用地表溫度計算公式計算得到地表溫度。本發(fā)明的有益效果是,考慮MODIS像元l公里分辨率的特征,利用MODI近紅外和紅波段計算的植被指數(shù)來計算熱紅外波段31和32的發(fā)射率,保證了發(fā)射率的同步獲取,而且簡單并保證了一定的精度。利用近紅外波段反演大氣水汽含量,并利用MODTRAN4建立大氣透過率與水汽的關系,從而實時的獲得大氣透過率參數(shù)。這樣解決了地表溫度反演中的兩個關鍵參數(shù)。提高了反演精度和計算時間,克服以往美國宇航局算法復雜,一般科研人員難以實現(xiàn)的缺點。為氣候變化研究,氣象預報,蒸散發(fā),農(nóng)情監(jiān)測以及災害監(jiān)測等提供了有效手段和技術支撐,這個方法已經(jīng)確定被農(nóng)業(yè)部農(nóng)情監(jiān)測系統(tǒng)采用。下面結合附圖和實施例對本發(fā)明進一步說明。圖1MODIS遙感器。圖2是本發(fā)明的主流程示意圖。圖3是環(huán)渤海地區(qū)的大氣水汽含量遙感反演結果。圖4是MODIS的第31波段大氣透過率圖。ii圖5是MODIS的第32波段大氣透過率圖。圖6是NDVI分布圖。圖7是MODIS第31波段的地表比輻射率分布圖。圖8是MODIS第32波段的地表比輻射率分布圖。圖9是反演而得的環(huán)渤海地區(qū)地表溫度分布圖。圖10對比驗證結果。具體實施方式這里用2003年8月上旬中國環(huán)渤海地區(qū)的一景MODIS影像數(shù)據(jù)做具體反演,本方法主要包括四個步驟,如圖2。第一步,通過大氣水汽含量,計算MODIS第31和32波段的透過率1-1)利用MODIS第19波段和第2波段計算比值T,然后利用這個比值計算大氣水汽含量W,如圖3:A0.02-lnr20.6511-2)利用大氣水汽含量與MODIS第31和32波段的關系f3l=2.89798—1.88366e2'.22704t32=—3.59289+4.60414e32639將l-l)中計算得到的大氣水汽含量代入上面公式,計算得到第31波段和32波段的透過率(r31,r32),得到第31和32波段的透過率,如圖4和5。第二步、通過植被指數(shù)NDVI(NormalizedDifferenceVegetationIndex)計算MODIS傳感器第31和32波段的發(fā)射率。2-1)讀入MODIS第l波段(近紅外波段)和第2波段(紅光波段),計算NDVI;5,+及計算得到的結果如圖6。2-2)利用NDVI來判斷1平方公里分辨率像元的地物類型,然后根據(jù)美國噴汽推進實驗室(JPL)測試的ASTER波譜庫計算M0DIS傳感器第31波段和32波段的發(fā)射率(&,^分別表示第31和32波段發(fā)射率)。計算如下當NDVK0時,是水體和雪,&=0.992,&=0.988當0<NDVI<0.05,是裸土,=0.986,&2=0.991當0.05<NDVI<0.65,為植被和裸土的混合像元。PV=(NDVI-0.05)/0.6f31=0.986*(l-PV)+0.972*PV,f32=0.991*(1-PV)+0.976*PV;當NDVI〉0.65,為植被s31=0.972,s32=0.976計算后得到MODIS第31和32波段的發(fā)射率分別如圖7和圖8。第三步建立MODIS第31波段和32波段簡化的輻射傳輸方程組第四步利用下面公式計算MODIS第31波段和32波段的星上亮度溫度T31=14380/(11.03*ln(2*59500000/(B31*11.03A5)+1》T32=14380/(12.02*ln(2*59500000/(B32*12.02A5)+l》將第一步和第二步計算得到的發(fā)射率和透過率,計算得到的星上亮度溫度分別代入系數(shù)A"A;,A〃A尸O.13787*7^+31.65677*r31-31.65677G;=(l-z"31)*(l+(l-f31)*r31)*0.13787仏尸("31)*(1+(1,)、)*31.65677^32=0.11849*£32*r32532=0.11849*7^+26.50036*r32*£"32-26.50036C^=(l-r32)*(l+(l-s32)*r32)*0.11849A2=(l-r32)*(l+(l^2)*r32)*26.50036利用地表溫度計算公式計算得到地表溫度,如圖9。7>(C32(^+D3;)-/(c"廣c"》在應用中,對方法反演精度評價是非常重要的。因為MODIS的像元分辨率是1公里*1公里),地表測量通常是點測量。我們很難用幾個點測量的數(shù)據(jù)的平均值來代表一個大尺度的像元,而且是要在衛(wèi)星過境時同時獲取。在這里我們用美國通量網(wǎng)數(shù)據(jù)來對我們的方法進行評價。通量數(shù)據(jù)庫0lttp:〃public.ornl.gov/ameriflux/datahandler.cfrn)是一個全球微氣象觀測網(wǎng)絡,主要是用來監(jiān)測二氧化碳交換、水蒸汽、陸地生態(tài)系統(tǒng)和大氣能量交換。也被用來對MODIS地表溫度產(chǎn)品進行驗證[Wang,W.,S.Liang,ValidatingMODISlandsurfacetemperatureproduct,ispmsrs05,17-19,October,Beijing,China,2005.〗。這里我們也選擇了3個地表比較比較平坦且均一站點(Brookings,Audubon,F(xiàn)ortPeck)的地表溫度作為實測量數(shù)據(jù)。具體介紹請參考(http:〃public.oml.gov/ameriflux/site-selectcfin)。從MODIS反演地表溫度和通量站點數(shù)據(jù)比較如圖10所示。平均精度大約是1.2K,能滿足我們目前的應用需求。權利要求1、從MODIS數(shù)據(jù)反演地表溫度方法,其步驟為第一步,通過大氣水汽含量,計算MODIS第31和32波段的透過率1-1)利用MODIS第19波段和第2波段計算比值T,然后利用這個比值計算大氣水汽含量W。利用大氣輻射傳輸模型MODTRAN4模擬計算得到大氣水汽含量與MODIS第31和32波段透過率的關系,計算得到的大氣水汽含量代入下面公式,計算得到第31波段和32波段的透過率。<mathsid="math0001"num="0001"><math><![CDATA[<mrow><msub><mi>τ</mi><mn>31</mn></msub><mo>=</mo><mn>2.89798</mn><mo>-</mo><mn>1.88366</mn><msup><mi>e</mi><mfrac><mi>w</mi><mn>21.22704</mn></mfrac></msup></mrow>]]></math></maths><mathsid="math0002"num="0002"><math><![CDATA[<mrow><msub><mi>τ</mi><mn>32</mn></msub><mo>=</mo><mo>-</mo><mn>3.59289</mn><mo>+</mo><mn>4.60414</mn><msup><mi>e</mi><mfrac><mrow><mo>-</mo><mi>w</mi></mrow><mn>32.70639</mn></mfrac></msup></mrow>]]></math></maths>第二步、通過植被指數(shù)NDVI計算MODIS傳感器第31和32波段的發(fā)射率。2-1)讀入MODIS第1波段(近紅外波段)和第2波段(紅光波段),計算NDVI。利用NDVI來判斷1平方公里分辨率像元的地物類型,然后根據(jù)美國噴汽推進實驗室(JPL)測試的ASTER波譜庫計算MODIS傳感器第31波段和32波段的發(fā)射率(ε31,ε32分別表示第31和32波段發(fā)射率)。計算如下當NDVI<0時,是水體和雪,ε31=0.992,ε32=0.988當0<NDVI<0.05,是裸土,ε31=0.986,ε32=0.991當0.05<NDVI<0.65,為植被和裸土的混合像元。利用PV指數(shù)混合像元中植被的覆蓋的比率。計算公式PV=(NDVI-0.05)/0.6(式5)ε31=0.986*(1-PV)+0.972*PV,ε32=0.991*(1-PV)+0.976*PV;當NDVI>0.65,為植被ε31=0.972,ε32=0.976第三步簡化輻射傳輸方程3-1)簡化普朗克(Planck)普朗克(Planck)函數(shù)是一個非線性函數(shù),為了簡化方程組,需要對普朗克(Planck)函數(shù)進行線性簡化。分別對MODIS的第31波段(10.780~11.280μm)和32波段(11.77-12.27μm))的熱輻射與溫度在273K-322K區(qū)間內的變化關系進行計算,得到如圖所示的計算結果。從圖中可以看出,熱輻射強度隨溫度的變化接近于線性關系。因此,對散點圖建立線性回歸方程,得到對第31波段B31(T)=0.13787T31-31.65677,R2=0.9971對第32波段B32(T)=0.11849T32-26.50036,R2=0.99783-2)簡化輻射傳輸方程組,對于MODIS的第31和32波段,簡化的輻射傳輸方程組可以寫成如下B31(T31)=τ31(θ)ε31(θ)B31(Ts)+[1-τ31(θ)][1+(1-ε31(θ)τ31(θ)]B31(Ta)B32(T32)=τ32(θ)ε32(θ)B32(Ts)+[1-τ32(θ)][1+(1-ε32(θ)τ32(θ)]B32(Ta)把3-1)中第31和32波段的簡化方程B31(T)=0.13787T31-31.65677和B32(T)=0.11849T32-26.50036分別代入輻射傳輸方程組,得到0.13787ε31τ31Ts=0.13787T31+31.65677ε31τ31-(1-τ31)[1+(1-ε31)τ31](0.13787Ta-31.65677)-31.656770.11849ε32τ32Ts=0.11849T32+26.50036ε32τ32-(1-τ32)[1+(1-ε32)τ32](0.11849Ta-26.50036)-26.50036為了便于計算,將上面方程組中的系數(shù)分別記為A31=0.13787*ε31*τ31B31=0.13787*T31+31.65677*τ31*ε31-31.65677C31=(1-τ31)*(1+(1-ε31)*τ31)*0.13787D31=(1-τ31)*(1+(1-ε31)*τ31)*31.65677A32=0.11849*ε32*τ32B32=0.11849*T32+26.50036*τ32*ε32-26.50036C32=(1-τ32)*(1+(1-ε32)*τ32)*0.11849D32=(1-τ32)*(1+(1-ε32)*τ32)*26.50036輻射傳輸方程組可以寫成A31TS=B31-C31Ta+D31A32TS=B32-C32TA+D32解方程組,陸地表面溫度的計算公式為Ts=(C32(B31+D31)-C31(D32+B32))/(C32A31-C31A32)。全文摘要本發(fā)明涉及一種從MODIS數(shù)據(jù)反演地表溫度的方法,能夠應用在氣象、環(huán)境監(jiān)測、土地管理、農(nóng)情監(jiān)測、以及災害監(jiān)測等遙感應用部門。該方法,包含四個步驟第一步驟是通過大氣水汽含量,計算MODIS第31和32波段的透過率。第二個步驟是通過植被指數(shù)NDVI計算MODIS傳感器第31和32波段的發(fā)射率。第三步驟是簡化輻射傳輸方程。第四步是對MODIS數(shù)據(jù)進行反演計算,得到地表目標對象溫度和發(fā)射率分布情況,可以用于氣象預報、環(huán)境監(jiān)測、農(nóng)情監(jiān)測和災情監(jiān)測等部門。文檔編號G01J5/00GK101629850SQ200910091030公開日2010年1月20日申請日期2009年8月24日優(yōu)先權日2009年8月24日發(fā)明者周清波,張立新,李三妹,毛克彪,王道龍,覃志豪申請人:中國農(nóng)業(yè)科學院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所;國家衛(wèi)星氣象中心