重磁數(shù)據(jù)的重加權(quán)正則化共軛梯度法約束反演

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

(成都理工大學(xué) a.地球物理學(xué)院，b.教育部地球探測與信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，成都 610059)

0 前言

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

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

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

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

1 正演理論

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

地面某一觀測點(diǎn)的響應(yīng)可根據(jù)地下每個剖分單元的重磁場計(jì)算結(jié)果進(jìn)行累加。剖分單元的重磁異常可采用直立長方體表達(dá)式計(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)為某觀測點(diǎn)地下長方體單元引起的重力異常;ΔT(x,y,z)為某觀測點(diǎn)地下長方體單元引起的磁異常;G為常數(shù);μ0為真空磁導(dǎo)率;σ為密度;M為磁化強(qiáng)度;I為磁傾角;D為磁偏角;(ξ;η,ζ)為剖分單元體積微元的中心坐標(biāo);(ξ1,ξ2,η1,η2,ζ1,ζ2)為地下長方體單元的角點(diǎn)坐標(biāo)位置;r為長方體單元到觀測點(diǎn)距離。

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

Gm=d

(3)

其中：

(4)

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

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

2 約束反演理論

2.1 深度加權(quán)

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

(5)

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

2.2 光滑約束

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

(6)

此時進(jìn)行深度加權(quán)以及光滑約束下的模型約束項(xiàng)為：

(7)

其中

(8)

令：

(9)

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

2.3 構(gòu)建正則化目標(biāo)函數(shù)

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

GZ-1Zm=d

(10)

令：

Gw=GZ-1

(11)

mw=Zm

(12)

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

(13)

(14)

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

(15)

2.4 正則化因子的選取

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

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

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

(16)

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

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

(17)

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

2.5 基于PRP公式的重加權(quán)正則化共軛梯度法

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

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

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

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

(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ù)下的磁反演結(jié)果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，磁化強(qiáng)度為0.8 A/m，向下延伸至200 m。右側(cè)異常體II剩余密度為1g / cm3，磁化強(qiáng)度為1 A / m，向下延伸至500 m，二者磁傾角皆為45°，磁偏角皆為45°，而反演時只考慮垂直磁化情況下的數(shù)據(jù)。將地下三維空間剖分為20×20×12=4 800個單元，其邊長為50 m。共有21×21=441個采集數(shù)據(jù)，數(shù)據(jù)采集點(diǎn)位于z=0處。在觀測數(shù)據(jù)中，加入了5%的高斯噪聲污染(圖3(c),圖3(f))，以測試反演針對實(shí)測數(shù)據(jù)的能力。反演的物性參數(shù)范圍控制為0 g/cm3～1 g / cm3與0 A/m～1 A/m。

迭代100次，重力反演花費(fèi)時間為8.593 s，磁反演花費(fèi)時間為9.135 s，其中因?yàn)榇女惓９綇?fù)雜，計(jì)算核矩陣的時間花費(fèi)較多。

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

圖6 不同迭代次數(shù)下的重力反演結(jié)果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ù)的相對擬合差變化表。當(dāng)?shù)螖?shù)較小時，F(xiàn)R公式下的相對擬合差略微小于PRP公式下的相對擬合差。隨著迭代次數(shù)的增加，PRP公式下的相對擬合差逐漸小于FR公式下的相對擬合差。反映該模型下PRP公式下的反演計(jì)算具有較高的收斂速度和計(jì)算效率。

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

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

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

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

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

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

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

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

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

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

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

4 實(shí)際資料應(yīng)用

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

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

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

反演結(jié)果表明測線下方飛仙關(guān)二段地層存在明顯的局部地質(zhì)體，其磁化強(qiáng)度最高可達(dá)3.5 A/m，這一磁化強(qiáng)度顯然不是該沉積地層的響應(yīng)。但我們注意到該測線的三個主要異常，有兩個異常的傾向于地層傾向角度近似(北西傾向45°～57°)，故初步判斷該異常可能為后期沿斷裂侵入的局部輝綠巖墻。但由于該區(qū)后期構(gòu)造改造復(fù)雜，使得該區(qū)地層出現(xiàn)大量的倒轉(zhuǎn)褶皺，并在倒轉(zhuǎn)翼上發(fā)生沖斷層或逆斷層，使得部分侵入巖體與地層傾向存在明顯差異，如圖8中的最南東部的異常。

5 結(jié)論

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

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

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

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