模擬退火算法在江蘇及鄰區地震定位中的應用

于悅穎 李正楷 周昱辰 楊 馳

（中國南京 210014 江蘇省地震局）

0 引言

作為一種全局優化算法，模擬退火算法因算法簡單、可對目標函數全局尋優且不易陷入局部極小值，最初被應用于計算機輔助設計、圖像處理、代碼設計等領域，幫助解決組合優化問題（Van Laarhoven et al，1987），繼而被廣泛用于地球物理反演。1994 年Billings將模擬退火算法引入地震定位，發現可以得到更可靠和準確的全局最小值解。此后，利用模擬退火算法，Gao 等（2002）對1999 年9 月21 日中國臺灣集集地震進行了重定位，竇立婷等（2018）、楊波等（2019）對模擬地震數據和少量大陸震例進行了地震定位，尹奇峰等（2019）對井中微地震事件進行了反演定位，均表明模擬退火算法能夠對實際地震事件進行地震定位，具有地球物理實用價值。與Geiger（1912）提出的經典線性方法相比，模擬退火算法在解決地震定位問題時避免了線性方法處理過程中容易出現的高階項省略、初始模型選定等問題，不會出現線性迭代中陷入局部極小值的情況（師學明等，2007；田玥等，2002）。此外，模擬退火算法不需要進行偏導數計算，可以使用更少的計算量定位到全局最小值附近，有效提高了計算效率和精度（Billings，1994）。

為驗證模擬退火算法在江蘇及鄰區地震事件定位的適用性，文中利用P 波走時數據構建全局尋優模型，以2019 年4 月6 日江蘇溧水ML3.3 地震為例，進行反演試算，確定合適的初始溫度、擾動函數和降溫策略，運用于2019 年江蘇省測震臺網中心記錄的105 次ML≥1.8 地震事件重定位，并對定位誤差進行分析。

1 算法原理

模擬退火算法（Simulated Annealing）是一種解決非線性問題的直接反演方法，其基本思路是，假設物質處于狀態m0（即模型參數），該模型對應的能態為E（m0），利用擾動函數，在一定范圍內對初始模型m0進行擾動，得到新模型mi，對應能態為E（mi）；計算2 個模型的能量差ΔE，即ΔE=E（mi）-E（m0），若ΔE＜ 0，則接受新模型mi，若ΔE＞ 0，則有一定概率接受該新模型，如果0—1 之間的隨機數R＜P[P=exp（-ΔE/T）（T為絕對溫度）]，則接受該新模型，并以之替代舊模型，否則重新產生新模型；當有新模型被接受時，即以一定策略進行降溫，不斷循環直至模型收斂，迭代結束（師學明等，2007；Sen et al，2013）。由于ΔE＞ 0 時，新模型有一定概率被選擇，因此的值尤為重要。若T很大，則P將接近1，那么每個新模型（ΔE＞ 0）均有類似概率被接受；若T很小，則P趨近于0，反演結果易陷入局部極小值。因此，T應當與ΔE具有相同數量級。本文使用2 種擾動函數，分別為非常快速模擬退火算法（VFSA）和局部收斂加強型方法。

（1）擾動函數1。Ingber（1989）提出非常快速模擬退火算法（VFSA），公式如下

式中，T為溫度，R為0—1 之間的隨機數，im′為新模型，mi為舊模型，mi的取值范圍是

（2）擾動函數2。Pishro-Nik（2014）提出局部收斂加強型方法，公式如下

式中，T為溫度，R為0—1 之間的隨機數，m′為新模型，m為舊模型，m的取值范圍是[mmin，mmax]。

2 反演模型構建

地震定位可以簡化為利用P 波到時來確定震源位置。設某次地震事件中共有m個臺站記錄到直達波，任一臺站位置為（xk，yk，hk），實際震源位置為（x0，y0，h0），則直達P波到達臺站的走時可以表示為

式中，vP為P 波速度。

利用式（5）進行實際反演，需將經緯度坐標轉化到平面坐標系下，反演結束后再轉化為經緯度坐標系。選用江蘇省測震臺網中心地震定位所用的速度結構，設定vP=6.01 km/s。定義模擬退火算法中的能態為，所有臺站記錄的直達P 波走時tk與P 波拾取震相走時ttrue之差的二范數，即

3 震例試算及參數分析

以江蘇省測震臺網中心記錄的2019 年4 月6 日江蘇溧水ML3.3 地震為例，進行反演試算，確定模擬退火算法中的關鍵參數。在反演過程中，使用20 個地震臺站記錄的P 波到時信息。為確定模擬退火算法中的相關參數，設定初始溫度T為1、25、50，降溫策略設置為T=0.9T和T=0.96T，分析擾動函數1、2 對反演結果的影響。反演參數見表1，定位結果見表2。

表1 反演參數Table1 Inversion parameters

表2 真實及反演定位結果Table 2 Real and inverted earthquake locations

3.1 初始溫度改變

將初始溫度設置為1、25 和50，分別代入降溫策略進行反演，對比相同降溫策略下（同為T=0.9T或T=0.96T）被接受模型的走時差二范數變化曲線，結果見圖1。

由圖1 可見，3 種初始溫度設置對模型收斂影響不明顯。同時，由反演參數可知，走時差二范數的值最終均收斂于Ecov=1.44，故認為1、25 和50 均為合理初始溫度。

圖1 被接受模型的走時差二范數變化曲線Fig.1 The curves of the 2-norm of travel time differences for selected models

由表2 可知，當初始溫度改變時，反演所得定位結果無明顯差別。

3.2 擾動函數改變

在每種溫度條件下，設置反演計算最多迭代200 次，以此選擇合適的反演模型，并按前述規則選擇接受或不接受，當模型被接受時就降溫到下一個溫度，依此循環，直至殘差收斂。由反演參數（表1）可知：①對于擾動函數1、2，當初始溫度升高時，被接受模型的總個數Nm也會增多；②對比總迭代次數i與迭代中被接受模型的總個數Nm可知，擾動函數1 的迭代次數明顯大于擾動函數2，但2 種情況下被接受模型的數量接近，故擾動函數2 比擾動函數1 更高效。

由表2 可知，擾動函數的改變對定位結果無明顯影響。

3.3 降溫策略改變

將降溫系數由0.9 改變為0.96，由反演參數（表1）可知，由于降溫速率變慢，相同初始溫度下的迭代次數增加，被接受模型的數量隨之增加。

由表2 可知，降溫策略的改變，對定位結果的影響并不顯著。

4 實際地震定位

2019 年，江蘇省測震臺網中心記錄到本省及行政區界線外30 km 以內ML≥1.8 地震事件105 次，采用模擬退火算法，對以上事件進行重新定位。根據對典型震例的試算分析，在地震重定位反演中，將初始溫度設置為25，使用效率較高的擾動函數2 產生隨機模型，降溫策略為T=0.96T，且設置低收斂閾值以保證模型收斂。將地震經緯度及深度反演結果與江蘇省測震臺網中心編目結果進行對比，并計算誤差。為清晰對比二者差別，剔除定位結果中經緯度誤差大于1°的4 次黃海海域地震，其余事件定位誤差見圖2、圖3。

圖2 經、緯度坐標誤差分布Fig.2 Error distributions of Longitude and Latitude

（1）震中經、緯度坐標位置誤差（圖2）。在105 次地震中，80 次地震事件的經緯度坐標定位誤差絕對值在0.1°以內，占比約76%；25 次地震事件定位誤差大于0.1°，其中海域地震19 次，陸地地震6 次。

江蘇東部臨海且地震臺站多建設在陸地上，在對海域地震進行定位時，可用臺站偏于一側，在反演過程中，模型范圍的設置存在較大誤差，所以海域地震定位結果較差。對于6 次誤差大于0.1°的陸地地震，可能是由于所使用含有效P 波到時信息的臺站數量少，且臺站空間分布較差，未能較好地包圍震中，導致反演的震中位置結果誤差較大。

（2）震源深度誤差（圖3）。在105 次地震中，31%的地震事件深度定位誤差絕對值大于10 km，由缺乏深度約束所致。

圖3 深度誤差分布Fig.3 Depth error distribution

5 結論

利用模擬退火算法，對2019 年江蘇及鄰區105 次地震重新定位，認為該算法對于江蘇及鄰區陸地地震的定位結果較為可信，具體結論如下：①選取典型震例進行反演試算，通過設置不同的初始溫度，分析不同降溫策略及在擾動函數1、2 下的反演結果，從而確定合適的反演參數；②使用模擬退火算法進行地震定位反演時，經度和緯度的反演結果較好，由于缺乏深度約束，定位深度有一定誤差；③模擬退火算法對P 波到時信息多且地震臺站分布較好的地震事件反演結果更優，對可用臺站偏向于震中一側的海域地震而言，因模型范圍存在明顯誤差，反演結果較差；④對于江蘇省測震臺網中心記錄的地震事件，使用模擬退火算法進行地震定位反演，大部分陸地地震事件重定位結果較為可信，但海域地震事件的重定位結果可信度不足。