楊何群，尹 球，周紅妹，葛偉強(qiáng)

(上海市衛(wèi)星遙感與測(cè)量應(yīng)用中心，上海 201199)

利用MATLAB實(shí)現(xiàn)FY－3/MERSI地表溫度反演及專題制圖

楊何群，尹 球，周紅妹，葛偉強(qiáng)

(上海市衛(wèi)星遙感與測(cè)量應(yīng)用中心，上海 201199)

我國(guó)新型自主研發(fā)的風(fēng)云三號(hào)衛(wèi)星MERSI(FY－3/MERSI)數(shù)據(jù)目前多見于試驗(yàn)研究，國(guó)家衛(wèi)星氣象中心(NSMC)發(fā)布的業(yè)務(wù)產(chǎn)品中尚缺乏地表溫度(land surface temperature，LST)產(chǎn)品。為此，在分析FY－3/MERSI L1級(jí)數(shù)據(jù)HDF5格式及通道特點(diǎn)的基礎(chǔ)上，采用Jiménez－Mun～oz和Sobrino普適性單通道地表溫度反演算法，借助MATLAB編程，實(shí)現(xiàn)了直接利用MERSI數(shù)據(jù)進(jìn)行250 m空間分辨率的LST反演與專題制圖。詳細(xì)介紹了LST反演算法的參數(shù)化流程，包括輻射定標(biāo)、云檢測(cè)處理以及地表比輻射率和大氣水汽含量的估計(jì)等，并在反演LST的基礎(chǔ)上，計(jì)算其他衍生指標(biāo)，建立了FY－3/MERSI數(shù)據(jù)反演LST及制圖流程。對(duì)上海市空間熱環(huán)境監(jiān)測(cè)分析的實(shí)驗(yàn)結(jié)果表明:該方法能根據(jù)設(shè)定范圍對(duì)FY－3/MERSI L1數(shù)據(jù)進(jìn)行快速、實(shí)時(shí)和自動(dòng)處理，適用于業(yè)務(wù)化產(chǎn)品制作與共享，節(jié)省人力資源;形式豐富的MERSI專題制圖產(chǎn)品可更直觀、精細(xì)地揭示出上海城市熱場(chǎng)空間格局與熱島效應(yīng)特征。

FY－3/MERSI;地表溫度;反演;專題制圖;MATLAB

0 引言

地表溫度(land surface temperature，LST)是描述地－氣系統(tǒng)相互作用中物質(zhì)和能量交換結(jié)果的重要物理參量，是地表熱量平衡的產(chǎn)物［1－2］。地球上的植物生長(zhǎng)與作物產(chǎn)量、地表水分蒸發(fā)與循環(huán)、城市生態(tài)與熱島效應(yīng)、氣候變遷與全球環(huán)境變化以及資源積累與人類生存等許多重要的自然、人文現(xiàn)象或過(guò)程，無(wú)不與其密切相關(guān)。地表溫度蘊(yùn)含了豐富的地學(xué)信息，是地面與空氣界面的能量計(jì)數(shù)器。由于其在時(shí)空上具有高度變異性［3］，區(qū)域乃至全球性的動(dòng)態(tài)監(jiān)測(cè)與應(yīng)用需完全依賴于衛(wèi)星遙感技術(shù)。

作為中國(guó)首顆高性能綜合探測(cè)衛(wèi)星，中國(guó)自主研發(fā)的第二代極軌氣象衛(wèi)星系列風(fēng)云三號(hào)(FY－3)是與美國(guó)氣象與環(huán)境監(jiān)測(cè)衛(wèi)星NPOESS、歐洲極軌氣象衛(wèi)星METOP等相當(dāng)?shù)男乱淮l(wèi)星，可實(shí)現(xiàn)對(duì)地球－大氣的全天候、全天時(shí)、三維、定量、多光譜的不間斷觀測(cè)［4］。其上搭載的中分辨率光譜成像儀(medium resolution spectral imager，MERSI)是最主要的光學(xué)成像載荷，具備從可見光—熱紅外20個(gè)通道的探測(cè)能力，可以探測(cè)大氣、陸地、海洋的可見光反射輻射和熱紅外發(fā)射輻射。區(qū)別于其他同類星載傳感器(如EOS/MODIS)的是，F(xiàn)Y－3/MERSI加強(qiáng)了對(duì)地物精細(xì)觀測(cè)的能力，將250 m分辨率通道從MODIS的2個(gè)增加到5個(gè)，特別是包含了一個(gè)中心波長(zhǎng)位于 11.25 μm，光譜帶寬達(dá) 2.5 μm 的熱紅外通道。由于空間分辨率和光譜響應(yīng)均得以提高［5］，該通道可以更為精細(xì)地捕捉地表、洋面及云頂?shù)臒釋傩孕畔?，同時(shí)也可為云檢測(cè)和其他遙感產(chǎn)品的生成提供高分辨率熱紅外圖像。

限于FY－3A(上午星)和FY－3B(下午星)是近年發(fā)射(前者2008年5月27日發(fā)射，后者2010年11月5日發(fā)射)，并各自經(jīng)歷長(zhǎng)達(dá)半年多的在軌測(cè)試，目前，F(xiàn)Y － 3/MERSI只見于試驗(yàn)研究［6－7］。國(guó)內(nèi)遙感研究機(jī)構(gòu)大都仍以EOS/MODIS數(shù)據(jù)為主要數(shù)據(jù)源進(jìn)行相關(guān)領(lǐng)域的監(jiān)測(cè)評(píng)估應(yīng)用。隨著衛(wèi)星系統(tǒng)穩(wěn)定、正常運(yùn)行及數(shù)據(jù)的積累，有必要發(fā)展區(qū)域性FY－3/MERSI地表溫度及衍生產(chǎn)品的業(yè)務(wù)化流程算法，這一方面有助于提升中國(guó)衛(wèi)星遙感自主應(yīng)用能力，另一方面可在業(yè)務(wù)系統(tǒng)中減少對(duì)國(guó)外遙感數(shù)據(jù)的依賴，確保熱力要素相關(guān)監(jiān)測(cè)的及時(shí)性和準(zhǔn)確性。針對(duì)FY－3/MERSI的數(shù)據(jù)特點(diǎn)，本文探索通過(guò)MATLAB語(yǔ)言編程實(shí)現(xiàn)FY－3/MERSI區(qū)域地表溫度反演及專題制圖。

1 FY－3/MERSI數(shù)據(jù)及MATLAB優(yōu)勢(shì)

FY －3/MERSI傳感器掃描角度為 ±55.4°，成像幅寬為2000 km，光譜范圍為0.41 ～12.5 μm，共有20個(gè)離散通道，其中250 m分辨率通道有5個(gè)，1000 m分辨率通道有15個(gè)，量化等級(jí)均為12 bit(表1)。

表1 FY－3/MERSI光譜通道及性能指標(biāo)Tab.1 FY －3/MERSI spectral channels and characteristics

在FY－3衛(wèi)星繞地球運(yùn)行過(guò)程中，MERSI探測(cè)的輻射信息在星上編排為MPT(medium resolution spectral picture transmission，中分辨率圖像傳輸)數(shù)據(jù)，通過(guò)X波段傳輸通道發(fā)送至地面接收站，解碼成為L(zhǎng)0級(jí)數(shù)據(jù)。中國(guó)氣象局DVBS(digital video broadcasting system，數(shù)字視頻廣播系統(tǒng))及風(fēng)云衛(wèi)星遙感數(shù)據(jù)服務(wù)網(wǎng)(fy3.satellite.cma.gov.cn)實(shí)時(shí)發(fā)布MERSI L1級(jí)產(chǎn)品，它是原始L0級(jí)數(shù)據(jù)經(jīng)過(guò)多站接收去重復(fù)、質(zhì)量檢驗(yàn)后進(jìn)入定位定標(biāo)預(yù)處理后生成的，存貯格式為HDF，使用的是HDF5版本。

MERSI數(shù)據(jù)開發(fā)應(yīng)用的前提是實(shí)現(xiàn)對(duì)MERSI L1級(jí)數(shù)據(jù)HDF5文件的提取。目前常用的遙感處理軟件(如ENVI)只支持MODIS 1B HDF4和HDFEOS格式的讀取，還沒有提供HDF5數(shù)據(jù)支持接口或插件，更不用說(shuō)提供直接的MERSI LST反演操作，應(yīng)用受到軟件自身功能限制。C，F(xiàn)ortran等傳統(tǒng)語(yǔ)言工具雖可讀取HDF5文件，但程序代碼較為復(fù)雜，且配置HDF庫(kù)函數(shù)繁瑣。

MATLAB作為美國(guó)MathWorks公司出品的商業(yè)數(shù)學(xué)軟件，用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語(yǔ)言和交互式環(huán)境，可以很容易地實(shí)現(xiàn)C或Fortran語(yǔ)言的幾乎全部功能，而且它提供了一組HDF應(yīng)用程序接口，可以方便訪問HDF5文件，并且以矩陣計(jì)算為基礎(chǔ)，可使圖像處理編程大大簡(jiǎn)化。此外，MATLAB具有出色的圖形處理功能及功能強(qiáng)大的模塊集/工具箱。

2 LST反演算法

隨著對(duì)熱紅外遙感基礎(chǔ)理論問題研究的深入，LST物理反演算法有了很大發(fā)展。結(jié)合傳感器成像特性，LST反演方法有單通道法、劈窗法、多通道法、多角度法以及溫度與比輻射率分離算法，等等。由于FY－3/MERSI繼承了FY－1的特點(diǎn)——單窗，即只有一個(gè)熱紅外通道，因而只適宜采用單通道反演法。經(jīng)調(diào)研對(duì)比，發(fā)現(xiàn) Jiménez－Mun～oz和 Sobrino［8］提出的普適性單通道算法相對(duì)簡(jiǎn)單實(shí)用。

2.1 Jiménez－Mun～oz和 Sobrino 單通道法

LST反演以地表熱輻射傳導(dǎo)方程為基礎(chǔ)，即在晴空大氣、局地?zé)崃ζ胶鈼l件下，假定地球表面為朗伯發(fā)射體，并且忽略大氣分子和氣溶膠的散射，認(rèn)為衛(wèi)星熱紅外輻射傳感器所接收的能量主要包括3個(gè)部分:經(jīng)大氣削弱后被傳感器接收的地表熱輻射能量、大氣下行輻射經(jīng)地表反射后再被大氣削弱而最終被傳感器接收的那部分能量和大氣上行輻射能量。該方程的核心是Planck函數(shù)，數(shù)學(xué)描述為

式中:Bi(Ti)和Bi(Ts)分別為第i通道的星上輻射強(qiáng)度和地面輻射強(qiáng)度;Ti為第i通道的星上亮度溫度;Ts為地表溫度;τi為第i通道的大氣透過(guò)率;εi為第i通道的比輻射率;θ為觀測(cè)視角;分別為第i通道的大氣下、上行輻射能量。

為了獲得真實(shí)的地表溫度，必須剔除掉衛(wèi)星信號(hào)中由于大氣影響而產(chǎn)生的噪聲，即需要進(jìn)行大氣校正。Jiménez－ Mun～oz和 Sobrino以單個(gè)熱紅外通道為對(duì)象，在對(duì)Planck方程作一階泰勒級(jí)數(shù)展開的基礎(chǔ)上，推導(dǎo)并簡(jiǎn)化了輻射傳輸方程，從而提出了把大氣校正包含在地表溫度反演過(guò)程中的普適性單通道算法，即

其中，

式中:Ts為地表溫度，K;Tb為亮度溫度，K;λ為有效波長(zhǎng)，μm;為傳感器接收到的熱輻射亮度，Wm－2sr－1μm－1;c1，c2為系數(shù); ε 為地表比輻射率(無(wú)量綱); ψ1，ψ2，ψ3分別為大氣水汽含量 w(g·cm－2)及有效波長(zhǎng)λ的函數(shù)。

可以看出，Jiménez－ Mun～oz和 Sobrino 單通道法只需通道有效波長(zhǎng)、地表比輻射率和整層大氣水汽含量3個(gè)中間參數(shù)，而不象其他單通道算法［9］需要探空資料或地面觀測(cè)的大氣有效平均作用溫度等，而且它可將相同的方程和系數(shù)應(yīng)用于不同的傳感器。基于該算法，應(yīng)用Landsat TM數(shù)據(jù)的LST反演已得到廣泛認(rèn)可，均方根誤差只在0.6～1.5 K左右［10］，說(shuō)明將其應(yīng)用于 FY－3A/MERSI數(shù)據(jù)也應(yīng)有較高精度。

2.2 地表比輻射率估計(jì)

地表比輻射率ε作為先驗(yàn)知識(shí)參與地表溫度反演，但由于其影響因素很多(包括地表物質(zhì)結(jié)構(gòu)、表面粗糙度、介電常數(shù)、含水量、波長(zhǎng)及觀測(cè)角度等等)，對(duì)其測(cè)量與計(jì)算至今尚未找到完善的數(shù)值解法和模型。Valor和 Casselles［11］提出了利用植被覆蓋度Pv建立估算比輻射率的模型，同時(shí)考慮了混合像元的影響，適用于類似FY－3/MERSI中分辨率混合像元的傳感器，計(jì)算方法為

利用遙感數(shù)據(jù)估算植被覆蓋度Pv的常用方法有三波段梯度差法［12］和歸一化植被指數(shù)NDVI像元二分法［13］，后者被證明具有較高精度，計(jì)算方法為

式中:NDVI=(ρ4－ ρ3)/(ρ4+ ρ3)，ρ4和 ρ3分別為經(jīng)過(guò)大氣校正的FY－3/MERSI近紅外波段(0.815～0.915 μm)及紅波段(0.600 ～0.700 μm)的反射率;NDVIs，NDVIv分別為裸露土壤或無(wú)植被建筑表面、完全植被覆蓋像元的NDVI值，取值分別為0.05和 0.75。

2.3 大氣水汽含量測(cè)算

大氣水汽含量是指單位截面大氣柱內(nèi)所含水汽總量。目前水汽測(cè)量具有多種多樣的形式，包括探空、紅外、微波探測(cè)，以及近年興建的GPS/MET站網(wǎng)等，而最成熟且最簡(jiǎn)便的仍是應(yīng)用衛(wèi)星遙感近紅外波段的方法，能精確到像元尺度［14］。其中，Kaufman和Gao［15］通過(guò)大量實(shí)驗(yàn)研究，從差分吸收概念得出了利用近紅外大氣水汽吸收帶及其邊沿大氣窗口區(qū)2通道比值法估算MODIS大氣柱水汽含量的公式;胡秀清等［16］利用查找表方法，開發(fā)了基于FY －3/MERSI 0.905 μm，0.940 μm 及 0.980 μm 3個(gè)近紅外水汽吸收通道加權(quán)反演水汽總量的算法，通過(guò)反演試驗(yàn)和結(jié)果檢驗(yàn)也證明了FY－3/MERSI具有較強(qiáng)的近紅外水汽探測(cè)能力，并能探測(cè)空中水汽空間分布的精細(xì)特征。但是，由于此算法還在優(yōu)化過(guò)程中，本文將Kaufman和Gao算法，即

其中，

應(yīng)用于FY－3/MERSI對(duì)應(yīng)通道，測(cè)算水汽總量。式中: τa是大氣透射率; ρ16，ρ18分別是 FY －3/MERSI第16通道(中心波長(zhǎng)0.865 μm)和第18通道(中心波長(zhǎng) 0.980 μm)的反射率。

對(duì)于系數(shù)α和β，可根據(jù)下墊面狀況進(jìn)行調(diào)整:完全植被覆蓋的地表(NDVI≥0.75)，α =0.012，β =0.651; 混合型地表(0.05≤NDVI＜ 0.75)，α =0.020，β=0.651;地表為裸土或建筑用地等不透水面(NDVI＜0.05)，α = －0.040，β =0.651。

3LST反演及專題制圖的MATLAB實(shí)現(xiàn)

3.1 MERSI HDF5數(shù)據(jù)的讀取

目前LST的業(yè)務(wù)產(chǎn)品主要為1000 m分辨率的MODIS LST，而FY－3/MERSI因具有1個(gè)獨(dú)具特色的250 m分辨率的熱紅外通道，為每日生成更精細(xì)的百米級(jí)地表溫度產(chǎn)品提供了新數(shù)據(jù)源。上文所述的 Jiménez－ Mun～oz和 Sobrino 普適性單通道算法需涉及 FY －3/MERSI第 3，4，5，16，18 共 5 個(gè)通道，但只有用于計(jì)算NDVI與亮度溫度Tb的第3，4，5通道為250 m分辨率，而用于反演大氣水汽含量的16和18通道分辨率僅為1000 m。因此，F(xiàn)Y－3/MERSI 250 m和1000 m二種分辨率數(shù)據(jù)均需作為輸入數(shù)據(jù)源。為了空間匹配，對(duì)后者還需進(jìn)行插值。另外，考慮到FY－3/MERSI采用的是多探元跨軌并掃成像模式，即250 m分辨率通道(通道1～5)是40個(gè)探元并掃，1000 m分辨率通道(通道6～20)是10個(gè)探元并掃，探元的非均勻性會(huì)使原始圖像存在條帶，即使采用星上校正裝置根據(jù)發(fā)射前確定的校正系數(shù)也不能很好地消除，故而需進(jìn)行均勻性校正和圖像消旋［17］處理。

輸入數(shù)據(jù)即為經(jīng)過(guò)均勻性校正和圖像消旋的FY－3/MERSI 250 m和1000 m L1產(chǎn)品。由于MERSI的光譜通道較多，再加上較寬的掃描范圍，因而每隔5 min獲取的1景數(shù)據(jù)量常常高達(dá)1.0 GB，如果要完全讀入，必將導(dǎo)致MATLAB內(nèi)存溢出，所以，這里只針對(duì)感興趣研究區(qū)域進(jìn)行相關(guān)數(shù)據(jù)讀取，其技術(shù)路線如圖1所示。

圖1 FY－3/MERSI HDF5數(shù)據(jù)讀取技術(shù)路線示意圖Fig.1 Technical route diagram of reading FY －3/MERSI HDF5 data

通過(guò)圖1所示步驟，形成了維數(shù)一致的5個(gè)通道DN計(jì)數(shù)值和2個(gè)地理位置的矩陣數(shù)據(jù)集，以及一些標(biāo)量信息，便于后續(xù)反演處理。

3.2 數(shù)據(jù)預(yù)處理

3.2.1 輻射定標(biāo)及反射率、亮溫計(jì)算

采用非線性多項(xiàng)式對(duì)FY－3/MERSI數(shù)據(jù)進(jìn)行輻射定標(biāo)(紅外通道是三次項(xiàng)定標(biāo)，可見光通道是二次項(xiàng)定標(biāo))。由于MERSI數(shù)據(jù)存貯的是DN計(jì)數(shù)值，需先通過(guò)DN*=Slope(DN－Intercept)轉(zhuǎn)換后再參與定標(biāo)計(jì)算。但是，紅外定標(biāo)過(guò)程比較復(fù)雜，包括逐探元定標(biāo)、多探元輻射響應(yīng)差異光譜補(bǔ)償和圖像消旋處理。所以，MERSI存貯的是已經(jīng)過(guò)定標(biāo)的輻亮度值，只需通過(guò)Planck定律逆求解得到等效黑體亮溫Teff，再進(jìn)一步計(jì)算，獲得亮溫溫度Tb。計(jì)算式為

式中:Lλ為經(jīng)過(guò)定標(biāo)的輻亮度，mW·m－2·sr－1·cm;v為通道中心波數(shù)，為 875.1379 cm－1;系數(shù) A=1.0103，B= － 1.8521，c1=1.191066 × 10－5mW·m－2·sr－1·cm，c2=1.438833 K·cm。

式中:h為普朗克常數(shù)，h=6.626 ×10－34J·s;c為光速，c=2.99793 ×108m·s－1;k 為玻爾茲曼常數(shù)，k=1.3806 × 10－23J·K－1。

反射率則通過(guò)二次項(xiàng)定標(biāo)計(jì)算，即

式中:ρTOAi為第i通道大氣層頂(top of the atmos-phere，TOA)反射率;θ為太陽(yáng)天頂角;為對(duì)應(yīng)i通道的3個(gè)定標(biāo)系數(shù)。參考權(quán)維?。?8］等利用MODTRAN輻射傳輸模型模擬和利用最小二乘回歸方法建立的FY－3/MERSI第3，4通道參數(shù)化方程

3.2.2 云檢測(cè)處理

需注意的是，在有云的情況下，利用熱紅外波段數(shù)據(jù)反演的不是真實(shí)的地表溫度，而是云頂溫度。為有效識(shí)別云像元，本文引進(jìn)歸一化云檢測(cè)指數(shù)(cloud detection index，CDI)進(jìn)行處理，計(jì)算式為

式中: ρ3，ρ18分別為 FY －3/MERSI第 3，18 通道的地物反射率。由于云在0.65 μm 和 0.94 μm 處(與FY－3/MERSI第3，18通道相對(duì)應(yīng))的光譜特征與其他地物的光譜特征形成明顯反差:云在0.65 μm處具有高反射率，能良好區(qū)分出陸地和云的邊界;水汽在0.94 μm處具有吸收谷。將它們歸一化處理，便可以突出云的信息。根據(jù)CDI＞0或CDI≤0的閾值劃分可以生成一個(gè)二值數(shù)據(jù)集，用0表示非云，1表示云，將其與地表溫度數(shù)據(jù)集進(jìn)行掩模計(jì)算，即可獲得晴空地域地表溫度信息。

但實(shí)際上，云頂溫度作為云微物理屬性的重要參量，對(duì)于天氣分析與預(yù)報(bào)非常有用，對(duì)于業(yè)務(wù)應(yīng)用實(shí)踐并不需進(jìn)行地表溫度的去云運(yùn)算，只需將云信息另行存儲(chǔ)或標(biāo)識(shí)。

3.3 LST反演及衍生指標(biāo)提取

文中采用“植被指數(shù)→植被覆蓋度→地表比輻射率”，“大氣透射率→大氣水汽含量”逐步進(jìn)行估算，借助于MERSI數(shù)據(jù)自身對(duì)地表反射/發(fā)射屬性的測(cè)量，以及MERSI近紅外波段對(duì)水汽參數(shù)敏感的特性，以確保MERSI各像元能較精確地獲得對(duì)應(yīng)的地表比輻射率與大氣透過(guò)率。這種參數(shù)化的方法克服了過(guò)去某種地表類型只用一個(gè)地表比輻射率以及同一景圖像只用一個(gè)大氣透過(guò)率的缺點(diǎn)，將對(duì)參數(shù)的估計(jì)從一個(gè)點(diǎn)擴(kuò)大到對(duì)整個(gè)圖像的各個(gè)像元的估計(jì)上，具有連續(xù)性變化，從而使LST反演更加符合實(shí)際情況。參照2.1 節(jié)中Jiménez－Mun～oz和Sobrino普適性單通道算法所涉及的公式，通過(guò)“地表比輻射率+大氣水汽含量+熱輻亮度→地表溫度”的矩陣運(yùn)算，即可快速實(shí)現(xiàn)FY－3/MERSI LST反演。

LST是熱環(huán)境最基礎(chǔ)的指標(biāo)，尤其是城市空間熱環(huán)境的主載體和量化基本。然而，僅憑LST的絕對(duì)值并不能完全反映熱環(huán)境的空間格局和動(dòng)態(tài)變化。因此，可在遙感反演地表溫度的基礎(chǔ)上，通過(guò)MATLAB設(shè)計(jì)相關(guān)程序搜索標(biāo)識(shí)熱島中心，提取重點(diǎn)地標(biāo)、地塊溫度，計(jì)算熱場(chǎng)溫度距平、熱場(chǎng)強(qiáng)度歸一化指數(shù)、熱島強(qiáng)度分級(jí)以及熱島比例指數(shù)等由LST衍生的指標(biāo)，勾畫熱場(chǎng)等值線及熱場(chǎng)橫縱向剖面，分析與下墊面介質(zhì)的對(duì)應(yīng)關(guān)系，實(shí)現(xiàn)高溫區(qū)分級(jí)、面積統(tǒng)計(jì)以及區(qū)域最高、最低和平均溫度的統(tǒng)計(jì)等，從而多角度展示熱環(huán)境空間形態(tài)及統(tǒng)計(jì)信息。這些指標(biāo)在MATLAB中均可方便計(jì)算統(tǒng)計(jì)。

3.4 專題制圖和結(jié)果保存

MATLAB具有強(qiáng)大的制圖功能，特別是MATLAB提供的地圖工具箱(mapping toolbox)，便于調(diào)用相關(guān)函數(shù)及圖形用戶界面(GUIs)，從而實(shí)現(xiàn)地理空間數(shù)據(jù)的分析及地圖顯示創(chuàng)建，而不需依賴于大家通常采用的通過(guò)ENVI等遙感處理軟件將結(jié)果格式轉(zhuǎn)換導(dǎo)入ArcGIS中，以人機(jī)交互方式統(tǒng)計(jì)面積、類型等相關(guān)信息，手動(dòng)生成專題圖。

經(jīng)實(shí)驗(yàn)，在MATLAB中，只要事先設(shè)計(jì)好圖形顯示模板，調(diào)用 axesm，surfm，shaperead，geoloc2grid，mapprofile及contourm等函數(shù)處理顯示柵格與矢量數(shù)據(jù)，設(shè)置投影方式、坐標(biāo)軸及圖形窗戶屬性等，添加標(biāo)題、指北針、顏色條、文字注釋、圖標(biāo)圖例等，就可自動(dòng)完成專題圖的繪制與保存。

此外，為了方便地表溫度及衍生指標(biāo)數(shù)據(jù)的使用和共享，可遵循標(biāo)準(zhǔn)HDF5格式，利用MATLAB的hdf5write函數(shù)將研究區(qū)域數(shù)據(jù)集進(jìn)行存貯輸出，但事先需設(shè)計(jì)好該類HDF5文件的組織架構(gòu)。通過(guò)創(chuàng)建一個(gè)新的HDF5文件，在指定位置(根組root group或子組subgroup)寫入地表溫度Ts數(shù)據(jù)集，同時(shí)可以添加屬性描述信息。另外，在此已存在的HDF5文件基礎(chǔ)上，再依次添加熱場(chǎng)溫度距平Ts_Ano、熱場(chǎng)強(qiáng)度歸一化指數(shù)HFII、經(jīng)度longitude、緯度latitude、歸一化云檢測(cè)指數(shù)CDI、二值化云掩模Cloud Mask等數(shù)據(jù)集。

綜上所述，采用FY－3/MERSI數(shù)據(jù)反演 LST算法及專題制圖流程如圖2所示。將各步驟的MATLAB代碼串聯(lián)編譯為可執(zhí)行程序。在此過(guò)程中，采取內(nèi)存動(dòng)態(tài)分配管理方法，隨著流程的進(jìn)行臨時(shí)釋放和加載數(shù)據(jù)。

圖2 FY－3/MERSI數(shù)據(jù)LST反演及制圖流程Fig.2Flow chart for FY －3/MERSI LST retrieval and mapping

4 應(yīng)用實(shí)例

以上海及其周邊地區(qū)LST反演及相關(guān)產(chǎn)品制圖為例，所用數(shù)據(jù)為2010年9月21日 FY－3A/MERSI數(shù)據(jù)。針對(duì)一個(gè)5 min的 FY－3/MERSI HDF Granule文件，進(jìn)行上海及周邊地區(qū)LST計(jì)算及相關(guān)產(chǎn)品制圖。從數(shù)據(jù)檢測(cè)到文件讀取、參數(shù)反演、產(chǎn)品制圖直至結(jié)果保存，整個(gè)流程耗時(shí)僅120 s左右，運(yùn)算速度快，無(wú)需人工干預(yù)即自動(dòng)完成。圖3為利用HDFView軟件展示的FY－3/MERSI LST產(chǎn)品HDF5格式結(jié)構(gòu)樣例圖。

圖3 FY－3/MERSI LST產(chǎn)品HDF5格式結(jié)構(gòu)樣例Fig.3 Sample of HDF5 format for FY－3/MERSI LST product

從LST反演結(jié)果(圖4)來(lái)看，上海市250 m LST分布產(chǎn)品的空間特征精細(xì)，區(qū)域內(nèi)LST分布差異由

圖4 上海市250m FY－3/MERSI LST分布圖Fig.4 LST Spatial distribution map in Shanghai，retrieved from 250 m FY－3/MERSI data

外圍到中心層次遞進(jìn)，城市熱島效應(yīng)明顯。具體地，城區(qū)地表溫度相對(duì)較高，高溫區(qū)域主要分布在人口密集的中心城區(qū)及北部寶山、嘉定等能耗大、熱源強(qiáng)度高的工業(yè)區(qū)內(nèi);中心城區(qū)周邊的近郊地區(qū)如閔行、浦東東北部、松江東部目前發(fā)展較快，已成為新的熱力中心基質(zhì)區(qū)域;其余周圍零散分布一些熱力斑塊。相反地，郊區(qū)地表溫度相對(duì)較低，東南部沿海和崇明縣溫度整體處于較低水平，比中心城區(qū)低約1.0～3.2℃;具有較大水域面積和較密植被的城中各大公園形成多個(gè)冷島。水域的溫度最低，可以看到在淀山湖區(qū)、長(zhǎng)江口與沿海溫度明顯偏低，另外黃浦江也形成了一條穿越城市的長(zhǎng)條形相對(duì)低溫帶。

LST統(tǒng)計(jì)直方圖(圖5)表明，2010年9月21日LST為23～34℃，且主要集中在27～30℃區(qū)間，像元點(diǎn)數(shù)多達(dá)1×104以上;30℃以上的高溫區(qū)面積比例約占23%，說(shuō)明城市整體溫度較高。

鑒于2010年9月21日上海天氣晴朗無(wú)風(fēng)，城市熱島效應(yīng)主要取決于下墊面介質(zhì)和結(jié)構(gòu)，以及城市建筑密度和高度［19］。以市中心 N 31.1°，E121.4°為原點(diǎn)，繪制城市熱場(chǎng)橫(EW向)、縱(SN向)剖面圖(圖6a和6b);對(duì)照上海市區(qū)土地覆蓋/利用分布圖(圖略)，分析上海城市熱場(chǎng)細(xì)部結(jié)構(gòu)。

圖5 上海市LST統(tǒng)計(jì)直方圖Fig.5 Histogram of LST bins in Shanghai

圖6 上海城市熱場(chǎng)剖面Fig.6 Profiles of thermal variation in Shanghai

從圖6可以看出，從遠(yuǎn)郊至近郊再到市區(qū)，地表溫度基本呈逐漸升高趨勢(shì)，不同下墊面的地表溫度差異明顯。高溫區(qū)主要由工業(yè)用地、商業(yè)區(qū)以及交通等混合的密集建筑用地、居民區(qū)用地組成，這主要是由于這些區(qū)域工業(yè)生產(chǎn)和人為生活排放熱較多，建筑物密集，不宜通風(fēng)，氣流擴(kuò)散慢，導(dǎo)致溫度很高;低溫區(qū)主要由農(nóng)田耕地、林地綠地、水體等介質(zhì)組成，由于水體的熱容量大，熱傳導(dǎo)率小，水體占低溫區(qū)的比例最多，降溫效果明顯，而農(nóng)田、耕地和林地綠地因覆蓋林草植被，可通過(guò)蒸騰作用減少土壤或地表中的熱量，地表溫度也較低，因此也表明植被和水體在城市降溫方面有著重要的作用，同時(shí)也說(shuō)明可通過(guò)增加植被、水體等景觀廊道的覆蓋面積來(lái)降低城市熱島強(qiáng)度，優(yōu)化城市空間布局。

以上分析說(shuō)明，F(xiàn)Y－3/MERSI數(shù)據(jù)獲得的250 m地表溫度能夠精確反映出城市熱場(chǎng)的分布格局以及城市熱島的特征和效應(yīng)。通過(guò)累積每日數(shù)據(jù)和產(chǎn)品，可用于分析不同天氣形式下城市熱場(chǎng)形態(tài)特征，摸清城市熱島效應(yīng)的日夜、季節(jié)、年代演變機(jī)理，為城市人居熱環(huán)境的模擬調(diào)控、優(yōu)化布局提供科學(xué)依據(jù)。

5 結(jié)語(yǔ)

本文在調(diào)研FY－3/MERSI HDF5數(shù)據(jù)結(jié)構(gòu)及普適性單通道LST反演算法參數(shù)化過(guò)程的基礎(chǔ)上，詳細(xì)介紹了利用MATLAB語(yǔ)言編程實(shí)現(xiàn)FY－3/MERSI L1級(jí)數(shù)據(jù)LST反演及專題制圖方法流程，并給出了應(yīng)用分析實(shí)例。結(jié)果表明:

1)利用MATLAB建立的FY－3/MERSI 250 m LST反演及專題制圖方法，不需依賴于GIS或遙感處理平臺(tái)就能從海量的FY－3/MERSI數(shù)據(jù)中直接提取有效數(shù)據(jù);處理過(guò)程中采取的內(nèi)存動(dòng)態(tài)分配管理方法，能快速實(shí)現(xiàn)LST反演及衍生產(chǎn)品制圖，便于我們建立自動(dòng)化、業(yè)務(wù)化工作流程，用于空間熱環(huán)境監(jiān)測(cè)服務(wù)。

2)基于FY－3/MERSI單窗特點(diǎn)，采用Jiménez－Mun～oz和Sobrino的普適性單通道算法，可獲得較合理的LST反演結(jié)果。應(yīng)用實(shí)例分析表明，F(xiàn)Y－3/MERSI數(shù)據(jù)質(zhì)量較好，能精細(xì)反映城市熱場(chǎng)的分布格局及城市熱島效應(yīng)，成為與MODIS匹敵的LST業(yè)務(wù)化產(chǎn)品，以滿足日常應(yīng)用需求。

3)需要注意的是，本文建立的流程目前只適用于白天地表溫度反演，對(duì)于夜間地表溫度的反演及產(chǎn)品制作還需另行設(shè)計(jì)修改。

針對(duì)FY－3/MERSI地表溫度反演結(jié)果的定量精度驗(yàn)證、內(nèi)部參數(shù)調(diào)整以及如何應(yīng)用于城市熱場(chǎng)特征機(jī)理研究是后續(xù)工作。

［1］Wan Z M，Dozier J.A Generalized Split－ window Algorithm for Retrieving Land Surface Temperature from Space［J］.IEEE Transactions on Geoscience and Remote Sensing，1996，34(4):892 －905.

［2］Sobrino J A，El Kharraz J，Li Z L.Suface Temperature and Water Vapour Retrieval from MODIS Data［J］.International Journal of Remote Sensing，2003，24:5161 －5182.

［3］Olesen F S，Kind O，Reutter H.High Resolution Time Series of IR Data from a Combination of AVHRR and METEOSAT［J］.Advances in Space Research，1995，16(10):141 －146.

［4］楊 軍，董超華，盧乃錳，等.中國(guó)新一代極軌氣象衛(wèi)星——風(fēng)云三號(hào)［J］.氣象學(xué)報(bào)，2009，67(4):501 －509.Yang J，Dong C H，Lu N M，et al.FY － 3A:The New Generation Polar－ orbiting Meteorological Satellite of China［J］.Acta Meteorological Sinica，2009，67(4):501 － 509(in Chinese with English Abstract).

［5］畢研盟，楊忠東，陸其峰，等.FY－3氣象衛(wèi)星中分辨率光譜成像儀和掃描輻射計(jì)11μm紅外窗區(qū)通道的比較［J］.紅外與毫米波學(xué)報(bào)，2009，28(5):330 －334.Bi Y M，Yang Z D，Lu Q F，et al.Comparison of Infrared 11μm Window Channels for Medium Resolution Spectral Imager and Visible and Infrared Radiometer on FY － 3 Meteorology Satellite［J］.Journal of Infrared and Millimeter Waves，2009，28(5):330 － 334(in Chinese with English Abstract).

［6］向大享，劉良明，韓 濤.FY－3A MERSI數(shù)據(jù)干旱監(jiān)測(cè)能力評(píng)價(jià)［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2010，35(3):334 －338.Xiang D X，Liu L M，Han T.Estimation of Drought Monitoring Ability of FY － 3A MERSI Data［J］.Geomatics and Information Science of Wuhan University，2010，35(3):334 －338(in Chinese with English Abstract).

［7］韓秀珍，鄭 偉，劉 誠(chéng)，等.基于MERSI和MODIS的太湖水體葉綠素 a含量反演［J］.地理研究，2011，30(2):291－300.Han X Z，Zhen W，Liu C，et al.Estimation of Chlorophyll－ a Using MERSI and MODIS Images in Taihu Lake，China［J］.Geographical Research，2011，30(2):291 － 300(in Chinese with English Abstract).

［8］Jiménez－ Mun～oz J C，Sobrino J A.A Generalized Single － channel Method for Retrieving Land Surface Temperature from Remote Sensing Data［J］.Journal of Geophysical Research，2003，108(D22):4688.

［9］覃志豪，Zhang M H，Arnon K，等.用陸地衛(wèi)星TM 6數(shù)據(jù)演算地表溫度的單窗算法［J］.地理學(xué)報(bào)，2001，56(4):456 －466.Qin Z H，Zhang M H，Arnon K，et al.Mono － window Algorithm for Retrieving Land Surface Temperature from Landsat TM6 data［J］.Acta Geographica Sinica，2001，56(4):456 － 466(in Chinese with English Abstract).

［10］Sobrino J A，Jiménez－ Mun～oz J C，Paolini L.Land Surface Temperature Retrieval from Landsat TM5［J］.Remote Sensing of Environment，2004，90(4):434 －440.

［11］Valor E，Casselles V.Mapping of Land Surface Emissivity from NDVI:Application to European，African，and South American Areas［J］.Remote Sensing of Environment，1996，57(3):167 －184.

［12］唐世浩，朱啟疆，周宇宇，等.一種簡(jiǎn)單的估算植被覆蓋度和恢復(fù)背景信息的方法［J］.中國(guó)圖象圖形學(xué)報(bào)，2003，8A(11):1304－1308.Tang S H，Zhu Q J，Zhou Y Y，et al.A Simple Method to Estimate Crown Cover Fraction and Rebuild the Background Information［J］.Journal of Image and Graphics，2003，8A(11):1304 － 1308(in Chinese with English Abstract).

［13］李苗苗，吳炳方，顏長(zhǎng)珍，等.密云水庫(kù)上游植被覆蓋度的遙感估算［J］.資源科學(xué)，2004，26(4):153 －159.Li M M，Wu B F，Yan C Z，et al.Estimation of Vegetation Fraction in the Upper Basin of Miyun Reservoir by Remote Sensing［J］.Resources Science，2004，26(4):153 － 159(in Chinese with English Abstract).

［14］唐仁茂，陳英英，葉建元.探空、地面及衛(wèi)星資料反演水汽含量的比較［J］.氣象科學(xué)，2010，30(3):373 －377.Tang R M，Chen Y Y，Ye J Y.The Comparision of Water Vapor Content Retrieved by Radiosonde，Ground Station and Satellite Data［J］.Scientia Meteorologia Sinica，2010，30(3):373 － 377(in Chinese with English Abstract).

［15］Kaufman Y J，Gao B C.Remote Sensing of Water Vapor in the Near IR from EOS/MODIS［J］.IEEE Trans Geosci Remote Sens，1992，30(5):871 －884.

［16］胡秀清，黃意玢，陸其峰，等.利用FY－3A近紅外資料反演水汽總量［J］.應(yīng)用氣象學(xué)報(bào)，2011，22(1):46 －56.Hu X Q，Huang Y B，Lu Q F，et al.Retrieving Precipitable Water Vapor Based on the Near－ infrared Data of FY －3A Satellite［J］.Journal of Applied Meterological Science，2011，22(1):46 －56(in Chinese with English Abstract).

［17］關(guān) 敏，郭 強(qiáng).光學(xué)像消旋系統(tǒng)在FY－3 MERSI圖像定位中的應(yīng)用［J］.應(yīng)用氣象學(xué)報(bào)，2008，19(4):420 －427.Guan M，Guo Q.Offsetting Image Rotation System in FY －3 MERSI’s Geolocation［J］.Journal of Applied Meterological Science，2008，19(4):420 －427(in Chinese with English Abstract).

［18］權(quán)維俊，葉彩華，劉勇洪.基于FY－3A MERSI數(shù)據(jù)的地表反射率反演參數(shù)化方法［C］//中國(guó)氣象學(xué)會(huì).第七屆全國(guó)優(yōu)秀青年氣象科技工作者學(xué)術(shù)研討會(huì)論文集.北京:［s.n.］，2010.Quan W J，Ye C H，Liu Y H.A Parameterization Method of Surface Reflectance Based on FY －3A/MERSI Data［C］//Chinese Meteorological Society.The Proceedings of Seventh National Outstanding Young Meteorological Science and Technology Workers Symposium.Beijing:［s.n.］，2010(in Chinese).

［19］周紅妹，周成虎，葛偉強(qiáng)，等.基于遙感和GIS的城市熱場(chǎng)分布規(guī)律研究［J］.地理學(xué)報(bào)，2001，56(2):189 －197.Zhou H M，Zhou C H，Ge W Q，et al.The Surveying on Thermal Distribution in Urban Based on GIS and Remote Sensing［J］.Acta Geographica Sinica，2001，56(2):189 －197(in Chinese with English Abstract).

Utilization of MATLAB to Realize LST Retrieval and Thematic Mapping from FY－3/MERSI Data

YANG He－qun，YIN Qiu，ZHOU Hong－mei，GE Wei－qiang
(Shanghai Center for Satellite Remote Sensing and Measurement Application，Shanghai 201199，China)

Currently，application－oriented researches on the data of Medium Resolution Spectral Imager(MERSI)，which is on board China’s new generation polar orbit meteorological satellite FY － 3，are very insufficient，due to the reason that the data as a new source have been delivered only since 2008.With the normal operation of FY－3 satellite system，it is necessary to develop an operational module for FY－3/MERSI regional land surface temperature(LST)retrieval and its post－processing，since LST is required for a wide variety of scientific studies but FY －3/MERSI’s operational LST products have not yet been provided by National Satellite Meteorological Center(NSMC).Based on an analysis of FY －3/MERSI L1 data’s HDF5 format and its channel characteristics，the authors selected the generalized single － channel algorithm developed by Jiménez － Mun～oz ＆Sobrino to directly realize the LST retrieval at 250 m spatial resolution with MATLAB programming and the thematic mapping of LST derivative products.This paper describes the parametric processes of LST retrieval algorithm in detail，which include radiometric calibration，cloud detection，estimation of two intermediate parameters－ surface emissivity and atmospheric water vapor，and calculation of thermal indexes from LST.On these bases，an automatic flowchart for FY －3/MERSI LST retrieval and thematic mapping was established.Experimental results of this flowchart applied in Shanghai thermal environmental monitoring show that it can process FY－3/MERSI L1 data in a fast，real－ time and automatic way，thus suitable for operational products producing and sharing，with the saving of human resources.It is also proved that FY －3/MERSI data and various forms of LST products can reveal the spatial pattern of Shanghai thermal field and the urban heat island effect more finely and intuitively.

FY－3/MERSI;land surface temperature(LST);retrieval;thematic mapping;MATLAB

TP 79

A

1001－070X(2012)04－0062－09

2011－10－20;

2011－12－22

上海市科委重點(diǎn)項(xiàng)目“上海城市熱環(huán)境災(zāi)害遙感監(jiān)測(cè)評(píng)估預(yù)警系統(tǒng)研究”(編號(hào):10511501700)、上海市科委重點(diǎn)海洋科技(臨港新城)專項(xiàng)“上海海洋區(qū)域海霧災(zāi)害預(yù)警預(yù)報(bào)系統(tǒng)研究”(編號(hào):10DZ1210802)共同資助。

10.6046/gtzyyg.2012.04.11

楊何群(1981－)，女，碩士，工程師，主要從事城市熱島、臺(tái)風(fēng)衛(wèi)星氣象應(yīng)用研究。E－mail:hequn_yang@163.com。

(責(zé)任編輯:刁淑娟)