模擬退火法在地震定位中的應(yīng)用研究＊

李小晗 李光科 焦大偉

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

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

3) 山東省投融資擔(dān)保集團(tuán)有限公司，濟(jì)南 250000

引言

地震定位是地震學(xué)研究中最經(jīng)典、最基本的問題之一[1-2]，是進(jìn)行其他地震相關(guān)研究工作的基礎(chǔ)。對(duì)于研究震源機(jī)制、地殼速度結(jié)構(gòu)、地球內(nèi)部動(dòng)力學(xué)和地質(zhì)構(gòu)造等基本地震學(xué)問題具有深遠(yuǎn)而重要的意義。此外，如何實(shí)現(xiàn)地震定位工作的高效率、高精確度，對(duì)于地震預(yù)警工程的建立健全和開展，震后的減災(zāi)、救災(zāi)工作都具有重要的意義[3]，從而對(duì)防震減災(zāi)工作發(fā)揮重要的影響和作用。因此，地震學(xué)工作者們一直在通過各種努力不斷改進(jìn)和更新地震定位的方法，積極促進(jìn)地震定位精確度地進(jìn)一步提升。

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

本文基于地震的一般傳播定義和相關(guān)原理公式，通過模擬退火算法反演地震的發(fā)震位置，對(duì)于在地殼中傳播的淺源地震忽略介質(zhì)差異對(duì)傳播速度的影響[4]，同時(shí)認(rèn)為臺(tái)站接收的到時(shí)數(shù)據(jù)所含噪聲均滿足標(biāo)準(zhǔn)正態(tài)分布。

1 模擬退火法的步驟

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

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

（1）確定模型參數(shù)的變化范圍，設(shè)定一個(gè)初始模型m0，計(jì)算相應(yīng)的走時(shí)殘差的二范數(shù)的平方值，即誤差方程E（m0）；

（2）對(duì)當(dāng)前模型的每一個(gè)參數(shù)進(jìn)行擾動(dòng)，產(chǎn)生一個(gè)新模型m1，計(jì)算相應(yīng)的誤差方程E（m1），并計(jì)算得到△E=E（m1）-E（m0）；

（3）若△E＜0，則接受新模型并賦值給舊模型，即m0=m1；若△E≥0，則新模型m1按概率P=exp（-△E/T）進(jìn)行賦值，exp表示自然指數(shù)，T表示溫度。因?yàn)镻（ △E）的函數(shù)取值范圍是（0，1），在0—1之間隨機(jī)產(chǎn)生一個(gè)數(shù)R，當(dāng)P（ △E）＞R時(shí)，接受修改，即m0=m1[6-7]；

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

（5）緩慢降低溫度T；

（6）重復(fù)步驟（2—5），直到△E滿足確定的收斂條件為止。為了避免最優(yōu)的解在重復(fù)過程中被篩掉，我們?cè)谡麄€(gè)搜索過程中隨時(shí)記錄最優(yōu)解?；谀M退火法的特性，如果收斂條件設(shè)置過大，隨著迭代次數(shù)的增加，可能會(huì)跳出全局最優(yōu)解，得到局部最優(yōu)解，此時(shí)，我們可以通過之前的隨時(shí)記錄排除偽全局最優(yōu)解[4]。

2 反演計(jì)算

2.1 建立模型

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

其中，ti為到時(shí)（i=1，2，3，4，…，n），τ為發(fā)震時(shí)刻，為便于計(jì)算，此處設(shè)τ=0。

2.2 數(shù)據(jù)準(zhǔn)備

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

表1 山東地區(qū)速度模型參數(shù)表Table 1 Parameter setting of seismic wave velocity model in Shandong

表2 重慶地區(qū)速度模型參數(shù)表Table 2 Parameter setting of seismic wave velocity model in Chongqing

表3 河北地區(qū)速度模型參數(shù)表Table 3 Parameter setting of seismic wave velocity model in Hebei

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

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

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

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

根據(jù)震中坐標(biāo)及地震觀測(cè)臺(tái)站坐標(biāo)匯總表可以得到震中位置及地震臺(tái)站分布圖（圖1—3）。

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

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

根據(jù)已知震中位置（為方便計(jì)算，假定發(fā)震時(shí)刻為0），利用已知模型公式，計(jì)算得到臺(tái)站到時(shí)。需要注意的是，已知臺(tái)站位置和震中位置均為地理坐標(biāo)，模型公式為平面坐標(biāo)，需要把地理坐標(biāo)轉(zhuǎn)化為平面坐標(biāo)后代入計(jì)算。

2.3 利用模擬退火法進(jìn)行計(jì)算

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

表7 反演結(jié)果對(duì)比（萊州ML3.3）Table 7 The real data and inverse data（Laizhou ML3.3）

表8 反演結(jié)果對(duì)比（武隆ML5.0）Table 8 The real data and inverse data（Wulong ML5.0）

表9 反演結(jié)果對(duì)比（唐山ML4.3）Table 9 The real data and inverse data（Tangshan ML4.3）

圖4—6中，綠色上三角形表示觀測(cè)臺(tái)站，紅色圓圈表示實(shí)際震中位置，藍(lán)色星號(hào)表示反演得到的震中位置，通過以上反演結(jié)果可以看出，模擬退火法反演的震中位置與實(shí)際震中位置較為一致。下面通過改變模擬退火的相關(guān)參數(shù)進(jìn)行對(duì)比分析。

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

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

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

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

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

表10 改變初始溫度T反演結(jié)果對(duì)比分析Table 10 The inverse results after changing T

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

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

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

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

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

由表11和圖8可知，降溫速度直接影響迭代次數(shù)，進(jìn)而影響反演結(jié)果精度。

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

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

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

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

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

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

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

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

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

①全局?jǐn)_動(dòng)

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

②與溫度有關(guān)

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

③和溫度相關(guān)的局部收斂加強(qiáng)型

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

④與溫度有關(guān)

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

表14 改變擾動(dòng)方式反演結(jié)果對(duì)比分析Table 14 The inverse results after changing disturbance mode

3 總結(jié)與建議

通過對(duì)試驗(yàn)結(jié)果的綜合分析可以得出如下結(jié)論：

（1）對(duì)于溫度：①初始溫度的選擇直接影響模擬退火算法全局搜索的整體性能，溫度的初始值越高，全局搜索的范圍會(huì)擴(kuò)大，找到最優(yōu)解的概率也會(huì)大大增加，但是溫度初始值越高需要花費(fèi)搜索的時(shí)間就越多，全局搜索的速率會(huì)下降；反之，溫度的初始值越低，全局搜索的速率會(huì)提高，但找到全局最優(yōu)解的概率可能會(huì)因?yàn)樗阉鞣秶蛔愣档汀＂跍囟乳撝翟礁?，精確度越低，閾值越低，精確度越高，但迭代次數(shù)的增加會(huì)導(dǎo)致計(jì)算時(shí)間的增加。③退火速度，即迭代次數(shù)，退火速度越快，迭代次數(shù)越少，3個(gè)方向的誤差會(huì)越大。通常來說，在相同溫度下進(jìn)行全面的搜索是很有意義的，但進(jìn)行全面充分的搜索意味著需要更多的計(jì)算時(shí)間，循環(huán)次數(shù)增加必定導(dǎo)致計(jì)算時(shí)間的增大。

（2）對(duì)模型參數(shù)進(jìn)行擾動(dòng)后，反演結(jié)果對(duì)3個(gè)參數(shù)的影響均衡，導(dǎo)致在垂直方向上的相對(duì)誤差較大。此方法對(duì)于震源深度的確定誤差相對(duì)較大。

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

（4）模型空間越小越好，這樣可以提高搜索效率。當(dāng)用模擬退火方法進(jìn)行反演，出現(xiàn)收斂結(jié)果總是在模型邊界時(shí)，考慮可能是模型邊界太小所致。

（5）擾動(dòng)函數(shù)有多種形式，可以根據(jù)實(shí)際情況選擇最優(yōu)的擾動(dòng)函數(shù)。

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