路 茜，張鐵寶

（四川省地震局，成都610041）

0 引言

在地震前地?zé)岙惓５脑缙谘芯恐校梢园l(fā)現(xiàn)不少關(guān)于地震、地質(zhì)構(gòu)造與地表及淺層地溫關(guān)系研究的記錄[1-3]。近年來，衛(wèi)星遙感技術(shù)在現(xiàn)代光學(xué)、航天技術(shù)和計(jì)算機(jī)技術(shù)的基礎(chǔ)上迅速地發(fā)展了起來。在地震預(yù)報(bào)和構(gòu)造分析上的應(yīng)用也取得了初步的進(jìn)展，由于其具有覆蓋范圍廣、實(shí)時監(jiān)測及較高的空間分辨率等優(yōu)點(diǎn)，成為了觀測、分析地表熱輻射的新技術(shù)手段。地表溫度是地表與大氣相互作用過程中的一個重要的物理參數(shù)。陸面溫度遙感反演一直是許多學(xué)者力圖解決的一個難題。其難點(diǎn)在于存在其他大氣成分、地表層熱狀況以及熱存儲釋放過程的復(fù)雜性和熱傳遞方式等因素的影響。盡管如此，隨著遙感技術(shù)的不斷發(fā)展，應(yīng)用衛(wèi)星紅外資料對陸面溫度的研究已取得了很大的發(fā)展，陸面 “分裂窗”方法使得陸面溫度反演由難變易。

1 MODIS數(shù)據(jù)前期處理

衛(wèi)星遙感熱紅外亮溫信息用于地震前兆異常分析研究的難點(diǎn)之一來自云層的干擾來自地表的熱輻射信息難以穿透云層被衛(wèi)星傳感器接收到.若在溫度反演中引入有云數(shù)據(jù)，則得到的反演結(jié)果存在干擾，其中必然存在云頂溫度，即比實(shí)際的地表溫度低很多。MODIS數(shù)據(jù)有36個光譜通道，多光譜聯(lián)合云檢測是MODIS數(shù)據(jù)的應(yīng)用研究特色之一。 例如， 波段 0.65 μm、 0.85 μm、 1.38 μm、1.6 μm、 8.6 μm 和 11 μm 的組合應(yīng)用就可以非常準(zhǔn)確地得到不同云層的特性。多光譜的研究應(yīng)用可以將一景數(shù)據(jù)的云和晴空分為多個不同的類別。云和雪在0.65 μm波段看起來極其相似，而在1.6 μm波段的圖像中卻大不相同；位于對流層之上的1.38 μm波段圖象對于高空的冰晶云特別敏感；8.6 μm和11 μm的組合應(yīng)用可以區(qū)分卷云和低空的可降水云[4]。根據(jù)有云數(shù)據(jù)特征，利用通道1及通道31閾值計(jì)算進(jìn)行云檢測，本文中的數(shù)據(jù)范圍為中國大陸西部地區(qū)，采用了通道1反射率大于0.2與通道31亮溫值小于273 K，該閾值會隨著季節(jié)的變化作調(diào)整。

2 溫度數(shù)據(jù)分裂窗法

為從6個熱紅外地表亮溫波段中獲取更為統(tǒng)一的數(shù)據(jù)源，并探尋多種MODIS數(shù)據(jù)應(yīng)用方法，現(xiàn)對熱紅外亮溫?cái)?shù)據(jù)反演為溫度數(shù)據(jù)的算法實(shí)現(xiàn)做初步研究。本文的方法是分裂窗地表溫度普遍性算法。我們已知地表光譜反射率，根據(jù)多通道海面溫度算法修正大氣因素就能得到陸地表面溫度[5]。這種方法跟其他方法的不同之處在于：①不要求精確的大氣廓線；②不需要坐標(biāo)系中模擬輻射轉(zhuǎn)換；③其精確性取決于表面輻射率。分裂窗法地表溫度普遍性算法采用一個簡單模式，即使在處理過程中模擬輻射轉(zhuǎn)換不能節(jié)省處理時間，處理時間也不長。該方法被作為MODIS地表溫度反演算法的首要算法，是基于以下關(guān)于大氣、地表溫度和地表輻射的幾點(diǎn)考慮的：① 由于大氣（特別是水汽廓線）在垂直、水平方向上的變化幅度大，要計(jì)算其他相對濕氣廓線，使之精確度達(dá)到誤差10%以內(nèi)是很困難的；②地表溫度隨著地區(qū)和時間變化，且與地面空氣溫度也聯(lián)系較小，不同地表覆蓋類型白天和夜晚溫度可能相差10°C以上[6]；③ 多數(shù)地表覆蓋類型的輻射率在MODIS 31、32通道中是相對穩(wěn)定的，如圖1所示。

圖1 各種不同的地表覆蓋類型在MODIS 31、32波段的平均輻射率（ε31、ε32分別為31、32波段的輻射率）Fig.1 Average radiation rate of MODIS band 31， 32 with various land cover types （ε31、 ε32are radiation rate of band31， 32）

3 NDVI植被歸一化指數(shù)理論

論經(jīng)驗(yàn)性植被指數(shù)是根據(jù)葉子的典型光譜反射率特征得到的。由于葉綠素吸收在藍(lán)色（470 nm）和紅色（670 nm）波段最敏感，可見光波段的反射能量很低。幾乎所有的近紅外（NIR）輻射都被散射，很少吸收，且散射程度因葉冠光學(xué)特性、結(jié)構(gòu)特性而異。因此紅色波段和近紅外波段的反差是對植被量很敏感的度量波段。無植被或少植被區(qū)反差最小，中等植被區(qū)的反差是紅色和近紅外波段的變化結(jié)果，高植被區(qū)則只有近紅外波段對反差有貢獻(xiàn)，紅色波段趨于飽和，不再變化。這種反差對比可以用比值（NIR/red）， 差分（NIR-red）， 線性組合（x1×red+x2×NIR）或上述三者的組合來增強(qiáng)。 一般來說，植被指數(shù)就是對這種對比的度量，也是葉冠結(jié)構(gòu) （LAI，LAD）和生物參數(shù) （FAPAR）的函數(shù)。Deering[7]提出 “歸一化差分植被指數(shù)”（NDVI），將比值限定在 [-1，1]范圍內(nèi)：

經(jīng)過了輻射校正和幾何初校正的MODIS數(shù)據(jù)源，再進(jìn)一步對每日、每景圖像進(jìn)行幾何精校正、除條紋、去云等處理，進(jìn)而進(jìn)行NDVI計(jì)算及合成。計(jì)算公式為：

其中b1、b2為MODIS的第1、2波段。

4 溫度計(jì)算的實(shí)現(xiàn)

4.1 分裂窗法公式參數(shù)選取

經(jīng)過對NDVI原理的初步了解以后，下面介紹運(yùn)用NDVI極限方法來進(jìn)行分裂窗法地表溫度反演的方法。

分裂窗法的經(jīng)典公式為[8]：

式（3）中， ε=0.5（ε31+ε32）， Δε=ε31-ε32。

利用NDVI極限方法來求解ε31、ε32的值，其關(guān)系為： 如果： NDVI＜NDVImin， 則 ε31=0.979 5-0.056 5 ρ1，ε32=0.981 5-0.027 5 ρ1， ρ1為第 1 波段反射率； 如果： NDVImin≤NDVI≤NDVImax， 則： ε31=0.968+0.021 fV，ε32=0.974-0.015 fV。fV=如果NDVI＞NDVImax,， 則 ε31=ε32=0.989 。

NDVImin、NDVImax是通過大氣修正的裸露土壤和葉面積指數(shù) LAI（Leaf Area Index）的 大小極NDVI值，取決于所研究的數(shù)據(jù)圖像。

從前面談到的植被指數(shù)理論可知，植被指數(shù)是對地表植被覆蓋率的簡單、有效和經(jīng)驗(yàn)的衡量。它會隨時間、空間變化。對于每個MODIS數(shù)據(jù)點(diǎn)來說，不同時間的植被指數(shù)就有所不同，因?yàn)橐^準(zhǔn)確的計(jì)算地表溫度，需要根據(jù)NDVI的變化情況來正確地選擇參數(shù)ε31、ε32。這里我們所采用的方法是建立NDVI最值數(shù)據(jù)庫來存儲每個數(shù)據(jù)點(diǎn)的NDVImin、NDVImax，在讀取數(shù)據(jù)文件時獲取每個MODIS數(shù)據(jù)點(diǎn)的NDVI值與當(dāng)月NDVImin、NDVImax進(jìn)行比較并替換，在此基礎(chǔ)上查找上月該點(diǎn)所有 數(shù)據(jù)的最大值和最小值和當(dāng)前值進(jìn)行比較從而選擇合適的參數(shù)ε31、ε32，這樣就可以保證計(jì)算溫度數(shù)據(jù)文件的實(shí)時準(zhǔn)確性。

由于接收到的各數(shù)據(jù)點(diǎn)的觀測角度都不盡相同，因此分裂窗法地表溫度普遍性算法中公式的系數(shù)也不相同（如表1）。

MODIS掃描觀測寬度達(dá)2 330 km，觀測角度達(dá)到±58°，接收到的數(shù)據(jù)與經(jīng)度線，衛(wèi)星運(yùn)行軌道方向與經(jīng)度存在一定的夾角。接收到的MODIS數(shù)據(jù)經(jīng)過定標(biāo)定位及等經(jīng)緯度投影以后，再按照所需要的區(qū)域截取包含所需熱紅外地表波段的數(shù)據(jù)文件，由此得到的亮溫?cái)?shù)據(jù)文件中的每個像元（Pixels，注下文用p表示）數(shù)據(jù)的觀測角度均不同，多年數(shù)據(jù)的海量像元參與計(jì)算就需要一個龐大的觀測角度數(shù)據(jù)庫。為了便于計(jì)算，本文采用了各觀測角度的系數(shù)取平均值的辦法，由此得到的平均值系數(shù)如上表1中所示。

4.2 溫度計(jì)算結(jié)果

綜上所述步驟， 各參數(shù)值 （ε31、ε32、 A1、 A2、A3、B1、B2、B3）帶入分裂窗法公式中進(jìn)行計(jì)算，由此得到各點(diǎn)亮溫?cái)?shù)據(jù)的地表溫度LST（Land Surface Temperature）。例如：下圖2為由2004年3月29日北緯30°、東經(jīng)100°為中心，長寬各為2 000 km的去云亮溫文件（圖 2a）與NDVI文件（圖2b）而得到的LST數(shù)據(jù)文件（圖 2c）。

4.3 與MODIS亮溫產(chǎn)品比較

經(jīng)過分裂窗算法計(jì)算出的溫度文件，首先與MODIS地表亮溫產(chǎn)品進(jìn)行比較。比較時段為2004年3月至2006年11月。在MODIS地表亮溫產(chǎn)品包含的6個熱紅外地表亮溫波段，通道21、22、23的數(shù)據(jù)形態(tài)較一致，通道31、32的數(shù)據(jù)形態(tài)較一致，而通道20則與他們都不盡相同。在此給出幾個像元點(diǎn)的反演溫度與代表性的通道20、通道21、通道31的比較結(jié)果（圖3）。

表1 不同觀測角度與分裂窗算法系數(shù)的對應(yīng)表Table 1 Corresponding relation between different observational angle and split window algorithm

圖2 溫度文件的反演過程示例圖Fig.2 The example of assimilation procedure of the temperature files

圖3 反演溫度與MODIS地表溫度產(chǎn)品的曲線圖（4p）Fig.3 The sequence diagram of assimilation of LST and MODIS thermal infrared band（4p）

圖3顯示，較高緯度地區(qū)（a點(diǎn)）的反演溫度及各數(shù)據(jù)產(chǎn)品之間相關(guān)性較中緯度地區(qū)（b、c、d點(diǎn)）高，可能因?yàn)閍點(diǎn)位于新疆、青海交界地區(qū)，植被覆蓋少，干擾較少，年變清晰。而b、c、d點(diǎn)分別位于緬甸熱帶雨林地區(qū)、四川甘孜州理塘地區(qū)、四川和貴州交界地區(qū)，植被覆蓋較高，或云多，地表溫度實(shí)測存在干擾因素較多。

經(jīng)過多點(diǎn)的反復(fù)比較（表2）可以看出，由分裂窗算法計(jì)算出的反演溫度與21和31通道的MODIS數(shù)據(jù)存在較好的相關(guān)性，其中與31通道相關(guān)系數(shù)接近1。而20通道的數(shù)據(jù)則存在較大的波動性，相關(guān)性較差。

表2 反演溫度與MODIS亮溫產(chǎn)品之間的相關(guān)系數(shù)Table 2 The correlation coefficient between assimilation of LST and MODIS thermal infrared band

5 MODIS反演溫度數(shù)據(jù)與地溫?cái)?shù)據(jù)的初步相關(guān)性分析

5.1 氣象地溫?cái)?shù)據(jù)和反演溫度數(shù)據(jù)收集情況

本文所采用的氣象數(shù)據(jù)地溫資料為地面氣象站觀測所得，這些觀測站是位于我省鮮水河、安寧河、龍門山等主要斷裂帶區(qū)域內(nèi)的24個氣象臺站觀測記錄到的地表溫度（0 cm）以及地下淺層（包括5～320 cm）溫度，由于與反演所得到的地表層溫度做比較，因此本文采用地表溫度（0 cm）和地下淺層（40 cm）的日均值數(shù)據(jù)及其平均所得到的月均值數(shù)據(jù)。大部分臺站的資料起止時間為2004年1月至2006年7月。

為了配合氣象臺站觀測點(diǎn)地表溫度數(shù)據(jù)的時間， MODIS反演溫度數(shù)據(jù)也選定了上述相同時間數(shù)據(jù)作為分析數(shù)據(jù)。MODIS熱紅外波段數(shù)據(jù)的空間分辨率為1 km，而地面氣象臺觀測地溫是用以代表一個區(qū)域的地溫。因此本文根據(jù)氣象觀測臺的經(jīng)緯度分別選取了MODIS反演溫度數(shù)據(jù)的不同大小區(qū)域的均值進(jìn)行相關(guān)性對比分析。反演溫度數(shù)據(jù)經(jīng)過多天數(shù)據(jù)拼接減少云數(shù)據(jù)覆蓋和去云處理。

5.2 反演溫度（20×20 p）與各深度層（0 cm、40 cm）地溫日均值的相關(guān)性

利用前面分裂窗法計(jì)算所得的多年MODIS反演溫度與地溫的日均值數(shù)據(jù)進(jìn)行對比分析。在此給出甘孜、宜賓、馬爾康、康定四個地溫觀測點(diǎn)與反演溫度的對比分析。本研究利用線性相關(guān)系數(shù)計(jì)算公式：

計(jì)算反演溫度與地溫的相關(guān)性。圖4給出了反演溫度曲線及四個地溫觀測點(diǎn)日均值曲線。

圖4可以看出，地溫日均值曲線較反演溫度曲線穩(wěn)定，波動小。這可能是由衛(wèi)星傳感器接收信號受到大氣顆粒、太陽輻射等多方面因素的影響造成的。反演溫度與地溫都出現(xiàn)了在夏季時波動較大、冬季時較為穩(wěn)定的現(xiàn)象。0 cm地溫與反演溫度的相似性比40 cm地溫與反演溫度的相似性好。表3給出了反演溫度與4個地區(qū)地溫的相關(guān)系數(shù)。

圖4 反演溫度與甘孜（a）、宜賓（b）、馬爾康（c）、康定（d）0 m、 0.4 m 地溫的日均值曲線圖Fig.4 The daily mean profiles of assimilation of land surface temperature （LST） and geotemperature of Ganzi， Yibin， Maerkang，Kangding area

表3 反演溫度（20×20p）與地溫（日均值）之間的相關(guān)系數(shù)Table 3 The correlation coefficient between assimilation of LST （20×20p） and daily mean geotemperature

5.4 反演溫度（20×20 p）與各深度層（0 cm、40 cm）地溫月均值的相關(guān)性

用日均值數(shù)據(jù)進(jìn)行月平均處理，可得到月均值時序圖。由于篇幅的原因，在此僅給出月均值相關(guān)性分析的結(jié)果。從月均值的相關(guān)系數(shù)來看，比日均值有了明顯的提高（圖4）。

5.3 不同大小區(qū)域范圍反演溫度與地溫日均值相關(guān)性的比較分析

地溫?cái)?shù)據(jù)來自與氣象臺站的測定，而反演溫度則是依據(jù)氣象臺站的準(zhǔn)確經(jīng)緯度來確定的數(shù)據(jù)點(diǎn)位的單個像元或多個像元的均值。衛(wèi)星數(shù)據(jù)受到大氣、輻射、地表形態(tài)等多種因素的影響，單個像元的數(shù)據(jù)可用率較低，可靠性也較低，通過多像元的數(shù)據(jù)求均值可以弱化干擾，更能反映區(qū)域的實(shí)際情況。本文就這個問題，采用不同大小區(qū)域范圍反演溫度與地溫進(jìn)行對比、相關(guān)性分析。表5給出了結(jié)果，由此可以看出除了馬爾康地區(qū)以外，其他三個地區(qū)經(jīng)過多像元數(shù)據(jù)平均后，一定程度地弱化了干擾因素，相關(guān)系數(shù)有所提高。甘孜、馬爾康、康定屬于高原地區(qū)，植被覆蓋率較低、地表覆蓋類型較單一，這可能是其實(shí)測地溫與各大小區(qū)域范圍反演溫度相關(guān)系數(shù)相差不大的原因，而宜賓位于四川盆地，長江流域，地貌以中低山地和丘陵為主，從表5中可以看出最大區(qū)域范圍比單個像元與反演溫度的相關(guān)系數(shù)相差0.2。

表4 反演溫度（20×20 p）與地溫（月均值）之間的相關(guān)系數(shù)Table 4 The correlation coefficient between assimilation of LST （20×20 p） and monthly mean

表5 不同大小區(qū)域范圍反演溫度與地溫（日均值）之間的相關(guān)系數(shù)Table 5 The correlation coefficient between assimilation of LST of different regional ranges and monthly mean geotemperature

6 討論

（1）分裂窗法地表溫度算法實(shí)現(xiàn)中存在的誤差因素：公式中的7個系數(shù)的確定存在一定的誤差，本文在確定的時候是采取的取均值的方法，但實(shí)際觀測角度不同時，該計(jì)算系數(shù)是有一定誤差的；還有的選取也不是非常準(zhǔn)確，由于的季節(jié)性變化，其值也是不穩(wěn)定的，值域范圍也不容易確定，該如何來確定其一定時間范圍的最大、最小值，使得參數(shù)更加準(zhǔn)確適合，還需要反復(fù)的試驗(yàn)。

（2）據(jù)上文分析得出，地表溫度與遙感數(shù)據(jù)反演溫度的日值相關(guān)性系數(shù)僅在0.6左右，而月值相關(guān)性卻在0.8左右。這可能受到兩方面原因的影響：①地溫?cái)?shù)據(jù)是由氣象臺站在每天多個時段測量所得溫度值的平均值，而遙感衛(wèi)星數(shù)據(jù)則是每天某個特定的時刻接收到的，本文所用的TERRA衛(wèi)星數(shù)據(jù)是在上午10點(diǎn)至12點(diǎn)的實(shí)測數(shù)據(jù)。這兩者之間必定存在差異。② 衛(wèi)星遙感屬于 “面”觀測，而地面溫度測點(diǎn)屬于 “點(diǎn)”觀測，通過衛(wèi)星遙感獲得的地表溫度代表一定空間區(qū)域內(nèi)的平均溫度，通過地表實(shí)測的方式獲得的地表溫度代表測點(diǎn)的溫度。遙感與實(shí)測之間，空間尺度存在不確定性[9]，這可能為二者帶來差異。③ 遙感熱紅外數(shù)據(jù)的每日數(shù)值隨機(jī)波動遠(yuǎn)大于地表測定溫度，可能是因?yàn)檫b感數(shù)據(jù)受到大氣吸收、植被覆蓋、人類活動等干擾引起的，而進(jìn)行月均值處理以后這些干擾明顯被弱化了。將遙感衛(wèi)星熱紅外數(shù)據(jù)應(yīng)用于地震監(jiān)測預(yù)報(bào)工作時，應(yīng)盡量采取一些去噪聲的處理，如小波變換、滑動濾波等，以便于在應(yīng)用中提取地殼內(nèi)部變化引起的熱紅外信息。

致謝：本文所采用的地表及淺層地表溫度數(shù)據(jù)來自四川省氣象局，在此對他們的大力支持表示感謝。

參考文獻(xiàn)：

[1]尤傳俠.淺層地溫在滇西南地震預(yù)報(bào)中的應(yīng)用 [J].地震研究，1989，12（3）：254-259.

[2]張治洮.淺層地溫變化與強(qiáng)震、旱澇的關(guān)系[J].中國減災(zāi),1992，2（3）：50-50.

[3]劉小鳳.青藏高原北部淺層地溫異常特征及中短期地震預(yù)測[J].高原地震，2002，14（2）：1-7.

[4]劉玉潔，楊忠東.MODIS遙感信息處理原理與算法[M].北京：科學(xué)出版社，2001.

[5]Wan Z.J.Dozier.A.Generalized split-window algorithm for retrieving land surface temperature measurement from space[J].IEEE Trans Geosci Remote sensing，1996（34）：892-905.

[6] Betts A K，Ball J H，Beljaars A C M，et al.The land surface-atmosphere interaction：A review based on observation and global modeling perspectives[J].J.Geophys.Res，1996，101（D3）：7209-7225.

[7] Deering D W.Rangeland Reflectance Characteristics measured by aircraft and Spacecraft Sensor[D].Ph.D.Diss,Texas A&M University，1978.

[8]Zheng M W.MODIS Land-Surface Temperature Algorithm Theoretical Basis Document （LST ATBD）.Institute for Computational Earth System Science [D].University of California，Santa Barbara，1999.

[9] 陳順云，劉培洵，劉力強(qiáng)，等.遙感與實(shí)測地表溫度的對比分析及在地震研究中的意義[J].地球物理學(xué)報(bào)，2011，54（3）：747-755.