兩種優化算法在反演斷層滑動速率中的對比研究＊

李 鵬

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

0 引 言

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

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

1 反演理論與算法

在大地測量反演問題中，觀測值與模型參數的關系常可表示為

式中：d為觀測值；m為模型參數；G為聯系觀測值和模型參數的函數；Δ為觀測誤差。

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

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

1.1 蒙特卡洛法

將反演過程中任何一個階段，用隨機（或偽隨機）發生器產生模型、以實現模型全空間搜索的方法統稱為蒙特卡洛反演法（MC）.

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

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

1.2 蟻群算法

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

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

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

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

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

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

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

2 算例及結果分析

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

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

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

表1 正演模擬斷層參數

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

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

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

表3 待反演參數取值范圍

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

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

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

表4 MC和ACA的反演結果

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

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

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

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

3 結 論

通過理論分析和模擬GPS數據的反演計算分析，可得出以下結論：

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

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］ 王文萍，王慶良.利用遺傳算法和最小二乘聯合反演共和地震位錯參數［J］.地震學報，1999，21（3）：285－290.

［6］ 王衛民，趙連鋒，李 娟，等.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.