基于數據空間和稀疏約束的三維重力和重力梯度數據聯合反演

張镕哲， 李桐林， 劉財*， 李福元， 鄧馨卉， 石會彥

1 吉林大學地球探測科學與技術學院， 長春 130026 2 中國地質調查局廣州海洋地質調查局， 廣州 510760 3 長春工程學院， 長春 130021

0 引言

重力探勘具有輕便、快捷、成本低等優點，已經廣泛應用于研究地殼內部結構、探明固體礦產和油氣資源分布等領域(Last and Kubik，1983；Guillen and Menichetti，1984；Li and Oldenburg，1998；Portniaguine and Zhdanov，1999；Nabighian et al.，2005；Silva and Barbosa，2006；Zhdanov，2009；Commer，2011). 隨著重力梯度全張量儀(FTG)的問世，人們不僅可以測量重力數據，還可以測量重力梯度數據. 重力數據是通過測量重力位的垂直一階導數獲得，而重力梯度數據是通過測量重力位在三個方向的二階導數獲得. 重力梯度數據包含多個分量信息，每個梯度分量包含的信息量不同，綜合利用各個分量信息有助于提高地質解釋的準確性(Vasco and Taylor，1991；Zhdanov et al.，2004；Martinez et al.，2013；陳閆等，2014；耿美霞等，2016；Zhang and Li，2019).

反演技術是地球物理勘探最為重要的定量解釋手段之一，它是通過地球物理異常場反推地下場源的空間分布情況，提供目標地質構造的物性和幾何形態等特征信息.在頻率域中重力數據包含更多的低頻信息，對反映深層異常體細節方面具有更高的分辨率，而重力梯度數據包含更多的高頻信息，對反映淺層異常體細節方面具有更高的分辨率.根據兩種數據各自的特點，結合兩種數據共同進行反演計算，勢必會降低反演多解性和提高地質解釋的可信性.Wu等(2013)將重力和重力梯度數據結合起來，通過自適應權值估計物體的質量、方向和距離.Capriotti和Li(2014)針對重力和重力梯度數據核函數隨深度衰減速率不同的情況，在模型平滑約束項中加入權重平衡兩種衰減速率，實現了重力和重力梯度數據的平滑聯合反演.秦朋波和黃大年(2016)采用一種空間深度加權函數，用來克服重力和重力梯度數據的趨膚效應，并通過有限內存BFGS擬牛頓法求解重力和重力梯度數據的聯合反演.Qin等(2016)采用非線性共軛梯度法對重力和重力梯度數據開展聯合反演研究，并將該算法應用到美國墨西哥灣海岸鹽丘調查中.李紅蕾等(2017)提出了一種數據加權矩陣，實現了最小二乘迭代法的重力和重力梯度數據聯合反演，并應用到青藏高原及鄰區巖石圈三維密度結構預測中.高秀鶴等(2019)采用閾值約束的協克里金法對重力和重力梯度數據進行了聯合反演研究.根據上述研究發現，目前重力和重力梯度聯合反演算法仍存在以下幾點問題，首先，采用平滑約束模型矩陣得到的反演結果過于模糊發散，不能很好地恢復出真實模型的尖銳邊界；然后，引入依賴先驗信息的深度加權函數，需要人為經驗來確定參數變量的大小，具有一定的不確定性；其次，采用傳統的最小二乘迭代法需要計算和存儲靈敏度矩陣，尤其在三維大數據情況下，靈敏度矩陣會占用大量的內存空間，而共軛梯度法會出現連續搜索小步長的現象，增加反演的迭代次數，從而降低反演的計算效率.最后，在模型空間下求解三維聯合反演時，求解的線性方程組維度較大，勢必會增加大量的計算時間和內存消耗量，并對計算機的性能提出挑戰.

針對上述問題，本文提出了基于數據空間和稀疏約束的重力和重力梯度數據聯合反演算法.首先，我們從基于幾何格架的重力快速正演算法出發，推導出多分量重力梯度數據的快速正演算法.然后，構建了重力和重力梯度數據聯合反演目標函數，目標函數中包含了重力和重力梯度數據擬合項、L0范數正則化模型約束項，模型約束項中摒棄了依賴先驗信息的深度加權函數，引入了自適應模型積分靈敏度矩陣來克服趨膚效應.其次，將模型計算從模型空間通過恒等式轉換到數據空間下，采用一種改進的共軛梯度算法進行反演求解，可以降低反演求解方程的維度和避免存儲大型的靈敏度矩陣.最后，通過模型試算和實測數據驗證本文提出算法的準確性和有效性.

1 快速正演計算理論

將地下目標空間劃分為一系列大小相同的長方體網格單元，且每個長方體單元的密度均勻, 此時地表重力數據或重力梯度張量數據與地下密度分布可以表示為如下線性關系：

d=A·m,

(1)

其中：d為重力或重力梯度數據；A為重力數據或重力梯度數據的正演核函數，是一個Nd×Nm的矩陣，其中Nd為觀測點個數，Nm為地下模型網格剖分個數；m為網格單元的剩余密度.

重力或重力梯度正演計算需要完成的積分次數為Nd×Nm，當地下網格剖分過大時，正演計算將面臨巨大的挑戰.為了提高正演的計算效率，本文采用姚長利等(2003)提出的幾何格架的快速計算策略，通過分析地下模型單元與觀測點之間的位置關系，可以發現同一層各模型單元與觀測點之間的相對位置關系存在等價性，利用這個等價關系，只對每一層第一個模型單元進行幾何格架核函數計算，得到該單元的核函數矩陣，而該層其他模型單元的核函數矩陣可以通過幾何格架的等效性進行查找獲得.通過上述方法，原本反演過程中需要多次重復進行的正演計算就變成了物性參數與對應的幾何格架之間簡單的乘積運算，每層模型單元只計算第一個模型單元，從而大大的提高了正演的計算效率.具體思想如下：

地下模型均勻剖分，觀測點位于模型單元中心正上方的地表處，如圖1所示.假設地表觀測面上任意一個觀測點p(u,v)(u=1,2,3,…,Nx;v=1,2,3,…,Ny，Nx和Ny分別為x軸和y軸方向上的觀測點個數)與地下第kk層模型單元Qkk(u,v)對應(kk=1,2,3,…,L，L為地下模型的層數).首先需要求取第kk層第一個模型單元Qkk(1,1)在觀測點p(u,v)處的幾何格架核矩陣Akk(1,1,u,v)，然后根據幾何格架的平移等效性、對稱互換性可以得到第kk層其他模型單元Qkk(i,j)(i=1,2,3,…,Nx;j=1,2,3,…,Ny)在觀測點p(u,v)處的幾何格架核函數矩陣Akk(i,j,u,v).

tt=|u-i|+1,pp=|v-j|+1.

(2)

Akk(i,j,u,v)=Akk(1,1,tt,pp).

(3)

為了簡約起見，將公式(3)改寫為一維列向量：

(4)

圖1 三維模型網格剖分圖Fig.1 3D model meshing

上述幾何格架計算策略僅適用于重力的核函數計算，而將上述方法應用到重力梯度多分量時，非對角線重力梯度分量(Vxy,Vzy,Vzx)將得不到正確的核函數矩陣.通過實驗模擬分析，模型單元的位置與觀測點的位置處于特殊關系時，公式(4)將不再成立.我們將幾何格架等效關系式進行修改, 重力梯度非對角線分量的幾何格架關系式可表示為：

對于分量Vzx：

ifu

end

對于分量Vzy:

ifv

end

對于分量Vxy:

ifu

ifu

end

end

2 反演計算理論

2.1 目標函數

為了避免因求解病態反問題所引起的不穩定性和多解性等問題，本文采用Tikhonov和Arsenin(1977)正則化方法求解反演線性方程.首先，我們構建重力和重力梯度數據聯合反演目標函數，其表達式如下：

Φ(m)=Φd(m)+α·Φm(m),

(5)

其中，Φd(m)為數據擬合項，Φm(m)為穩定器構成目標函數中的模型約束項，α為正則化因子.

目標函數中數據擬合項包括兩部分，

(6)

將數據擬合項進行簡化，

Φd(m)=‖WqWd(dobs-Am)‖2,

(7)

(8)

L0范數正則化模型約束項的表達式如下：

(9)

求解公式(9)可以將L0范數最小化近似為迭代重加權L2范數最小化問題(Wolke and Schwetlick，1988).

Φm(m)=‖Wmm‖2,

(10)

(11)

其中，ε是很小的常數，該參數決定了反演結果的尖銳程度.與傳統的平滑反演方法相比，基于稀疏約束的L0范數反演方法可以獲得邊界更清晰、對比度更強的模型.將L0范數作為模型約束項，目標函數表達式可以表示為：

Φ(m)=‖WqWd(dobs-Am)‖2

+α·‖Wm(m-mref)‖2,

(12)

其中，mref為參考模型.利用公式(12)求解的反演結果會出現異常體集中在地表的現象，是由于重力和重力梯度數據的核函數隨深度增加而衰減，因此會引起反演結果的趨膚效應.為了改善趨膚效應的影響，傳統的方法是在模型約束項中加入深度加權函數(Li and Oldenburg，1996；Commer，2011)，但是重力和重力梯度核函數衰減速率不同，采用相同的深度加權函數勢必會得到不準確的反演結果，況且深度加權函數中存在未知的變量，需要人為經驗來定義變量的大小，具有一定的不確定性和不穩定性.針對上述問題，Zhdanov (2002)提出了模型積分靈敏度矩陣，它可以消除不同模型參數對觀測數據的貢獻是不同的，觀測數據對不同模型參數的積分靈敏度也是不同所帶來的影響，使得觀測數據對不同模型參數的積分靈敏度變化一致.目前該方法已經應用在重力梯度三維反演中(陳閆，2014；Capriotti and Li，2014；高秀鶴和黃大年，2017).本文也采用模型積分靈敏度矩陣，用來消除因重力和重力梯度核函數衰減速率不同引起的趨膚效應.因此，仍然是形成一個權重，平衡重力和重力梯度兩個核函數衰減速率不同的問題.

首先分析模型參數mk的微擾對數據的影響.觀測數據與模型參數變化的關系可表示為：

δdi=Aikδmk,

(13)

數據靈敏度對模型參數mk的積分可表示為：

(14)

故模型積分靈敏度矩陣可表示為：

(15)

將式(15)加入到目標函數中，最終得到重力梯度和重力數據聯合反演的目標函數表達式如下：

Φ(m)=‖WqWd(dobs-Am)‖2

+α·‖WmS(m-mref)‖2.

(16)

2.2 目標函數的優化

通常目標函數最優化方法是采用最小二乘法，即目標函數對模型參數的求導為零，得到反演模型結果.但最小二乘法存在一定的缺陷，很容易陷入局部極小，求解不準確，所以一些非線性優化方法經常被使用，例如，牛頓法、梯度法、高斯牛頓法、共軛梯度法等.共軛梯度法對初始模型要求較少，具有存儲量小、收斂速度快等優勢.傳統的共軛梯度法迭代求解過程是在同一個循環內完成，而本文采用內外雙循環的迭代過程進行目標函數的求解，該方法可以減少迭代次數，提高計算效率.

2.2.1 傳統共軛梯度法反演

對目標函數公式(16)求梯度，得到第k次迭代的梯度表達式：

(17)

準備觀測數據dobs，數據加權矩陣Wd，初始模型m0，模型積分靈敏度矩陣S和數據加權項Wq，初始正則化參數α0.

計算目標函數初始梯度g0和RR0；

Whilek

k=k+1；

更新模型參數：mk=mk-1-tk-1dk-1；

計算數據擬合差：

計算目標函數對模型mk的梯度gk和RRk；

end While

輸出mk.

2.2.2 改進的共軛梯度法反演

對目標函數公式(16)求極值，得到第k次迭代的高斯-牛頓法模型改變量表達式：

Δmk=mk-mk-1=(Rk)-1bk,

(18)

準備觀測數據dobs，數據加權矩陣Wd，初始模型m0，模型積分靈敏度矩陣S和數據加權項Wq，初始正則化參數α0.

Whilekσ(外部循環)

計算Rk和bk；

令x0=0,r0=d0=bk；

i=i+1；

更新參數：xi=xi-1+ti-1di-1；

ri=ri-1-ti-1Rkdi-1；

end While

k=k+1；

mk=mk-1+xi；

計算數據擬合差：

end While

輸出mk.

2.3 數據空間共軛梯度法

數據空間法(Siripunvaraporn and Egbert, 2000; Siripunvaraporn et al., 2005)是將模型計算從模型空間通過恒等式轉換到數據空間下進行求解.由于模型個數Nm遠遠大于數據個數Nd，反演計算屬于欠定問題，采用模型空間法進行反演計算需要求解一個Nm×Nm的方程組，而在數據空間中只需要求解一個Nd×Nd的方程組，因此可以有效地降低反演計算的內存占用量和計算時間(Pilkington, 2009; 李澤林, 2015; Zhang and Li, 2019; Zhang et al., 2020).本文將模型空間下的模型改變量表達式(18)轉換為數據空間下的模型改變量表達式：

(19)

對高斯-牛頓法推導出的模型改變量表達式(19)直接求解，雖然可以降低反演計算的維度，但是需要存儲靈敏度矩陣A，在三維情況下勢必也會產生大量的計算內存和時間.本文采用改進的共軛梯度反演方法求解數據空間下的模型改變量，令Rk=

正則化因子的選取與Portniaguine和Zhdanov(2002)采用的方法相同，初始正則化因子α0=0.1. 隨著迭代的增加，模型約束項可能會增加，為了確保目標函數的收斂在全區最小，此時需要修正正則化因子，表達式如下：

(20)

(21)

為了保障聯合反演迭代收斂，我們會在擬合差增大或擬合差變化緩慢時，按比例減小正則化因子.

αk+1=αkq,

(22)

其中，q為縮放因子，設q= 0.6.

3 模型試算

3.1 單獨和聯合反演對比

為了證明重力與重力梯度數據聯合反演優于單獨重力數據反演，我們設計了一個不同高度相同大小的雙異常體組合模型，如圖2所示.異常體的大小都為200 m×200 m×150 m，剩余密度為1 g·cm-3.地下模型空間尺度x,y,z為1000 m×1000 m×600 m，將地下模型剖分為20×20×12個緊密排列的規則立方體單元，每個單元大小為50 m×50 m×50 m，背景剩余密度為0 g·cm-3.地面觀測點個數為20×20，觀測點位于網格中心，正演計算的觀測重力和重力梯度數據如圖2所示.

表1 反演模型和理論模型的均方根誤差Table 1 The root mean square error of each inversion model and the theoretical model

圖2 模型一的理論模型和正演響應等值線圖(a) 重力gz; (b) 重力梯度分量Vxx; (c) 重力梯度分量Vxy; (d) 重力梯度分量Vzx; (e) 重力梯度分量Vzy; (f) 重力梯度分量Vzz.Fig.2 Theoretical model and forward response contour map of model 1(a) Gravity gz; (b) Gravity gradient component Vxx; (c) Gravity gradient component Vxy ; (d) Gravity gradient component Vzx; (e) Gravity gradient component Vzy ; (f) Gravity gradient component Vzz .

圖3 不同觀測數據反演結果對比圖(a) 理論模型; (b) gz反演結果; (c) gz和Vxx聯合反演結果; (d) gz和Vxy聯合反演結果; (e) gz和Vzx聯合反演結果; (f) gz和Vzy聯合反演結果; (g) gz和Vzz聯合反演結果; (h) gz和全張量聯合反演結果; 垂向切片位于y=550 m處.Fig.3 Comparison of different observation data inversion results(a) Theoretical model; (b) gz inversion results; (c) gz and Vxx joint inversion results; (d) gz and Vxy joint inversion results; (e) gz and Vzx joint inversion results; (f) gz and Vzy joint inversion results; (g) gz and Vzz joint inversion results; (h) gz and full tensor joint inversion results; The cross section is located at y=550 m.

3.2 傳統和改進的共軛梯度法反演對比

為了驗證改進的共軛梯度法更適合應用于重力和重力梯度數據的聯合反演，我們建立了一個雙異常體組合模型，如圖4所示.模型中包含了兩個大小相同的異常體，兩個異常體的大小為200 m×200 m×200 m，剩余密度為1 g·cm-3，埋深均為150 m，相距300 m.地下空間網格剖分和觀測方式與模型一完全相同，正演計算的觀測重力和重力梯度數據如圖4所示.

圖4 模型二的理論模型和正演響應等值線圖(a) 重力gz; (b) 重力梯度分量Vxx; (c) 重力梯度分量Vxy; (d) 重力梯度分量Vzx; (e) 重力梯度分量Vzy; (f) 重力梯度分量Vzz.Fig.4 Theoretical model and forward response contour map of model 2(a) Gravity gz; (b) Gravity gradient component Vxx; (c) Gravity gradient component Vxy; (d) Gravity gradient component Vzx; (e) Gravity gradient component Vzy; (f) Gravity gradient component Vzz.

分別對傳統和改進共軛梯度算法進行重力和重力梯度數據聯合反演計算，兩種方法都選擇模型積分靈敏度矩陣，初始模型均采用半空間為0 g·cm-3的剩余密度模型.為了獲得更穩定的反演結果，傳統共軛梯度法反演需要在加權模型參數域下求解，反演經過48次迭代達到擬合差閾值1以下，耗時290 s，獲得的反演結果切片如圖5(b,e,h)所示.而改進共軛梯度法反演經過4次迭代達到擬合差閾值1以下,耗時70 s，獲得的反演結果切片如圖5(c，f，i)所示.可以發現，兩種方法都可以恢復出異常體的中心位置和幾何大小，但是傳統共軛梯度法反演結果邊界更模糊，而改進的共軛梯度法反演結果聚焦程度更高，異常體邊界更清晰.無論是在計算時間方面，還是在反演結果分辨率方面，改進的共軛梯度法反演都具有明顯的優勢，更適合于三維重力和重力梯度聯合反演計算.

圖5 傳統和改進的共軛梯度法反演結果對比圖(a,d,g) 理論模型; (b,e,h) 傳統的共軛梯度法反演結果; (c,f,i) 改進的共軛梯度法反演結果; (a,b,c) z=300 m處的橫向切片; (d,e,f) y=500 m處的垂向切片; (g,h,i) x=250 m處的垂向切片圖.Fig.5 Comparison of traditional and improved conjugate gradient inversion results(a,d,g) Theoretical model; (b,e,h) Traditional conjugate gradient inversion results; (c,f,i) Improved conjugate gradient inversion results; (a,b,c) In the z=300 m depth section; (d,e,f) In the y=500 m cross section; (g,h,i) In the x=250 m cross section.

3.3 模型積分靈敏度矩陣和深度加權矩陣對比

為了驗證模型積分靈敏度矩陣相比于深度加權矩陣，可以更有效的消除因重力和重力梯度核函數衰減速率不同引起的趨膚效應，我們繼續選用模型二進行重力和重力梯度數據聯合反演試算.將公式(16)中的S矩陣選取深度加權矩陣(Li and Oldenburg，1996)和模型積分靈敏度矩陣分別進行平滑和稀疏約束反演.所有反演方法最終的擬合差均收斂到閾值0.1以下，反演結果如圖6所示.圖6(a，e)和圖6(b，f)分別為加入深度加權矩陣后的平滑和稀疏約束反演結果切片圖，可以發現，反演異常體的中心位置與真實異常體存在一定的差異，由于趨膚效應，中心位置有上移的趨勢，并且異常體的物性參數值要遠小于真實值.圖6(c，g)和圖6(d，h)分別為加入模型積分靈敏度矩陣后的平滑和稀疏約束反演結果切片圖，由于模型積分靈敏度矩陣的加入，反演異常體物性參數值得到了明顯地改善，并且異常體中心下移，有效地克服了核函數衰減引起的趨膚效應.模型積分靈敏度矩陣與稀疏約束反演相結合，反演得到的結果無論是幾何形態還是物性參數值方面，都與真實模型基本吻合，尤其在縱向分辨率方面得到了明顯地改善.稀疏約束反演相比于平滑反演能有效地捕捉反演解的小尺度細節，增強稀疏性以保持尖銳的邊界.模擬試算驗證了采用模型積分靈敏度矩陣和稀疏約束算法，可以有效地消除因重力和重力梯度核函數衰減速率不同引起的趨膚效應，將該方法應用到重力和重力梯度數據的聯合反演中可以提高縱向分辨率，降低反演多解性，獲得更加準確的反演結果.

圖6 不同S矩陣的平滑和稀疏約束反演結果的對比圖(a,b,c,d) z=250 m處的橫向切片; (e,f,g,h) y=500 m處的垂向切片; (a,e) 深度加權平滑反演結果; (b,f) 深度加權稀疏約束反演結果; (c,g) 模型積分靈敏度平滑反演結果; (d,h) 模型積分靈敏度稀疏約束反演結果.Fig.6 Comparison of smooth and sparse constraint inversion results of different S matrices(a,b,c,d) In the z=250 m depth section; (e,f,g,h) In the y=500 m cross section; (a,e) The depth weighted smooth inversion results; (b,f) The depth weighted sparse inversion results; (c,g) The model integral sensitivity smooth inversion results; (d,h) The model integral sensitivity sparse inversion results.

3.4 模型空間和數據空間法對比

為了證明本文提出的算法在計算時間和內存占用量方面的優勢，我們首先設計了一個地下網格剖分數量較大的簡單模型，如圖7所示.模型中包含了兩個大小不同的異常體，異常體的大小分別為300 m×200 m×150 m和300 m×500 m×200 m，剩余密度均為1 g·cm-3，頂部埋深分別為150 m和200 m，相距350 m.地下模型空間尺度x,y,z為1500 m×1500 m×1000 m，將地下模型剖分為30×30×20個緊密排列的規則立方體單元，每個單元大小為50 m×50 m×50 m，背景剩余密度為0 g·cm-3.地面觀測點個數為30×30，觀測點位于網格中心，正演計算的觀測重力和重力梯度數據如圖7所示.

圖7 模型三的理論模型和正演響應等值線圖(a) 重力gz; (b) 重力梯度分量Vxx; (c) 重力梯度分量Vxy; (d) 重力梯度分量Vzx; (e) 重力梯度分量Vzy; (f) 重力梯度分量Vzz.Fig.7 Theoretical model and forward response contour map of model 3(a) Gravity gz; (b) Gravity gradient component Vxx; (c) Gravity gradient component Vxy; (d) Gravity gradient component Vzx; (e) Gravity gradient component Vzy; (f) Gravity gradient component Vzz .

本文在改進的共軛梯度法反演基礎上，分別采用模型空間和數據空間法對重力和重力梯度數據進行聯合反演計算(計算機配置：Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz 2.21 GHz, 內存16 GB).所有反演的初始模型均采用均勻半空間剩余密度為0 g·cm-3的模型.模型空間法反演經過6次迭代，最終獲得的擬合差為0.99，反演結果如圖8(b，e，h)所示；數據空間法反演經過6次迭代，最終獲得的擬合差為1.03，反演結果如圖8(c，f，i)所示.可以發現，兩種算法得到的反演結果相似，都可以較好地恢復出異常體的邊界位置、幾何大小和物性參數值，從而證明本文提出的算法具有一定的準確性.在反演計算時間方面，模型空間法耗時2061.85 s，而數據空間法耗時527.25 s，可以說明改進的共軛梯度法結合數據空間法有效地提高了計算效率；在內存消耗方面，模型空間法占用的最大內存約為12 GB，而數據空間法占用的最大內存只需要約0.03 GB.在這個例子中，數據空間法內存最大占用量相比于傳統模型空間法減小了約百分之幾.

圖8 模型空間和數據空間聯合反演結果對比圖(a,d,g) 理論模型; (b,e,h) 模型空間法反演結果; (c,f,i) 數據空間法反演結果; (a,b,c) z=300 m處的橫向切片; (d,e,f) y=750 m處的垂向切片; (g,h,i) x=1150 m處的垂向切片圖.Fig.8 Comparison of model space and data space joint inversion results(a,d,g) Theoretical model; (b,e,h) The model space inversion results; (c,f,i) The data space inversion results; (a,b,c) In the z=300 m depth section; (d,e,f) In the y=750 m cross section; (g,h,i) In the x=1150 m cross section.

為了進一步驗證本文提出算法的準確性、抗噪性和有效性，我們又設計了一個多異常體組合復雜模型，如圖9所示，由四個不同大小的長方體組成，鑲嵌在剩余密度為0 g·cm-3的地下均勻半空間中.地下模型空間大小為2000 m×2000 m×1000 m，將地下剖分為40×40 ×20個緊密排列的直立立方體單元，每個單元大小為50 m×50 m×50 m.地面觀測點個數為40×40，觀測點位于網格中心，加入5%高斯隨機噪聲的理論模型正演響應如圖9所示.異常體的物性值、幾何大小和頂部埋深情況如表2所示.

圖9 模型四的理論模型和正演響應等值線圖(a) 重力gz; (b) 重力梯度分量Vxx; (c) 重力梯度分量Vxy; (d) 重力梯度分量Vzx; (e) 重力梯度分量Vzy; (f) 重力梯度分量Vzz.Fig.9 Theoretical model and forward response contour map of model 4(a) Gravity gz; (b) Gravity gradient component Vxx; (c) Gravity gradient component Vxy; (d) Gravity gradient component Vzx; (e) Gravity gradient component Vzy; (f) Gravity gradient component Vzz.

表2 理論模型異常體屬性情況Table 2 Theoretical model anomaly attributes

由于網格剖分個數的增加，采用模型空間法進行反演計算所需要的內存空間已經超過計算機的承受范圍，而數據空間法反演經過17次迭代，最終獲得的擬合差為0.97，反演結果如圖10所示.可以發現，加入噪聲后本文提出的算法仍然能夠準確地恢復出地下異常體的幾何形態和物性參數值大小，與真實模型基本吻合，表現出一定的準確性和抗噪性.反演計算耗時2206.26 s，占用的最大內存約為0.039 GB.在這個例子中，數據空間法反演內存最大占用量相比于傳統模型空間法減小了約百分之幾，可以證明改進的共軛梯度法結合數據空間法有效地降低了計算內存消耗量，表現出一定的有效性，更適合應用于大數據量的重力和重力梯度數據聯合反演計算.

圖10 數據空間共軛梯度法聯合反演結果圖(a,b,c,d) 理論模型; (e,f,g,h) 聯合反演結果; (a,e) z=300 m處的橫向切片; (b,f) x=600 m處的垂向切片; (c,g) y=550 m處的垂向切片; (d,h) x=1600 m處的垂向切片.Fig.10 The data space conjugate gradient method joint inversion results(a,b,c,d) Theoretical model; (e,f,g,h) Joint inversion results; (a,e) In the z=300 m depth section; (b,f) In the x=600 m cross section; (c,g) In the y=550 m cross section; (d,h) In the x=1600 m cross section.

4 實測數據應用

為了進一步說明本文提出算法的有效性和實用性，我們將反演算法應用于黑龍江省大興安嶺呼中區碧水鎮鉛鋅多金屬礦區的重力和重力梯度數據解釋中.在碧水鎮鉛鋅多金屬礦區內共采集了27條測線，測線間距40 m，每條測線上分布16個觀測點，測點間距40 m，將實測布格剩余重力異常數據網格化，網格化后的重力異常數據等值線如圖11a所示.利用傅里葉變換計算得到布格剩余重力異常的重力梯度數據，網格化后的重力梯度異常數據如圖11b所示.方便起見，實測數據應用中只考慮重力和重力梯度Vzz分量數據進行聯合反演計算.

圖11 實測重力和重力梯度數據等值線圖(a) 重力數據; (b) 重力梯度數據Vzz.Fig.11 Contour plots of measured gravity and gravity gradient data(a) The gravity data; (b)The gravity gradient data Vzz.

觀測區域大小為640 m×1040 m，將反演目標區域的地下空間劃分為16×27×10個緊密排列的直立長方體單元，每個單元的大小為 40 m×40 m×40 m.我們對重力和重力梯度數據進行聯合反演，反演初始模型采用剩余密度為0 g·cm-3的均勻半空間模型.反演經過5次迭代收斂到擬合差閾值以下，反演結果沿測線方向的密度分布切片如圖12所示，不同深度的密度分布切片如圖13所示.我們發現，本文提出的反演方法可以快速地分辨出地下不同密度的分布情況，密度異常分布從南到北逐漸升高，高密度異常區域主要出現在研究區的東北方向，由于鉛鋅多金屬礦床具有較高的密度，因此可以通過反演高密度分布情況有效地圈定多金屬礦床的區域位置和劃分礦體與圍巖的邊界.同時，該反演算法獲得最終模型結果耗時約434.5 s，占用內存約0.028 GB.進一步說明了所提出的算法適用于實測數據處理，可以在普通計算機下快速獲得地下密度分布模型，為精準礦產勘探提供重要依據.

圖12 重力和重力梯度數據聯合反演獲得的密度分布垂向切片圖Fig.12 Cross section of density distribution obtained by joint inversion of gravity and gravity gradient data

圖13 重力和重力梯度數據聯合反演獲得的密度分布橫向切片圖(a) z=100 m處的橫向切片; (b) z=200 m處的橫向切片; (c) z=300 m處的橫向切片.Fig.13 Depth section of density distribution obtained by joint inversion of gravity and gravity gradient data(a) In the z=100 m depth section; (b) In the z=200 m depth section; (c) In the z=300 m depth section.

5 結論

本文將數據空間法和改進的共軛梯度算法結合應用到重力和重力梯度數據處理中，開發出計算速度更快、占用內存更小的高分辨率三維重力和重力梯度數據聯合反演算法.理論模型試算表明相比于傳統的共軛梯度反演算法，改進的算法可以有效地降低反演的迭代次數，提高反演的收斂速度，獲得更接近于真實模型的反演結果；聯合反演采用自適應模型積分靈敏度矩陣，可以有效地消除因重力和重力梯度核函數衰減速率不同引起的趨膚效應，相比于依賴先驗信息的深度加權方法，可以自適應調解矩陣大小，有效提高反演結果的縱向分辨率；稀疏約束反演方法相比于傳統平滑反演方法能有效地捕捉反演解的小尺度細節，增強稀疏性以保持尖銳的邊界；將數據空間法和改進的共軛梯度算法結合，可以更好地降低聯合反演求解方程的維度，避免存儲靈敏度矩陣，有效地提高了計算效率和大幅度的降低內存占用空間.野外實例表明，本文提出的算法可以在普通計算機下快速地獲得地下密度分布模型，有效地圈定多金屬礦床的區域位置和劃分礦體與圍巖的邊界，表現出較強的穩定性和實用性.

致謝感謝匿名評審專家為本文提出的寶貴意見.