實時演化地震定位方法及其不確定性分析

盧建旗，李山有，謝志南，馬 強

(中國地震局工程力學研究所 中國地震局地震工程與工程振動重點實驗室，黑龍江 哈爾濱 150080)

地震預警作為一種有效的減災手段，已經在重大工程地震預警(如高速鐵路、燃氣管網、核電)以及城市地震災害預警中得到了廣泛應用[1,15-16]。地震預警的有效性在于其能夠在地震發生后最短的時間內為預警用戶提供反映地震破壞強弱或地震動大小的預警信息，使用戶根據該信息采取相應的緊急處置策略。為了給預警用戶提供盡可能多的預警時間，研究者希望地震預警系統能夠在地震發生后最短的時間內(也就意味著利用盡可能少的信息)對該次地震的最終規模(震級及斷層破裂尺度等)及其影響范圍(地震動場或烈度場)做出準確估計。然而，利用少量信息獲得的地震最終規模及其影響范圍往往帶有很大的不確定性[15]，例如：預警信息提供的估計工程場址的烈度為Ⅵ度，而真實場址的烈度可能為Ⅳ度，也可能為Ⅷ度。這種不確定性來源于場點烈度估計的每一個環節，如P波撿拾的誤差造成震級估算的誤差以及地震定位的誤差，震級估算結果的不確定性，地殼速度結構的不確定性造成的定位方結果的不確定性，場點地震動參數估計的不確定性等等。這使得重大工程難以根據預警信息采取合理的緊急處置策略，極大地降低了地震預警的適用性。

與傳統的地震定位方法不同，地震預警為了提高時效性的要求，定位方法主要依靠臺站的P波走時信息進行定位。例如主站法[10], 該方法不需要地震發生的絕對時間，只要知道兩個臺站的P波到時就可以將震源限定在一個空間的雙曲面內，隨著觸發臺站的增多，形成了多個雙曲面相交的情況，將震源位置限定在多個雙曲面的交線或交點上，從而完成了地震定位。除此之外，還有雙差定位法[6,17,18,21]，基于P波走時差分的最大似然估計定位方法[4,5]等。

利用已觸發臺站的走時信息進行定位只利用到了已觸發臺站的P波走時信息，而沒有用到未觸發臺站的位置信息。而未觸發臺站可以將震中限定在一個更合理的范圍內，不但可以提高定位的精度，還可以提高定位的效率。Horiuchi等[7,9,11,19-20]將未觸發臺站的信息成功運用到地震預警實時定位中，形成了“著未著”地震預警定位方法，從而很好地將未觸發臺站的信息運用到地震預警定位中。Rosenberger等[2]提出了基于高階圖以及P波到時順序的定位方法(AOL)。金星等[12]在現有定位方法的基礎上，提出了一種基于地震臺網資料的網格化地震定位方法以及從單臺觸發到四臺觸發過渡的一套完整的地震預警連續定位方法[13-14,22]。Satriano等[3]在主站法及“著未著”方法的基礎上，提出了一種等時差面與未觸發臺站信息相結合的概率定位方法。該方法很好地將兩種方法進行了結合，并且，通過該方法能夠提高定位概率。

目前，對于定位方法的不確定性概率模型還缺乏研究，在地震預警信息的不確定性分析中，缺乏關于地震定位的不確定性概率模型。

1 基于等時差面及未觸發臺站信息的實時演化定位方法

假設一次地震發生時，各個臺站接收到P波到時均不存在誤差，當有2個臺站被地震觸發后，利用2個臺站的理論到時差與實際到時差可以將震源位置限定在空間的一個雙曲面內，該雙曲面即為等時差面(EDT)。當第3個臺站觸發后，利用3個臺站兩兩到時差與理論到時差就可以將震源位置限定在3個雙曲面的交線或交點上。當有更多臺站P波到時信息后，震源位置就可以限定在空間一點(如圖1(a))。這種定位方法稱為等時差面法[10]。等時差面法的優點在于其不依賴于地震發生的絕對時刻，僅僅利用P波到時就可以進行地震定位。因此，在地震預警定位中得到了廣泛應用。然而，各個臺站的P波到時由于地殼速度結構的不確定性和P波到時識別結果的不確定性，會造成等時差面不相交的情況。為了增加更多的限定信息，人們想到了利用未觸發臺站的信息。假設某一局部平面內有如圖1b所示的臺站分布，該平面由均勻、各向同性介質構成，連接每兩個臺站構成線段的垂直平分線便可形成維諾圖，維諾圖將每一個臺站劃分在一個封閉的多邊形區域內。在該平面內發生一次地震，如果有一個臺站首先觸發，則震源位于包含該臺站的維諾圖多邊形內；如果兩個臺站同時觸發，則震源位于包含兩個臺站的多邊形公共邊上，如果有三個臺站同時觸發則震源位于三個多邊形的公共頂點上。Horiuchi將未觸發臺站的信息成功運用到地震預警實時定位中，形成了“著未著”地震預警定位方法，從而很好地將未觸發臺站的信息運用到地震預警定位中[7,9,11,19,20]

圖1 等時差面與臺站維諾圖示意圖Fig.1 Sketch map for EDT surface and Voronoi cells of seismic stations

基于等時差面及未觸發臺站信息的實時演化定位方法則是綜合利用了已觸發臺站的到時信息和未觸發臺站的到時信息，建立了一種概率定位方法。

1.1 等時差面的概率表述

在不考慮走時誤差和P波撿拾誤差的條件下，當有兩個臺站被觸發，則震源到兩個臺站的理論走時差與實際到時差應該相等(公式1)，則滿足公式(1)的空間點組成了一個空間雙曲面。每兩個已觸發臺站的到時差都可以確定一個雙曲面，隨著觸發臺站的增多，震源位置被限定在多個雙曲面形成的交線或交點上。

(ttm-ttn)i,j,k=tm-tn.

(1)

其中：m和n表示兩個已觸發臺站的編號；tm和tn分別為臺站m和臺站n的P波實際到時；ttm和ttn分別表示震源位于空間網格點(i,j,k)處時，根據地殼速度結構計算得到臺站m和臺站n的P波理論到時。

然而，地殼速度結構具有不確定性，我們根據地殼速度結構計算的P波首波理論到時與實際觀測到的P波首波實際到時差呈以均值為0，標準差為σ1的正態分布。因此，公式(1)可以寫為：

(ttm-tm)±σ1-(ttn-tn)?σ1=0.

(2)

同時，由于臺站的實際到時是通過P波初至自動撿拾方法得到，P波真實到時與自動撿拾到時的殘差呈均值為0，標準差為σ2的正態分布。因此，公式(2)可以進一步改寫為：

(ttm-tm)±σ1±σ2-(ttn-tn)?σ1?σ2=0.

(3)

則根據獨立同分布的規律，等時差面上的理論到時與實際到時差也滿足正態分布：

(4)

其概率密度函數的形式為：

(5)

(6)

由公式(6)可知，當實際到時差與理論到時差相差為0時，定位概率最大，隨著差值變大，概率變小。對所有已觸發臺站進行兩兩組合，并在每個格點對所有臺站對計算的概率求和可得等時差面總概率：

(7)

從概率理論的角度分析，在多個臺站觸發條件下，每個等時差面的概率密度函數在相互獨立的條件下應為求積的關系而不是求和的關系。但在定位過程中，如果存在被噪聲觸發的臺站，采用求積的形式會導致真實等時差面的概率被誤觸發臺站在該處的低概率值所稀釋，從而丟失真實震源位置，而采用求和的形式則不會出現這一現象。因此，這里每個空間網格點的總概率只表示了該點為震源位置的相對可能性大小，而不是絕對意義上的概率。

假設一共有N個臺站參與定位，有nT個臺站已觸發，N-nT個臺站未觸發。則空間網格點處的等時差面概率最大值為：

Qmax=nT(nT-1)/2.

(8)

即，當只有一個臺站觸發時，沒有等時差面，當有2個臺站觸發時有1個等時差面，以此類推。當每個觸發臺站的理論到時與實際到時差在等時差面上正好為0且多個等時差面能夠恰好交于同一網格點，則該網格點的概率值恰好為等時差面概率的最大值Qmax。因此，等時差面概率值的大小能夠體現定位過程中等時差面相交情況。

1.2 未觸發臺站控制范圍的概率表述

如果利用臺站分布建立的維諾圖以及首個觸發臺站位置進行震源位置的空間限定，則存在維諾圖邊界的實時識別、誤觸發臺站的處理等一系列技術問題。Horiuchi首次將未觸發臺站的限制思想引入到地震定位方法中，建立了“著未著”定位方法[7]。在“著未著”算法中，將未觸發臺站的走時與當前絕對時間之間建立不等式關系：

(ttl-ttn)i,j,k≥tnow-tn.

(9)

其中：l表示未觸發臺站的編號；n表示已觸發臺站的編號；tnow表示定位過程中的當前絕對時間；ttl和ttn分別表示震中位于空間網格點(i,j,k)處時計算得到未觸發臺站l和已觸發臺站n的理論到時；tn為已觸發臺站的實際到時。

未觸發臺站的約束不等式定義為概率的形式有：

(10)

即，當震源位于空間格點(i,j,k)處時，如果一個已觸發臺站n和一個未觸發臺站l的理論到時滿足不等式(10)，則該格點為可能的震源位置，其概率值設定為1，否則設定為0。如圖2(a)，圖中有1個已觸發臺站(圖中紅色三角)，6個未觸發臺站(圖中藍色三角)，其余臺站暫不考慮(灰色三角)。在首臺觸發時刻，利用1個已觸發臺站與6個未觸發臺站計算的可能為震源的區域是6條邊界線(概率0、1邊界)合圍的區域，每個未觸發臺站與已觸發臺站計算的概率邊界線正好在維諾圖的邊界上。隨著時間的推移，邊界線合圍區域也在逐漸縮小(如圖2(b))，即可能的震源區域在逐漸縮小。

圖2 隨時間演化的維諾圖示意圖Fig.2 Sketch map for the evolutionary Voronoi cells

同樣，為了避免誤觸發臺站造成真實震源位置被誤觸發臺站排除在外的可能，當有多個未觸發臺站和多個已觸發臺站時，利用每一對已觸發臺站和未觸發臺站計算每個空間網格點的概率Pl,n(i,j,k)，然后采用概率求和的形式計算該點的維諾圖總概率PVoro(i,j,k)：

(11)

假設一共有N個臺站參與定位，有nT個臺站已觸發，剩余N-nT個臺站未觸發，在無誤觸發臺站的條件下，位于震源附近格點處維諾圖總概率的最大值應為：

Pmax=(N-nT)nT.

(12)

即，在無誤觸發臺站的條件下，由維諾圖限定的震源可能區域內網格點處的維諾圖概率最大值應為Pmax，其余區域格點上的維諾圖概率均小于該值。

1.3 聯合定位概率表述

為了讓未觸發臺站形成的實時維諾圖和已觸發臺站形成的實時等時差面在定位過程中同時起到限定震源位置的作用，將等時差面概率與維諾圖概率在對應格點處相加得到每一個格點的聯合定位概率R，并進行歸一化：

(13)

通過搜索聯合定位概率中的最大值即可找到震源所在的位置。歸一化后的聯合定位概率限定其值介于0到1之間，聯合定位概率的大小反映了所有P波到時信息的吻合程度，概率值越大說明吻合程度越高，概率值越小說明吻合程度越低。在利用該值確定的震源點處，聯合定位概率值為1，說明：(1)該點位于滿足所有未觸發臺站限定的可能震源位置區域內；(2)該點位于所有已觸發臺站形成的等時差面相交的公共面、公共線或公共點上。聯合定位概率值小于1說明：(1)該點可能位于未觸發臺站限定的具有相矛盾的區域，即部分未觸發臺站將該點確定為不可能震源位置點，而大部分未觸發臺站則將該點限定為可能的震源位置點；(2)所有已觸發臺站形成的等時差面可能并未完全相交于同一條線或者同一點。因此，聯合定位概率值的大小客觀反映了P波到時數據的不確定性大小，從而在一定程度上也反映了定位精度。但是，該概率值并非真正意義上的概率，而是一種表示震源位置可能性的相對值。

1.4 快速搜索方法

為了提高定位算法的效率，減少一次計算中網格格點的數量，本文采用了固定網格點數的網格間距縮小搜索方法，網格搜索的步驟如下：

(1)在定位算法中事先設定網格點數，當首臺觸發后，將網格中心點置于首臺位置，網格大小可以設置到能夠包含可能的震源區域，計算震源位于每個格點時的聯合定位概率；

(2)搜索網格中最大概率值對應的網格格點，按照固定縮小倍數縮小網格間距，網格點數不變，然后將縮小后的網格中心置于概率值最大的網格點上，再次計算每個網格格點的概率；

(3)重復最大值搜索與網格間距縮小過程，直到網格間距縮小到可接受的程度。

以網格格點數5×5×5為例(如圖3)，網格間距縮小倍數設為0.5。為了便于顯示，本文僅展示了平面內網格搜索過程：

圖3 網格搜索方法示意圖Fig.3 Sketch map for grid searching

1)初次搜索時，將網格中心點置于第一個觸發臺站的位置，計算每個格點的概率。根據圖3(a)中的概率等值線分布情況，最大概率可能位于空間坐標(50,25)處，也可能位于(75,50)處，假設自動搜索的最大值在(50,25)處；

2)將網格間距縮小到原來的一半，且將中心點置于(50,25)處(如圖3(b))，再次計算縮小后網格點處的概率(即，紅色網格格點處的概率)。根據圖3b可知，本次搜索的最大概率位于網格點(62.5,32.5)處；

3)再次縮小網格間距，并將網格中心置于(62.5,32.5)處(如圖3c)，計算縮小后網格點處的概率(即，藍色網格格點處的概率)。依次循環，直到網格間距達到預先設定的最小值，定位結束。

這種搜索方法的優點在于其數組開銷量小，不會在計算機內存中頻繁改變內存地址，只需改變地址中的值，有助于提高效率。在實際定位過程中，如果網格縮小倍數太小，也會導致陷入局部極值的可能。因此，在實際計算過程中，可以適當增大網格間距縮小的倍數，使得定位結果更加可靠。本文采用了9×9×9的網格數，初始網格間距40 km，網格間距縮減倍數為0.8。

2 數據

為了讓本研究中的定位不確定性能夠更真實地反映實際臺網分布以及實際地殼速度模型不確定造成定位結果的不確定性，本文選取了真實地震記錄進行重新定位。由于中國強震觀測臺網數據目前沒有絕對到時，為了滿足這一要求，本文從日本K-NET強震觀測臺網中選取了能夠涵蓋2-9級的地震數據樣本。其中網內地震204次，網外地震22次，震中分布見圖4。由于定位的可靠性與臺網的分布以及臺網密度密切相關，本文計算了日本K-NET臺網間距等值線圖(如圖5 所示)。結果顯示，本文選取的日本K-NET臺站的平均間距約25 km，最小臺站間距可達2 km。該臺站平均間距并不代表日本預警系統的臺站間距，僅代表本文定位分析所用到的臺站的平均間距。在本次定位不確定性分析中，沒有考慮P波撿拾方法造成的P波到時拾取的不確定性，也沒有考慮地震預警系統延時。定位過程中的P波到時使用了各個臺站P波的實際到時，該實際到時通過人工判別方法進行拾取。P波的理論到時采用Crust1.0公布的速度結構進行首波到時計算。

圖4 震中分布圖 圖5 選取的日本K-NET臺站間距等值線圖

3 算例及結果

為了展示地震定位過程中等時差面概率以及未觸發臺站約束概率在定位過程中的變化情況，以及網內地震及網外地震定位概率的變化情況，本文選擇了一次網內地震和一次網外地震。分別利用3臺定位、4臺、5臺至8臺對這次地震進行了定位計算。

選擇的網內地震為2014年6月11日19時52分發生在日本東京附近的一次4.0級地震，震中位于135.685°E,35.085°N，震源深度10 km，共觸發臺站73個。定位過程中的定位概率分布圖如圖6所示。由圖6可以看出，網內地震由于臺站處于震中的四周，每兩個已觸發臺站生成的等時差面與其余等時差面的相交角度較大，定位精度也高。同時，由于未觸發臺站也分布在震中的四周，能夠有效將震源位置限定在一個很小的區域內，并且隨著首播觸發后時間的增加，未觸發臺站限定的震中位置也在不斷縮小。在臺站到時無誤的情況下，3臺就可以獲得很高精度的定位結果。

圖6 網內地震定位概率分布圖Fig.6 Location probability distribution map for an inner-net event

選擇的網外地震為2011年3月11日15時15分發生在日本東部海域的一次7.7級地震，震中位于141.265°E,36.108°N，震源深度43 km，共觸發臺站412個。定位過程中的定位概率分布圖如圖7所示。和網內地震相比，網外地震由于臺站位于震中的一側，而震中的另一側無臺站約束，形成的等時差面的交角也比較小，形成的相交區域大，定位精度低。同時，由于所有臺站都分布在震中一側，未觸發臺站不能將震中約束在一個封閉的區域內，而是只能約束在一個半封閉的區域內。因此，未觸發臺站的約束對定位所起到的效果有限。

注：圖中綠色五角星表示實際震中；藍色三角形表示已觸發臺站；白色三角形表示未觸發臺站；第一列表示等時差面概率，第二列表示未觸發臺站約束概率，第三列表示聯合概率。

利用該方法對本文搜集到的所有地震進行了重新定位，從定位殘差與震級、網內地震、網外地震之間的分布情況(如圖8)可以看出，定位殘差的分布與震級的相關性較小，而與網內地震或網外地震的相關性較大。

圖8 定位殘差隨震級分布圖Fig.8 Residual versus magnitude distribution map

本文對定位殘差的統計分成了網內地震和網外地震兩組，定位殘差如圖9所示。每次地震分別采用了3臺、4臺、5臺及8臺定位。由圖9可知，網內地震的定位殘差普遍較小，網外地震的定位殘差普遍較大，定位殘差符合對數正態分布：

圖9 定位殘差分布圖Fig.9 Distribution map of locating residuals

(14)

其中：x表示定位殘差，取正值；μ和σ分別表示lnx的均值和標準差；f(x)表示對數正態分布的概率密度函數。

4 結論

地震預警定位結果的不確定性是地震預警信息不確定性的主要來源之一。為了分析地震預警定位的不確定性，本文選取日本K-NET臺網記錄到的204次網內地震及22次網外地震，利用地震預警系統中常用的地震定位方法——“實時演化地震定位方法”對這些地震進行了重新定位。結果表明，利用該定位方法獲得的定位殘差符合對數正態分布。網內地震的定位殘差普遍較小，殘差的均值在1～4 km之間；而網外地震由于臺站分布在震中的單側，造成等時差面相交角度較小，未觸發臺站不能將震中約束在一個封閉的區域內，從而導致定位殘差普遍較大。由此可知，網外地震定位精度低也是預警參數不確定性分析中的一個不可忽視的因素。如果將網外地震的震中可能位置放在一個非常大的區域內進行考慮，則會因為震源位置變化太大而導致預警信息的不確定性增加，造成地震預警系統的實用性降低。我國東南沿海也存在大面積網外地震的發震構造。因此，如何提高網外地震的定位精度是地震預警應該考慮的問題。同時，如果地震預警系統能夠識別網外地震，在此基礎上，利用網外地震的定位方法[8]可能是解決這一問題的另一個途徑。

本文選取日本K-NET臺站的臺站平均間距約25 km，與我國地震預警臺網的密度較為接近。因此，該不確定性分析結果對我國地震預警系統定位不確定性有參考價值。

致謝：感謝日本K-NET強震觀測臺網為本文提供數據支持。