柳旭峰 許才軍

（中國武漢430079武漢大學(xué)測繪學(xué)院）

引言

視震源時間函數(shù)（apparent source time function，簡寫為ASTF）的提取對于研究震源破裂機制和反演震源破裂過程非常重要.臺站記錄到的地震波形相當(dāng)于格林函數(shù)、視震源時間函數(shù)（對于點源地震可視為震源時間函數(shù)）與儀器響應(yīng)三者的褶積，通過反褶積過程，可以得到視震源時間函數(shù)（張勇，2008；許力生，陳運泰，2002）.

在實際求取視震源時間函數(shù)時，往往要限制一些條件，一般包括：非負(fù)性、因果性和持續(xù)時間有限性（Bertero et al，1997；張勇，2008）.反演視震源時間函數(shù)，可以在頻率域中進(jìn)行，常用的方法有水平線方法（water-level method）（Mori，Hartzell，1990）和光滑頻譜方法（Radulian，Popa，1996）等.其優(yōu)點是反演效率較高，缺點是反演結(jié)果誤差較大，而且無法添加約束條件；也可以在時間域中進(jìn)行，如采用共軛梯度法（許力生，陳運泰，2004）和Tikhonov正則化方法（Kraeva，2004）等，這些方法能夠顧及約束條件，但反演結(jié)果隨反演方法的不同而不同，并有可能陷入局部極值.Bertero等（1997）提出了在頻率域中迭代的映射Landweber反褶積（projected landweber deconvolution，簡寫為PLD）方法，在迭代過程中將結(jié)果變換到時間域來添加約束條件，計算效率比較高，但缺點是當(dāng)假定的持續(xù)時間與實際持續(xù)時間不符時，反演結(jié)果誤差較大，需要在一個時間區(qū)間內(nèi)對持續(xù)時間搜索來確定實際持續(xù)時間.張勇（2008）用遺傳算法（generic algorithm，簡寫為GA）對視震源時間函數(shù)進(jìn)行提取，得到了比PLD方法更好的結(jié)果，但是計算效率偏低.本文旨在利用改進(jìn)的粒子群算法（particle swarm optimization，簡寫為PSO）反演視震源時間函數(shù)，以克服遺傳算法效率偏低的缺點并獲取較好的反演結(jié)果.

1 改進(jìn)的PSO算法

PSO算法是由Kennedy和Eberhart（1995）提出的一種仿生學(xué)進(jìn)化算法.在這個算法中，解空間中的每個解稱之為“粒子”，每個粒子都有一個判別函數(shù)來量度到最優(yōu)解的距離，同時還有一個速度決定粒子更新的方向和大小，在迭代的過程中，每個粒子都追隨當(dāng)前空間的最優(yōu)粒子（Poli et al，2007）.

設(shè)PSO算法解空間有m個粒子，初始值為x0＝（x01，x02，…，x0m），初始速度為v0＝（v01，v02，…，v0m）.在每次迭代中，每個粒子通過跟蹤兩個極值來更新自己，第一個極值為粒子本身在迭代過程中的最優(yōu)解，稱為個體極值（pbest）；另一個極值是整個種群當(dāng)前的最優(yōu)解，稱為全局極值（gbest）.設(shè)在第n次迭代后，解空間為xn＝（xn1，xn2，…，xnm），對應(yīng)的速度為vn＝（vn1，vn2，…，vnm），則第n＋1次解空間xn＋1可以表達(dá)為

式中，w為慣性因子，表示前一時刻粒子速度對當(dāng)前時刻粒子速度的影響，取值范圍為［0，1］；c1是粒子自身速度學(xué)習(xí)因子，表示粒子向自身最優(yōu)值學(xué)習(xí)；c2是群體速度學(xué)習(xí)因子，表示粒子向群體最優(yōu)值學(xué)習(xí)；r1和r2為0與1之間的隨機值，是調(diào)節(jié)粒子向自身學(xué)習(xí)和群體學(xué)習(xí)的權(quán)重系數(shù).每次迭代中，計算每個個體的判別函數(shù)值F，并更新個體極值和全局極值.判別函數(shù)值F是當(dāng)前解與最優(yōu)解之間差值的量度，一般在［0，1］之間，當(dāng)F為0時，表示達(dá)到最優(yōu)解.同時為了限制PSO算法的更新速度過大或者過小，往往根據(jù)經(jīng)驗或者試算結(jié)果假定粒子最大速度vmax和最小速度vmin，超過則取最大值，不足則取最小值.在迭代過程中，還可以隨機對空間中的粒子的速度進(jìn)行變異操作，以達(dá)到全局優(yōu)化的目的.

PSO算法雖然屬于全局優(yōu)化算法，但在迭代過程中，所有解都趨向于最優(yōu)解進(jìn)行更新，有可能陷入局部解.另外，相對于遺傳算法，PSO算法中引入速度概念，加快了算法的收斂過程.但如果慣性因子不變，則算法在后期容易來回震蕩，導(dǎo)致效率降低.

圖1 改進(jìn)的PSO算法流程圖Fig.1 Flow-chart of the improved PSO algorithm

考慮到PSO算法存在的問題，本文對PSO算法進(jìn)行了改進(jìn)（改進(jìn)方法的流程如圖1所示）.具體內(nèi)容包括以下4點：

1）初值的選取.為了提高PSO算法的收斂速度，并且克服陷入局部解的傾向，本文利用水平線方法得到的初值添加約束之后的結(jié)果作為PSO算法初始全局最優(yōu)解.水平線方法原理如下：

式中，ε稱為閾值因子；max［ε，｜G（w）｜2］表示取ε和｜G（w）｜2的較大值.得到的S（w）經(jīng)過反傅里葉變換，可以得到在時間域中表示的視震源時間函數(shù).根據(jù)非負(fù)性和持續(xù)時間有限性約束，令結(jié)果中負(fù)值為0，并且在假定持續(xù)時間之外的值為0，將得到的結(jié)果作為PSO算法的初始全局最優(yōu)解.在初始全局最優(yōu)解周圍一定范圍內(nèi)隨機產(chǎn)生m個解，形成PSO算法的初始解空間.

2）判別函數(shù)的選擇.判別函數(shù)表示得到的結(jié)果與真實結(jié)果的距離大小的量度.本文中采用直觀的方式.

設(shè)U（t）表示真實主震記錄，Usyn（t）表示通過解空間中的解得到的合成記錄.定義數(shù)據(jù)誤差Vd為

式中，t表示持續(xù)時間.令F＝Vd表示判別函數(shù).

3）對慣性因子的改進(jìn).為了使算法迭代后期不至于在最優(yōu)值附近來回震蕩，要求慣性因子在判別函數(shù)變小時減小，因此設(shè)計

式中，c0和c1為常數(shù)；F（n）min表示第n次迭代解空間判別函數(shù)的最小值；F（n）i表示第n次迭代解空間中第i個解的判別函數(shù)值；后面的指數(shù)函數(shù)隨著迭代次數(shù)n的增大而減小.c1的數(shù)值不能取太大或者太小.太大會使得這個函數(shù)沒有意義，起不到相應(yīng)的減小速度的效果；太小則會使迭代速度減小過于迅速，算法的效率下降.實際計算中，c0和c1的值根據(jù)經(jīng)驗取得.

4）為了提高收斂速度，去掉了PSO算法的自身學(xué)習(xí)成分，即令PSO算法的自身學(xué)習(xí)權(quán)重c1＝0.

改進(jìn)的PSO算法的迭代公式為

2 利用改進(jìn)的PSO算法反演模擬視震源時間函數(shù)

采用一個包含兩個峰值的模擬視震源時間函數(shù)作為真值，一個經(jīng)過低通濾波的地震波信號作為經(jīng)驗格林函數(shù)（empirical Green’s function，簡寫為EGF），合成模擬主震記錄.同時，對合成的主震記錄添加噪聲，在不同的噪聲水平下，采用改進(jìn)的PSO算法、PLD方法以及GA算法分別對視震源時間函數(shù)進(jìn)行反演，對比說明改進(jìn)方法的優(yōu)缺點.

在合成的模擬主震記錄中加入噪聲，關(guān)系式為

圖2 合成主震記錄（a）、0.25Hz濾波后的EGF（b）和模擬 ASTF（c）圖Fig.2 Synthetic main shock record（a），empirical Green function（EGF）filtered by 0.25Hz low pass（b）and simulated ASTF（c）

式中，＊表示卷積；α為噪聲水平；r（t）為［0，1］上均值為0、方差為1的高斯白噪聲序列；max［Ge（t）＊Sreal（t）］表示取模擬主震記錄振幅的最大值.

模擬視震源時間函數(shù)總共持續(xù)時間40s，包含兩段過程，每段采用高斯函數(shù)生成.經(jīng)驗格林函數(shù)采用GNI臺站觀測到的2005年克什米爾MW7.6地震的一次余震記錄的P波數(shù)據(jù)，經(jīng)過8階Butterworth低通濾波器進(jìn)行濾波，濾波上限為0.25Hz，通過與模擬的視震源時間函數(shù)褶積得到模擬的主震記錄（圖2）.

定義模型誤差Vm為

式中，S（t）為反演得到的視震源時間函數(shù)，t為持續(xù)時間.

采用0，0.1和0.3的噪聲水平，當(dāng)模型誤差Vm＝0.001時終止迭代.假定最大迭代次數(shù)為2 000（試驗表明，PLD方法迭代2 000次之后模型誤差基本不變），并且假定視震源時間函數(shù)持續(xù)時間長度為100s，反演結(jié)果如圖3—5所示.

反演資料時間間隔為0.25s，反演過程中水平線方法的閾值因子采用0.3，三種方法的初始值都選擇了水平線方法得到的結(jié)果并加入了約束條件.在迭代過程中構(gòu)造拉普拉斯方程進(jìn)行時間光滑.令▽2s（t）＝0，即s（t－1）＋s（t＋1）－2s（t）＝0.PLD方法的松弛因子選擇為1／｜Gmax（w）｜2；PSO算法種群個數(shù)為200，式（4）中c0取值為2，式（5）中c2取值為2；GA算法中種群個數(shù)為200，種群變異因子為0.4，交叉因子為0.5，選擇因子為0.1，在計算過程中采用了退化方法：設(shè)第i次迭代的結(jié)果為Si（t），m為Si（t）的中位數(shù)，那么第i次迭代的實際結(jié)果為si（t）－mi.

從圖3—5反演的結(jié)果可以看出，改進(jìn)的PSO算法與GA算法結(jié)果相似，但都好于PLD方法.在持續(xù)時間（40s）之外，隨著噪聲的增大，PLD方法的誤差變大，而改進(jìn)的PSO方法和GA算法誤差都比PLD方法誤差小.

從模型誤差可以看出，改進(jìn)的PSO算法能夠用較少的迭代次數(shù)得到較優(yōu)結(jié)果.PLD方法得到的模型誤差在一定的迭代次數(shù)之后趨于穩(wěn)定，反演結(jié)果無法繼續(xù)收斂，表現(xiàn)為最大幅值一直處于0.9以下.而且PLD方法在各種噪聲水平下都無法收斂到終止迭代誤差（即模型誤差為0.001）.其原因在于，由于添加了拉普拉斯平滑，使得PLD迭代結(jié)果失真，從而導(dǎo)致最終結(jié)果失真，而對不添加拉普拉斯平滑的PLD方法進(jìn)行的試驗結(jié)果表明，噪聲存在情況下，PLD方法的反演結(jié)果不平滑，模型誤差較大.在0.3噪聲水平下，3種方法都無法收斂.從圖5可以看出，改進(jìn)的PSO算法得到的結(jié)果與實際視震源時間函數(shù)最為吻合.同樣迭代2 000次，改進(jìn)的PSO算法的計算時間比GA算法要少一倍，而且PSO算法在迭代200次左右模型誤差就不發(fā)生變化.考慮到這一點，在0.3噪聲水平下3種方法所用時間比例與低噪聲水平時基本不變.本研究中多次計算統(tǒng)計表明，改進(jìn)的PSO算法相比GA算法計算時間減少了80%以上，效率提高了至少5倍.

圖6 不同噪聲水平下改進(jìn)的PSO算法、PLD方法和GA算法計算得到的模型誤差與數(shù)據(jù)誤差（a），（b），（c）分別對應(yīng)0.05，0.1和0.3的噪聲水平.圓圈表示PLD方法所得結(jié)果，星號為改進(jìn)的PSO算法所得結(jié)果，十字為GA算法所得結(jié)果Fig.6 Model errors and data errors inversed from improved PSO，PLD and GA under different noise levels（a），（b）and（c）correspond to noise levels 0.05，0.1and 0.3，respectively.Circles，stars and crosses stand for inversion results from PLD，the improved PSO and GA，respectively

圖7 PLD、改進(jìn)的PSO和GA三種方法的平均運行時間Fig.7 Average computing time for using PLD，the improved PSO and GA

為了進(jìn)一步說明改進(jìn)的PSO算法、GA算法以及PLD方法在反演結(jié)果上的優(yōu)劣，分別在0.05，0.1和0.3噪聲水平下采用3種方法進(jìn)行了計算.改進(jìn)的PSO算法粒子個數(shù)為200，迭代次數(shù)200；PLD方法迭代次數(shù)為2 000；GA算法種群個數(shù)采用1 000，迭代次數(shù)500，其它參數(shù)不變.以模型誤差和數(shù)據(jù)誤差作為指標(biāo)，此時3種方法得到的反演結(jié)果基本不變（圖6）.計算圖6所需的3種方法的平均時間如圖7所示.

從圖6和圖7可以看到，采用改進(jìn)的PSO算法和GA算法得到的結(jié)果相對于PLD方法更為穩(wěn)定，表現(xiàn)為分布更為集中.改進(jìn)的PSO算法在種群個數(shù)和迭代次數(shù)均較少的情況下能夠得到與GA算法相當(dāng)?shù)男Ч?可見改進(jìn)的PSO算法一方面能夠克服GA算法效率低的特點，另外還能夠得到相比PLD方法更為準(zhǔn)確的結(jié)果.

3 2005年巴基斯坦克什米爾MW7.6地震P波視震源時間函數(shù)反演

2005年10月8日巴基斯坦克什米爾地區(qū)發(fā)生了MW7.6地震.該地震位于印度板塊與歐亞板塊碰撞帶上喜馬拉雅山地震帶西部.此處印度板塊向北運動，與歐亞板塊發(fā)生碰撞，印度洋板塊向下俯沖，引發(fā)地震.繼主震發(fā)生之后，又發(fā)生了多次余震（表1）.因此，我們可以利用經(jīng)驗格林函數(shù)的方法提取視震源時間函數(shù).

表1 2005年克什米爾MW7.6地震主震及余震序列Table 1 Information on the 2005Kashmir MW7.6main shock and its aftershocks

張勇（2008）采用GA算法和PLD方法，利用31個臺站觀測到的主震和余震數(shù)據(jù)對視震源時間函數(shù)進(jìn)行了反演.結(jié)果表明，此次地震的視震源時間函數(shù)持續(xù)時間在25s以內(nèi)，P波視震源時間函數(shù)相比S波的要短.

張勇等（2009）利用PLD方法對此次地震的P波、S波、瑞雷波和勒夫波視震源時間函數(shù)進(jìn)行了提取.結(jié)果表明，此次地震的視震源時間函數(shù)平均持續(xù)時間為25s，不同震相的結(jié)果形狀大致相似，細(xì)節(jié)略有不同.

本文利用改進(jìn)的PSO算法對此次地震的P波視震源時間函數(shù)進(jìn)行提取.從美國地震學(xué)聯(lián)合研究會（IRIS）數(shù)據(jù)中心下載了64個臺站記錄的克什米爾地震的主震及11次余震的Z分量數(shù)據(jù)，并選擇信噪比大于4的臺站數(shù)據(jù).通過去趨勢項、去平均項，并進(jìn)行降采樣處理，使資料采樣間隔統(tǒng)一變?yōu)?.5s.采用8階Butterworth濾波器進(jìn)行低通濾波，濾波上限為0.25Hz.參考P波理論到時，截取P波初到時與S波初到時之間的數(shù)據(jù).假設(shè)視震源時間函數(shù)持續(xù)時間為50s，利用經(jīng)驗格林函數(shù)的方法，對P波視震源時間函數(shù)進(jìn)行提取，并且對每個臺站得到的多個視震源時間函數(shù)歸一化后進(jìn)行疊加.剔除掉反演較差的結(jié)果，總共得到40個臺站的P波視震源時間函數(shù)，如圖8、圖9所示.

通過對臺站的視震源時間函數(shù)進(jìn)行對比可以看出，得到的P波視震源時間函數(shù)持續(xù)時間在25s以內(nèi)，各個臺站得到的視震源時間函數(shù)相似，說明利用改進(jìn)的PSO方法能夠有效地提取視震源時間函數(shù).對各個不同方位角的臺站的視震源時間函數(shù)的對比顯示，位于震中西北方向的視震源時間函數(shù)相對于東南方向較快地出現(xiàn)峰值，由此可見，震源的破裂方向應(yīng)該是西北向的.這些結(jié)果與張勇（2008）和張勇等（2009）的結(jié)果一致.

研究中還發(fā)現(xiàn)，臺站資料的信噪比對于提取視震源時間函數(shù)非常重要，臺站記錄到的信噪比越高越好.

4 討論與結(jié)論

視震源時間函數(shù)的提取對于研究震源破裂過程十分重要.傳統(tǒng)的頻率域反演方法，如水平線方法，往往無法加入約束；而傳統(tǒng)的時間域方法，如共軛梯度法，則反演效率低，并且易于陷入局部極值.本文研究表明，PLD方法同時顧及約束和效率，得到的結(jié)果好于傳統(tǒng)方法.而GA算法能夠得到比PLD更好的結(jié)果，但GA算法計算效率相對較低.考慮到PSO算法在GA算法基礎(chǔ)上引入了速度的概念，能夠極大提高計算效率.因此本文提出利用改進(jìn)的PSO算法對視震源時間函數(shù)進(jìn)行反演研究.

利用模擬數(shù)據(jù)，對改進(jìn)的PSO算法、GA算法和PLD方法的計算效率和計算結(jié)果的準(zhǔn)確性進(jìn)行了對比.結(jié)果表明，改進(jìn)的PSO算法比GA的計算效率提高了5倍以上，而相比于PLD算法，改進(jìn)的PSO算法計算所得結(jié)果的準(zhǔn)確性則有很大提高.說明了改進(jìn)的PSO算法能夠有效地反演視震源時間函數(shù).

利用改進(jìn)的PSO算法對2005年巴基斯坦克什米爾MW7.6地震的P波視震源時間函數(shù)進(jìn)行了提取，得到了40個臺站的視震源時間函數(shù).結(jié)果表明，此次地震的P波視震源時間函數(shù)持續(xù)時間在25s以內(nèi)，破裂方向為西北方向.該結(jié)果與張勇（2008）和張勇等（2009）的結(jié)果一致.

此外需要指出的是，臺站數(shù)據(jù)質(zhì)量對于提取視震源時間函數(shù)也是非常重要的.

許力生，陳運泰.2002.震源時間函數(shù)與震源破裂過程［J］.地震地磁觀測與研究，23（6）：1－－7.

許力生，陳運泰.2004.從全球長周期波形資料反演2001年11月14日昆侖山口地震時空破裂過程［J］.中國科學(xué)：D輯，34（3）：256－－264.

張勇.2008.震源破裂過程反演方法研究［D］.北京：北京大學(xué)地球物理學(xué)系：1－－70.

張勇，許力生，陳運泰.2009.提取視震源時間函數(shù)的PLD方法及其對2005年克什米爾MW7.6地震的應(yīng)用［J］.地球物理學(xué)報，52（3）：672－－680.

Bertero M，Bindi D，Boccacci P，Cattaneo M，Eva C，Lanza V.1997.Application of the projected Landweber method to the estimation of the source time function in seismology［J］.Inverse Probl，13（2）：465－－486.

Kennedy J，Eberhart R.1995.Particle swarm optimization［C］∥Proceedings of IEEE International Conference on Neural Networks，4：1942－－1948.

Kraeva N.2004.Tikhonov’s regularization for deconvolution in the empirical Green function method and vertical directivity effect［J］.Tectonophysics，383（1）：29－－44.

Mori J，Hartzell S H.1990.Source pulse enhancement by deconvolution of an empirical Green’s function［J］.Geophys Res Lett，12（1）：33－－36.

Poli R，Kennedy J，Blackwell T.2007.Particle swarm optimization［J］.Swarm Intell，1（1）：33－－57.

Radulian M，Popa M.1996.Scaling of source parameters for Vrancea（Romania）intermediate depth earthquake［J］.Tectonophysics，261（1）：67－－81.