重力梯度張量種植反演

曹書錦，朱自強，魯光銀

?

重力梯度張量種植反演

曹書錦1, 2, 3，朱自強2，魯光銀2

(1. 湖南科技大學(xué)頁巖氣資源利用湖南重點實驗室，湘潭 411201；2. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083；3. 湖南科技大學(xué)資源環(huán)境與安全工程學(xué)院，湘潭 411201)

針對在重力梯度張量正演中計算耗時過長和核矩陣內(nèi)存消耗過大等制約反演實施的瓶頸問題，在L1范數(shù)的基礎(chǔ)上，引入種植反演，用累加求和分析替換迭代求解，避免計算或存儲反演核矩陣，以減少內(nèi)存占用和加快反演迭代；針對種植反演容易導(dǎo)致相鄰異常源相互侵入的問題，引入一個基于位場水平衰減特性加權(quán)函數(shù)來限制密度吸引子的作用范圍，以期使密度吸引子忽略較遠(yuǎn)的異常源，抑制相鄰異常源相互干擾。反演結(jié)果及分析表明重力及重力梯度張量種植反演所需計算機內(nèi)存小和水平衰減特性加權(quán)函數(shù)能有效的抑制相鄰異常源的 侵入。

種植反演；水平加權(quán)特性函數(shù)；重力梯度張量；L1范數(shù)

對于位場反演而言，正演計算是其重要組成部分，正演計算精確性和計算效率是主要制約反演充分發(fā)揮作用的瓶頸問題[1?4]，尤其在大尺度/海量位場數(shù)據(jù)的反演解釋中，這種瓶頸效應(yīng)更為明顯[5]。針對在重力場中正演計算耗時過長、核矩陣內(nèi)存消耗過大等制約重磁反演的瓶頸問題，很多學(xué)者自正反演兩個方面開展富有成效的研究。在正演研究中，1) 在不改變計算量的前提下，開展高性能計算[1]，但由于不儲存核矩陣，在實際應(yīng)用中，每次反演迭代都涉及到多次正演計算，導(dǎo)致其應(yīng)用于反演計算的效率相對較低，同時其硬件和軟件的準(zhǔn)入門檻較高，導(dǎo)致研究成果不多；2) 采用壓縮技術(shù)，以降低核矩陣的存儲成本，如等效構(gòu)架技術(shù)[3]和小波壓縮技術(shù)[6]等，其在減少核矩陣的所需存儲空間的同時，能極大的提高正演計算速度，但等效構(gòu)架技術(shù)不適用帶地形觀測數(shù)據(jù)，小波壓縮技術(shù)在一定程度上損失了信號精度[1]。在反演研究中，1) 引入預(yù)條件技術(shù)，以優(yōu)化系數(shù)矩陣的條件數(shù)，加快反演迭代的進(jìn)行[7]；2) 引入密度上下界函數(shù)，以加快重構(gòu)密度向密度上下界附近靠攏，加快反演收斂速 度[8?10]；3) 引入新的反演方法，以規(guī)避反演核矩陣的存儲或計算，如逆時偏移技術(shù)[11?12]、小波變換技術(shù)[6]和足印反演[13]等。

基于以上考慮，在L1范數(shù)的基礎(chǔ)上[14]，引入種植反演，以累加求和分析替換最優(yōu)化問題中的迭代求解，減少每次參與反演迭代的物性網(wǎng)格的個數(shù)，避免計算或存儲反演核矩陣，以減少內(nèi)存占用和加快反演迭代；針對在L1范數(shù)下的重力及重力梯度張量種植反演易于導(dǎo)致相鄰異常源相互侵入的問題，通過引入加權(quán)函數(shù)，使密度吸引子或種子忽略較遠(yuǎn)的異常源，減弱相鄰異常源相互侵入的趨勢。反演結(jié)果及分析表明重力及重力梯度張量種植反演所需計算機內(nèi)存小和水平衰減特性加權(quán)函數(shù)能有效的抑制相鄰異常源的侵入。

1 重力梯度張量種植反演

為引入種植反演，首先假設(shè)為實際觀測數(shù)據(jù)，為相對應(yīng)的預(yù)測模型通過正演得到的預(yù)測數(shù)據(jù)，此處構(gòu)建一個殘差向量：

=?

由基于最小二乘估計定義的數(shù)據(jù)泛函可知：

式中：下標(biāo)對應(yīng)于各觀測點；為觀測點個數(shù)。

相比于L2范數(shù)，L1范數(shù)不敏感于含大誤差的數(shù)據(jù)，能提供穩(wěn)健的估計策略[15]。

傳統(tǒng)反演應(yīng)用正則化方法降低地球物理反問題中解的不穩(wěn)定性和非唯一性，而在種植反演中，為使原有病態(tài)問題轉(zhuǎn)化為良態(tài)，引入如下約束[16?18]：

1) 緊致或聚焦，以確保在反演所得密度模型中沒有空洞；

2) 僅以種植反演算法中所提供的密度吸引子為中心逐步選擇最優(yōu)網(wǎng)格，且最優(yōu)網(wǎng)格的密度值與相鄰的密度吸引子密度相同，即=；

3) 對于其他未涉及的網(wǎng)格，強制密度值等于零。

基于以上約束策略，可以構(gòu)建種植反演的目標(biāo) 函數(shù)：

(3)

針對種植反演中密度吸引子的水平位置不易確定及鄰域網(wǎng)格沒有空間范圍的限制，導(dǎo)致相鄰重構(gòu)異常源易于相互侵入等問題，UIEDA等[18]提出利用位場數(shù)據(jù)中的水平形態(tài)特性(“Shape-of-anomaly” data misfit)來確定異常源的水平空間分布，以確保在異常響應(yīng)幅值大的區(qū)域，優(yōu)先更新物性網(wǎng)格，在式(1)的基礎(chǔ)上，數(shù)據(jù)泛函可重新定義為

式中：為縮放因子，對于任意給定預(yù)測數(shù)據(jù)，可通過變換上式計算：

但式(4)主要關(guān)注預(yù)測數(shù)據(jù)和觀測數(shù)據(jù)峰值間的差異，而并沒有考慮異常響應(yīng)幅值較小的區(qū)域如：深部/弱異常源等，同時在重力梯度張量多分量聯(lián)合反演中，雖然可以使用多個縮放因子進(jìn)行聯(lián)合反演，但沒有考慮重力梯度張量分量在深度方向及水平方向的衰減特性，即重力梯度張量各分量的縮放因子趨勢并不一致，當(dāng)g分量處于峰值時，而g分量卻趨于零。因而在對深部/弱等異常源反演時存在一定的缺陷。

圖1 種植反演算法的4個階段

與傳統(tǒng)反演基于迭代更新不同，種植反演算法主要通過選取與密度吸引子相鄰的最優(yōu)網(wǎng)格單元，并更新該網(wǎng)格單元密度，以實現(xiàn)反演迭代。以圖1為例，來闡述種植算法的計算流程。

圖1所示為種植反演算法的4個階段。1) 初始化，首先設(shè)置兩個密度吸引子(黑色的網(wǎng)格單元)，并確定其空間位置及密度值，密度吸引子的空間位置可根據(jù)歐拉反褶積給出，密度值依賴于人工經(jīng)驗選擇；2) 構(gòu)建鄰域網(wǎng)格，確定與兩個密度吸引子相鄰的物性網(wǎng)格單元(淺綠色)；3) 遴選最優(yōu)網(wǎng)格，確保逐步變小和確保P()泛函逐步減小的前提下，確定最優(yōu)網(wǎng)格單元(藍(lán)色)，并賦予該網(wǎng)格與密度吸引子相同的密度值，將其自鄰域網(wǎng)格選集中剔除，將這一過程定義為最優(yōu)網(wǎng)格選擇策略。4) 進(jìn)入下一次迭代，選擇其他密度吸引子，在其相應(yīng)的鄰域網(wǎng)格中，利用最優(yōu)網(wǎng)格選擇策略選擇最優(yōu)網(wǎng)格單元，并賦予該網(wǎng)格與密度吸引子相同的密度值(即紅色網(wǎng)格)。種植算法遍歷所有密度吸引子及所有的網(wǎng)格以滿足反演終止條件。

2 算法驗證

圖2所示為正演模型Ⅰ及其不同方法的反演效果對比。如圖2(a)所示，構(gòu)建含有2個同深度異常體的模型對上述方法原理進(jìn)行對比分析與驗證。2個異常體的長×寬×高均為200 m×200 m×200 m，異常體頂板埋深均為200 m，底板埋深均為400 m，剩余密度均為1.0 g/cm3。將場源空間劃分為20×20× 10=4000個單元格，每個單元格沿、和軸方向的長度均為50 m；觀測點高程為地面上25 m，測點間距為50 m，共有20×20=400個采集數(shù)據(jù)。

反演數(shù)據(jù)為添加3%的高斯白噪聲的正演模擬結(jié)果。在Marquardt反演和Occam反演中，=4，=1，=1，=1，=0.0005，=1.8，如圖2(b)~(c)[24?26]所示；在圖2(d)中，聚焦算子為=diag[(2+2)?1/2]，=1.2×10?6。

在圖2(b)和(c)中兩者效果相當(dāng)，兩者在深度上的分辨率均較低，異常底板分界面非常模糊，難以辨識；而聚焦反演取得了較清晰的邊界，由于聚焦作用的影響，存在過度聚焦的現(xiàn)象(見圖2(d))，即反演得出異常的規(guī)模遠(yuǎn)比實際的要小，而異常幅值要比實際正演模型的大得多。

圖3所示為種植反演結(jié)果及其密度模型切片。此處需要注意的是：在反演結(jié)果中，所有已賦值的網(wǎng)格單元的密度均為0.30 g/cm3。在圖3(a)中，出現(xiàn)部分淺色且低密度網(wǎng)格，表明切片未能切割到該網(wǎng)格，而其密度值僅為切片兩側(cè)網(wǎng)格密度值和的平均值，即 0.15 g/cm3。

由于在種植反演所有網(wǎng)格僅涉及到至多一次正演計算，因而其反演速度較快、及內(nèi)存消耗較小。相比于圖2中多種反演方法相比，本研究中種植反演的重構(gòu)輪廓與正演模型的相當(dāng)，且更為規(guī)整，這表明了本研究中算法的可靠性。

圖2 正演模型Ⅰ及其不同方法的反演效果對比

圖3 重力梯度張量種植反演結(jié)果

3 計算性能分析

對于傳統(tǒng)點元法而言，反演核系數(shù)矩陣所需存儲核函數(shù)的個數(shù)N，為觀測點數(shù)N與物性網(wǎng)格個數(shù)N的乘積；對于種植反演而言，反演核系數(shù)矩陣所需存儲核函數(shù)的個數(shù)是變化的。隨著迭代次數(shù)iter和種子個數(shù)seeds增大，種植反演搜索空間呈現(xiàn)出急劇增大的趨勢，即鄰域網(wǎng)格的個數(shù)。如圖 2所示，對于單一種子的種植反演而言，迭代一次將增加了5個派生網(wǎng)格。在不考慮某些最優(yōu)網(wǎng)格與其它最優(yōu)網(wǎng)格存在相同鄰域網(wǎng)格的情況下，第iter次迭代時，由等差數(shù)列關(guān)系可以計算所需存儲核函數(shù)的個數(shù)：

(6)

因而隨著種子個數(shù)的增加，能有效的降低種植反演的迭代次數(shù)；進(jìn)而，可進(jìn)一步降低核系數(shù)矩陣所需存儲核函數(shù)的個數(shù)。

圖4所示為不同方法的反演參與計算的網(wǎng)格數(shù)的對比圖，曲線1為種植反演的迭代次數(shù)與涉及反演計算網(wǎng)格數(shù)的關(guān)系曲線，曲線2為類似Marquardt 反演、Occam反演和聚焦反演等傳統(tǒng)反演方法的迭代次數(shù)與涉及反演計算網(wǎng)格數(shù)的關(guān)系曲線，曲線3為基于式(6)預(yù)測給出的迭代次數(shù)與涉及反演計算網(wǎng)格數(shù)的關(guān)系曲線。從圖4可以看出，隨著迭代次數(shù)的增加，曲線3能大體預(yù)測出涉及種植反演的計算網(wǎng)格數(shù)但其仍然遠(yuǎn)比實際涉及種植反演實際涉及的計算網(wǎng)格數(shù)要大；然而對比種植反演(曲線1和曲線3)和傳統(tǒng)反演方法(曲線2)，可以發(fā)現(xiàn)傳統(tǒng)反演方法所涉及計算網(wǎng)格數(shù)不隨迭代次數(shù)的變化，且一直在比較高的水平。因而重力及重力梯度張量種植反演所需計算機內(nèi)存小。

圖4 對比不同方法的反演參與計算的網(wǎng)格數(shù)

4 種子密度值對種植反演的影響

針對二進(jìn)制反演中，所設(shè)定密度值對反演結(jié)果密度分布有非常大的影響。在圖1示例的基礎(chǔ)上，不改變原有種子的空間位置，而改變其(1，2)密度值，通過設(shè)置多組模型，研究種子密度值對種植反演結(jié)果的影響。

如圖5所示，共測試了6組模型。由種植反演方法的計算量或反演迭代次數(shù)與反演中所賦予密度的網(wǎng)格數(shù)目以及這些網(wǎng)格的鄰域網(wǎng)格成正比?？梢园l(fā)現(xiàn)種子的密度值能極大的影響種植反演方法重構(gòu)的結(jié)果。如圖5(a)~(d)所示，當(dāng)種子的密度值過低時，易于導(dǎo)致相鄰異常的侵入，導(dǎo)致最終重構(gòu)模型與實際情況差別很大；當(dāng)種子的密度值與實際情況較為一致時(見圖5(e))，種植反演方法重構(gòu)模型與實際情況較為一致；當(dāng)種子的密度值小于實際情況時(見圖5(c)和(d))，易于導(dǎo)致反演重構(gòu)密度分布范圍遠(yuǎn)大于實際情況，這將極大的增加反演迭代次數(shù)，降低反演效率；與之相對的是，當(dāng)種子的密度值遠(yuǎn)大于實際情況時(見圖5(a)和(f))，將使得重構(gòu)模型極度聚焦。

圖5 不同密度值的種子種植反演結(jié)果

5 基于位場水平衰減特性加權(quán)函數(shù)

在種植反演中，密度吸引子的位置決定了重構(gòu)密度模型的空間分布，隨著重構(gòu)模型增大，鄰域網(wǎng)格的數(shù)量也在逐步增大，密度吸引子到鄰域網(wǎng)格間的間距也在增大，使得某一個方向鄰域網(wǎng)格始終被確定為最優(yōu)網(wǎng)格，最終導(dǎo)致相鄰異常源間出現(xiàn)相互侵入的現(xiàn)象。為此，引入一個基于位場水平衰減特性加權(quán)函數(shù)w限制密度吸引子的作用范圍，以期使密度吸引子忽略較遠(yuǎn)的異常源，減弱相鄰異常源相互侵入的趨勢：

式中：s和s為密度吸引子的水平位置，和為觀測點水平位置，L為密度吸引子的作用范圍，為水平衰減因子，對于重力梯度張量數(shù)據(jù)的反演，一般在4~6之間。需要注意的是，w僅作用于種植反演的差殘估計部分的計算，而不直接作用于目標(biāo)函數(shù)。當(dāng)差殘估計過大時，最優(yōu)網(wǎng)格將始終不被賦于密度值，從而使得密度吸引子僅在一定的區(qū)域內(nèi)起作用。

首先，以一孤立異常源為例，分別對比水平加權(quán)和非加權(quán)狀態(tài)下的效果。

圖6所示分別為孤立異常源不帶水平加權(quán)函數(shù)和帶水平加權(quán)函數(shù)種植反演g數(shù)據(jù)的結(jié)果，在圖6(b)中=450平面上二維圖形為水平加權(quán)函數(shù)的權(quán)重分布圖，由w定義可知，w為0~1分布的衰減函數(shù)。通過以上結(jié)果的對比分析，可以發(fā)現(xiàn)帶水平加權(quán)函數(shù)的種植反演在水平方向上更為規(guī)整，而不帶水平的反演結(jié)果在各方向上均較為散亂。

圖6 gzz分量種植反演結(jié)果

對于多異常源而言，水平加權(quán)函數(shù)的作用最為關(guān)鍵的是：權(quán)函數(shù)的加載方式和最優(yōu)網(wǎng)格的判斷方式及多分量聯(lián)合反演的方式，對于前者主要有如下兩種方式：其一、分別加載各單一密度吸引子的水平加權(quán)函數(shù)w，其二、將各密度吸引子的水平加權(quán)函數(shù)w進(jìn)行疊加；對于后者主要對應(yīng)如下兩種方式：其一、以式(1)或式(4)單獨計算各密度吸引子的數(shù)據(jù)泛函和目標(biāo)函數(shù)P()，其二、以式(1)或式(4)計算數(shù)據(jù)泛函和目標(biāo)函數(shù)P()；對于多分量反演而言，基于式(1)的種植反演于分量的組合形式無關(guān)，而基于式(4)的種植反演，則差別相對較大，主要是由于縮放因子的計算方式于分量組合形式相關(guān)。

以埋深相同密度不同的組合模型為例，分析、研究水平加權(quán)函數(shù)的有效加載方式。特構(gòu)建如下所示模型：在地下空間存在埋深為250 m的兩個大小均為200 m×200 m×200 m異常體，剩余密度分別為0.1 g/cm3(綠色線框)和0.3 g/cm3(藍(lán)色線框)，觀測點高程為地面上25 m，觀測點間距為50 m，共有40×40=1600個觀測點。并將含3%的高斯噪聲的正演結(jié)果中作為反演數(shù)據(jù)。

圖7所示分別為各密度吸引子水平加權(quán)函數(shù)疊加及單獨作用的g分量種植反演結(jié)果。由于在同一次種植反演迭代中，遍歷不同密度吸引子的鄰域網(wǎng)格和更新其最優(yōu)網(wǎng)格后的重構(gòu)密度模型的差異極小，僅有一個新的最優(yōu)網(wǎng)格產(chǎn)生，故只顯示迭代終止時的重構(gòu)密度模型及與最后一個密度吸引子相對應(yīng)的水平加權(quán)函數(shù)。通過以上種植反演結(jié)果的對比，兩種權(quán)函數(shù)加載模式的反演效果相似。

圖8所示分別為不同水平加權(quán)函數(shù)加載方式及不同多分量處理方式的種植反演結(jié)果。由于各重力梯度張量分量核函數(shù)幅值差異極大，故在圖8(c)和(d)中各分量在合并為單一分量之前已經(jīng)進(jìn)行了歸一化。通過以上反演結(jié)果的對比分析，可以看出圖8(a)在重構(gòu)密度差異上部最不規(guī)整，導(dǎo)致上述問題的主要原因是相鄰異常源相互干擾；而圖8(b)和(d)在深度方向則遠(yuǎn)超出真實模型，導(dǎo)致上述問題的主要原因是采用式(4)計算多個分量并將各個結(jié)果求和，將最終放大數(shù)據(jù)泛函，但隨著深度的增加，物性網(wǎng)格的響應(yīng)變小，此時深部的物性網(wǎng)格的數(shù)據(jù)泛函變小，因而最優(yōu)網(wǎng)格將逐步沿深度方向增加，并最終造成重構(gòu)密度模型在深度方向上遠(yuǎn)遠(yuǎn)超出真實模型的范圍。

圖7 gzz分量種植反演結(jié)果

圖8 種植反演結(jié)果

6 結(jié)論

1) 在種植反演中，當(dāng)密度吸引子的密度值過低時，易于導(dǎo)致相鄰的異常侵入，導(dǎo)致最終重構(gòu)模型發(fā)生畸變，且與實際情況差別很大；當(dāng)密度吸引子的密度值小于實際情況時，易于導(dǎo)致重構(gòu)密度分布范圍遠(yuǎn)大于實際情況，這將極大的增加反演迭代次數(shù)和降低反演效率；當(dāng)密度吸引子的密度值遠(yuǎn)大于實際情況時，將使得重構(gòu)模型過度聚焦。

2) 當(dāng)密度吸引子相對較淺時，易于導(dǎo)致“趨膚”效應(yīng)的產(chǎn)生；而當(dāng)密度吸引子相對較深時，重構(gòu)密度模型整體下移。

3) 對于基于位場水平衰減特性加權(quán)函數(shù)的加載方式，水平加權(quán)函數(shù)單獨加載的全張量聯(lián)合反演的效果最優(yōu)。

4) 當(dāng)密度吸引子水平位置偏離異常源中心時，種植反演通過最優(yōu)網(wǎng)格選取策略，修正因密度吸引子空間位置偏離異常源的問題，使得新近賦值的網(wǎng)格單元逐步接近于實際異常體中心，但無法改變其前期已經(jīng)錯誤賦值的情況，而這些錯誤賦值的網(wǎng)格單元將最終導(dǎo)致重構(gòu)密度模型出現(xiàn)一定地畸變。故此，在本文中，不就密度吸引子水平位置偏離異常源中心的情況進(jìn)行討論，將另撰文詳述并引入其他方法處理該問題。

REFERENCES

[1] 陳召曦, 孟小紅, 郭良輝, 劉國峰. 基于GPU并行的重力, 重力梯度三維正演快速計算及反演策略[J]. 地球物理學(xué)報, 2012,55(12): 4069?4077. CHEN Zhao-yi, MENG Xiao-hong, GUO Liang-hui, LIU Guo-feng. Three-dimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data based on GPU[J]. Chinese Journal of Geophysics, 2012, 55(12): 4069?4077.

[2] 吳文鸝, 高艷芳, 顧觀文. 起伏地形重磁三維快速正演計算[J]. 物探化探計算技術(shù), 2009,31(3): 179?182. WU Wen-peng, GAO Yang-fang, GU Guan-wen. Gravity and magnetic 3-D fast forward computing with varying Topography[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2009, 31(3): 179?182.

[3] 姚長利, 郝天珧, 管志寧, 張聿文. 重磁遺傳算法三維反演中高速計算及有效存儲方法技術(shù)[J]. 地球物理學(xué)報, 2003,46(2): 2252?2258. YAO Chang-li, HAO Tian-yao, GUAN Zhi-ning, ZHANG Yu-wen. High speed computation and efficient storage in 3D gravity and magnetic inversion based on genetic algorithms[J]. Chinese Journal of Geophysics, 2003, 46(52): 2252?2258.

[4] 曹書錦, 朱自強, 魯光銀, 曾思紅, 郭文斌. 基于單元細(xì)分 H-自適應(yīng)有限元全張量重力梯度正演[J]. 地球物理學(xué)進(jìn)展, 2010,25(3): 1015?1023. CAO Shu-jin, ZHU Zi-qiang, LU Guang-yin, ZENG Si-hong, GUO Wen-bin. Forward modelling of full gravity gradient tensors based H-Adaptive mesh refinement[J]. Progress in Geophysics, 2010, 25(3): 1015?1023.

[5] WANG Jun, MENG Xiao-hong, LI Fang. A computationally efficient scheme for the inversion of large scale potential field data: Application to synthetic and real data[J]. Computers & Geosciences, 2015, 85 (Part A): 102?111.

[6] TALWANI M. Nonlinear inversion of gravity gradients and the GGI gradiometer[J]. Central European Journal of Geosciences, 2011, 3(4): 424?434.

[7] WU L, CHEN L. Fourier forward modeling of vector and tensor gravity fields due to prismatic bodies with variable density contrast[J]. Geophysics, 2016, 81(1): G13?G26.

[8] CHEN Zhao-xi, MENG Xiao-hong, GUO Liang-hui, LIU Guo-feng. GICUDA: A parallel program for 3D correlation imaging of large scale gravity and gravity gradiometry data on graphics processing units with CUDA[J]. Computers & Geosciences, 2012, 46: 119?128.

[9] UMA MARTIN, ZHDANOV M S. Massively parallel regularized 3D inversion of potential fields on CPUs and GPUs[J]. Computers & Geosciences, 2014, 62: 80?87.

[10] WILSON G, ?UMA M, ZHDANOV M S. Massively parallel 3D inversion of gravity and gravity gradiometry data[J]. Preview, 2011, 2011(152): 29?34.

[11] LI Y, OLDENBURG D W. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method[J]. Geophysical Journal International, 2003,152(2): 251?265.

[12] 吳小平, 徐果明. 利用共軛梯度法的電阻率三維反演研究[J]. 地球物理學(xué)報, 2000,43(3): 420?427. WU Xiao-ping, XU Guo-ming. Study on 3D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics, 2000, 43(3): 420?427.

[13] 劉銀萍, 王祝文, 杜曉娟, 劉菁華, 許家姝. 基于Extrapolation Tikhonov正則化算法的重力數(shù)據(jù)三維約束反演[J]. 地球物理學(xué)報, 2013,56(5): 1650?1659. LIU Yin-ping, WANG Zhu-wen, DU Xiao-juan,LIU Jing-hua, XU Jia-shu. 3D constrained inversion of gravity data based on Extrapolation Tikhonov regularization algorithm[J]. Chinese Journal of Geophysics, 2013, 56(5): 1650?1659.

[14] COMMER M, NEWMAN G A. New advances in three-dimensional controlled-source electromagnetic inversion[J]. Geophysical Journal International, 2008,172(2): 513?535.

[15] CARDARELLI E, FISCHANGER F. 2D data modelling by electrical resistivity tomography for complex subsurface geology[J]. Geophysical Prospecting, 2006,54(2): 121?133.

[16] ZHDANOV M S, LIU X, WILSON G A, WAN L. 3D migration for rapid imaging of total-magnetic-intensity data[J]. Geophysics, 2012,77(2): J1?J5.

[17] ZHDANOV M S, LIU X, WILSON G A, WAN L. Potential field migration for rapid imaging of gravity gradiometry data[J]. Geophysical Prospecting, 2011,59: 1052?1071.

[18] COX L H, WILSON G A, ZHDANOV M S. 3D inversion of airborne electromagnetic data using a moving footprint[J]. Exploration Geophysics, 2010,41(4): 250?259.

[19] LORIS I, NOLET G, DAUBECHIES I, Dahlen F. Tomographic inversion using L1-norm regularization of wavelet coefficients[J]. Geophysical Journal International, 2007,170(1): 359?370.

[20] CLAERBOUT J F, MUIR F. Robust modeling with erratic data[J]. Geophysics, 1973,38(5): 826?844.

[21] LAST B, KUBIK K. Compact gravity inversion[J]. Geophysics, 1983,48(6): 713?721.

[22] STOCCO S, GODIO A, SAMBUELLI L. Modelling and compact inversion of magnetic data: A matlab code[J]. Computers and Geosciences, 2009,35(10): 2111?2118.

[23] UIEDA L, BARBOSA V C. Robust 3D gravity gradient inversion by planting anomalous densities[J]. Geophysics, 2012,77(4): G55?G66.

[24] THANASSOULAS C, TSELENTIS G A, DIMITRIADIS K. Gravity inversion of a fault by marquardt’s method[J]. Computers & Geosciences, 1987, 13(4): 399?404.

[25] CONSTABLE S C, PARKER R L, CONSTABLE C G. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics, 1987, 52(3): 289?300.

[26] ZHDANOV M S, ELLIS R, MUKHERJEE S. Three- dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004, 69(4): 925?937.

(編輯 何學(xué)鋒)

Planting inversion of tensor gravity data

CAO Shu-jin1, 2, 3, ZHU Zi-qiang2, LU Guang-yin2

(1. Hunan Provincial Key Laboratory of Shale Gas Resource Utilization,Hunan University of Science and Technology, Xiangtan 411201, China;2. School of Geosciences and Info-physics, Central South University, Changsha 410083, China;3. School of Resource Environment and Safety Engineering,Hunan University of Science and Technology, Xiangtan 411201, China)

Large-scale inversion of gravity gradient tensor data is a time-consuming problem with high demands on computational and physical memory usage. To avoid extraordinary matrix-vector multiplications in each inverse iteration and to speed up the forward of geophysical models, planting inversion is introduced and conjugate gradient iteration replaced by accumulation summary based on L1 norm. The planting inversion easily leads to adjacent anomalies mutually invasive, a horizontal weighted function is proposed to suppress mutual interference between the adjacent anomaly sources. These results of the inversions and analysis results show that planting inversion with horizontal weighted function obtain a meaningful geophysical model. And these methods require little memory and high efficiency.

planting inversion; horizontal weighted function; gravity gradient tensor; L1 norm

Project(E51651) supported by the Doctoral Foundation of Hunan University of Science and Technology; Project(16K031) supported by the Scientific Research Fund of Hunan Provincial Education Department of China; Project(2017JJ3069) supported by Natural Science Foundation of Hunan Province of China; Projects(41174061, 41374120) supported by the National Natural Science Foundation of China

2014-07-13; Accepted date:2017-03-31

ZHU Zi-qiang; Tel: +86-13507319431; E-mail: 13507319431@csu.edu.cn

10.19476/j.ysxb.1004.0609.2017.05.017

1004-0609(2017)-05-0997-09

P631.1

A

湖南科技大學(xué)博士科研啟動基金資助項目(E51651)；湖南省教育廳科研項目(16K031)；湖南省自然科學(xué)基金資助項目(2017JJ3069)；國家自然科學(xué)基金資助項目(41174061，41374120)

2014-07-13；

2017-03-31

朱自強，教授，博士；電話：13507319431；E-mail：13507319431@139.com