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

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

0 引言

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

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

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

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

2 原理

2.1 輻射定標(biāo)

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

2.2 輻射傳輸方程法

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

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

其中，L↑、L↓可通過(guò)影像的中心緯度與成像時(shí)間在網(wǎng)站（http://atmcorr.gsfc.nasa.gov）獲取（張菊等，2020）。

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

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

2.3 地表比輻射率及植被覆蓋度計(jì)算

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

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

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

2.4 最大似然分類

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

可導(dǎo)出貝葉斯定理公式

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

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

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

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

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

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

為反演結(jié)果的準(zhǔn)確性，本文使用監(jiān)督分類方法中的最大似然分類方法將影像進(jì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）。本文使用輻射傳輸方程法對(duì)研究區(qū)域影像進(jìn)行地表溫度反演，反演結(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)計(jì)了6 種地表物所在區(qū)域內(nèi)的地表溫度及植被覆蓋度的平均值與標(biāo)準(zhǔn)差，統(tǒng)計(jì)結(jié)果如表1、2 所示。

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

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

表2 植被覆蓋度Table 2 Vegetation coverage

4 結(jié)論與討論

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

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