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

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

0 引言

基于強震動臺網的烈度速報對于推動我國防震減災事業(yè)的發(fā)展，促進地震監(jiān)測臺網資源更好地服務我國經濟社會進步將發(fā)揮重要作用［1－2］。遼寧省自2007年開展強震動觀測以來，至目前已擁有了82個實時傳輸的數據強震動臺站，但利用現有臺站進行烈度速報，不但臺站數量不足，且分布并不均勻，需要進行空間插值才能滿足烈度速報的需要［3－4］。因此，研究適當的插值方法，完成烈度等值線圖的繪制成為關鍵。

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

1 插值方法及驗證

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

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

距離反比加權平均插值是根據估算的網格中心點附近臺站實際烈度值插值網格中心點烈度值，計算公式如下:

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

圖1 FIDW插值法網格示意圖Fig.1 Sketch map of FIDW interpolation method

1.2 克里格法

(1)克里格法原理

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

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

(2)克里格法變異分析

變異函數通過其自身的結構及其各項參數從不同的角度反映空間變異性。當空間點x在x軸上變化，區(qū)域變化量xi和xi + h在空間位置的實測值是Z(xi)和Z(xi)+ h，其中i =1，2，3，...，N(h)，變異函數的離散公式為:

h為兩樣本點空間分隔距離。

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

變異函數理論模型，是當計算得到變異函數值之后，選取合適的模型進行參數的擬合，即所謂的“結構分析”。變異函數模型常用的有球狀模型、指數模型、高斯模型等。

球狀模型公式為:

指數模型公式為:

高斯模型公式為:

球狀模型、指數模型和高斯模型對應的圖形如圖2所示。

(4)克里格法估算網格點烈度值

利用克里格法估算網格點烈度值主要通過以下步驟:

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

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

1.3 空間插值方法精度評價——交叉驗證法

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

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

2 兩種插值方法應用分析

2.1 實際震例

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

2.2 烈度圖繪制

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

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

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

(3)利用距離反比加權插值法及克里格插值法構建仿真臺站的烈度網格;

(4)根據烈度網格點繪制整個區(qū)域內的地震烈度分布等值線。

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

2.3 與實際調查結果的比較

如圖4所示宏觀調查烈度圖，是根據遼寧省地震局震害宏觀調查結果以及咨詢強震觀測臺站當值人員核實后繪制而成，烈度圖由6 度～3度四個橢圓形烈度等值線組成，長軸方向與震中附近NE 向佟二堡斷裂［10］走向基本一致。從圖3與圖4比較而言，克里格法繪制的烈度圖更貼近于實際宏觀調查結果。

表1 強震動記錄情況Tab.1 Strong motion records

距離反比加權插值法是相近相似的原理，即離插值點越近的樣本點賦予的權重越大，對鄰近采樣點具有依賴性，因此對臺站密度及分布均勻性要求較高。而克里格法強調變量的空間自相關性，考慮了各個鄰近樣點彼此之間的相互關系，所以對臺站密度及均勻性要求相對較低。

由圖3、圖4可見，圖中臺站分布不均勻且密度不足。因此，根據上述分析，克里格插值法相比距離反比加權插值法所繪烈度圖更接近于實際調結果。但鑒于克里格法計算相對復雜，如果臺網足夠密集且均勻，采用計算相對簡單的距離反比加權插值法更利于烈度速報需求。

由圖4與圖3對比可知，宏觀調烈度區(qū)域范圍明顯大于插值方法所繪制烈度區(qū)域，這可能與插值方法所得烈度圖，臺站記錄均為基巖臺站記錄，且沒有經過場地校正，所以烈度區(qū)域偏小也是可以理解的。

圖3 烈度等值線圖Fig.3 Intensity isogram

2.4 兩種插值方法的精度對比

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

表2 交叉驗證誤差統(tǒng)計結果Tab.2 Error statistics of interpolation cross validation

圖4 宏觀烈度調查結果Fig.4 Survey results of macro intensity

3 結論

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

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

［2］張敏政.地震烈度及其評定［J］.防災科技學院學報，2010(1):1－6.

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

［4］李永振，翟文杰，于沈平，等.遼寧省數字強震臺網建設與運行［J］，東北地震研究，2009，25(3):18－22.

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

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

［7］陳晶晶，胡蓓蓓，王軍，等.天津降水數據的空間插值分析［J］，安徽師范大學學報，2010，33(4):15－21.

［8］靳國棟，劉衍聰，牛文杰，等.距離加權反比插值法與克里金插值法的比較［J］，長春工業(yè)大學學報，2003，24(3):53－57.

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

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