吳昊天,李 軍,b,王堃鵬,b,李 霞

(成都理工大學 a.地球物理學院，b.教育部地球探測與信息技術重點實驗室，成都 610059)

0 前言

重磁勘探作為地球物理勘探中的一種常用方法，它是根據(jù)地下目標體與圍巖的密度或磁性差異引起的重磁場的空間變化，解決某些地質(zhì)問題的一種勘探方法。重磁反演是重磁勘探解釋工作中的重要內(nèi)容，但由于重磁位場的“趨膚效應”和體積效應明顯，其間接觀測手段也決定了重磁反演問題多是不適定問題，導致了重磁反演存在現(xiàn)實困難和技術瓶頸。

重磁物性反演過程本質(zhì)是求解大型欠定型線性方程組，計算效率和求解過程的穩(wěn)定性和解的非唯一性是重磁物性反演的主要技術瓶頸。計算效率目前主要通過大型并行計算和改進矩陣計算和儲存方式加以克服，而求解的穩(wěn)定性和非唯一性問題則主要是通過先驗約束和算法實現(xiàn)。B.J.Last[1]提出的緊密反演算法方法，就是在反演迭代的過程中加入最小體積這一先驗信息，作為一約束條件，這一方法能夠讓結果異常體邊界更加尖銳，物性分布更加集中，對整塊地質(zhì)體有較好的反映；Valeria Cristina F.Barbosa等[2]也提出了先加入先驗信息的反演算法，這種方法對于軸線性的目標體有較強的分辨能力；Li等[3]所提出了最光滑反演方法，通過在反演算法中加入模型的最小梯度這一先驗信息，作為約束條件，能夠讓異常體的邊界緩慢過渡，也可以較好地反映目標體的空間位置；Oleg Portniaguine等[4]提出了聚焦反演算法，可以讓異常體邊界更加明顯。Li 和 Oldenburg以及Pilkington[3,5-6]在反演算法中引入物性深度函數(shù)加權和光滑矩陣的約束，一定程度克服了反演目標的趨膚效應。Fedi等[7]對三維位場數(shù)據(jù)進行反演計算，解決了因為核函數(shù)隨深度增加而遞減造成的趨膚效應；Fedi等[8]針對通過SVD和GSVD方法對反演過程中的深度分辨率的問題進行分析；Fedi[9]中提出了DEXP方法，通過極值點預測異常體的深度，這是一種快速且穩(wěn)定的方法；Cella[10]應用構造指數(shù)作為加權函數(shù)的參數(shù)來抵消核函數(shù)隨深度的衰減。

顯然，先驗約束和算法改進能有效克服密度反演的不穩(wěn)定性和非唯一性。然而，有效建立先驗約束，并在算法上尋求優(yōu)化協(xié)作，保證兩者能貫穿整個反演求解過程，充分改進其計算的穩(wěn)定性和結果的可靠性，仍是一個難點。

共軛梯度反演方法具有計算速度快、內(nèi)存需求少，能夠避免大型矩陣的乘積計算的特點，已成為大規(guī)模計算的實用反演方法。近年來預條件共軛梯度法和重加權共軛梯度法也被應用在位場數(shù)據(jù)反演中[11-15]。筆者試圖在目標函數(shù)中加入約束項以及深度加權函數(shù)，利用基于PRP公式的重加權正則化共軛梯度法進行求解，期望能有效改進大型欠定方程的求解效率與精度和重磁物性反演結果的可靠性。

1 正演理論

重磁正演效率決定了重磁異常反演效率及其實用性。多面體構建的復雜重磁模型在重磁反演過程人機交互過高，其可操作性和實用性較差。這里設計基于網(wǎng)格剖分的模型構建和重磁異常正演計算，即對地下半空間進行規(guī)則剖分，剖分的每個規(guī)則六面體單元都賦予以密度值或磁化強度值(圖1)。

地面某一觀測點的響應可根據(jù)地下每個剖分單元的重磁場計算結果進行累加。剖分單元的重磁異常可采用直立長方體表達式計算：

Δg(x,y,z)=-Gσ(ξ-x){Ln[(η-y)+r]+(η-y)Ln[(ξ-x)+r]

(1)

(2)

若不考慮剩磁的作用，則：

k1=2MN,k2=2NL,k3=2ML,

k4=L2,k5=M2,k6=-N2，

L=cosIcosD,M=cosIsinD,N=sinI

式中：Δg(x,y,z)為某觀測點地下長方體單元引起的重力異常;ΔT(x,y,z)為某觀測點地下長方體單元引起的磁異常;G為常數(shù);μ0為真空磁導率;σ為密度;M為磁化強度;I為磁傾角;D為磁偏角;(ξ;η,ζ)為剖分單元體積微元的中心坐標;(ξ1,ξ2,η1,η2,ζ1,ζ2)為地下長方體單元的角點坐標位置;r為長方體單元到觀測點距離。

地表某點的觀測數(shù)值與地下各剖分單元其物性參數(shù)所產(chǎn)生的場值呈線性關系：

Gm=d

(3)

其中：

(4)

圖1 長方體單元模型Fig.1 Rectangular unit model

式中：G為N×M階的核矩陣;其元素gi,j表示第i個測點由第j個剖分單元產(chǎn)生的場值;m為求解模型向量，共有M個模型;d為觀測數(shù)據(jù)向量，共有N個觀測數(shù)據(jù)。

2 約束反演理論

2.1 深度加權

由于重磁場具有“趨膚效應”特征，即趨于地表的地質(zhì)體對觀測點處所觀測的重力異常權重較大，而離地表較遠的地質(zhì)體對觀測結果影響較小，故反演物性會反映到更敏感的淺部，結果基本趨于地表，這會導致反演結果與實際情況相差較大。Li和Oldenburg[3,5]提出在反演過程中引入一個深度加權函數(shù)進行約束，能有效克服了重磁場隨深度呈指數(shù)衰減的特征。這里采用這一深度加權函數(shù)：

(5)

式中：z為剖分單元的中心點埋深；z0取決于觀測面位置；β為深度加權因子。根據(jù)重磁場隨深度的衰減規(guī)律，一般在重力反演中取2.0，磁異常反演中取3.0。

2.2 光滑約束

在三維重磁線性反演，模型變量參數(shù)遠遠大于采集數(shù)據(jù)，所以反演往往是欠定的問題，則必須要在目標函數(shù)中加入模型約束項，以約束所求得解的范圍。筆者引入光滑約束函數(shù)，其意義在于使得相鄰的模型之差最小，使模型在各個方向上是光滑過渡的。定義模型粗糙度R為模型m在X、Y和Z方向的一階偏微分的平方和，即：

(6)

此時進行深度加權以及光滑約束下的模型約束項為：

(7)

其中

(8)

令：

(9)

式(7)中：Z為深度加權函數(shù)組成的對角矩陣；D為模型約束矩陣；as、ax、ay、az為光滑權重。

2.3 構建正則化目標函數(shù)

根據(jù)Tikhonov正則化理論，可以建立正則化目標函數(shù)來避免反演中出現(xiàn)的多解性問題。正則化目標函數(shù)包括兩部分：①約束穩(wěn)定函數(shù);②數(shù)據(jù)擬合函數(shù)。為了抵消核函數(shù)的“趨膚效應”，增大其垂向分辨率，用深度加權函數(shù)對核矩陣進行改造，則有以下線性關系：

GZ-1Zm=d

(10)

令：

Gw=GZ-1

(11)

mw=Zm

(12)

可得到目標函數(shù)的數(shù)據(jù)擬合函數(shù)式(13)與約束穩(wěn)定函數(shù)式(14)。

(13)

(14)

因此，正則化目標函數(shù)表達式為式(15)。

(15)

2.4 正則化因子的選取

在Tikhonov正則化反演過程中，正則化因子體現(xiàn)了數(shù)據(jù)擬合與模型約束之間的取舍與平衡。當正則化因子較大時，會優(yōu)先進行模型約束，正則化因子過大可能造成無法擬合數(shù)據(jù)。當正則化因子較小時，會優(yōu)先擬合數(shù)據(jù)，正則化因子過小可能會造成無法對模型約束，使得反演結果無法可靠反映地下的模型分布。

一般采用經(jīng)驗選取法、廣義交叉驗證法和L-curve法選擇正則化因子。經(jīng)驗選取法憑借經(jīng)驗選取正則化因子具有較強的主觀性。廣義交叉驗證法選取正則化因子包含了大規(guī)模矩陣的多次運算，L-curve法需要進行多次的反演運算才能選擇合適的正則化因子。后兩者的計算大，增加了計算成本。

筆者選取了一種自適應的正則化因子的選擇方法，選擇隨一定步長下降的正則化因子[11]。令μ0=0，計算出經(jīng)過初次迭代后的模型量，再使用式(16)計算出首個的正則化因子。

(16)

每次迭代后根據(jù)式(17)計算正則化因子。

μk+1=qμk，k=1,2,…,n

(17)

式中，q通常取0.5～0.9。根據(jù)式(16)與式(17)計算正則化因子。隨著迭代次數(shù)的增加，正則化因子減小，能夠優(yōu)先擬合模型而后對數(shù)據(jù)進行擬合，保證了反演的穩(wěn)定。該選取正則化因子的方法與廣義交叉驗證法以及L-curve法比較，有效地減少了計算量，達到了正則化因子自適應的目的。

2.5 基于PRP公式的重加權正則化共軛梯度法

共軛梯度法用迭代替代了大型矩陣的乘積與矩陣求逆的過程，加快了反演求解速度。在反演中使用傳統(tǒng)的共軛梯度算法，進行界限約束或改變正則化因子時都需要重構殘差，對算法進行“重啟”。該情況下的共軛梯度算法便會退化為梯度下降法，降低了收斂速度與穩(wěn)定性。筆者采用的是重加權正則化共軛梯度法，可以實現(xiàn)每次迭代都改變權重，以配合自適應的正則化因子，實現(xiàn)了數(shù)據(jù)擬合與模型約束的權重的動態(tài)變化。保證了約束質(zhì)量的情況下加快收斂速度。

圖2 基于PRP公式的重加權正則化共軛梯度算法流程圖Fig.2 RRCG based on PRP formula

2.6 稀疏矩陣的優(yōu)化儲存

每個光滑約束矩陣中僅存在±C、0三種數(shù)值(C取值見本文2.2部分)，在剖分數(shù)量大時，若將該矩陣儲存下來將會占用大量內(nèi)存，加上核矩陣后，將會儲存(N+4M)×M大小的矩陣，對如此大的矩陣進行乘積運算將耗費大量的計算時間。

(18)

圖3 正演模型以及噪聲污染的正演數(shù)據(jù)Fig.3 The forward model and the forward data contaminated by noise(a)三維模型示意圖；(b)重力異常正演數(shù)據(jù)；(c)噪聲污染的重力異常正演數(shù)據(jù)；(d)磁異常正演數(shù)據(jù)；(e)磁異常垂直磁化數(shù)據(jù)；(f)噪聲污染的磁異常數(shù)據(jù)

圖4 不同公式下的重力反演與磁反演相對擬合差變化圖Fig.4 The relative misfit of different formulas in gravity inversion and magnetic inversion(a)重力反演下不同公式的相對擬合差變化圖；(b)磁反演下不同公式的相對擬合差變化圖

(19)

圖5 不同迭代次數(shù)下的磁反演結果Fig.5 Magnetic inversion results under different iteration times(a)迭代60次y=525 m處剖面；(b)迭代60次z=125 m處截面；(c)迭代80次y=525 m處剖面；(d)迭代80次z=125 m處截面；(e)迭代100次y=525 m處剖面；(f)迭代100次z=125 m處截面

3 模型反演

筆者引用了姚長利等[13,17]使用過的Y型巖脈模型。如圖3所示，左側(cè)異常體I剩余密度為0.8 g/cm3，磁化強度為0.8 A/m，向下延伸至200 m。右側(cè)異常體II剩余密度為1g / cm3，磁化強度為1 A / m，向下延伸至500 m，二者磁傾角皆為45°，磁偏角皆為45°，而反演時只考慮垂直磁化情況下的數(shù)據(jù)。將地下三維空間剖分為20×20×12=4 800個單元，其邊長為50 m。共有21×21=441個采集數(shù)據(jù)，數(shù)據(jù)采集點位于z=0處。在觀測數(shù)據(jù)中，加入了5%的高斯噪聲污染(圖3(c),圖3(f))，以測試反演針對實測數(shù)據(jù)的能力。反演的物性參數(shù)范圍控制為0 g/cm3～1 g / cm3與0 A/m～1 A/m。

迭代100次，重力反演花費時間為8.593 s，磁反演花費時間為9.135 s，其中因為磁異常公式復雜，計算核矩陣的時間花費較多。

重力反演與磁反演在PRP公式與FR公式下的相對擬合差隨迭代次數(shù)的變化曲線如圖4所示。

圖6 不同迭代次數(shù)下的重力反演結果Fig.6 Gravity inversion results under different iteration times(a)迭代60次y=525 m處剖面；(b)迭代60次z=125 m處截面；(c)迭代80次y=525 m處剖面；(d)迭代80次z=125 m處截面；(e)迭代100次y=525 m處剖面；(f)迭代100次z=125 m處截面

表1 重力反演的相對擬合差Tab.1 The relative misfit of different formulas in gravity inversion

表2 磁反演的擬合差Tab.2 The relative misfit of different formulas in magnetic inversion

由圖4可以看到，在重力反演與磁反演中，分別迭代40次與16次之后，PRP公式開始低于FR公式下的相對擬合差，且相對擬合差保持穩(wěn)定下降的趨勢。

表1和表2為重、磁反演于PRP公式與FR公式下不同迭代次數(shù)的相對擬合差變化表。當?shù)螖?shù)較小時，F(xiàn)R公式下的相對擬合差略微小于PRP公式下的相對擬合差。隨著迭代次數(shù)的增加，PRP公式下的相對擬合差逐漸小于FR公式下的相對擬合差。反映該模型下PRP公式下的反演計算具有較高的收斂速度和計算效率。

如圖5與圖6所示為加入噪聲的重力異常數(shù)據(jù)和磁異常數(shù)據(jù)在PRP公式和FR公式下分別進行60次、80次和100次的反演迭代結果?？傮w而言，重力反演以及磁異常反演在加入噪聲的情況下能夠在一定程度上的反映異常體的位置以及形態(tài)，沒有出現(xiàn)過擬合的現(xiàn)象。

由于深度加權對核函數(shù)進行了改造，反演計算中的模型量是由深部向淺部移動的，所以在收斂慢的情況下，模型量在深部聚集的更多。

對于重力反演，與FR公式相比，PRP公式下的反演剖面淺部位置的模型量能夠更快迭代得到，異常體形態(tài)相對更加緊湊，模型底界面深度較為可靠。對于地下截面，PRP公式下的反演結果，分辨率相對較高，邊界刻畫更加明顯，且模型量更接近實際模型的值。同時，PRP公式下的數(shù)據(jù)擬合結果更加接近原始數(shù)據(jù)。所以PRP公式在該模型下的重力反演有著較快的收斂速度和更高的計算精度。

圖7 馬角壩3號地質(zhì)剖面磁異常圖Fig.7 Magnetic anomaly map of geological section no.3 of Majiaoba

圖8 馬角壩3號地質(zhì)剖面磁異常磁化強度反演圖Fig.8 Inversion map of magnetic anomaly magnetization in geological section No.3 of Majiaoba(a)數(shù)據(jù)擬合結果；(b)反演結果

在磁反演中，與FR公式相比，PRP公式下的反演剖面，其淺部與深部模型分布更加合理。FR公式中的結果中，其異常體I下方聚集的模型量較大且恢復的模型在垂向上的形態(tài)寬大。對于地下截面，在迭代60次之下，PRP公式的結果中反映的模型量稍大，更加接近真實值。而在迭代80次與與100次時，二者的結果并沒有很明顯的差異，且出現(xiàn)了少量因為噪音而產(chǎn)生的模型分布，但可以忽略不計。同時，PRP公式下的數(shù)據(jù)擬合與FR公式下的數(shù)據(jù)擬合結果，稍微更接近原始數(shù)據(jù)曲線，但差別較小。則PRP公式在該模型下的磁反演同樣有著較快的收斂速度和較高的計算精度。

對比重力反演與磁反演結果，對于y=525 m處剖面，由于異常體I較為寬厚，且異常體II在深部延伸至異常體I之下，由于等值現(xiàn)象，在重力反演以及磁異常反演的結果中，均有左側(cè)模型量較右側(cè)模型量大的結果出現(xiàn)，且在磁異常反演中因為異常體II延伸至異常體I之下，抵消了部分異常體I的負異常，導致反演出的磁化強度偏大且聚集。二者在深部反演的結果也有一定程度向下延伸，呈現(xiàn)光滑過渡的形態(tài)。對于z=125 m處截面，由于磁異常數(shù)據(jù)本身的較高分辨率，且因磁反演的深度權重較重力反演大的影響，磁異常反演的結果更加緊湊，分辨率較好。

從二者整體反演結果來看，反演結果基本與Y型巖脈一致。在淺部的形態(tài)相對于深部的清晰，但反演結果不可避免地隨著深度下降而降低分辨率。經(jīng)以上討論表明：

1)基于PRP公式的重加權正則化共軛梯度的約束反演能夠在一定噪聲干擾下從一定程度上反映地下物性分布情況，相對FR公式有較高的計算效率和計算精度。

2)在重磁反演結果對比中，由于磁異常數(shù)據(jù)較為復雜，相鄰緊密的模型異常等效作用更明顯，會存在南北極耦合及正負異性抵消的情況，而重力數(shù)據(jù)這種情況表現(xiàn)相對較弱。

3)由于磁異常數(shù)據(jù)本身的分辨率更高，且其深度權重較大反演能夠得以較快地收斂，其反演結果的分辨率更高。

4 實際資料應用

成都理工大學馬角壩實習基地地質(zhì)3號線自火焰包到灌林包，長約1 200 m，地處龍門山系前山帶的邊緣。在三疊紀末的印支運動中，受西北造山帶造山作用的影響，全區(qū)褶皺上升成陸，結束海侵。在由西北向東南推擠的應力作用下，本區(qū)地層形成軸向北東的褶皺和一系列北東向延伸、向西北傾斜走向的逆沖斷層和與之配套的近南北向和北西西方向的平移斷層。同時伴有小規(guī)模的巖漿沿斷裂或裂隙侵入，形成輝綠巖墻。

高精度磁測異常表明地質(zhì)線路上在飛仙關二段有明顯的磁異常，標本測量揭示該段地層粉砂質(zhì)泥巖的磁化率最大可達100×10-5(SI)，相對其他地層有較明顯的磁性，可引起50 nT～100 nT以上的異常。但結合本次測量，可以注意到在飛仙關二段的整體異常高達600 nT～800 nT，這說明該地質(zhì)線路的高磁異常并非完全由飛仙關二段地層引起。

筆者利用該條磁測剖面進行了三維磁化強度反演，設計正演網(wǎng)格模型為240×41×40個單元，其網(wǎng)格單元邊長為5 m。共有400個采集數(shù)據(jù)。采用重加權正則化共軛梯度法約束反演，迭代150次，耗時26 min 35 s，得到如下異常(圖8)的反演結果與擬合結果。

反演結果表明測線下方飛仙關二段地層存在明顯的局部地質(zhì)體，其磁化強度最高可達3.5 A/m，這一磁化強度顯然不是該沉積地層的響應。但我們注意到該測線的三個主要異常，有兩個異常的傾向于地層傾向角度近似(北西傾向45°～57°)，故初步判斷該異?？赡転楹笃谘財嗔亚秩氲木植枯x綠巖墻。但由于該區(qū)后期構造改造復雜，使得該區(qū)地層出現(xiàn)大量的倒轉(zhuǎn)褶皺，并在倒轉(zhuǎn)翼上發(fā)生沖斷層或逆斷層，使得部分侵入巖體與地層傾向存在明顯差異，如圖8中的最南東部的異常。

5 結論

針對重磁異常反演的趨膚效應的問題，采用深度加權函數(shù)約束；針對重磁物性參數(shù)反演的不適定問題，建立了基于光滑約束的正則化目標函數(shù)。利用基于PRP公式的重加權正則化共軛梯度法約束反演，完成了復雜模型的重磁數(shù)據(jù)的反演效果評價，并通過實際資料的應用，論證了算法理論的有效性和可靠性，得出以下結論：

1)使用了基于PRP公式的重加權正則化共軛梯度法，改善了重磁反演的計算效率和計算精度，動態(tài)調(diào)整正則化因子，降低了運算成本，達到自適應正則化因子的目的；

2)根據(jù)模型反演，PRP公式下反演得到的剖面、截面、數(shù)據(jù)擬合曲線和相對擬合差曲線相對于FR公式下的反演結果有相對優(yōu)勢。

3)對比重磁反演結果，磁異常的分辨率略優(yōu)于重力異常，這一方面與磁異常模型數(shù)據(jù)本身具有較高分辨率有關外，還可能與反演中磁異常深度權重的相對較快收斂作用有關。