梁生賢

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

0 引言

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

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

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

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

1 反演理論

1.1 正則化反演

重力正演公式可寫為

d=Am。

(1)

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

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

(2)

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

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

1/(zM)β/2。

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

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

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

(3)

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

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

(4)

式中，I為單位矩陣。

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

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

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

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

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

(5)

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

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

(6)

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

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

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

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

注：it為L(zhǎng)SQR算法的迭代次數(shù)，一般而言it<

圖1 理論模型Fig.1 Theory model

a.深度加權(quán)；b.互相關(guān)系數(shù)加權(quán)；c.深度與初始互相關(guān)系數(shù)加權(quán)；d.反演迭代兩次后的深度與互相關(guān)系數(shù)加權(quán)；e.反演迭代三次后的深度與互相關(guān)系數(shù)加權(quán)。黑色線框?yàn)槔碚撃Ｐ?。圖2 理論模型不同加權(quán)結(jié)果Fig.2 Synthetic model of different weighted results

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

2 分塊矩陣LSQR方法

2.1 阻尼LSQR算法

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

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

(7)

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

方程A*×Δm*=Δd關(guān)于阻尼因子λ最小二乘問(wèn)題的LSQR算法見(jiàn)表2[28]。

表2 阻尼LSQR算法

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

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

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

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

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

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

(8)

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

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

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

3 反演實(shí)例

3.1 理論模型

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

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

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

3.2 實(shí)例

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

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

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

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

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

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

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

4 結(jié)果與結(jié)論

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

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

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

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