包寶小 席文雅 王志勇 李騰宇 賈寶金 王 慧

（中國內(nèi)蒙古自治區(qū)137400 內(nèi)蒙古自治區(qū)地震局烏蘭浩特地震臺）

0 引言

自1988 年Gomy 等發(fā)現(xiàn)1984 年加茲利地震震前衛(wèi)星熱紅外異常后（張帆等，2016），有關(guān)衛(wèi)星熱紅外的研究逐漸開展起來，研究結(jié)果已在震后緊急救援、災(zāi)情監(jiān)測、災(zāi)后重建等方面得到應(yīng)用（葉秀薇等，2010）。

20 世紀80 年代末，馬瑾在國內(nèi)最先開展了衛(wèi)星熱紅外與地震間相關(guān)性的研究（張治廣等，2017）；在此基礎(chǔ)上，強祖基與徐秀登等合作開展了熱紅外遙感研究（徐好民，2002）。隨著遙感技術(shù)的發(fā)展，衛(wèi)星數(shù)據(jù)質(zhì)量不斷提高，利用氣象衛(wèi)星資料，獲取海面溫度的技術(shù)趨于成熟，海面溫度遙感反演技術(shù)的成功使得地表溫度遙感反演成為遙感領(lǐng)域又一研究熱點（黃啟廳，2010）。近年來，地表溫度異常、植被覆蓋度變化與地震異常間相關(guān)性的研究逐漸興起。胡光庭等（2017）分析了陽江市地表溫度變化特征，認為植被覆蓋度的增加能很好地降低地表溫度；王蕾（2014）等基于遙感影像發(fā)現(xiàn)，植被覆蓋度對緩解城市熱島效應(yīng)具有一定作用，說明地表溫度與植被覆蓋度間具有較好的的負相關(guān)性；龐光輝等（2016）對沈陽市植被覆蓋度降溫效應(yīng)進行了分析，發(fā)現(xiàn)植被覆蓋度每增加10%，地表溫度可降低0.71℃。前人研究結(jié)果大多適用于城市生態(tài)環(huán)境規(guī)劃和建設(shè)，且反演方法較簡單，對地表物未進行詳細分類而直接提取地表溫度數(shù)據(jù)，且未統(tǒng)計地表溫度平均值及標準差。本文采用輻射傳輸方程算法對研究區(qū)域地表溫度進行了反演，為提高反演結(jié)果的準確性，使用最大似然分類算法對地表物進行了耕地、林地、草地、水體、建設(shè)用地、未利用土地的詳細分類，依據(jù)分類結(jié)果及反演結(jié)果，分別計算了6 種地表物所在區(qū)域內(nèi)的地表溫度、植被覆蓋度的平均值、標準差，分析總結(jié)了研究區(qū)域地表溫度與植被覆蓋度間的相關(guān)性。

1 研究區(qū)及數(shù)據(jù)

烏蘭浩特市地理坐標為45°22′—46°18′N、121°51′—122°20′E，總面積2 353.5 km2，屬中溫帶半干旱季風氣候，四季分明，氣溫的年差、日差都較大，年平均氣溫5℃，年降水量為440 mm 左右。烏蘭浩特地震臺位于內(nèi)蒙古自治區(qū)興安盟烏蘭浩特，是興安盟唯一有人值守的測震數(shù)字化臺站，承擔著興安盟區(qū)域內(nèi)地震監(jiān)測預(yù)測任務(wù)（包寶小等，2019）。本文選取過境時間為2010、2015、2020 年3 年的9 月11 日2 時烏蘭浩特地震臺監(jiān)測區(qū)域的Landsat 5、Landsat 8 遙感影像數(shù)據(jù)（云量較少，符合使用需求），進行輻射定標、大氣校正預(yù)處理，使用烏蘭浩特矢量圖進行矢量裁剪。利用影像數(shù)據(jù)第1、2、3、4、5、7 波段提取NDVI（歸一化植被指數(shù)）及土地利用分類，利用Landsat 5 的第6 波段及Landsat 8的第10 波段反演地表溫度數(shù)據(jù)。

2 原理

2.1 輻射定標

輻射定標是將傳感器記錄的無量綱DN 值轉(zhuǎn)換成具有實際物理意義的大氣頂層輻射亮度或反射率。其原理是建立數(shù)字量化值與對應(yīng)視場中輻射亮度之間的定量關(guān)系，以消除傳感器本身產(chǎn)生的誤差（央金卓嘎等，2019）。

2.2 輻射傳輸方程法

衛(wèi)星傳感器接收到的熱紅外輻射亮度Lλ由3 部分組成：大氣向上輻射亮度L↑、地面的真實輻射亮度經(jīng)過大氣層之后到達衛(wèi)星傳感器的能量、大氣向下輻射到達地面后反射的能量L↓。衛(wèi)星傳感器接收到的熱紅外輻射亮度Lλ的表達式可寫為（輻射傳輸方程）

式中，ε為地表比輻射率；Ts為地表真實溫度；B(Ts)為通過反映物體溫度和輻射強度之間關(guān)系的普特朗克公式推算所得的溫度T相對應(yīng)的黑體輻射亮度；τ為大氣在熱紅外波段的透射率；L↑、L↓為大氣的向上、向下輻射亮度（崔耀輝，2018）。由此可得B(Ts)的計算公式如下

其中，L↑、L↓可通過影像的中心緯度與成像時間在網(wǎng)站（http://atmcorr.gsfc.nasa.gov）獲?。◤埦盏龋?020）。

運用普朗克公式的反函數(shù)，可推導出亮度溫度TR

其中，K1、K2的值可在影像MLT.txt 文件中獲得。

2.3 地表比輻射率及植被覆蓋度計算

地表比輻射率又稱發(fā)射率，是指在同一溫度下地表發(fā)射的輻射量與一黑體發(fā)射的輻射量間的比值，是利用熱紅外遙感獲取地表溫度的必不可少的參數(shù)。針對Landsat 影像的第10 段熱紅外波段，使用Sobrino 提出的NDVI 閾值法計算地表比輻射率（齊廣慧等，2020）

其中，Pv為植被覆蓋度，其公計算式如下

其中，NDVI 為歸一化植被指數(shù)；NDVIs為完全是裸土或無植被覆蓋區(qū)域的NDVI；NDVIv為完全植被覆蓋的像元的NDVI。一般取經(jīng)驗值，NDVIs=0.05，NDVIv=0.70（蔡賢等，2020）。

2.4 最大似然分類

為提高地表溫度反演的準確性，使用最大似然分類方法對影像進行地表物的分類。在2 類或多類判別中，假定各類分布函數(shù)為正態(tài)分布，并選擇訓練區(qū)，計算各待分類樣區(qū)的歸屬概率而進行分類的圖像分類方法稱為貝葉斯分類方法，本文根據(jù)貝葉斯準則對遙感影像進行分類。貝葉斯定理：假設(shè)B1、B2，…互斥且構(gòu)成一個完全時間，A伴隨它們出現(xiàn)，已知它們分別發(fā)生的先驗概率P(Bi)（i=1，2，…）及A的條件概率P(Bi)，則可以得到事件A的后驗概率P(Bi|A)。由概率乘法公式

可導出貝葉斯定理公式

設(shè)有s個類別，用ω1、ω2、…、ωs來表示，每個類別發(fā)生的概率（先驗概率）分別為P(ω1)、P(ω2)、…、P(ωs)。設(shè)有未知類別的樣本X，其類條件概率分別為

則根據(jù)貝葉斯定理可以得到樣本X出現(xiàn)的后驗概率為

此時，以樣本X 出現(xiàn)的后驗概率作為判別函數(shù)來確定樣本X的所屬類別，其分類準則為：如果P(ωi|X)=max1

3 計算結(jié)果與分析

利用影像數(shù)據(jù)的紅光波段及近紅外波段提取了2010、2015、2020 年的歸一化植被指數(shù)（NDVI）（圖1）。NDVI 的取值范圍-1—1，當NDVI 值＜0 時，表示水域，即烏蘭浩特及科爾沁右翼前旗區(qū)域間的歸流河；當NDVI 值為0 時，表示地表被裸土及巖石覆蓋，為建筑用地，即烏蘭浩特主城區(qū)區(qū)域；當NDVI 值為1 時，表示植被覆蓋豐富，即草地、林地及耕地、主城區(qū)四周區(qū)域。由圖1 可見，2010—2020 年NDVI 低值區(qū)域沒有太大變化，高值區(qū)域明顯向四周擴散。

圖1 2010 年（a）、2015 年（b）、2020 年（c）NDVI 值Fig.1 NDVI value of the year of 2010,2015,and 2020

為反演結(jié)果的準確性，本文使用監(jiān)督分類方法中的最大似然分類方法將影像進行6 種地表物分類，分為耕地、林地、草地、水體、建設(shè)用地及未利用地，分類結(jié)果如圖2所示。

圖2 2010 年（a）、2015 年（b）、2020 年（c）土地利用結(jié)果Fig.2 Results of land use in 2010,2015 and 2020

遙感影像Landsat TM/ETM+的第6、10、11波段屬于熱紅外波段，可用于地表溫度反演，地表溫度反演算法主要有輻射傳輸方程法、單窗算法、單通道算法、基于影像的反演算法（張菊等，2020）。本文使用輻射傳輸方程法對研究區(qū)域影像進行地表溫度反演，反演結(jié)果如圖3 所示。

圖3 2010 年（a）、2015 年（b）、2020 年（c）輻射傳輸法地表溫度反演結(jié)果Fig.3 Retrieval results of land surface temperature by radiative transfer method in 2010,2015 and 2020

依據(jù)分類結(jié)果及反演結(jié)果影像，分別統(tǒng)計了6 種地表物所在區(qū)域內(nèi)的地表溫度及植被覆蓋度的平均值與標準差，統(tǒng)計結(jié)果如表1、2 所示。

表1 地表溫度（單位：℃）Table 1 Surface temperature

由表1、2 可見，水體因自身的物理特征，溫度最低；其次是耕地，溫度均值為29.5℃，標準差也相對較低；再次為耕地及林地；草地地表溫度均值最高，說明不同地表物之間地表溫度差異較大。建設(shè)用地與未利用土地相比，未利用土地地表溫度較高，建設(shè)用地上雖然生長的植物較少，但烏蘭浩特市綠化較好，所以建設(shè)用地地表溫度比未利用土地低。因建設(shè)用地的綠化結(jié)構(gòu)較復(fù)雜，由表1、表2 可見，建設(shè)用地地表溫度及植被覆蓋度的標準差較大。地表溫度較低的地表物類型，其植被覆蓋度相對較高；而地表溫度較高的地表物類型，植被覆蓋度則較低，這說明植被有很好的降溫作用。所以，進行地表溫度映震異常分析時要多方面考慮植被覆蓋度的影響。

表2 植被覆蓋度Table 2 Vegetation coverage

4 結(jié)論與討論

選取烏蘭浩特地震臺監(jiān)測區(qū)域2010、2015、2020 年3 年的9 月11 日2 時的Landsat 5、Landsat 8 遙感影像數(shù)據(jù)，經(jīng)過輻射定標、大氣校正、裁剪、拼接等預(yù)處理，提取NDVI（歸一化植被指數(shù)），計算植被覆蓋度，使用最大似然分類法進行土地利用分類，使用Landsat 5的第6波段及Landsat 8的第10波段反演了地表溫度。研究得到如下結(jié)論：①將林地、草地、耕地分為一類，按照地表溫度進行排序，從低到高為耕地、林地、草地，建設(shè)用地及未利用土地分為一類，未利用土地地表溫度高于建設(shè)用地；②按照植被覆蓋度從低到高排序結(jié)果為草地、林地、耕地；③因水體自身特征，地表溫度及植被覆蓋度都很低，除水體之外，建設(shè)用地地表溫度及植被覆蓋度標準差高于其他地表物；④由不同地表物反演的地表溫度有明顯差異，不同地表物植被覆蓋度亦各異；⑤地表溫度較低的地表物類型植被覆蓋度相對較高，地表溫度較高的相對較低，說明對于不同地表物類型植被均有降溫作用。

綜上所述，地表溫度異常不一定為地震短臨異常，需綜合考慮如下干擾因素：當?shù)乇頊囟瘸霈F(xiàn)增溫現(xiàn)象時，首先應(yīng)考慮是否植被覆蓋度降低所致，如表1 中2015 年與2020 年建設(shè)用地地表溫度均值出現(xiàn)上升現(xiàn)象，但并沒出現(xiàn)地震事件，而是市政規(guī)劃將草地變?yōu)榻ㄔO(shè)用地所致。當?shù)乇頊囟瘸霈F(xiàn)降溫現(xiàn)象時，需實地取證，以明確是否增添植被而導致地表溫度下降的，如表1 中草地、林地、耕地地表溫度逐年下降，但并沒出現(xiàn)地震事件，而是因為植被覆蓋度增加了。當?shù)乇頊囟劝邏K性異常時，需根據(jù)市政規(guī)劃文件，實地考察是否建設(shè)用地進行了重新規(guī)劃改造所致。本文給出的2010、2015、2020 年3 年的9 月地表溫度、植被覆蓋度的均值及標準差，所得結(jié)果可為研究地震熱紅外異常工作提供有效的數(shù)據(jù)支持及科學依據(jù)。