吳瑋瑩 單新建 屈春燕 李新艷

1)中國地震局地質(zhì)研究所，地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100029 2)寧夏回族自治區(qū)地震局，銀川 750001

0 引言

地震是造成財(cái)產(chǎn)損失和人員傷亡最為嚴(yán)重的自然災(zāi)害之一。預(yù)測(cè)未來地震發(fā)生的時(shí)間、 空間、 強(qiáng)度信息是預(yù)防地震災(zāi)害最為有效的手段(Mukhopadhyayetal.，2019)。隨著遙感技術(shù)的發(fā)展，大量利用熱紅外衛(wèi)星數(shù)據(jù)的地震異常研究表明： 地震前存在地表溫度(文翔等，2015； Ahamdetal.，2019； Khalilietal.，2019； Chenetal.，2020)、 長波輻射(Xiongetal.，2010； Kongetal.，2018； Zhongetal.，2020)等多種熱參數(shù)的異常變化現(xiàn)象，熱參數(shù)的時(shí)空分布特征對(duì)未來地震的震中、 發(fā)震斷層具有一定的指示作用。

由于熱參數(shù)受其他因素(如局部水熱條件、 季節(jié)、 氣候等)影響，其本身的變化具有周期性、 偶發(fā)性的特點(diǎn)。因此，采用合適的地震熱參數(shù)異常提取方法抑制其他因素對(duì)熱參數(shù)的影響、 突出地震因素的作用，對(duì)于地震熱參數(shù)異常變化研究是十分重要的。目前，ZS法(Z-score，也稱偏移指數(shù)法)和RST法(Robust Satellite Technology)是應(yīng)用較為廣泛的地震熱參數(shù)異常提取方法(Filizzolaetal.，2004； Tramutolietal.，2005； Ouzounovetal.，2011； Genzanoetal.，2015； Khalilietal.，2019)。ZS法和RST法均利用衛(wèi)星多年觀測(cè)的熱參數(shù)值的平均值作為熱參數(shù)正常變化的期望值，以多年觀測(cè)值的標(biāo)準(zhǔn)差描述觀測(cè)值的可變性以構(gòu)建背景場(chǎng)，以當(dāng)前觀測(cè)值和背景場(chǎng)的差作為分母部分，多年的標(biāo)準(zhǔn)差作為分子部分，利用該比值定義熱參數(shù)變化。不同的是，RST法在構(gòu)建背景場(chǎng)時(shí)，通過劃分地物覆蓋類型并求取不同地物類型的平均值，以降低地物覆蓋類型對(duì)最終結(jié)果的影響。ZS法和RST法都是基于統(tǒng)計(jì)分析的方法，原理清晰明確，操作較為簡便。

盡管大量研究成果表明地震前存在熱參數(shù)時(shí)空分布特征的異常變化，但由于不同地震的熱參數(shù)異常提取方法對(duì)同一參數(shù)提取的結(jié)果存在差異，使得熱參數(shù)異常的時(shí)空特征及其與地震活動(dòng)的關(guān)聯(lián)受到質(zhì)疑。因此，需要對(duì)不同地震熱參數(shù)異常提取方法的適用性和優(yōu)缺點(diǎn)進(jìn)行比較分析，但目前尚缺乏相關(guān)研究(Jiaoetal.，2018)。 本研究以2014年MW6.9 于田地震為典型震例，以2008年汶川MW7.9 地震為驗(yàn)證震例，同時(shí)利用ZS法和RST法提取典型震例震前3個(gè)月的地表溫度(Surface temperature，ST)和晴空條件下向外出射長波輻射(Outgoing longwave radiation，OLR)數(shù)據(jù)異常變化的時(shí)空特征，定性分析提取結(jié)果之間的差異。隨后，利用地震平靜年(無震級(jí)>5.0級(jí)地震發(fā)生的年份)數(shù)據(jù)對(duì)2種提取方法在異常變化敏感性、 對(duì)背景信息的抑制能力和對(duì)地震信息的指示性3個(gè)方面進(jìn)行了定量比較，并利用驗(yàn)證震例驗(yàn)證定量化比較的結(jié)果，詳細(xì)探討了造成差異的原因，總結(jié)了各方法的優(yōu)缺點(diǎn)和適用性，為后續(xù)進(jìn)一步探究地震震前熱參數(shù)異常的時(shí)空分布及其可能的物理機(jī)制奠定基礎(chǔ)。

1 數(shù)據(jù)和方法

1.1 ERA5再分析數(shù)據(jù)

本文主要對(duì)ZS法和RST法對(duì)同一震例、 不同熱參數(shù)在震前時(shí)空變化特征提取方面的適用性和可靠性進(jìn)行對(duì)比研究，選取地表溫度(Surface temperature，ST)和晴空條件下向外長波輻射(Clear-sky Outgoing longwave radiation，OLR)作為研究參數(shù)。地表溫度是地表和大氣邊界層之間能量交換的重要影響因素，能夠更直觀地反映地震孕育過程中的熱輻射變化。晴空條件下向外長波輻射是地球表面(陸地和海洋)、 大氣反射、 吸收和發(fā)射的長波輻射最終離開地球進(jìn)入太空的輻射能量，是整個(gè)地面-大氣系統(tǒng)熱輻射的綜合反映。

熱參數(shù)數(shù)據(jù)來自歐洲中期天氣預(yù)測(cè)中心(ECMWF，European Centre for Medium-Range Weather Forecasts)第5代再分析數(shù)據(jù)集ERA5(1)https： ∥cds.climate.copernicus.eu/cdsapp#！/dataset/reanalysis-era5-single-levels？tab=overview。。該數(shù)據(jù)集是通過數(shù)據(jù)同化和再分析技術(shù)將先進(jìn)的預(yù)報(bào)模型得到的預(yù)報(bào)結(jié)果與多種觀測(cè)手段得到的大量觀測(cè)數(shù)據(jù)結(jié)合為精度更高、 覆蓋更廣的數(shù)據(jù)集。與傳統(tǒng)的衛(wèi)星遙感數(shù)據(jù)(如MODIS(Moderate-resolution Imaging Spectroradiometer)數(shù)據(jù))相比，再分析數(shù)據(jù)具有以下優(yōu)勢(shì)： 1)時(shí)間序列完整，目前能夠提供從1979年至今的完整數(shù)據(jù)集； 2)空間上無空值，能夠完整反映熱參數(shù)的空間分布特征。與其他再分析數(shù)據(jù)集(如NCEP(National Centers for Environmental Prediction)數(shù)據(jù))相比，ERA5數(shù)據(jù)集的空間分辨率(0.25°×0.25°)更高。

由于地震震前熱異常往往發(fā)生在震中和斷裂帶附近的大尺度區(qū)域內(nèi)，因此，我們采用Dobrovolsky等(1979)提出的地震影響區(qū)范圍估計(jì)公式(式(1))劃定研究區(qū)域的范圍：

r=100.43M

(1)

其中，r為地震影響區(qū)域的半徑(單位為km)，M為地震震級(jí)。根據(jù)計(jì)算結(jié)果，2014年于田地震研究區(qū)域的范圍為30.91°～40.91°N，77.59°～87.59°E，2008年汶川地震的研究區(qū)范圍為26°～36°N，98.3°～108.3°E。 2次地震的研究時(shí)段為地震前3個(gè)月—地震當(dāng)月，共獲取13a的觀測(cè)數(shù)據(jù)。除地震年外，其余年份內(nèi)均無震級(jí)>5.0級(jí)的地震發(fā)生。熱參數(shù)數(shù)據(jù)為夜間平均值(北京時(shí)間00：00—03：00，UTC： 16：00—19：00)，以避免太陽輻射和人為活動(dòng)的影響。

1.2 地震熱參數(shù)異常的提取方法

ZS法和RST法是目前地震熱參數(shù)異常時(shí)空變化特征研究中較為常用的2種地震熱參數(shù)異常提取方法。ZS法在圖像差分的基礎(chǔ)上引入了信噪比的概念，將地震活動(dòng)區(qū)當(dāng)前觀測(cè)值與多年觀測(cè)數(shù)據(jù)平均值之間的差值作為信號(hào)部分，并將多年觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差作為噪聲部分。該方法可表示為(Ouzounovetal.，2011)

(2)

式中，v(x，y，t)是地表溫度的當(dāng)前觀測(cè)值，μ(x，y)為背景場(chǎng)，即同一位置多年同期觀測(cè)數(shù)據(jù)的平均值，δ(x，y)為相應(yīng)的標(biāo)準(zhǔn)差。ZS指數(shù)的大小不僅反映了當(dāng)前觀測(cè)值偏離多年背景值的程度，同時(shí)還反映了這種偏離程度相對(duì)于多年噪聲水平的強(qiáng)弱。

RST法利用在類似觀測(cè)條件(如同月、 同時(shí)、 同傳感器)下獲得的多年觀測(cè)數(shù)據(jù)集的預(yù)期值(即多年平均值)及可變性(多年標(biāo)準(zhǔn)差)定義被測(cè)信號(hào)的異常變化。RST法將鄰域差分與ZS法中的背景場(chǎng)歸一化相結(jié)合，明確定義了統(tǒng)計(jì)學(xué)中地震前異常的數(shù)學(xué)表達(dá)式，熱紅外異常變化通過特定的指數(shù)： REIRA(Robust Estimator of TIR Anomalies)定義(Genzanoetal.，2015； Ouzounovetal.，2018； Khalilietal.，2019)。RST法可表示為

(3)

本文同時(shí)采用ZS法和RST法對(duì)2014年于田地震和2008年汶川地震震前的地表溫度和長波輻射異常時(shí)空變化特征進(jìn)行提取分析。具體步驟包括(圖 1)： 1)數(shù)據(jù)下載和預(yù)處理(包括數(shù)據(jù)格式轉(zhuǎn)換、 重投影等)； 2)對(duì)用于構(gòu)建背景場(chǎng)的數(shù)據(jù)進(jìn)行滑動(dòng)平均以剔除偶發(fā)因素對(duì)背景場(chǎng)的影響； 3)針對(duì)研究區(qū)的地物類型，在RST法中，大致將于田地震和汶川地震研究區(qū)分為山區(qū)和盆地2種地物類型，并分別計(jì)算這2種地物類型各自的平均值； 4)按照ZS法和RST法的算法原理分別計(jì)算多年數(shù)據(jù)的平均值和標(biāo)準(zhǔn)差作為背景場(chǎng)； 5)按照式(2)和式(3)計(jì)算得到ZS指數(shù)和REIRA指數(shù)； 6)為完整展示研究時(shí)段內(nèi)熱參數(shù)異常變化的過程，我們參照MODIS等數(shù)據(jù)的形式(MODIS產(chǎn)品提供8d平均值數(shù)據(jù))，以8d為一個(gè)階段繪制三維ZS指數(shù)和REIRA指數(shù)的平均值圖像。

圖 1 地震異常提取的流程圖Fig. 1 Flow chart of seismic anomalies extraction.

在ZS法和RST法的實(shí)際應(yīng)用中，通常假設(shè)數(shù)據(jù)集具有近似高斯分布的概率分布。在這個(gè)假設(shè)下，68.3%、 95.4%和99.7%的值分布在1、 2、 3個(gè)標(biāo)準(zhǔn)差的范圍內(nèi)(Jiaoetal.，2018)。為在統(tǒng)一異常閾值下得到可對(duì)比的異常提取結(jié)果并討論不同異常強(qiáng)度閾值對(duì)異常提取結(jié)果的影響，我們分別將異常強(qiáng)度閾值設(shè)定為2和3，此時(shí)的置信區(qū)間分別為95%和99%。

2 于田地震的概況及區(qū)域構(gòu)造背景

本文以2014年于田地震為典型震例，因此本節(jié)將著重介紹2014年于田地震的概況及區(qū)域構(gòu)造背景。據(jù)中國地震臺(tái)網(wǎng)測(cè)定，北京時(shí)間2014年2月12日17時(shí)19分49秒，中國新疆維吾爾自治區(qū)和田地區(qū)于田縣發(fā)生MW6.9 地震，震中位于無人區(qū)內(nèi)，坐標(biāo)為(36.1°N，82.5°E)，震源深度12km，震源機(jī)制為左旋走滑。截至北京時(shí)間2014年2月13日14：00：00，共記錄到余震1097個(gè)，其中MS5.0 ～5.9地震1個(gè)，MS4.0 ～4.9地震12個(gè)，MS3.0 ～3.9地震22個(gè)。

2014年于田地震發(fā)生在阿爾金斷裂帶的西南段(圖 2)。阿爾金斷裂帶位于青藏高原西北邊緣，是亞洲大陸內(nèi)最長的走滑型活動(dòng)斷裂帶，呈NEE走向，長約2000km。阿爾金斷裂帶分割了青藏高原和塔里木盆地，并劃分了青藏高原北部不同的構(gòu)造單元，對(duì)吸收印度-歐亞大陸碰撞而產(chǎn)生的地殼縮短具有重要作用。根據(jù)斷裂帶屬性，阿爾金斷裂帶可分為3個(gè)部分： NNE—E向逆沖的喀拉喀什斷裂、 E—NEE走向的阿什庫勒斷裂和NE—NEE走向的硝爾庫勒斷裂(Yinetal.，2002； 徐錫偉等，2011； Lietal.，2015； 宋春燕等，2015； Lietal.，2016)。該地區(qū)歷史上共發(fā)生過2次震級(jí)>7級(jí)的地震，分別為1924年阿爾金斷裂帶西段7.2級(jí)地震，震中位于2014年于田地震震中東北約250km處； 2008年于田MW7.2 地震，震中位于2014年于田地震震中西南150km處。

圖 2 研究區(qū)的區(qū)域構(gòu)造背景圖Fig. 2 Tectonic map of the study area.

圖 3 2014于田地震震前3個(gè)月的地表溫度和長波輻射ZS 8d平均值的三維時(shí)空特征分布圖Fig. 3 Surface temperature(ST) and outgoing longwave radiation(OLR) three-dimensional maps based on 8-day averaged Z-score index values for 3 months before the 2014 Yutian earthquake.圖m為包含發(fā)震時(shí)段的結(jié)果

圖 4 2014年于田地震震前3個(gè)月的地表溫度和長波輻射REIRA 8d平均值的三維時(shí)空特征分布圖Fig. 4 Surface temperature(ST) and outgoing longwave radiation(OLR) three-dimensional maps based on 8-day averaged REIRA index values for 3 months before the 2014 Yutian earthquake.圖m為包含發(fā)震時(shí)段的結(jié)果

圖 5 地表溫度(a)和長波輻射(b)敏感性測(cè)試數(shù)據(jù)時(shí)間序列圖Fig. 5 Time series plot of sensitivity test data for ST(a) and OLR(b).

3 2014年于田地震震前熱異常的時(shí)空變化特征

圖 3 展示了研究時(shí)段內(nèi)地表溫度(ST)和晴空向外長波輻射(OLR)ZS指數(shù)8d平均值三維圖像，其中紅色部分表示ZS指數(shù)>2的區(qū)域，紅色矩形框代表地震的發(fā)震時(shí)段。從圖 3 中可以看出，在震前3個(gè)月，ST和OLR都出現(xiàn)多次高值異常。從時(shí)間序列上來看，隨著地震的臨近，ST和OLR異常變化的頻次逐漸增加； 高值異常最早出現(xiàn)在地震前2個(gè)月(圖3f)，而地震前1個(gè)月出現(xiàn)持續(xù)的高值異常(圖3j，3k，3l)。從空間分布上來看，隨著地震的臨近，ST和OLR的高值異常區(qū)域由研究區(qū)的北部和西南部(圖3f，3j)逐漸向地震震中和阿爾金斷裂帶周圍集中(圖3k，3l)，在包含發(fā)震時(shí)刻的時(shí)段內(nèi)，ST無明顯高值異常，而OLR的高值異常區(qū)則明顯分布于震中周圍。

圖 4 展示了研究時(shí)段內(nèi)ST和OLR的REIRA指數(shù)8d平均值的三維圖像，其中紅色部分表示REIRA指數(shù)>2的區(qū)域，紅色矩形框代表地震的發(fā)震時(shí)段。對(duì)比ZS法和RST法的提取結(jié)果可發(fā)現(xiàn)： 1)對(duì)比圖3f、 4f，3k、 4k，3l、 4l可知，RST法與ZS法的提取結(jié)果在地理空間分布上具有相似性，基本分布于研究區(qū)的同一位置； 2)在地震震前3個(gè)月，2種方法的提取結(jié)果差異較大(圖3a—e，圖4a—e)，ZS法得到的結(jié)果無明顯異常變化，而RST法得到的結(jié)果則有異常變化且分散分布于整個(gè)研究區(qū)內(nèi)； 此外，2種方法的提取結(jié)果在異常出現(xiàn)的頻率和幅度上明顯不同，RST法得到的異常頻次和幅度都高于ZS法的提取結(jié)果。為何使用同樣的數(shù)據(jù)會(huì)得到不盡相同的時(shí)空分布特征？哪種方法得到的時(shí)空分布特征更為可靠？為了回答以上問題，我們?cè)谙鹿?jié)中針對(duì)2種方法的異常變化敏感性、 對(duì)背景信息抑制能力和對(duì)地震信息指示性3個(gè)方面進(jìn)行探究。

4 地震熱參數(shù)異常提取方法的定量化對(duì)比

4.1 對(duì)異常變化敏感性的對(duì)比

熱參數(shù)時(shí)空特征的異常變化受到地形、 氣候變化、 人為活動(dòng)等多種因素的影響。相較于其他影響因素，地震孕震過程對(duì)熱參數(shù)異常變化的影響十分微弱。地震熱參數(shù)異常提取方法的敏感性代表了其識(shí)別地震孕震過程中引起的熱參數(shù)微弱異常變化的能力。本文選取2014年于田地震前(2009年)和地震后(2015年)2個(gè)無震級(jí)>5級(jí)地震發(fā)生年份的震中所在像元的時(shí)間序列數(shù)據(jù)作為測(cè)試數(shù)據(jù)(圖5a 為ST數(shù)據(jù)，圖5b 為OLR數(shù)據(jù))，在時(shí)間序列中隨機(jī)選取20個(gè)數(shù)據(jù)值，人為設(shè)定1～15(單位為K或W/m2)的異常變化，隨后利用2種方法提取異常，并統(tǒng)計(jì)在不同異常閾值和異常變化幅度下檢測(cè)出的設(shè)定異常個(gè)數(shù)及其準(zhǔn)確率(即檢測(cè)出的設(shè)定異常個(gè)數(shù)與總異常個(gè)數(shù)之比)。

圖 6 為利用2009年地表溫度(a、 c)和長波輻射(b、 d)數(shù)據(jù)通過ZS法和RST法提取得到設(shè)定異常提取天數(shù)(a、 b)和準(zhǔn)確率(c、 d)。圖 7 為利用2015年地表溫度(a、 c)和長波輻射(b、 d)數(shù)據(jù)通過ZS法和RST法提取得到的設(shè)定異常提取天數(shù)(a、 b)和準(zhǔn)確率(c、 d)。異常閾值的設(shè)定對(duì)異常變化檢測(cè)具有較大影響，當(dāng)設(shè)定的異常閾值較高時(shí)，雖然不能檢測(cè)出全部人為設(shè)定的異常變化，但異常檢測(cè)的準(zhǔn)確率較高。對(duì)比2種方法得到的結(jié)果可知，ZS法和RST法都對(duì)小幅度的異常變化具有一定的敏感性(異常變化幅度為1～5K(W/m2)時(shí))，但RST法的異常檢測(cè)準(zhǔn)確率優(yōu)于ZS法。當(dāng)異常變化幅度較大時(shí)(12～15K(W/m2))，利用2009年的數(shù)據(jù)，RST法的提取效果更好； 利用2015年的數(shù)據(jù)，當(dāng)異常閾值為2.5時(shí)，ZS法的提取效果也較好，這可能是由于提取參數(shù)不同所致。綜合上述研究成果，我們認(rèn)為ZS法和RST法都對(duì)于小幅度的異常變化具有一定的敏感性，但RST法的提取效果優(yōu)于ZS法。

圖 6 利用2009年地表溫度(a、 c)和長波輻射(b、 d)數(shù)據(jù)通過ZS法、 RST法提取得到設(shè)定異常提取天數(shù)(a、 b)和準(zhǔn)確率(c、 d)Fig. 6 Artificial anomaly extraction days(a and b)and accuracy rate(c and d)obtained from 2009 surface temperature(a and c)and outgoing longwave radiation(b and d)data extracted by ZS and RST methods.

圖 7 利用2015年地表溫度(a、 c)和長波輻射(b、 d)數(shù)據(jù)通過ZS法、 RST法提取得到設(shè)定異常提取天數(shù)(a、 b)和準(zhǔn)確率(c、 d)Fig. 7 Artificial anomaly extraction days(a and b)and accuracy rate(c and d)obtained from 2015 surface temperature(a and c)and outgoing longwave radiation(b and d)data extracted by ZS and RST methods.

4.2 對(duì)背景信息抑制能力的對(duì)比

地震熱參數(shù)異常提取方法對(duì)背景信息的抑制能力決定了提取結(jié)果的準(zhǔn)確性。本文提取了2個(gè)無震年(2009年和2015年)的地表溫度和長波輻射異常時(shí)空變化特征。無震年的數(shù)據(jù)沒有受到地震活動(dòng)的影響，因此理想的異常提取結(jié)果應(yīng)無明顯異常變化。然而，由于異常提取方法無法完全剔除其他因素造成的影響，故提取結(jié)果中可能有少量由其他因素引起的熱異常。本研究以無震年的異常像元密度，即異常像元個(gè)數(shù)與總像元個(gè)數(shù)的比值定量描述異常提取方法對(duì)背景信息的抑制能力。我們統(tǒng)計(jì)了利用ZS法提取無震年異常像元密度大于RST法提取的天數(shù)和利用RST法提取無震年異常像元密度大于ZS法提取的天數(shù)，分別記為ρZR和ρRZ。無震年異常像元密度越低，則該方法對(duì)其他因素引起的熱異常的抑制能力越強(qiáng)，相應(yīng)地，其在地震年提取出的熱異常時(shí)空變化特征可能越可靠。

圖 8 為2009年(a、 c)和2015年(b、 d)地表溫度(a、 b)和長波輻射(c、 d)異常像元密度在不同異常標(biāo)準(zhǔn)下的統(tǒng)計(jì)結(jié)果。無論異常標(biāo)準(zhǔn)設(shè)定為多少，ρZR的值遠(yuǎn)小于ρRZ，說明ZS法在無震年得到的異常像元明顯少于RST法，基于2009年和2015年數(shù)據(jù)得到的結(jié)果呈現(xiàn)類似的特征。由此可知，ZS法相較于RST法對(duì)其他因素引起的熱參數(shù)異常變化的抑制作用更強(qiáng)。

圖 8 2009年(a、 c)和2015年(b、 d)地表溫度(a、 b)和長波輻射(c、 d)異常像元密度在不同異常標(biāo)準(zhǔn)下的統(tǒng)計(jì)結(jié)果Fig. 8 Statistical results of surface temperature(a and b)and outgoing longwave radiation(c and d)anomaly pixel densities under different anomaly standards in 2009(a and c)and 2015(b and d).

4.3 對(duì)地震信息指示性的對(duì)比

前人研究表明，地震震前熱異常的空間分布與震中位置具有一定的相關(guān)性(Jiaoetal.，2018)。因此，地震震前熱異常的空間分布是預(yù)測(cè)地震信息(發(fā)震時(shí)段和震中位置)的重要指標(biāo)。本文在第3節(jié)中說明了2種方法提取的異常結(jié)果存在時(shí)空分布差異。那么，哪種方法提取得到的熱異常空間分布能夠更好地指示地震信息？為進(jìn)一步探究該問題，我們根據(jù)地理學(xué)第一定律，即“所有事物都與其他事物相關(guān)，但近處的事物比遠(yuǎn)處的事物更相關(guān)”(Tobler，1970)，通過歸一化距離指數(shù)(Qietal.，2020)定量評(píng)估熱參數(shù)異常時(shí)空分布對(duì)地震信息(發(fā)震時(shí)段和震中位置)的指示性作用(式4)：

(4)

式中，(xi，yi)為異常像元坐標(biāo)，(xepi，yepi)為震中所在的像元坐標(biāo)，n為整幅圖像異常像元總數(shù)。NWI值越大，說明得到的空間分布越集中于發(fā)震震中，指示性越好。

在對(duì)地震信息指示性進(jìn)行對(duì)比時(shí)，需要同時(shí)關(guān)注2個(gè)方面： 1)在都具有異常像元分布時(shí)，哪種方法計(jì)算得到的NWI值更大； 2)哪種方法計(jì)算得到NWI最大值的時(shí)間離發(fā)震時(shí)段更近。因此，我們首先統(tǒng)計(jì)了當(dāng)ZS法和RST法都具有異常像元分布時(shí)，其中一種方法計(jì)算得到的NWI值大于另一種方法的天數(shù)，然后對(duì)地震震前2種方法計(jì)算出的NWI最大值距地震發(fā)震時(shí)段的天數(shù)進(jìn)行了統(tǒng)計(jì)。圖 9 為基于地表溫度(圖9a)和長波輻射(圖9b)數(shù)據(jù)的ZS法和RST法對(duì)地震震中指示性的比較結(jié)果圖。從圖中可以看出，ZS法和RST法基于不同數(shù)據(jù)得到的結(jié)果存在不同。對(duì)于ST數(shù)據(jù)，當(dāng)同時(shí)具有異常像元分布時(shí)，ZS法計(jì)算得到的NWI值略高于RST法計(jì)算得到的NWI值的天數(shù)更多，特別是當(dāng)異常標(biāo)準(zhǔn)設(shè)定為3時(shí)。然而，對(duì)于OLR數(shù)據(jù)，RST法計(jì)算得到的NWI值高于ZS法計(jì)算出的NWI值的天數(shù)更多。這一現(xiàn)象說明對(duì)于ST數(shù)據(jù)，ZS法對(duì)于地震震中分布的指示性略好于RST法； 而對(duì)于OLR數(shù)據(jù)，RST法得到的結(jié)果對(duì)地震震中空間分布的指示性更好。

表1為在震前時(shí)段內(nèi)使用ZS法、 RST法計(jì)算得到的NWI最大值出現(xiàn)時(shí)間距地震發(fā)生時(shí)的天數(shù)。從表中可以看出，相較于RST法，ZS法得到的NWI最大值出現(xiàn)時(shí)間更靠近地震的發(fā)震時(shí)段，對(duì)地震發(fā)震時(shí)段的指示性更強(qiáng)。綜合上述研究成果表明： 1)針對(duì)ST數(shù)據(jù)，ZS法對(duì)于地震震中位置的指示性略好于RST法，對(duì)于OLR數(shù)據(jù)則是RST法更優(yōu)； 2)ZS法得到的NWI最大值出現(xiàn)時(shí)段距離發(fā)震時(shí)刻更近。

圖 9 基于地表溫度(a)和長波輻射(b)數(shù)據(jù)的ZS法和RST法對(duì)地震震中指示性的比較結(jié)果圖Fig. 9 Comparison of epicentre indication ability of surface temperature(a)and long wave radiation(b)obtained by ZS and RST methods.

表 1 ZS法、 RST法計(jì)算得到最大NWI值時(shí)段距離地震發(fā)震的天數(shù)(數(shù)字表示震前天數(shù))Table1 The days from the maximum NWI value period to the occurrence of earthquake calculated by ZS and RST methods (The number represents the days before the earthquake)

表 2 ZS法、 RST法針對(duì)驗(yàn)證震例計(jì)算得到的最大NWI值時(shí)段距離地震發(fā)震天數(shù)(數(shù)字表示震前天數(shù))Table2 The days between the maximum NWI value period and the earthquake occurrence calculated by ZS and RST methods from the verification earthquake case(The number represents the days before the earthquake)

4.4 驗(yàn)證震例對(duì)比結(jié)果

據(jù)中國地震臺(tái)網(wǎng)測(cè)定，2008年5月12日14點(diǎn)28分(北京時(shí)間)四川省阿壩藏族自治州汶川縣發(fā)生MW7.9 地震，震中坐標(biāo)為(31.0°N，103.4°E)，震源深度14km，是21世紀(jì)以來發(fā)生在中國大陸境內(nèi)破壞性最強(qiáng)的地震。我們以此次地震作為驗(yàn)證地震，利用汶川地震震前平靜年的數(shù)據(jù)(2007年)作為測(cè)試數(shù)據(jù)，采用與于田地震相同的方法，對(duì)ZS法和RST法在異常變化敏感性、 對(duì)背景信息抑制能力和對(duì)地震信息的指示性3個(gè)方面進(jìn)行對(duì)比，驗(yàn)證以于田地震作為典型震例的對(duì)比結(jié)果的可靠性。ZS法和RST法的對(duì)比策略如4.1、 4.2和4.3節(jié)所述。

圖 9 為驗(yàn)證震例異常變化敏感性(a、 b)、 對(duì)背景信息抑制能力(c、 d)和對(duì)地震信息指示性(e、 f)的測(cè)試結(jié)果圖。表2 為NWI最大值出現(xiàn)時(shí)段距地震發(fā)震的天數(shù)統(tǒng)計(jì)。利用汶川地震數(shù)據(jù)得到的異常變化敏感性結(jié)果與于田地震略有不同，當(dāng)異常閾值為2時(shí)ZS法對(duì)于微弱異常變化敏感性最好。2種方法對(duì)背景信息抑制能力的測(cè)試結(jié)果與于田地震類似，ZS法對(duì)于背景信息的抑制能力優(yōu)于RST法。在對(duì)地震信息指示性的測(cè)試中，總體而言，RST法提取得到的異常空間分布對(duì)地震震中的指示性更好，而ZS法對(duì)于地震發(fā)震時(shí)段的指示性更好，NWI最大值出現(xiàn)的時(shí)間更靠近地震發(fā)震時(shí)段。

圖 10 針對(duì)驗(yàn)證震例的異常變化敏感性(a、 b)、 對(duì)背景信息抑制能力(c、 d)和對(duì)地震信息指示性(e、 f)的測(cè)試結(jié)果圖Fig. 10 The test results of sensitivity to anomaly change(a and b)，ability to suppress background information(c and d)and indicate seismic information(e and f)of the verification earthquake case.

5 討論和結(jié)論

本研究以2008年于田地震為典型震例，以2008年汶川地震為驗(yàn)證震例，對(duì)ZS法和RST法在對(duì)異常變化的敏感性、 對(duì)背景信息的抑制能力和對(duì)地震信息的指示性3個(gè)方面的提取效果進(jìn)行了定性和定量化比較。除異常變化敏感性外，在對(duì)背景信息的抑制能力和對(duì)地震信息的指示性上，2個(gè)震例得到的對(duì)比結(jié)果類似； ZS法對(duì)背景信息的抑制能力更好； RST法對(duì)地震震中的指示性更好； ZS法對(duì)地震發(fā)震時(shí)段的指示性更好。

圖 11 ZS法和RST法信號(hào)部分和噪聲部分在整個(gè)研究時(shí)段內(nèi)每天最大值與最小值的差異時(shí)間序列曲線Fig. 11 The time series curve of difference values between the maximum and minimum values of the signal part and the noise part of the ZS method and the RST method for each day during the whole study period.a、 b 信號(hào)部分的時(shí)間序列曲線； c、 d 噪聲部分的時(shí)間序列曲線； a、 c 利用ST數(shù)據(jù)得到的結(jié)果； b、 d 利用OLR數(shù)據(jù)得到的結(jié)果

為何 2 種方法在這3個(gè)方面具有顯著的差異？如何說明在第3節(jié)中2種方法得到的異常變化時(shí)空特征的不同？這些問題仍需要進(jìn)一步分析和討論。ZS法和RST法都是以指數(shù)形式表示的異常提取方法。在計(jì)算過程中，2種方法的不同之處在于RST法針對(duì)區(qū)域進(jìn)行了地物劃分，求取了不同地物類型的平均值。在構(gòu)建背景場(chǎng)時(shí)，首先計(jì)算多年數(shù)據(jù)中每個(gè)像元與其所在的地物類型平均值的差值，之后再計(jì)算多年平均值和標(biāo)準(zhǔn)差以構(gòu)建背景場(chǎng)。ZS法則直接計(jì)算多年平均值和標(biāo)準(zhǔn)差作為背景場(chǎng)。為確定求取每個(gè)像元與其所在地物類型平均值的差值對(duì)最終信噪比結(jié)果的影響，我們統(tǒng)計(jì)了ZS法和RST法信號(hào)部分(即地震年與多年平均值的差值，如圖11a 和11b所示)和噪聲部分(即多年標(biāo)準(zhǔn)差，如圖11c 和11d所示)在整個(gè)研究時(shí)段內(nèi)每天的最大值與最小值的差異。從圖 11 中可以看出，求取每個(gè)像元與其所在地物類型平均值的差值對(duì)信號(hào)部分影響不大，對(duì)噪聲部分的影響較大。整體上，ZS法噪聲部分的最大值、 最小值的差異大于RST法。因此，在異常變化幅度較小時(shí)，由于RST法計(jì)算得到的信噪比大于ZS法，使得前者對(duì)于異常變化更為敏感，在相同的異常標(biāo)準(zhǔn)下能夠提取出更多的異常變化，且在地震平靜年也能提取出大量異常。當(dāng)ZS法提取結(jié)果無明顯變化時(shí)，RST法也能提取出分布于研究區(qū)內(nèi)的異常； 當(dāng)2種方法提取結(jié)果同時(shí)存在異常變化時(shí)，RST法提取結(jié)果的空間分布面積更大。ZS法和RST法的異常變化敏感性在不同震例中得到了不同的比較結(jié)果，是由研究區(qū)地物分布類型不同所造成的。對(duì)于田地震而言，區(qū)域內(nèi)兩大地物類型——山區(qū)和沙漠的分布較為均一，不同地物類型的平均值差異較大。而汶川地震研究區(qū)內(nèi)的地物類型則更為復(fù)雜，不同地物類型的平均值差異較小，影響了RST法對(duì)于微弱變化的敏感性。由于熱參數(shù)影像的空間分辨率較大，一個(gè)像元內(nèi)包含了大量地物，無法利用光學(xué)影像進(jìn)行精確的地物類型分類，這也降低了RST法提取結(jié)果的可靠性。此外，異常標(biāo)準(zhǔn)的設(shè)定對(duì)于提取結(jié)果的影響也很大，但目前并沒有有效定義異常標(biāo)準(zhǔn)的方法，大多數(shù)研究都依據(jù)經(jīng)驗(yàn)定義異常標(biāo)準(zhǔn)，因此需要進(jìn)一步研究定義異常標(biāo)準(zhǔn)的有效方法，以提高熱異常提取的可靠性。

綜上所述，我們認(rèn)為ZS法和RST法對(duì)于異常變化都具有一定的提取能力，但相比較而言，由于RST法對(duì)于背景信息抑制能力弱，且在對(duì)地物類型進(jìn)行分類過程中存在較大的人為因素影響，因此，ZS法是較為合適且簡便的地震異常變化提取方法。