直接反演多層密度界面的方法研究*

王 筍 申重陽

1)中國地震局地震研究所(地震大地測量重點實驗室)，武漢 430071)

2)廈門地震勘測研究中心，廈門 361021

3)中國地震局地殼應力研究所武漢創新基地，武漢430071

1 引言

利用重力異常資料反演地下介質密度差異明顯的分界面起伏形態稱為界面反演，對盆地基底研究、區域構造和深部構造起伏研究都具有重要作用［1］。實際應用中往往需要對地下的多個密度界面進行反演，由于重力異常疊加在一起，導致反演的難度很大。目前反演多層密度界面主要用場分離方法，對重力異常的空間頻譜特征進行分析，使特定頻率成分對應于單個界面或深度范圍，把多層密度界面的重力場逐層分離，將剝離后的剩余場按單層界面的方法逐層反演。但場的分離很難做到與各個界面完全對應。

將地下場源區域規則劃分成很多多邊形單元，通過反演確定這些單元的密度變化來約束異常形態，其反演條件自由度高，易操作，模擬的地質模型范圍更廣［2］。但如果想精細地描述場源的物性分布，就要增加單元數量，這樣必然擴大解的空間范圍［3］，而且使計算量急劇增加和出現無意義的解。雖然計算機技術的發展解決了前一個問題，但后一個問題只能通過改進反演方法來解決。

本文選用經典的二度體模型模擬各層密度界面，用地震等資料給出各層界面的深度控制，運用參數范圍約束和光滑約束壓縮解空間，并使用初始模型要求不高的非線性全局優化方法——遺傳算法進行多層密度界面反演計算。

2 算法原理

2.1 模型的建立和參數化

考慮二維情形，若地下空間被呈層狀分布的數種介質充填，不同介質具有不同的密度(圖1)，其間的分界面的起伏將產生附近空間重力位場的變化。根據位場的疊加原理，地面觀測到的重力擾動等于各層分界面產生的重力異常之和;要求取單個界面產生的重力異常，可以其固定的上下界(圖1 上紅色虛線)將附近剖面劃分出來，計算此部分剖面場源產生的重力效應;為簡化計算，記界面兩側密度差值為σ，將界面一側視為密度為0 的空間，另一側用密度σ 充填。

圖1 地下密度界面的等效場源Fig.1 Equivalent mass distribution of density interface

如圖2，把界面起伏視作相對深度D 的一系列水平長方柱體高度Δh 的變化，用一維數組dp0［N］來存放離散化的柱體高度值，以其作為初始模型模型參數。確定這些柱體的剩余密度值σ 后，再根據二度矩形截面體公式計算地面測點處的重力效應。

初始模型參數發生擾動，生成演進模型，隨之發生的重力場變化等于以初始界面為底邊、演進界面為頂邊的所有長方柱體產生的引力之和(圖3)。

用一維數組dpa［N］存放離散化的柱體高度值即演進界面偏離初始界面的程度，以其作為反演對象參數，dpa［N］的取值范圍即界面深度約束范圍，該約束值可以根據鉆孔、地震、電法及地質等已知資料獲得。

2.2 目標函數和光滑約束

用于反演的目標函數包含由模型參數計算得到的重力場值與觀測數據的誤差:考慮M 個不規則分布的重力臺站的觀測數據，對于第i 個重力臺站Pi(xi，zi)，設重力臺站點Pi的平面坐標、重力異常值和重力計算值分別為Pi(xi，zi)、gi，obs(i=1，…，M)和gi，cal(i=1，…，M)，則其均方差為:

在沒有約束的情況下，重力反演是病態的。為避免出現一些明顯不合理的結果，引入緊鄰參數估值相似，使被反演的眾多參數中空間相鄰參數之間的緊密聯系、數值上盡可能地接近。由于這種約束傾向于產生物性光滑過渡的場源，所以這種約束又稱為光滑約束［4］。其形式為:

為使反演問題穩定，采用Tikhonov 正則化方法［5］，使包括兩項的、關于模型擬合差和光滑約束項的目標函數最小，這里選擇的目標函數為:

式中b 為正則化參數。

2.3 遺傳算法

遺傳算法是一種以概率論為基礎的求解多極值非線性問題的最優化方法。

應用遺傳算法進行界面反演，關鍵在于如何對模型參數進行編碼。以往使用的二進制編碼表示一個參數就需要幾位，要求的分辨率越高，編碼越長，當參數眾多時，使用起來很不方便，故不能適應較多變量的反演計算［6］。為適應物性模型場源精細劃分的需要，采用十進制整數編碼方式［7］。

參數編碼采用十進制整數編碼，即把每個參數的編碼用一個整數來表示。N 個參數的L 組模型共有N×L 個編碼。參數群體中的一個編碼代表一個基因，表示一個模型的N 個基因構成一條染色體，群體中共有L 條染色體。對每條染色體譯碼，然后計算出這個模型的目標函數。首先計算初始群體的L 個目標函數，以目標函數作為判據進行搜索。

3 理論模型反演

給定一個試驗模型:三層密度界面，密度差分別為0.08、0.11 和0.43，每個界面各有一個斜坡或凹槽(圖4)。應用非均勻分布的223 個觀測點，其觀測精度相同。不添加任何區域特征的正演計算得到的重力異常(圖5)其特征為鞍形隆起。為驗證反演算法對誤差的敏感度，對異常加入8%隨機誤差。

初始模型設置為3 層水平界面，深度取平均深度。應用本文所用方法搜索解空間，求得擬合差達到設計要求的解，以確定密度界面的形態(圖6)。

由圖6 可見，從只有界面平均深度信息的初始模型出發，通過在光滑約束下對重力場特征的擬合，反演結果具有了和驗證模型相同的特征，結構基本一致，說明本反演程序分離重力異常的效果是可靠的，添加光滑度約束對多層密度界面進行同時反演的方法是可行的。

4 實際模型反演

類烏齊-玉樹-瑪多剖面位于青藏高原東部，剖面所在區域存在4 個主要密度界面:沉積層底界，上地殼底界，中地殼底界，下地殼底界(Moho 面)①王夫運，等.青海玉樹7.1級地震科學考察項目研究報告［R］.2010.。利用鉆井和地震資料可探明結晶基底面的形態，而上、中、下地殼的底界在地震資料上不夠清晰，需要利用重力資料［8］進行計算。若用頻率域場分離的方法反演結晶基底面下方的3 個密度界面，最淺的上地殼底界的重力異常只保留高頻成分，長波長變化信息喪失，界面只在平均深度附近跳動，與地震資料給出的此剖面上地殼底界埋深變化較平緩，總體呈南高北低的趨勢的結果偏離較大。本文根據P波速度-深度剖面，以及P 波速度和密度的經驗公式，建立巖石圈初始界面-密度模型，確定界面深度變化范圍(圖7)。先將沉積層的異常剝離，再運用前述方法同時反演這3 個密度界面的形態，重力異常擬合效果如圖8，均方差約為2 ×10－5ms－2。輸出模型如圖9 所示，Moho 面形態與文獻［7，8］得出的結果基本相同，并展現了上、中、下地殼分界面的起伏變化。

圖9 類烏齊-玉樹-瑪多密度界面模型Fig.9 Final density model of Leiwuqi-Yushu-Maduo profile

5 結論

位場異常的分離不可簡單視為某個頻率段對應于特定深度。一個可行的辦法是對模型使用盡量多的已知資料作控制約束，用鉆井、地震及地質資料等給出深度控制，添加光滑約束消除包含虛假高頻成分的解，達到壓縮解空間的目的，再運用全局優化方法直接反演多層密度界面。

以往用場分離的方法反演多層密度界面，是對重力異常進行濾波，把通過高通濾波器的重力異常歸為淺處的界面的，忽略了較淺界面可能的長波長變化(如整體向一側傾的)，這與很多區域的地下情況不符，無法得到理想的反演結果;本文提出應用光滑約束，其效果與對各層密度界面進行低通濾波類似，同樣使密度界面與一定頻率范圍的重力異常對應，但保留了淺部界面的低頻異常，與實際地下情形符合得更好。另外，觀測到的重力場值不可避免地混有誤差和干擾，混染重力異常的頻譜，嚴重影響場分離的效果;本文所述方法從界面的形態特征出發，有對重力異常的干擾成分不敏感的優點。

1 陳軍，等.應用改進的遺傳算法反演多層密度界面［J］.地球科學(中國地質大學學報)，2000，25(6):651－655.(Chen Jun，et al.Application of improved genetic algorithm to inversion of multi-layer density interface［J］.Earth Science-Journal of China University of Geosciences，2000，25(6):651－655)

2 姚長利，等.重磁遺傳算法三維反演中高速計算及有效存儲方法技術［J］.地球物理學報，2003，46(2):252－258.(Yao Changli，et al.High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms［J］.Chinese Journal of Geophysics，2003，46(2):252－258)

3 Chen Shi，Zhang Jian and Shi Yaolin.Gravity inversion using the frequency characteristics of the density distribution［J］.Applied Geophysics，2008，5(2):99－106.

4 姚長利，郝天珧，管志寧.重磁反演約束條件及三維物性反演技術策略［J］.物探與化探，2002，26(4):253－257.(Yao Changli，Hao Tianyao and Guan Zhining.Restrictions in gravity and magnetic inversions and technical strategy of 3D properties inversion［J］.Geophysical ＆ Geochemical Exploration，2002，26(4):253－257)

5 金其年.Tikhonov 正則化的飽和性與逆結果及后驗參數選取［J］.中國科學(A 輯)，1999，29(8):715－723.(Jin Qinian.Saturation and inverse-results of Tikhonov regularization and selection of posteriori parameter［J］.Science In China(Series A)，1999，29(8):715－723)

6 楊光亮，等.類烏齊-玉樹-瑪多剖面重力異常研究［J］.大地測量與地球動力學，2011，(5):1－4.(Yang Guangliang，et al.Study on gravity anomaly of profile riwoqe-yushumaduo［J］.Journal of Geodesy and Geodynamics，2011，(5):1－4)

7 Yang Guangliang，et al.Joint inversion of gravity and seimic data along a profile across the seimogenic fault of 2010 Yushu Ms7.1 earthquake［J］.Geodesy and Geodynamics，2011，2(4):21－27.Doi:10.3724/SP.J.1246.2011.00021.

8 王謙身，安玉林.青藏高原東部瑪多-沙馬地區的重力場與深部構造［J］.地球物理學進展，2001，16(4):4－10.(Wang Qianshen and An Yulin.Gravity field and deep structure of Maduo-Shama region in eastern Qinghai-Xizang(Tibetan)Plateau［J］.Progress of Geophysics，2001，16(4):4－10.)