基于交叉梯度函數的重磁三維聯合反演

徐紫陽，朱自強

（中南大學地球科學與信息物理學院，湖南長沙 410083）

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

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

本文通過引入交叉梯度算子，基于重磁多約束單一反演開展三維聯合反演工作，反演過程使用預條件共軛梯度法進行求解。通過建立不同三維地質理論模型，分別開展單一反演和聯合反演研究工作，并且將單一反演結果和聯合反演結果進行了對比分析。

1 重磁正反演理論

1.1 正則化約束

在重力和磁法勘探過程中，通過由地下密度或者磁性異常體在地表所引起的重力場和磁場的變化來計算獲得地下地質異常密度和磁化率分布情況，從而確定地下地質異常體的賦存情況。

重磁正反演問題可以用下式來表示：

其中：d 為長度為Nd=nx×ny的觀測數據向量，可以是磁法或重力數據，nx×ny分別表示測區范圍內的測線與測點數；m 為長度為Nm=nx×ny×nz的模型參數向量，nz是地下網格劃分層數；G 是反演中靈敏度矩陣或核矩陣。正演問題是已知G 和m 的值，求取觀測數據d的過程；反演問題則是通過地表觀測數據d，推導計算地下模型參數m。

地表觀測數據是有限的離散的數據，且采集過程存在觀測誤差和噪聲干擾，同時在反演過程中，地下場源參數個數遠遠大于觀測數據個數，這導致反演問題的結果具有多解性和非唯一性。為了提高反演精確度和穩定性，本文將Tikhonov 正則化思想引入到反演計算中進化求解：

其中：φ（m）為數據擬合泛函；S（m）為模型穩定泛函；α 代表正則化參數，也稱阻尼參數，它通過調節數據擬合泛函與模型目標函數之間的平衡，以防止數據過度擬合或模型過度平滑。

其中數據擬合函數可表示為：

1.2 預條件共軛梯度法

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

假設觀測數據為b，系數矩陣為A，模型參數為m，初始模型參數為m0。則對于共軛梯度法中的方程ATAm=ATb 可改進為

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

2 基于交叉梯度法重磁聯合反演

2.1 交叉梯度函數基本形式

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

其中：mg和mm為聯合反演中的物性參數。在重磁聯合反演中，mg代表密度參數、mm代表磁化率參數。式（6）為上式的具體展開形式：

實測數據是離散的且地下網格剖分也是離散的。因此，在應用交叉梯度函數時，還需要對其進行離散化處理。本文采用泰勒展開處理，只保留一次項，簡化可得式（7）：

式（7）在零點處展開，其形式可進一步簡化為：

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

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

2.2 聯合反演目標函數構建

多約束條件的單一反演算法以及交叉梯度函數，進行重磁聯合反演研究工作，構建其聯合反演目標函數。

交叉梯度項可以表示為：

綜合可得聯合反演目標函數如下：

其中：Φg為重力反演項；Φm為磁法反演項；β 為重力反演項的加權因子；γ 為磁法反演項的加權因子；μ表示交叉梯度項的加權因子。

3 理論模型試算

本文采用MATLAB 語言開發了重力、磁法聯合反演算法，為了檢驗算法的有效性，通過建立長方體模型進行數值模擬反演研究。

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

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

根據所設置單一模型體，分別進行重力和磁法正演計算，得到異常體在地表所引起的重力異常響應和磁異常響應，用所得重磁異常數據分別開展重磁單一反演工作和重磁聯合反演工作，從而得到地下介質密度和磁化率分布情況。

通過重磁兩種方法的單一反演結果與聯合反演結果進行對比（見圖2），可以看出無論是對于密度反演結果還是磁化率反演結果，聯合反演相比于單一反演，在邊界的分辨能力上都具有較大的提升，對異常體形態大小的描述更加接近真實情況。

圖2 單一模型單獨反演與聯合反演結果對比圖Fig.2 The comparison between separate and joint inversion of the single prism model

4 結論

（1）本文對交叉梯度函數進行討論，基于重磁單一反演，開展了基于預條件共軛梯度法的重磁交叉梯度三維聯合反演工作，該算法收斂速度快，穩定性好，可較好的應用于三維反演計算。

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