李 鵬

（四川西南交大鐵路發(fā)展有限公司，四川 成都610032）

0 引 言

通過地面大地測量數(shù)據(jù)反演大斷裂的滑動(dòng)速率等動(dòng)態(tài)參數(shù)，從而通過地面觀測到的地表變形來認(rèn)識(shí)斷層滑動(dòng)的動(dòng)力過程，是大地測量研究的主要問題之一。一般，在大地測量反演研究中，首先要根據(jù)實(shí)際問題選擇一個(gè)近似地球物理模式，例如在震源參數(shù)反演中，一般采用彈性位錯(cuò)模式［1－2］。若選擇的模型在理論上過于“真實(shí)”，由于數(shù)學(xué)處理困難，常使反演工作無法進(jìn)行。因此，實(shí)際應(yīng)用中必須在“模型”與“算法”之間取折衷方案。從數(shù)學(xué)的角度來看，反演是一個(gè)優(yōu)化過程，就是要尋找一個(gè)能最好地解釋觀測數(shù)據(jù)的物理模型，使實(shí)際觀測值與理論計(jì)算值的差異最小。

大地測量反演問題的算法，國內(nèi)外的有關(guān)學(xué)者作了大量有意義的嘗試，數(shù)值算法，尤其是優(yōu)化反演方法的研究最為普遍，并取得了相當(dāng)豐富的研究成果。如用Monte Carlo法或改進(jìn)的Monte Carlo法反演斷層參數(shù)［3－4］；利用水準(zhǔn)測量資料結(jié)合遺傳算法和最小二乘聯(lián)合反演了共和地震的同震位錯(cuò)參數(shù)［5］；采用模擬退火法結(jié)合GPS和遠(yuǎn)場地震波資料對(duì)1999年臺(tái)灣地震震源破裂過程進(jìn)行了反演計(jì)算研究［6］。為了豐富大地測量數(shù)據(jù)反演問題的求解方法，將一種新的啟發(fā)式算法——蟻群算法進(jìn)行斷層的參數(shù)反演，并基于位錯(cuò)模式，采用模擬的地面位移場，進(jìn)行了斷層三維滑動(dòng)速率的反演，與傳統(tǒng)的優(yōu)化算法蒙特卡洛法反演的結(jié)果進(jìn)行了比較，得出了蟻群算法的可靠性和穩(wěn)定性優(yōu)于蒙特卡洛法等結(jié)論。

1 反演理論與算法

在大地測量反演問題中，觀測值與模型參數(shù)的關(guān)系?？杀硎緸?/p>

式中：d為觀測值；m為模型參數(shù)；G為聯(lián)系觀測值和模型參數(shù)的函數(shù)；Δ為觀測誤差。

在大地測量反演問題中，G通常是非線性的，所以反演是非線性反演。目標(biāo)函數(shù)可表示為

式中：‖·‖2表示2－范數(shù)；E是目標(biāo)函數(shù)值。

1.1 蒙特卡洛法

將反演過程中任何一個(gè)階段，用隨機(jī)（或偽隨機(jī)）發(fā)生器產(chǎn)生模型、以實(shí)現(xiàn)模型全空間搜索的方法統(tǒng)稱為蒙特卡洛反演法（MC）.

假如，已知待求模型參數(shù)的上下界限

有兩種方法對(duì)模型空間進(jìn)行搜索，一種是徹底搜索法，把模型空間允許的范圍都搜索到，看哪一個(gè)模型，或哪一組模型的計(jì)算值和觀測數(shù)據(jù)擬合最好。這種方法也叫窮舉法；另一種搜索法是在模型空間允許的范圍內(nèi)隨機(jī)地搜索，對(duì)每一個(gè)隨機(jī)產(chǎn)生的模型計(jì)算其理論值并把它與觀測值進(jìn)行比較，看其是否可以接受，這就是傳統(tǒng)的蒙特卡洛法。

1.2 蟻群算法

蟻群算法［8］（ACA）是一種嶄新的仿生模擬進(jìn)化算法，由Dorigo等人于20世紀(jì)90年代首先提出。ACA的思想是模擬螞蟻尋食行為，使用大量的螞蟻在搜索空間中隨機(jī)搜索，并且用信息素來加強(qiáng)搜索路線，引導(dǎo)其他螞蟻的搜索，同時(shí)引入信息素的揮發(fā)機(jī)制來避免陷入局部最優(yōu)，這種引入揮發(fā)機(jī)制的正反饋使得該算法能夠找到全局的多個(gè)最優(yōu)解。

Dorigo等人充分利用蟻群搜索食物的過程與著名的旅行商問題 （TSP）之間的相似性，提出了蟻群算法，用人工螞蟻模擬自然螞蟻，通過模擬螞蟻搜索食物的過程求解復(fù)雜的組合優(yōu)化問題。一般地，可以用ACA在TSP的求解來說明ACA算法思想。

TSP是指對(duì)于給定的一組城市，設(shè)定其每兩個(gè)城市之間的距離為已知，要求出一條總長度最小的封閉路線，該路線剛好只通過每個(gè)城市一次。記城市數(shù)為n，城市i與城市j之間的距離為dij（i，j＝1，2，…，n），螞蟻數(shù)目為m，城市i與城市j之間路徑上在t時(shí)刻的信息素量為τij（t），初始時(shí)刻每條路徑上的信息素為τij（0）。螞蟻k（k＝1，2，…，m）在面臨路徑選擇的時(shí)候按照如下的轉(zhuǎn)移概率來決定下一個(gè)去的城市

假設(shè)tabuk表示螞蟻k已經(jīng)訪問過的城市列表，訪問規(guī)則為每個(gè)城市只訪問一次，所以要從可訪問集合中刪除已訪問的城市，記allowedk＝｛N－tabuk｝為螞蟻k還能去的城市集合。信息素更新規(guī)則為：螞蟻在兩個(gè)城市之間移動(dòng)一步之后，都會(huì)增加城市之間路徑上的信息素濃度，考慮到要防止快速陷入局部解，引入信息素?fù)]發(fā)機(jī)制，所以一步之后更新為

式中：ρ為信息素殘留系數(shù)（0≤ρ≤1）；Δτij（t）和分別為蟻群與螞蟻k在時(shí)間段t到（t＋1）內(nèi)，在邊（i，j）上留下的信息素濃度，Δτkij表示為

式中：Q為常量；Lk為螞蟻k在本次循環(huán)中所選擇路徑的總長度。

參數(shù)Q，α，β，ρ的最佳組合可由實(shí)驗(yàn)確定［9］，在蟻群算法中，當(dāng)路徑穩(wěn)定或者停止條件滿足后，搜索結(jié)束。

2 算例及結(jié)果分析

為了驗(yàn)證兩種優(yōu)化算法的有效性和穩(wěn)定性，基于位錯(cuò)模式，將位錯(cuò)理論模型［1－2］模擬計(jì)算的地面位移場作為觀測值，模擬計(jì)算采用的某斷層參數(shù)如表1所示。

表1中緯度、經(jīng)度、H為斷層起始端點(diǎn)坐標(biāo)；U1、U2、U3為斷層面上盤在走向、傾向和斷層面法線方向的滑動(dòng)量。L、W、D分別表示斷層面的長度、寬度和下底面深度，A和f表示斷層的走向和傾向。

表2示出了利用位錯(cuò)理論正演模擬計(jì)算得出的部分地面位移場。

表1 正演模擬斷層參數(shù)

表2 部分正演模擬的地面水平位移

利用表2所示的模擬GPS的觀測數(shù)據(jù)在斷層的其它參數(shù)不變的情況下，對(duì)斷層的三維滑動(dòng)速率進(jìn)行了反演計(jì)算分析，表2形式的模擬數(shù)據(jù)共527行，表2只示出了其中的一小部分；采用VB 6.0語言結(jié)合位錯(cuò)理論模型分別編制了蒙特卡洛法和蟻群算法反演的計(jì)算程序。

為了比較兩種算法的可靠性和穩(wěn)定性，在進(jìn)行反演時(shí)，采用相同的取值范圍，并且離散程度相同，計(jì)算相同的次數(shù)，共進(jìn)行計(jì)算35次，3個(gè)待反演參數(shù)的取值范圍如表3所示。

表3 待反演參數(shù)取值范圍

蒙特卡洛法反演的基本過程簡述如下：通過文件讀入模擬GPS觀測數(shù)據(jù)及觀測點(diǎn)的坐標(biāo)，確定每個(gè)參數(shù)的先驗(yàn)信息（即取值范圍見表3），并將每個(gè)參數(shù)空間離散成10 000份，然后采用隨機(jī)函數(shù)產(chǎn)生一組參數(shù)值。另外，由位錯(cuò)理論模型根據(jù)隨機(jī)斷層初始參數(shù)計(jì)算觀測點(diǎn)的位移場，再由隨機(jī)參數(shù)計(jì)算的位移場與模擬觀測的位移場求出目標(biāo)函數(shù)E，如果滿足式（2）的約束條件，則輸出反演結(jié)果。

蟻群算法反演的基本過程簡述如下：通過文件讀入斷層初始參數(shù)和模擬GPS觀測數(shù)據(jù)及觀測點(diǎn)的坐標(biāo)，將每個(gè)參數(shù)空間的取值在（表3）區(qū)間內(nèi)離散成10 000份，采用隨機(jī)函數(shù)產(chǎn)生一組參數(shù)值。另外，由位錯(cuò)理論模型根據(jù)隨機(jī)斷層初始參數(shù)計(jì)算觀測點(diǎn)的位移場，再由隨機(jī)參數(shù)計(jì)算的位移場與模擬觀測的位移場求出目標(biāo)函數(shù)E，并設(shè)螞蟻k在某次循環(huán)中所選擇路徑的長度Lk＝E，由公式（4）～（6）建立蟻群算法的遞推關(guān)系，最后由 （4）的大小來確定斷層參數(shù)的更新，在實(shí)際計(jì)算中，由于參數(shù)之間不存在距離的概念，能見度ηij（t）取1進(jìn)行計(jì)算。

表4分別示出了兩種算法在進(jìn)行35次計(jì)算中的一組最佳反演結(jié)果。

表4 MC和ACA的反演結(jié)果

圖1、圖2和圖3分別為在上述取值范圍和計(jì)算次數(shù)為35次的情況下，蒙特卡洛法和蟻群算法在斷層三維滑動(dòng)速率走滑、傾滑和張裂的反演計(jì)算結(jié)果對(duì)比圖。從表4的結(jié)果以及圖1、圖2和圖3所示的反演結(jié)果對(duì)比分析可以看出，蟻群算法得出的結(jié)果的穩(wěn)定性和可靠性優(yōu)于蒙特卡洛法。

圖1 MC和ACA走滑反演值對(duì)比

圖2 MC和ACA傾滑反演值對(duì)比

圖3 MC和ACA張裂反演值對(duì)比

3 結(jié) 論

通過理論分析和模擬GPS數(shù)據(jù)的反演計(jì)算分析，可得出以下結(jié)論：

1）蟻群算法在斷層滑動(dòng)速率反演結(jié)果中的穩(wěn)定性和可靠性優(yōu)于蒙特卡洛法，而且蒙特卡洛法在傾滑、張裂的反演值與理論值的偏離程度較大，明顯差于蟻群算法。

2）算例的斷層主要以走滑影響為主，因此，兩種算法走滑的反演值與理論值的吻合程度較傾滑、張裂好些。

［1］ OKADA Y.Surface deformation due to shear and tensile faults in a half－space［J］.BSSA，1985（82）：1018－1040.

［2］ OKADA Y.Internal deformation due to shear and tensile faults in a half－space［J］.BSSA，1992（82）：1018－1040.

［3］ MURRAY M H，MARSHALL G A，LISOWSKI M，et al.The 1992M＝7Cape Mendocino，California，earthquake：coseismic deformation at the south end of the Cascadia megathrust［J］.J.Geophys.Res.1996，101（B8）：17707－17725.

［4］ FTEYMULLER J.Kinematics of the pacific－north America plate boundary zone，northern California［J］.J.Geophys.Res.，1999，104（B4）：7419－7441.

［5］ 王文萍，王慶良.利用遺傳算法和最小二乘聯(lián)合反演共和地震位錯(cuò)參數(shù)［J］.地震學(xué)報(bào)，1999，21（3）：285－290.

［6］ 王衛(wèi)民，趙連鋒，李 娟，等.1999年臺(tái)灣集集地震震源破裂過程［J］.地球物理學(xué)報(bào)，2005，48（1）：132－147.

［7］ 王家映 .地球物理資料非線性反演方法講座（二）蒙特卡洛法［J］.工程地球物理學(xué)報(bào)，2007，4（2）：81－85.

［8］ DORIGO M，MANIEZZO V，COLOMI A.The ant system：optimization by ant colony of cooperating agents［J］.IEEE Trans.Syst.Man.Cybern－PartB，1996，26（1）：29－41.

［9］ DORIGO M，GAMBARDELL L M.Ant colonies for the traveling salesman problem，bioSystems［J］.1997（43）：73－78.