徐紫陽(yáng)，朱自強(qiáng)

（中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙 410083）

重磁勘探是一種重要勘探地球物理方法，在礦產(chǎn)資源勘察和地質(zhì)調(diào)查中有著重要的作用[1-3]，三維物性反演是重磁資料處理和解釋中至關(guān)重要的一個(gè)環(huán)節(jié)。但是，由于觀測(cè)數(shù)據(jù)量的局限性以及不可避免的噪聲干擾，使得反演結(jié)果常常是病態(tài)的。為了克服這一問(wèn)題，眾多學(xué)者通過(guò)引入各種約束條件來(lái)降低反演的多解性，使得病態(tài)問(wèn)題趨于良態(tài)。Tikhonov 等[4]提出在反演過(guò)程中通過(guò)加入正則化項(xiàng)來(lái)獲得穩(wěn)定的解；Li 等[5]針對(duì)位場(chǎng)反演過(guò)程中在深部方向上的分辨率低的問(wèn)題，通過(guò)引入深度加權(quán)函數(shù)來(lái)消除這種“趨膚效應(yīng)”；蘭學(xué)毅等[6]通過(guò)收集不同地質(zhì)先驗(yàn)信息進(jìn)行基于先驗(yàn)信息約束的重磁三維反演，梁生賢[7]利用反演迭代過(guò)程中的擬合殘差計(jì)算互相關(guān)系數(shù)，在此基礎(chǔ)上，采用互相關(guān)系數(shù)與深度加權(quán)對(duì)重力反演進(jìn)行約束求解。在進(jìn)行地球物理數(shù)據(jù)處理和解釋的過(guò)程中，不同的地球物理方法所對(duì)應(yīng)的物性參數(shù)不同，誤差和噪聲來(lái)源也不同；相較單一物性反演而言，多種地球物理方法的聯(lián)合反演能夠更好的綜合應(yīng)用不同類型的地球物理數(shù)據(jù)，增強(qiáng)反演穩(wěn)定性和精確度[8-11]。因此，多種勘探地球物理方法的聯(lián)合反演也越來(lái)越受到重視。

交叉梯度約束的結(jié)構(gòu)耦合方法自Gallardo 和Meju[12]提出以來(lái)，聯(lián)合反演應(yīng)用效果良好，迅速得到了廣泛的應(yīng)用。Emilia 和Fregoso 等[13]將交叉梯度函數(shù)引入到重力觀測(cè)數(shù)據(jù)和磁法觀測(cè)數(shù)據(jù)聯(lián)合反演工作中，數(shù)值計(jì)算結(jié)果表明其在橫向和深度分辨率方面得到改進(jìn)。Wang 等[14]于球坐標(biāo)下構(gòu)造交叉梯度算子，實(shí)現(xiàn)了大尺度重磁聯(lián)合反演。修春曉等[15]研究地質(zhì)體的結(jié)構(gòu)特性，實(shí)現(xiàn)重磁聯(lián)合反演，保證了反演結(jié)果與先驗(yàn)構(gòu)造信息具有較強(qiáng)的一致性。侯宇健等[16]利用交叉梯度約束實(shí)現(xiàn)了三維極化率/電阻率的聯(lián)合反演；閆政文等[17]對(duì)多種物性參數(shù)進(jìn)行了交叉梯度聯(lián)合反演工作。

本文通過(guò)引入交叉梯度算子，基于重磁多約束單一反演開(kāi)展三維聯(lián)合反演工作，反演過(guò)程使用預(yù)條件共軛梯度法進(jìn)行求解。通過(guò)建立不同三維地質(zhì)理論模型，分別開(kāi)展單一反演和聯(lián)合反演研究工作，并且將單一反演結(jié)果和聯(lián)合反演結(jié)果進(jìn)行了對(duì)比分析。

1 重磁正反演理論

1.1 正則化約束

在重力和磁法勘探過(guò)程中，通過(guò)由地下密度或者磁性異常體在地表所引起的重力場(chǎng)和磁場(chǎng)的變化來(lái)計(jì)算獲得地下地質(zhì)異常密度和磁化率分布情況，從而確定地下地質(zhì)異常體的賦存情況。

重磁正反演問(wèn)題可以用下式來(lái)表示：

其中：d 為長(zhǎng)度為Nd=nx×ny的觀測(cè)數(shù)據(jù)向量，可以是磁法或重力數(shù)據(jù)，nx×ny分別表示測(cè)區(qū)范圍內(nèi)的測(cè)線與測(cè)點(diǎn)數(shù)；m 為長(zhǎng)度為Nm=nx×ny×nz的模型參數(shù)向量，nz是地下網(wǎng)格劃分層數(shù)；G 是反演中靈敏度矩陣或核矩陣。正演問(wèn)題是已知G 和m 的值，求取觀測(cè)數(shù)據(jù)d的過(guò)程；反演問(wèn)題則是通過(guò)地表觀測(cè)數(shù)據(jù)d，推導(dǎo)計(jì)算地下模型參數(shù)m。

地表觀測(cè)數(shù)據(jù)是有限的離散的數(shù)據(jù)，且采集過(guò)程存在觀測(cè)誤差和噪聲干擾，同時(shí)在反演過(guò)程中，地下場(chǎng)源參數(shù)個(gè)數(shù)遠(yuǎn)遠(yuǎn)大于觀測(cè)數(shù)據(jù)個(gè)數(shù)，這導(dǎo)致反演問(wèn)題的結(jié)果具有多解性和非唯一性。為了提高反演精確度和穩(wěn)定性，本文將Tikhonov 正則化思想引入到反演計(jì)算中進(jìn)化求解：

其中：φ（m）為數(shù)據(jù)擬合泛函；S（m）為模型穩(wěn)定泛函；α 代表正則化參數(shù)，也稱阻尼參數(shù)，它通過(guò)調(diào)節(jié)數(shù)據(jù)擬合泛函與模型目標(biāo)函數(shù)之間的平衡，以防止數(shù)據(jù)過(guò)度擬合或模型過(guò)度平滑。

其中數(shù)據(jù)擬合函數(shù)可表示為：

1.2 預(yù)條件共軛梯度法

共軛梯度法的基本思想是把共軛性與最速下降法相結(jié)合，利用已知點(diǎn)處的梯度構(gòu)造一組共軛方向，沿著這個(gè)方向去搜索目標(biāo)函數(shù)極小值點(diǎn)。根據(jù)共軛方向的性質(zhì)，共軛梯度法具有二次終止性。但是，由于方程組本身是欠定的，模型參數(shù)規(guī)模較大，導(dǎo)致雅可比矩陣的條件數(shù)較大、求解過(guò)程不穩(wěn)定。本文通過(guò)引入預(yù)條件因子，對(duì)迭代法進(jìn)行改性，以改善反演收斂次數(shù)多，反演不穩(wěn)定等問(wèn)題。

假設(shè)觀測(cè)數(shù)據(jù)為b，系數(shù)矩陣為A，模型參數(shù)為m，初始模型參數(shù)為m0。則對(duì)于共軛梯度法中的方程ATAm=ATb 可改進(jìn)為

其中：S 為預(yù)條件因子，在理想條件下，S 是ATA逆的近似值，因此SATA 約等于單位矩陣。與原始對(duì)稱正定矩陣ATA 相比，SATA 的特征值更加聚類。因此，減少了條件數(shù)，并加快了CG 方法的收斂速度。但是，以這種方式求解S 計(jì)算量較大，求解較困難。在實(shí)際應(yīng)用中，用式（8）中的深度加權(quán)函數(shù)作為預(yù)條件因子進(jìn)行計(jì)算。

2 基于交叉梯度法重磁聯(lián)合反演

2.1 交叉梯度函數(shù)基本形式

交叉梯度函數(shù)最早由Gallardo 和Meju[12]提出，因其應(yīng)用限制較小，且應(yīng)用效果良好，迅速在地球物理學(xué)中得到廣泛應(yīng)用。兩種物性參數(shù)的三維形式交叉梯度函數(shù)如下：

其中：mg和mm為聯(lián)合反演中的物性參數(shù)。在重磁聯(lián)合反演中，mg代表密度參數(shù)、mm代表磁化率參數(shù)。式（6）為上式的具體展開(kāi)形式：

實(shí)測(cè)數(shù)據(jù)是離散的且地下網(wǎng)格剖分也是離散的。因此，在應(yīng)用交叉梯度函數(shù)時(shí)，還需要對(duì)其進(jìn)行離散化處理。本文采用泰勒展開(kāi)處理，只保留一次項(xiàng)，簡(jiǎn)化可得式（7）：

式（7）在零點(diǎn)處展開(kāi)，其形式可進(jìn)一步簡(jiǎn)化為：

以中心差分代替微分，Bg可表示為如下形式：

Bm的表示式也可由式（8）形式得到。

2.2 聯(lián)合反演目標(biāo)函數(shù)構(gòu)建

多約束條件的單一反演算法以及交叉梯度函數(shù)，進(jìn)行重磁聯(lián)合反演研究工作，構(gòu)建其聯(lián)合反演目標(biāo)函數(shù)。

交叉梯度項(xiàng)可以表示為：

綜合可得聯(lián)合反演目標(biāo)函數(shù)如下：

其中：Φg為重力反演項(xiàng)；Φm為磁法反演項(xiàng)；β 為重力反演項(xiàng)的加權(quán)因子；γ 為磁法反演項(xiàng)的加權(quán)因子；μ表示交叉梯度項(xiàng)的加權(quán)因子。

3 理論模型試算

本文采用MATLAB 語(yǔ)言開(kāi)發(fā)了重力、磁法聯(lián)合反演算法，為了檢驗(yàn)算法的有效性，通過(guò)建立長(zhǎng)方體模型進(jìn)行數(shù)值模擬反演研究。

設(shè)置一個(gè)理論的地下三維物性模型（見(jiàn)圖1），地下場(chǎng)源空間大小為200 m×200 m×100 m，將地下場(chǎng)源空間劃分為獨(dú)立的nx×ny×nz=20×20×20=8 000 個(gè)模型單元體，每個(gè)單元格規(guī)模為10 m×10 m×5 m。設(shè)置一個(gè)棱柱體異常模型，其中心坐標(biāo)為（0，0，30）。異常體x 軸方向，即東向規(guī)模為80 m，y 軸方向，即北向長(zhǎng)為80 m，棱柱體頂部埋深為20 m，底部埋深為40 m。地下物性異常體與圍巖間的剩余密度為ρ=1 g/cm3，地磁場(chǎng)傾角50°，磁偏角0°，無(wú)剩磁存在，物性異常體與圍巖間的相對(duì)磁化率設(shè)置為κ=1 SI。

圖1 單一棱柱體模型示意圖Fig.1 The single-prism model

根據(jù)所設(shè)置單一模型體，分別進(jìn)行重力和磁法正演計(jì)算，得到異常體在地表所引起的重力異常響應(yīng)和磁異常響應(yīng)，用所得重磁異常數(shù)據(jù)分別開(kāi)展重磁單一反演工作和重磁聯(lián)合反演工作，從而得到地下介質(zhì)密度和磁化率分布情況。

通過(guò)重磁兩種方法的單一反演結(jié)果與聯(lián)合反演結(jié)果進(jìn)行對(duì)比（見(jiàn)圖2），可以看出無(wú)論是對(duì)于密度反演結(jié)果還是磁化率反演結(jié)果，聯(lián)合反演相比于單一反演，在邊界的分辨能力上都具有較大的提升，對(duì)異常體形態(tài)大小的描述更加接近真實(shí)情況。

圖2 單一模型單獨(dú)反演與聯(lián)合反演結(jié)果對(duì)比圖Fig.2 The comparison between separate and joint inversion of the single prism model

4 結(jié)論

（1）本文對(duì)交叉梯度函數(shù)進(jìn)行討論，基于重磁單一反演，開(kāi)展了基于預(yù)條件共軛梯度法的重磁交叉梯度三維聯(lián)合反演工作，該算法收斂速度快，穩(wěn)定性好，可較好的應(yīng)用于三維反演計(jì)算。

（2）模型反演結(jié)果表明，交叉梯度聯(lián)合反演算法較單一反演而言，可以挖掘不同物性信息之間的聯(lián)系，更為準(zhǔn)確的刻畫出地下異常體位置信息和規(guī)模，降低反演問(wèn)題的多解性，使反演結(jié)果更加收斂。