專利名稱：：從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>&tau;</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>&tau;</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)星氣象中心