李 鵬

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

0 引 言

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

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

1 反演理論與算法

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

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

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

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

1.1 蒙特卡洛法

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

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

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

1.2 蟻群算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 MC和ACA走滑反演值對比

圖2 MC和ACA傾滑反演值對比

圖3 MC和ACA張裂反演值對比

3 結(jié) 論

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

1）蟻群算法在斷層滑動速率反演結(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)合反演共和地震位錯參數(shù)［J］.地震學報，1999，21（3）：285－290.

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

［7］ 王家映 .地球物理資料非線性反演方法講座（二）蒙特卡洛法［J］.工程地球物理學報，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.