互相關系數(shù)自約束的重力三維反演與高效求解

梁生賢

中國地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心，成都 610081

0 引言

重力反演就是根據(jù)觀測場值求解場源的信息，是資料定量解釋的重要環(huán)節(jié)之一。物性反演將模型空間離散化為若干個單元，只求解各單元相應的密度值，這種方法易于模擬復雜的地質(zhì)體[1]，逐漸成為重力三維反演的重要方向[2]。本文所討論的就是物性反演。它是線性離散不適定問題，加之三維反演的數(shù)據(jù)量比待求解的模型量更少、系數(shù)矩陣(核矩陣)嚴重病態(tài)性和觀測數(shù)據(jù)受誤差影響等，多解性和不穩(wěn)定性嚴重；因此，如何使得解模型更加符合實際情況是反演關注的主要問題之一。此外，當觀測數(shù)據(jù)量很大時，系數(shù)矩陣大型稠密，反演涉及到大規(guī)模數(shù)據(jù)的存儲與運算，在實際運算中須考慮到這一棘手問題。

在不適定問題求解方面，Tikhonov正則化方法是目前應用最為廣泛的算法之一，它在一定意義上有效減少了解的多解性和不穩(wěn)定性。已有研究表明，在正則化反演中盡可能地利用先驗信息對場源施加約束，可使解模型更加符合實際情況，如：Li等[3-4]在重磁反演中通過加入深度加權函數(shù)來克服位場的“趨膚效應”；Paoletti等[5]提出位場自約束反演概念，指出位場本身包含了大量的潛在信息，利用這些潛在的信息對模型自約束可提高反演結果的可靠性；Sun等[6]利用磁場兩個方向上的水平梯度數(shù)據(jù)計算互相關系數(shù)，得到空間加權函數(shù)，并將空間加權函數(shù)引入正則化反演的模型約束中，顯示這種模型自約束反演結果對場源的邊界刻畫更清晰。

在節(jié)省計算成本方面，主要有并行計算[7-8]、核矩陣壓縮、重構[9-13]以及快速正演計算[1-2,14]等方法，其中并行計算需要依靠計算機硬件設備，核矩陣壓縮則會導致部分信號失真等。對于Tikhonov正則化反演，最終都歸結為大型線性方程組問題的求解。Krylov子空間迭代法是求解大型線性方程組問題的一種有效途徑，在位場反演中應用較多的為共軛梯度法[9-10, 15]，近年來LSQR(最小二乘QR分解)算法[16]也越來越多地被應用于位場數(shù)據(jù)反演中[17-19]；對于非病態(tài)線性方程組兩者等效，而對于病態(tài)線性方程組后者求解更穩(wěn)定。

本文從提升重力反演結果的可靠性和節(jié)省計算成本兩個方面著手。其一，將快速掃描的互相關系數(shù)作為先驗信息，通過處理引入到重力正則化反演的目標函數(shù)中，以提高反演結果的可靠性；其二，利用LSQR算法求解線性方程組問題，并對其進行相應的改進，與快速正演計算方法結合，以節(jié)省計算成本。最后通過理論模型和實際數(shù)據(jù)來展示上述方法的應用效果。

1 反演理論

1.1 正則化反演

重力正演公式可寫為

d=Am。

(1)

式中：m為M階模型向量；d為N階數(shù)據(jù)向量；A為N×M階核矩陣，其元素A(i,j)為第j個模型單元在地表i處的重力響應核函數(shù)，在反演中保持不變。

物性反演的任務是根據(jù)觀測數(shù)據(jù)向量反求模型向量。由于反演是不適定的，Tikhonov正則化目標函數(shù)構制如下：

(2)

式中：Δd為觀測數(shù)據(jù)向量與模型響應向量之間的殘差；Δm為模型修改向量；λ為正則化因子或拉格朗日因子，體現(xiàn)了數(shù)據(jù)擬合與模型約束之間的某種“折衷”；D為模型約束矩陣。此外，為了克服重磁觀測幅值隨場源深度增加而迅速衰減的“趨膚效應”，通常加入一個深度加權函數(shù)，在僅考慮深度方向的情況下可寫為[3-4]

Ddepth=diag1/(z1)β/2,1/(z2)β/2,…,

1/(zM)β/2。

式中：zj為模型單元中心埋深，j∈[1,M]；β取1.5～2.0。

根據(jù)極值條件，求目標函數(shù)E(Δm,λ)關于模型修改向量Δm或ΔmT的偏導數(shù)，并令其等于0，可得模型空間迭代求解公式：

Δm=(ATA+λDTD)-1(ATΔd)。

(3)

Siripunvaraporn等[20]基于Occam反演策略提出數(shù)據(jù)空間算法，在數(shù)據(jù)量N遠小于模型量M的情況下，可大幅度提高計算效率，數(shù)據(jù)空間迭代求解公式為

Δm=(DTD)-1AT[λI+A(DTD)-1AT]-1Δd。

(4)

式中，I為單位矩陣。

對于大規(guī)模重力數(shù)據(jù)三維正則化反演而言(數(shù)據(jù)量可能為幾千甚至上萬，模型量可能達幾十萬甚至百萬)，無論模型空間還是數(shù)據(jù)空間，直接法(奇異值分解)求解都需要耗費巨大的計算量和存儲量(表1)。

1.2 利用互相關系數(shù)與深度加權的約束反演

在反演中，盡可能地利用已知地質(zhì)信息(包括地表地質(zhì)、鉆孔以及地震資料等)對場源施加約束是提高解的可靠性有效而實用的措施，但在有些情況下，地質(zhì)信息并不充足。

互相關系數(shù)[21-22]表征了實測重力異常與模型單元核函數(shù)的相關程度，其絕對值越接近1，表示該單元對重力異常的貢獻可能性越大，且具有計算簡單而快速的特點[23-24]。因而，可根據(jù)互相關系數(shù)的絕對值判斷模型單元的重要性，從而將相對重要的模型單元作為先驗信息對模型進行加權約束。并鑒于利用剩余異常計算所得的互相關系數(shù)成像結果可靠性較高[23-26]，利用每步反演迭代的擬合殘差計算互相關系數(shù)，取絕對值和歸一化后對模型進行約束：

Dω=diag1/ω1,1/ω2,…,1/ωM。

(5)

式中：Dω為互相關系數(shù)約束矩陣；ωj為互相關系數(shù)，其數(shù)學表達式[23-24]為

(j=1,2，…，M)。

(6)

此外，相鄰單元的互相關系數(shù)往往是漸變的，整體具有平緩變化的特點，因而無需再加入模型粗糙度約束矩陣?？紤]到位場“趨膚效應”，令D=Dω·Ddepth。

我們根據(jù)文獻[4]中的理論模型來展示不同算法的加權結果。該理論模型(圖1)由兩個傾斜方向相反的脈狀體異常體組成[4]，其中：向西傾斜的長脈狀體剩余密度為1.0 g/cm3，向東傾斜的短脈狀體剩余密度為0.8 g/cm3。正演數(shù)據(jù)密度：東西方向點距50 m，南北方向點距100 m，共861個數(shù)據(jù)。將數(shù)據(jù)加入4%的隨機噪聲，模型水平方向剖分與數(shù)據(jù)網(wǎng)格一一對應，垂向剖分24層，共20 664個模型單元。

圖2為理論模型計算的單元權重結果。由圖2可見：僅利用深度加權時(圖2a)，水平向各單元權重相同，加權值僅在垂向隨深度增加而變大；僅利用互相關系數(shù)加權時(圖2b)，模型單元權重依賴于互相關系數(shù)，而互相關系數(shù)對場源的準確位置反映較差，求解結果可能會出現(xiàn)偏差；同時利用深度與互相關系數(shù)加權時(圖2c、d、e)，由于一般情況下深度加權值的量級遠大于互相關系數(shù)的量級，因而不至讓反演結果過分地依賴于互相關系數(shù)，模型迭代求解從深部開始(圖2c)，隨著迭代進行，擬合殘差逐漸變小，根據(jù)擬合殘差計算所得的互相關系數(shù)逐漸趨于0附近，各單元權重對深度加權的依賴減小，相應的加權值也總體變小，其形態(tài)逐漸逼近理論模型(圖2d、e)?？梢?，同時利用互相關系數(shù)與深度加權的模型約束無疑使得先驗信息更加豐富，其本質(zhì)上屬于自約束反演[5]的一種，有助于提高反演結果的可靠性。

表1 模型空間、數(shù)據(jù)空間與分塊矩陣LSQR方法計算成本對比

注：it為LSQR算法的迭代次數(shù)，一般而言it<

圖1 理論模型Fig.1 Theory model

a.深度加權；b.互相關系數(shù)加權；c.深度與初始互相關系數(shù)加權；d.反演迭代兩次后的深度與互相關系數(shù)加權；e.反演迭代三次后的深度與互相關系數(shù)加權。黑色線框為理論模型。圖2 理論模型不同加權結果Fig.2 Synthetic model of different weighted results

在正則化反演過程中，拉格朗日因子(正則化因子)體現(xiàn)了模型約束與數(shù)據(jù)擬合之間的平衡：過大的拉格朗日因子往往偏重于模型約束中各單元體的重要性，而忽視了數(shù)據(jù)的擬合程度，會使得迭代收斂緩慢，甚至于出現(xiàn)不收斂的情況；過小的拉格朗日因子則偏重于數(shù)據(jù)擬合，而忽視了模型各單元應有權重的影響，會使得反演出現(xiàn)過擬合現(xiàn)象，無疑增加了解的非唯一性，使得求解結果不夠穩(wěn)定。目前常用的拉格朗日因子求取算法有廣義偏差準則、廣義交叉驗證準則和L曲線法，3種計算方法都包含了大規(guī)模矩陣的多次計算，無疑增加了計算成本。本文拉格朗日因子求取采用在一定步長下逐次遞減的方法[27]，當出現(xiàn)迭代發(fā)散時，拉格朗日因子相應地增大并重新求解。該方法求解過程較為穩(wěn)定，且避免了多次拉格朗日因子搜索所需的額外計算成本。

2 分塊矩陣LSQR方法

2.1 阻尼LSQR算法

由于核矩陣為大型稠密的，直接法求解式(3)或式(4)的計算成本巨大。Krylov子空間迭代法僅出現(xiàn)矩陣與向量的乘積，具有收斂速度快、對計算機內(nèi)存要求低等優(yōu)勢。CGLS(共軛梯度最小二乘法)和LSQR算法同屬于Krylov子空間投影方法，兩者的運算量與存儲量相當，但當執(zhí)行多次迭代或系數(shù)矩陣為病態(tài)時，LSQR的數(shù)值穩(wěn)定性更好[16]。

因此，我們引入LSQR算法求解反問題。令A*=AD-1，Δm*=DΔm；由于D為對角矩陣，有D-1D=I，則式(2)可寫為

(7)

此時式(7)為阻尼最小二乘法，可引入阻尼LSQR算法求解線性方程組問題。在求得Δm*后，根據(jù)Δm=D-1Δm*反求模型修改量。

方程A*×Δm*=Δd關于阻尼因子λ最小二乘問題的LSQR算法見表2[28]。

表2 阻尼LSQR算法

2.2 系數(shù)矩陣分塊運算

由上述阻尼LSQR算法可見，它僅涉及到矩陣與向量的乘積。若將A*按列劃分為k個子矩陣，且向量b(表2中的μi或νi)及子向量bi的階數(shù)與對應的矩陣及子矩陣列數(shù)一致，則有：

(b=[b1,b2,…,bk])，

假設面積數(shù)據(jù)是網(wǎng)格化水平分布的，并且模型剖分單元與數(shù)據(jù)網(wǎng)格采取一一對應的關系，則根據(jù)平移等效性和互換對稱性，在同一層模型各單元之間，核矩陣A的元素具有特定的規(guī)律，即等效幾何格架；實際中只需計算和存儲每一層的第一個單元在所有觀測點處的重力響應核函數(shù)，同一層的其他元素可根據(jù)下式進行索引查找[2]：

A(k,l),(i,j)=

A(k-i+1,l-j+1),(1,1)。

(8)

式中：i、j為任意模型單元在x方向和y方向上的排列號；k、l為任意數(shù)據(jù)觀測點在x方向和y方向上的排列號。

若子矩陣按照模型單元劃分，根據(jù)式(8)可快速獲取某個模型單元在所有觀測點處的重力響應核函數(shù)。大型稠密系數(shù)矩陣不再被完整地表示出來，且由于D為對角矩陣，不會明顯增加額外的計算量，實現(xiàn)分塊矩陣LSQR方法與等效幾何格架技術的結合。

在具體計算成本方面，分塊矩陣LSQR方法的空間復雜度僅為O(2M+2N)，時間復雜度為O(2itMN)，由于一般情況下迭代次數(shù)i遠小于數(shù)據(jù)量N和模型量M，相比直接法求解式(3)、(4)，分塊矩陣LSQR方法節(jié)省了大量存儲空間和計算時間(表1)。假設測區(qū)網(wǎng)格數(shù)據(jù)為100×100的規(guī)模，模型為100×100×50的三維網(wǎng)格，LSQR迭代次數(shù)為5 000次，則分塊矩陣LSQR方法的計算量不到數(shù)據(jù)空間的一半，反演速率至少提高了2 倍；假設數(shù)據(jù)以雙精度存儲，則模型空間求解需1.89 TB存儲量，數(shù)據(jù)空間求解需75 GB存儲量，分塊矩陣LSQR方法則僅需約8 MB存儲量，在普通計算機上就能實現(xiàn)較大規(guī)模的三維反演計算。

3 反演實例

3.1 理論模型

理論模型及數(shù)據(jù)密度、網(wǎng)格剖分等見2.2節(jié)。我們分別利用Ddepth和DωDdepth對模型進行約束反演，附加密度約束范圍為0～1 g/cm3，最終迭代反演的均方誤差分別為0.031、0.112 mGal。

圖3a、b分別為利用Ddepth和DωDdepth進行模型約束的反演。兩種模型約束的反演異常在淺部差別不大，但在深度為400 m以下則差別逐漸變大：利用Ddepth進行模型約束的反演異常呈較正的“v”字型(圖3a)，這與真實模型不一致；而利用DωDdepth進行模型約束的反演異常呈“y”字型(圖3b)，較清晰地反映了真實異常體的基本輪廓。

a.利用深度加權函數(shù)進行模型約束的反演結果；b.利用互相關系數(shù)與深度加權進行模型約束的反演結果。圖3 理論模型兩種約束方法的反演結果對比圖Fig.3 Comparison of the inversion results of the synthetic model with two model constraints

3.2 實例

蘆子園鐵鉛鋅銅多金屬礦位于保山—鎮(zhèn)康地塊的南端，成礦類型為巖漿熱液型，已有地質(zhì)、物探及化探工作均認為該區(qū)成礦作用與隱伏中酸性巖體密切相關[29]。但由于礦區(qū)地表未見巖體出露，且在礦區(qū)南東部出露的印支期木場花崗巖體與成礦作用無明顯關系，因而針對巖體的研究總體較少，譬如巖體的埋深、規(guī)模及其空間形態(tài)等均鮮有研究。區(qū)內(nèi)主要出露寒武系、奧陶系、志留系、泥盆系、石炭系、二疊系的碳酸鹽巖、碎屑巖，石炭系、三疊系的火山巖以及第四系。物性資料(圖4)表明：木廠出露的花崗巖體密度常見值為2.59 g/cm3，明信壩出露石英閃長玢巖密度常見值為2.66 g/cm3，兩者相對于區(qū)內(nèi)廣泛出露的碳酸鹽巖(密度通常在2.73～2.75 g/cm3)圍巖表現(xiàn)為低密度特征。在不同類型巖石中，泥巖、粉砂巖密度最小，其次為中酸性巖體，矽卡巖化、礦石以及火山巖密度最高，碳酸鹽巖密度居中；區(qū)內(nèi)第四系以及新近系、古近系、白堊系和侏羅系規(guī)模有限。綜合密度測試結果可知，重力測量在區(qū)內(nèi)尋找隱伏中酸性巖體具備較好的物性條件。我們以云南地質(zhì)調(diào)查院提供的蘆子園地區(qū)1∶5布格重力數(shù)據(jù)為例，進行三維反演以推測區(qū)內(nèi)隱伏中酸性巖體的空間分布特征。

圖4 云南蘆子園地區(qū)巖石密度測試統(tǒng)計結果Fig.4 Statistical results of rock density in Luziyuan, Yunnan

反演前采用矩形窗口滑動平均法求取剩余重力異常(圖5)。利用網(wǎng)格化的剩余重力異常數(shù)據(jù)進行三維反演，共4 148個數(shù)據(jù)，模型水平方向剖分與數(shù)據(jù)一一對應，垂向剖分34層，共141 032個網(wǎng)格單元。在普通Thinkpad臺式機電腦上(處理器：Intel i5-2400CPU，3.10 GH；內(nèi)存：4 GB)反演迭代3次，共耗時51 h 54 min，最終反演的均方誤差為0.104 mGal。

由于中酸性巖體表現(xiàn)為低密度，為方便起見，我們只提取三維反演結果的低密度體(根據(jù)物性測試統(tǒng)計結果，并結合MT測量與密度成像結果的形態(tài)特征[30-31]，低密度體取剩余密度小于0.1 g/cm3)，并將其套合在地質(zhì)圖上(圖6)。由于區(qū)內(nèi)礫巖、砂巖、粉砂巖等低密度體規(guī)模有限，且表現(xiàn)為低電阻率特征，而中酸性巖體表現(xiàn)為高電阻率特征，結合MT測量結果[30-31]以及其他地質(zhì)、物探、化探[29-32]證據(jù)，推測這種大規(guī)模低密度體主要反映了中酸性巖體。

圖5 云南蘆子園地區(qū)剩余重力異常圖Fig.5 Residual gravity anomaly contour of Luziyuan, Yunnan

1.第四系砂、礫、黏土；2.新近系南林組礫巖、砂巖、粉砂巖；3.古近系勐臘組礫巖、砂巖、泥巖；4.白堊系南新組砂礫巖；5.侏羅系中統(tǒng)粉砂巖、頁巖；6.三疊系粉砂巖、灰?guī)r、玄武巖；7.二疊系灰?guī)r；8.石炭系玄武巖、灰?guī)r；9.泥盆系灰?guī)r；10.志留系栗柴壩組灰?guī)r；11.奧陶系灰?guī)r；12.寒武系上統(tǒng)碳酸鹽巖；13.堿長花崗巖；14.鐵礦；15.鉛鋅礦；16.鉛鋅銅礦；17.錫礦；18.推測中酸性巖體。圖6 云南蘆子園地區(qū)地質(zhì)及重力三維反演結果圖Fig.6 Results of 3D gravity inversion and geology of Luziyuan, Yunnan

由圖6可見：研究區(qū)外圍南東部出露的印支期木場花崗巖體與反演結果的南東部低密度體在平面位置上一致，總體傾向南東，顯示與成礦作用關系不大，這與地質(zhì)、化探認識一致；而研究區(qū)北東的低密度體在總體趨勢上與鎮(zhèn)康復背斜一致(走向北東)，顯示復背斜對巖漿侵位的控制；已知礦床在平面位置上均位于巖體形成的凹入部位或轉(zhuǎn)折部位，為巖漿熱液成礦的最有利地段，其中，蘆子園、天生橋一帶的侵入巖枝隆起中心為天生橋，并向蘆子園礦段側(cè)伏，這與天生橋到蘆子園的“背形”隆起條帶狀矽卡巖型磁鐵礦(鐵礦體自天生橋向蘆子園呈隱伏延伸，埋深逐漸加大)形態(tài)一致。根據(jù)上述推斷的巖體空間分布特征與已知地質(zhì)信息[29-32]，顯示了本算法的實用性。

4 結果與結論

1)將擬合殘差計算所得的互相關系數(shù)作為先驗信息，與深度加權函數(shù)同時引入正則化反演模型約束中，既體現(xiàn)了互相關系數(shù)在模型約束中的作用，又不至讓反演結果過分依賴于互相關系數(shù)。理論模型反演結果表明，這種自約束反演方法對真實異常體的輪廓反映較為清晰，解模型更加符合實際情況。

2)引入阻尼LSQR算法求解反問題，由于一般情況下迭代次數(shù)遠小于數(shù)據(jù)量和模型量，反演速率可提高數(shù)倍。算法中只涉及到矩陣與向量的乘積，便于實現(xiàn)分塊運算，結合等效幾何格架技術，將原矩陣按照模型單元劃分為若干個子矩陣進行存儲與運算，極大地節(jié)省了反演對計算機存儲空間的需求。

3)將本文的重力三維反演方法應用于云南蘆子園隱伏花崗巖體的定位中，共4 148個數(shù)據(jù)，141 032個模型單元，在普通計算機上運算僅需約52 h。反演結果顯示已知礦床均位于推測巖體形成的凹入部位內(nèi)側(cè)，為巖漿熱液成礦的最有利地段，驗證了方法的有實用性。

致謝：王永華教授在程序編寫及論文撰寫期間提供了幫助，在此表示感謝。