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

梁生賢

中國地質調查局成都地質調查中心，成都 610081

0 引言

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

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

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

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

1 反演理論

1.1 正則化反演

重力正演公式可寫為

d=Am。

(1)

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

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

(2)

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

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

1/(zM)β/2。

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

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

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

(3)

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

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

(4)

式中，I為單位矩陣。

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

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

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

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

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

(5)

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

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

(6)

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

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

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

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

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

圖1 理論模型Fig.1 Theory model

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

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

2 分塊矩陣LSQR方法

2.1 阻尼LSQR算法

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

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

(7)

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

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

表2 阻尼LSQR算法

2.2 系數矩陣分塊運算

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

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

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

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

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

(8)

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

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

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

3 反演實例

3.1 理論模型

理論模型及數據密度、網格剖分等見2.2節。我們分別利用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.利用深度加權函數進行模型約束的反演結果；b.利用互相關系數與深度加權進行模型約束的反演結果。圖3 理論模型兩種約束方法的反演結果對比圖Fig.3 Comparison of the inversion results of the synthetic model with two model constraints

3.2 實例

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

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

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

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

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

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

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

4 結果與結論

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

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

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

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