王 筍 申重陽

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

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

3)中國地震局地殼應(yīng)力研究所武漢創(chuàng)新基地，武漢430071

1 引言

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

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

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

2 算法原理

2.1 模型的建立和參數(shù)化

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

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

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

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

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

2.2 目標(biāo)函數(shù)和光滑約束

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

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

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

式中b 為正則化參數(shù)。

2.3 遺傳算法

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

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

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

3 理論模型反演

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

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

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

4 實際模型反演

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

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

5 結(jié)論

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

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

1 陳軍，等.應(yīng)用改進的遺傳算法反演多層密度界面［J］.地球科學(xué)(中國地質(zhì)大學(xué)學(xué)報)，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 姚長利，等.重磁遺傳算法三維反演中高速計算及有效存儲方法技術(shù)［J］.地球物理學(xué)報，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 姚長利，郝天珧，管志寧.重磁反演約束條件及三維物性反演技術(shù)策略［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 正則化的飽和性與逆結(jié)果及后驗參數(shù)選?。跩］.中國科學(xué)(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］.大地測量與地球動力學(xué)，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 王謙身，安玉林.青藏高原東部瑪多-沙馬地區(qū)的重力場與深部構(gòu)造［J］.地球物理學(xué)進展，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.)