兩種空間插值方法在燈塔5.1級(jí)地震烈度分布中的對(duì)比研究

姜金征，梁永朵，劉琳婷，戴盈磊，于 浩

(遼寧省地震局，遼寧沈陽(yáng) 110034)

0 引言

基于強(qiáng)震動(dòng)臺(tái)網(wǎng)的烈度速報(bào)對(duì)于推動(dòng)我國(guó)防震減災(zāi)事業(yè)的發(fā)展，促進(jìn)地震監(jiān)測(cè)臺(tái)網(wǎng)資源更好地服務(wù)我國(guó)經(jīng)濟(jì)社會(huì)進(jìn)步將發(fā)揮重要作用［1－2］。遼寧省自2007年開(kāi)展強(qiáng)震動(dòng)觀測(cè)以來(lái)，至目前已擁有了82個(gè)實(shí)時(shí)傳輸?shù)臄?shù)據(jù)強(qiáng)震動(dòng)臺(tái)站，但利用現(xiàn)有臺(tái)站進(jìn)行烈度速報(bào)，不但臺(tái)站數(shù)量不足，且分布并不均勻，需要進(jìn)行空間插值才能滿(mǎn)足烈度速報(bào)的需要［3－4］。因此，研究適當(dāng)?shù)牟逯捣椒ǎ瓿闪叶鹊戎稻€(xiàn)圖的繪制成為關(guān)鍵。

所謂插值是通過(guò)對(duì)離散數(shù)據(jù)場(chǎng)插值后建立數(shù)學(xué)模型，進(jìn)而提取曲線(xiàn)或曲面等幾何信息，獲取數(shù)據(jù)場(chǎng)內(nèi)部信息并加以顯示［5－6］?？臻g插值常用于將離散點(diǎn)的測(cè)量數(shù)據(jù)轉(zhuǎn)換為連續(xù)的數(shù)據(jù)曲面，以便與其它空間現(xiàn)象的分布模式進(jìn)行比較。一般常用的插值方法有距離反比加權(quán)插值法、克里格插值法、最小曲率、改進(jìn)謝別德法、自然鄰點(diǎn)插值法、最近鄰點(diǎn)插值法、多元回歸法、徑向基函數(shù)法、線(xiàn)性插值三角網(wǎng)法、移動(dòng)平均法、局部多項(xiàng)式法等。本文選用距離反比加權(quán)插值法和克里格插值方法，在應(yīng)用燈塔5.1 級(jí)地震基礎(chǔ)上，對(duì)兩種方法計(jì)算結(jié)果進(jìn)行了比較研究。

1 插值方法及驗(yàn)證

1.1 距離反比加權(quán)插值法(FIDW，Inverse Distance Weighted Interpolation)

距離反比加權(quán)插值法是基于相近相似的原理，它以插值點(diǎn)與樣本點(diǎn)間的距離為權(quán)重進(jìn)行平均，離插值點(diǎn)越近的樣本點(diǎn)賦予的權(quán)重越大。距離反比加權(quán)插值法，可以進(jìn)行確切的或者圓滑的方式插值［7］。方次參數(shù)控制著權(quán)系數(shù)如何隨著離開(kāi)一個(gè)格網(wǎng)結(jié)點(diǎn)距離的增加而下降。對(duì)于一個(gè)較大的方次，較近的數(shù)據(jù)點(diǎn)被給定一個(gè)較高的權(quán)重份額，對(duì)于一個(gè)較小的方次，權(quán)重比較均勻地分配給各數(shù)據(jù)點(diǎn)。計(jì)算一個(gè)格網(wǎng)結(jié)點(diǎn)時(shí)給予一個(gè)特定數(shù)據(jù)點(diǎn)的權(quán)值與指定方次的從結(jié)點(diǎn)到觀測(cè)點(diǎn)的該結(jié)點(diǎn)被賦予距離倒數(shù)成比例。當(dāng)計(jì)算一個(gè)格網(wǎng)結(jié)點(diǎn)時(shí)，配給的權(quán)重是一個(gè)分?jǐn)?shù)，所有權(quán)重的總和等于1.0。當(dāng)一個(gè)觀測(cè)點(diǎn)與一個(gè)格網(wǎng)結(jié)點(diǎn)重合時(shí)，該觀測(cè)點(diǎn)被給予一個(gè)實(shí)際為1.0 的權(quán)重，所有其它觀測(cè)點(diǎn)被給予一個(gè)幾乎為0.0 的權(quán)重。換言之，該結(jié)點(diǎn)被賦給與觀測(cè)點(diǎn)一致的值，這就是一個(gè)準(zhǔn)確插值(圖1)。

距離反比加權(quán)平均插值是根據(jù)估算的網(wǎng)格中心點(diǎn)附近臺(tái)站實(shí)際烈度值插值網(wǎng)格中心點(diǎn)烈度值，計(jì)算公式如下:

式中:Zp是相鄰點(diǎn)的高程，d是插值點(diǎn)到p點(diǎn)的距離;n是參數(shù)，范圍從1.0 至6.0，通常用的值是2.0。－n表示越靠近被插值點(diǎn)越重要。

圖1 FIDW插值法網(wǎng)格示意圖Fig.1 Sketch map of FIDW interpolation method

1.2 克里格法

(1)克里格法原理

克里格法(Kriging)，也稱(chēng)空間局部估計(jì)或空間局部插值，是建立在變異函數(shù)理論及結(jié)構(gòu)分析基礎(chǔ)之上，在有限區(qū)域內(nèi)對(duì)區(qū)域化變量的取值進(jìn)行無(wú)偏最優(yōu)估計(jì)的一種方法。該方法將任一個(gè)點(diǎn)的估計(jì)值通過(guò)該點(diǎn)影響范圍內(nèi)的n個(gè)有效樣本值Z(xi)的線(xiàn)性組合得到，即

式中，λi是與樣點(diǎn)觀察Z(xi)有關(guān)的加權(quán)系數(shù)，用來(lái)表示各個(gè)樣點(diǎn)值Z(xi)對(duì)估計(jì)值的貢獻(xiàn)。對(duì)于任一給定的區(qū)域和數(shù)據(jù)信息Z(xi)，i=1，2，…，n，存在一組加權(quán)系數(shù) λi，λi由公式［8］確定。如果使得估計(jì)方差為最小，其區(qū)域內(nèi)的真值就能在最小的可能置信區(qū)間內(nèi)產(chǎn)生。

(2)克里格法變異分析

變異函數(shù)通過(guò)其自身的結(jié)構(gòu)及其各項(xiàng)參數(shù)從不同的角度反映空間變異性。當(dāng)空間點(diǎn)x在x軸上變化，區(qū)域變化量xi和xi + h在空間位置的實(shí)測(cè)值是Z(xi)和Z(xi)+ h，其中i =1，2，3，...，N(h)，變異函數(shù)的離散公式為:

h為兩樣本點(diǎn)空間分隔距離。

(3)克里格法變異函數(shù)理論模型

變異函數(shù)理論模型，是當(dāng)計(jì)算得到變異函數(shù)值之后，選取合適的模型進(jìn)行參數(shù)的擬合，即所謂的“結(jié)構(gòu)分析”。變異函數(shù)模型常用的有球狀模型、指數(shù)模型、高斯模型等。

球狀模型公式為:

指數(shù)模型公式為:

高斯模型公式為:

球狀模型、指數(shù)模型和高斯模型對(duì)應(yīng)的圖形如圖2所示。

(4)克里格法估算網(wǎng)格點(diǎn)烈度值

利用克里格法估算網(wǎng)格點(diǎn)烈度值主要通過(guò)以下步驟:

①先對(duì)需要插值的區(qū)域網(wǎng)格化，選定一個(gè)待估點(diǎn)，確定地理經(jīng)緯度位置;②其次確定待估點(diǎn)最近距離3～4個(gè)已知烈度的臺(tái)站;③再次根據(jù)變異函數(shù)的參數(shù)及各向異性情況，選用合適模型構(gòu)造方程組;④求解方程組，得到系數(shù)λi;⑤根據(jù)公式1，求出估計(jì)點(diǎn)烈度值;⑥重復(fù)1～5 步，求得所有網(wǎng)格點(diǎn)的烈度值。

圖2 三種模型示意圖Fig.2 Sketch maps of the three models

1.3 空間插值方法精度評(píng)價(jià)——交叉驗(yàn)證法

交叉驗(yàn)證法是一種常見(jiàn)的精度驗(yàn)證方法。常用的交叉驗(yàn)證統(tǒng)計(jì)指標(biāo)主要有平均誤差(ME)，平均絕對(duì)誤差(MAE)，平均相對(duì)誤差(MRE)、均方根誤差(RMSE)。ME、MAE、MRE和RMSE 的值越小，表明空間插值結(jié)果的精度越高。各指標(biāo)計(jì)算公式如下:

式中，z(xi)為預(yù)測(cè)值，本文中為兩種空間插值法的計(jì)算結(jié)果;z'(xi)為原始采樣值，本文中為宏觀烈度值。

2 兩種插值方法應(yīng)用分析

2.1 實(shí)際震例

2013年1月23日遼寧燈塔發(fā)生了5.1 級(jí)地震，遼寧省地震局強(qiáng)震動(dòng)臺(tái)網(wǎng)中心回收此次地震強(qiáng)震動(dòng)記錄，剔除了兩個(gè)記錄有問(wèn)題的臺(tái)站記錄，實(shí)際使用記錄20 組，合計(jì)60 條，取得強(qiáng)震記錄震中距范圍從16 至370km 不等，距離震中最近的遼陽(yáng)臺(tái)，強(qiáng)震記錄經(jīng)校正后，水平向加速度最大峰值為21gal。本文根據(jù)各臺(tái)站加速度峰值，利用美國(guó)ShakeMap 烈度計(jì)算方法［9］得到的烈度值如表1所示。

2.2 烈度圖繪制

本文利用以上2 種插值方法，應(yīng)用表1烈度計(jì)算結(jié)果，繪制了遼寧燈塔5.1 級(jí)地震儀器烈度分布圖。繪制烈度圖的具體作法如下:

(1)利用已有的臺(tái)站加速度值大小構(gòu)造一個(gè)粗略的、間隔統(tǒng)一的仿真臺(tái)站網(wǎng)格;

(2)根據(jù)臺(tái)站的加速度值，利用ShakeMap 方法確定臺(tái)站的烈度值;

(3)利用距離反比加權(quán)插值法及克里格插值法構(gòu)建仿真臺(tái)站的烈度網(wǎng)格;

(4)根據(jù)烈度網(wǎng)格點(diǎn)繪制整個(gè)區(qū)域內(nèi)的地震烈度分布等值線(xiàn)。

圖3為按上述方法繪制的烈度等值線(xiàn)圖，圖中實(shí)線(xiàn)烈度線(xiàn)為距離反比加權(quán)插值法繪制，虛線(xiàn)烈度線(xiàn)為克里格插值法繪制，5 度～3 度烈度等值線(xiàn)由震中向外圍呈環(huán)形展布。如圖3所示，兩種插值方法所繪烈度區(qū)域面積相當(dāng)，但距離反比加權(quán)插值法圖形接近于圓形，而克里格插值法接近于橢圓形。

2.3 與實(shí)際調(diào)查結(jié)果的比較

如圖4所示宏觀調(diào)查烈度圖，是根據(jù)遼寧省地震局震害宏觀調(diào)查結(jié)果以及咨詢(xún)強(qiáng)震觀測(cè)臺(tái)站當(dāng)值人員核實(shí)后繪制而成，烈度圖由6 度～3度四個(gè)橢圓形烈度等值線(xiàn)組成，長(zhǎng)軸方向與震中附近NE 向佟二堡斷裂［10］走向基本一致。從圖3與圖4比較而言，克里格法繪制的烈度圖更貼近于實(shí)際宏觀調(diào)查結(jié)果。

表1 強(qiáng)震動(dòng)記錄情況Tab.1 Strong motion records

距離反比加權(quán)插值法是相近相似的原理，即離插值點(diǎn)越近的樣本點(diǎn)賦予的權(quán)重越大，對(duì)鄰近采樣點(diǎn)具有依賴(lài)性，因此對(duì)臺(tái)站密度及分布均勻性要求較高。而克里格法強(qiáng)調(diào)變量的空間自相關(guān)性，考慮了各個(gè)鄰近樣點(diǎn)彼此之間的相互關(guān)系，所以對(duì)臺(tái)站密度及均勻性要求相對(duì)較低。

由圖3、圖4可見(jiàn)，圖中臺(tái)站分布不均勻且密度不足。因此，根據(jù)上述分析，克里格插值法相比距離反比加權(quán)插值法所繪烈度圖更接近于實(shí)際調(diào)結(jié)果。但鑒于克里格法計(jì)算相對(duì)復(fù)雜，如果臺(tái)網(wǎng)足夠密集且均勻，采用計(jì)算相對(duì)簡(jiǎn)單的距離反比加權(quán)插值法更利于烈度速報(bào)需求。

由圖4與圖3對(duì)比可知，宏觀調(diào)烈度區(qū)域范圍明顯大于插值方法所繪制烈度區(qū)域，這可能與插值方法所得烈度圖，臺(tái)站記錄均為基巖臺(tái)站記錄，且沒(méi)有經(jīng)過(guò)場(chǎng)地校正，所以烈度區(qū)域偏小也是可以理解的。

圖3 烈度等值線(xiàn)圖Fig.3 Intensity isogram

2.4 兩種插值方法的精度對(duì)比

表2為兩種插值方法的交叉驗(yàn)證結(jié)果，對(duì)比分析可知，兩種插值方法誤差差別不是太大，總體上看，克里格插值法的精度更高，但距離反比加權(quán)插值法的誤差也在可接受范圍內(nèi)。故兩種插值方法都可用于震后烈度估計(jì)。

表2 交叉驗(yàn)證誤差統(tǒng)計(jì)結(jié)果Tab.2 Error statistics of interpolation cross validation

圖4 宏觀烈度調(diào)查結(jié)果Fig.4 Survey results of macro intensity

3 結(jié)論

本文用兩種插值方法，即距離反比加權(quán)插值法和克里格插值法，繪制了燈塔5.1 級(jí)地震的烈度等值線(xiàn)圖，并與宏觀烈度調(diào)查結(jié)果進(jìn)行了比較。結(jié)果發(fā)現(xiàn):(1)當(dāng)臺(tái)站分布不均勻且密度相對(duì)不足時(shí)，采用克里格插值法繪制烈度圖更符合實(shí)際宏觀調(diào)查結(jié)果;(2)當(dāng)臺(tái)網(wǎng)密度較密集且分布相對(duì)均勻，本文推薦使用計(jì)算相對(duì)簡(jiǎn)單的距離反比加權(quán)插值法，此法更利于烈度速報(bào)。

［1］李山有，金星，陳先，等.地震動(dòng)強(qiáng)度與地震烈度速報(bào)研究［J］.地震工程與工程振動(dòng)，2002，22(6):1－7.

［2］張敏政.地震烈度及其評(píng)定［J］.防災(zāi)科技學(xué)院學(xué)報(bào)，2010(1):1－6.

［3］李永振，于沈平，金震，等.遼寧省強(qiáng)震臺(tái)網(wǎng)的地震應(yīng)急效果初探［J］，東北地震研究，2009，25(1):37－42.

［4］李永振，翟文杰，于沈平，等.遼寧省數(shù)字強(qiáng)震臺(tái)網(wǎng)建設(shè)與運(yùn)行［J］，東北地震研究，2009，25(3):18－22.

［5］周寶峰，溫瑞智，謝禮立.非因果濾波器在強(qiáng)震數(shù)據(jù)處理中的應(yīng)用［J］.地震工程與工程振動(dòng)，2012，32(2):25－34.

［6］胡進(jìn)軍，謝禮立.地震破裂的方向性效應(yīng)相關(guān)概念綜述［J］.地震工程與工程振動(dòng)，2011，31(4):1－8.

［7］陳晶晶，胡蓓蓓，王軍，等.天津降水?dāng)?shù)據(jù)的空間插值分析［J］，安徽師范大學(xué)學(xué)報(bào)，2010，33(4):15－21.

［8］靳國(guó)棟，劉衍聰，牛文杰，等.距離加權(quán)反比插值法與克里金插值法的比較［J］，長(zhǎng)春工業(yè)大學(xué)學(xué)報(bào)，2003，24(3):53－57.

［9］李俊，蘇楓，米宏亮，等.ShakeMap 及其在地震動(dòng)快速預(yù)估中的應(yīng)用［J］.中國(guó)地震，2010，26(1):103－111.

［10］單家增，張占文，孫紅軍，等.營(yíng)口－佟二堡斷裂帶成因機(jī)制的構(gòu)造物理模擬實(shí)驗(yàn)研究［J］.2004，石油勘探與開(kāi)發(fā)，2004，31(1):15－16.