模擬退火法在地震定位中的應用研究＊

李小晗 李光科 焦大偉

1) 山東省地震局，濟南 250102

2) 重慶市地震局，重慶 401147

3) 山東省投融資擔保集團有限公司，濟南 250000

引言

地震定位是地震學研究中最經典、最基本的問題之一[1-2]，是進行其他地震相關研究工作的基礎。對于研究震源機制、地殼速度結構、地球內部動力學和地質構造等基本地震學問題具有深遠而重要的意義。此外，如何實現地震定位工作的高效率、高精確度，對于地震預警工程的建立健全和開展，震后的減災、救災工作都具有重要的意義[3]，從而對防震減災工作發揮重要的影響和作用。因此，地震學工作者們一直在通過各種努力不斷改進和更新地震定位的方法，積極促進地震定位精確度地進一步提升。

地震定位指的是通過對地震波到達地震臺站的時間進行觀測，來確定地震震源的坐標以及發震時刻[1]。震源位置的確定，為我們進行地震活動性研究、社會地震預警、開展抗震救援等工作提供了科學基礎。目前最經典的地震定位方法是蓋格法，是一種線性反演方法，該方法的優點是計算速度快，但是對于非線性問題，它的缺陷也十分明顯，容易陷入局部極小值，反而會降低定位效率。而模擬退火算法作為一種非線性反演方法，避免了反演結果陷入局部最優解，并且隨著退火不斷選取當前最優解的過程，在概率上能夠更快地獲得全局最優解。綜合來看，該方法可以更好地提高地震定位的效率和精確度。

本文基于地震的一般傳播定義和相關原理公式，通過模擬退火算法反演地震的發震位置，對于在地殼中傳播的淺源地震忽略介質差異對傳播速度的影響[4]，同時認為臺站接收的到時數據所含噪聲均滿足標準正態分布。

1 模擬退火法的步驟

模擬退火算法（Simulated Annealing，SA）是一種非線性反演，可以避免反演結果陷入局部范圍內，不準確的局部范圍只會得到局部極大值而不是全局最優值。模擬退火本質上是現代蒙特卡洛方法，是持續找尋最優解的過程，在這個過程中，擬合度隨著迭代次數的增加跳躍起伏，不斷變化，但總體趨勢是擬合度越來越高，正是由于模擬退火算法允許擬合度變小或擬合誤差變大，從而使模型從局部最優解中跳出，得到全局最優解[4-5]。

本文使用了Metropolis algorithm模擬退火方法，其相應步驟如下：

（1）確定模型參數的變化范圍，設定一個初始模型m0，計算相應的走時殘差的二范數的平方值，即誤差方程E（m0）；

（2）對當前模型的每一個參數進行擾動，產生一個新模型m1，計算相應的誤差方程E（m1），并計算得到△E=E（m1）-E（m0）；

（3）若△E＜0，則接受新模型并賦值給舊模型，即m0=m1；若△E≥0，則新模型m1按概率P=exp（-△E/T）進行賦值，exp表示自然指數，T表示溫度。因為P（ △E）的函數取值范圍是（0，1），在0—1之間隨機產生一個數R，當P（ △E）＞R時，接受修改，即m0=m1[6-7]；

（4）在溫度T下，重復步驟（2）和（3）；

（5）緩慢降低溫度T；

（6）重復步驟（2—5），直到△E滿足確定的收斂條件為止。為了避免最優的解在重復過程中被篩掉，我們在整個搜索過程中隨時記錄最優解。基于模擬退火法的特性，如果收斂條件設置過大，隨著迭代次數的增加，可能會跳出全局最優解，得到局部最優解，此時，我們可以通過之前的隨時記錄排除偽全局最優解[4]。

2 反演計算

2.1 建立模型

在此，我們使用震源的三維坐標為（x0，y0，h0），其中（x0，y0）為震源投影到地球上的經度和緯度的平面坐標，h0為震源深度[2]。震源附近有n個地震臺觀測到了直達波P的走時，臺站坐標為（xi，yi，hi），根據上述參數和條件可以建立如下非線性模型：

其中，ti為到時（i=1，2，3，4，…，n），τ為發震時刻，為便于計算，此處設τ=0。

2.2 數據準備

采用山東、重慶和河北地區的一維速度模型，速度模型參數如表1—3[8-10]。

表1 山東地區速度模型參數表Table 1 Parameter setting of seismic wave velocity model in Shandong

表2 重慶地區速度模型參數表Table 2 Parameter setting of seismic wave velocity model in Chongqing

表3 河北地區速度模型參數表Table 3 Parameter setting of seismic wave velocity model in Hebei

使用2012年8月8日17時48分55秒山東省煙臺市萊州ML3.3地震、2017年11月23日17時43分33秒重慶市武隆區ML5.0地震和2016年9月10日18時09分37秒河北唐山ML4.3地震觀測到地震波的臺站地理坐標。接收到地震信號的臺站地理坐標如表4—6（本文所使用的地震經緯度精確到小數點后3位，考慮到臺站經緯度的保密性要求，信息表中只寫到小數點后一位）。

表4 地震臺站信息表（萊州ML3.3）Table 4 Seismic station information（Laizhou ML3.3）

表5 地震臺站信息表（武隆ML5.0）Table 5 Seismic station information（Wulong ML5.0）

表6 地震臺站信息表（唐山ML4.3）Table 6 Seismic station information（Tangshan ML4.3）

根據震中坐標及地震觀測臺站坐標匯總表可以得到震中位置及地震臺站分布圖（圖1—3）。

圖2 震中及臺站分布圖（武隆ML5.0）Fig. 2 Distribution of epicenter and station coordinates（Wulong ML5.0）

圖3 震中及臺站分布圖（唐山ML4.3）Fig. 3 Distribution of epicenter and station coordinates（Tangshan ML4.3）

根據已知震中位置（為方便計算，假定發震時刻為0），利用已知模型公式，計算得到臺站到時。需要注意的是，已知臺站位置和震中位置均為地理坐標，模型公式為平面坐標，需要把地理坐標轉化為平面坐標后代入計算。

2.3 利用模擬退火法進行計算

初始溫度T= 258，降溫速度T=T×0.989，收斂閾值為1×e-15時，進行模擬退火所得結果如圖4—6，得出的結果如表7—9。

表7 反演結果對比（萊州ML3.3）Table 7 The real data and inverse data（Laizhou ML3.3）

表8 反演結果對比（武隆ML5.0）Table 8 The real data and inverse data（Wulong ML5.0）

表9 反演結果對比（唐山ML4.3）Table 9 The real data and inverse data（Tangshan ML4.3）

圖4—6中，綠色上三角形表示觀測臺站，紅色圓圈表示實際震中位置，藍色星號表示反演得到的震中位置，通過以上反演結果可以看出，模擬退火法反演的震中位置與實際震中位置較為一致。下面通過改變模擬退火的相關參數進行對比分析。

圖4 震中位置與模擬退火反演震中位置對比（萊州ML3.3）Fig. 4 Coordinates of epicenter and inverse result（Laizhou ML3.3）

圖5 震中位置與模擬退火反演震中位置對比（武隆ML5.0）Fig. 5 Coordinates of epicenter and inverse result（Wulong ML5.0）

圖6 震中位置與模擬退火反演震中位置對比（唐山ML4.3）Fig. 6 Coordinates of epicenter and inverse result（Tangshan ML4.3）

由表10和圖7可知，溫度初始值設定越大，模擬退火的搜索范圍就越大，迭代次數會增加，同時，反演結果更接近真實結果。

（1）以山東萊州ML3.3地震為例，設降溫速度為T=T×0.98，Tm=0.01，收斂閾值為1×e-10，依次改變模擬退火的溫度初始值，進行反演結果（表10和圖7）的對比分析。

表10 改變初始溫度T反演結果對比分析Table 10 The inverse results after changing T

圖7 改變初始溫度T反演結果與實際結果差值圖Fig. 7 The difference between inverse result and real result after changing T

圖1 震中及臺站分布圖（萊州ML3.3）Fig. 1 Distribution of epicenter and station coordinates（Laizhou ML3.3）

（2）假設溫度T= 100，Tm= 0.01，收斂閾值為1×e-10，依次改變降溫速度，進行反演結果（表11和圖8）的對比分析。

表11 改變降溫速度反演結果對比分析Table 11 The inverse results after changing cooling rate

圖8 改變降溫速度反演結果與實際結果差值圖Fig. 8 The difference between inverse result and real result after changing cooling rate

由表11和圖8可知，降溫速度直接影響迭代次數，進而影響反演結果精度。

（3）假設溫度為T= 100，Tm= 0.01，降溫速度為T=T×0.98，依次改變收斂閾值，進行反演結果（表12和圖9）的對比分析。

表12 改變收斂閾值反演結果對比分析Table 12 The inverse results after changing convergence threshold

圖9 改變收斂閾值結果與實際結果差值圖Fig. 9 The difference between inverse result and real result after changing convergence threshold

由表12和圖9可知，收斂閾值大時，迭代次數少，精度相對差；收斂閾值小時，迭代次數多，精度好。

（4）假設溫度為T= 100，降溫速度為T=T×0.98，收斂閾值為1×e-10，依次改變溫度閾值Tm，進行反演結果（表13和圖10）的對比分析。

表13 改變溫度閾值反演結果對比分析Table 13 The inverse results after changing temperature threshold

圖10 改變溫度閾值結果與實際結果差值圖Fig. 10 The difference between inverse result and real result after changing temperature threshold

由表13和圖10可知，設置溫度的下限決定了最大迭代次數，進而影響模擬退火的迭代次數，溫度下限越小，迭代次數增加，反演結果相對較好。

（5）假設溫度為T= 100，降溫速度為T=T×0.98，Tm=0.01，收斂閾值為1×e-10，分別使用4種擾動方式：

①全局擾動

其中，m1為模型新值，r=U[0，1]。

②與溫度有關

其中，U為擾動因子，T為溫度，m1為模型新值，m0為模型舊值，r=U[0，1]。

③和溫度相關的局部收斂加強型

其中，U為擾動因子，T為溫度，m1為模型新值，m0為模型舊值，r=U[0，1]。

④與溫度有關

由表14可知，擾動方式②、③和④的結果較為一致，擾動方式①誤差過大，不建議采用。

表14 改變擾動方式反演結果對比分析Table 14 The inverse results after changing disturbance mode

3 總結與建議

通過對試驗結果的綜合分析可以得出如下結論：

（1）對于溫度：①初始溫度的選擇直接影響模擬退火算法全局搜索的整體性能，溫度的初始值越高，全局搜索的范圍會擴大，找到最優解的概率也會大大增加，但是溫度初始值越高需要花費搜索的時間就越多，全局搜索的速率會下降；反之，溫度的初始值越低，全局搜索的速率會提高，但找到全局最優解的概率可能會因為搜索范圍不足而降低。②溫度閾值越高，精確度越低，閾值越低，精確度越高，但迭代次數的增加會導致計算時間的增加。③退火速度，即迭代次數，退火速度越快，迭代次數越少，3個方向的誤差會越大。通常來說，在相同溫度下進行全面的搜索是很有意義的，但進行全面充分的搜索意味著需要更多的計算時間，循環次數增加必定導致計算時間的增大。

（2）對模型參數進行擾動后，反演結果對3個參數的影響均衡，導致在垂直方向上的相對誤差較大。此方法對于震源深度的確定誤差相對較大。

（3）收斂閾值越小，反演結果越好。

（4）模型空間越小越好，這樣可以提高搜索效率。當用模擬退火方法進行反演，出現收斂結果總是在模型邊界時，考慮可能是模型邊界太小所致。

（5）擾動函數有多種形式，可以根據實際情況選擇最優的擾動函數。

（6）在此問題中，反演的震中位置與實際位置總是存在相應的誤差，可能原因有：①采用簡單一維速度模型，換算過程中存在近似，影響了精確度；②計算過程中加入的隨機噪聲有影響；③本文使用的模擬退火方法需要完善，增加迭代次數，選擇更合適的初始溫度和退火方式等，從而獲得最優的結果。