陳 輝，陳煜烽，邱姜歆，鄧居智

(東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西 南昌 330013)

重力布格異常包含諸多信息，是地下半空間所有密度不均勻體所引起的重力異常綜合響應(yīng)(曾華霖，2005)。在實(shí)際重力數(shù)據(jù)處理中，人們往往只關(guān)注某些目標(biāo)地質(zhì)體(如斷裂構(gòu)造、基底起伏或成礦巖體等)的重力異常，由此必須從綜合布格重力異常中提取單個(gè)或若干目標(biāo)體引起的重力異常進(jìn)行反演解釋，即異常分離。然而，由于實(shí)際地質(zhì)體間存在錯(cuò)綜復(fù)雜的組合關(guān)系，提取某個(gè)地質(zhì)體或幾個(gè)目標(biāo)體的重力異常相對較困難(管志寧等，2002)。因此，重力異常分離成為當(dāng)前重力數(shù)據(jù)處理關(guān)鍵環(huán)節(jié)和難點(diǎn)問題。

經(jīng)過幾十年來的深入研究，目前具有代表性的重力異常分離方法有多項(xiàng)式擬合(Martín et al.，2011;Ojo et al.，1997)、濾波技術(shù)(Guo et al.，2013;郭良輝等，2012;Keating et al.，2011)、延拓技術(shù)(駱遙，2011;Meng et al.，2009)、趨勢分析(羊春華，2005;Abdelrahman et al.，1991;王彥國等，2017)、小波分解(Singh et al.，2015;Xu et al.，2009;楊文采等，2001)、經(jīng)驗(yàn)?zāi)B(tài)分解(曾琴琴等，2010)、最小曲率法(紀(jì)曉琳等，2015)、主成分分析(Pei et al.，2012;Zhang et al.，2009)和有限元法(Agarwal et al.，2010)等。其中，多項(xiàng)式擬合、濾波技術(shù)和延拓技術(shù)是較為成熟的方法。這些方法能夠利用特定理論模型和實(shí)測數(shù)據(jù)驗(yàn)證其有效性，但需要通過大量試驗(yàn)選擇合理參數(shù)組合，如小波分析中母函數(shù)、濾波技術(shù)中的濾波因子選擇，延拓技術(shù)的延拓高度選擇等。在實(shí)際數(shù)據(jù)處理中，有時(shí)可以采用多種方法在不同參數(shù)條件下進(jìn)行對比試驗(yàn)，并結(jié)合地質(zhì)資料進(jìn)行優(yōu)化組合(Martínez-Moreno et al.，2015;劉彥等，2012;Sedighi et al.，2009)。近年來，隨著重力三維反演技術(shù)逐漸走向成熟，在地質(zhì)和地球物理等相關(guān)資料約束下可獲取地下三維密度分布信息(Lelièvre et al.，2009;Williams，2008;Shin et al.，2006;陳輝等，2015;Li et al.，1998a)。Li等(1998b)提出基于三維反演的磁測數(shù)據(jù)異常分離技術(shù)，通過理論模型和實(shí)際資料驗(yàn)證有效性。該方法能夠比較直觀地提取不同深度的磁異常。

筆者介紹重力三維反演的基本原理，給出利用三維反演技術(shù)實(shí)現(xiàn)重力異常分離的方法和技術(shù)流程;進(jìn)而通過理論模型試驗(yàn)探討方法的有效性;最后將該方法應(yīng)用于江西九瑞礦集區(qū)實(shí)測布格重力異常處理中，提取了不同深度的局部重力異常，并結(jié)合地質(zhì)資料進(jìn)行分析和解釋。

1 理論與方法

1.1 三維重力反演基本原理

筆者采用加拿大哥倫比亞大學(xué)UBC-GIF小組研發(fā)的三維重力反演軟件Grav3D。其基本原理是基于Li等(1998a)提出的反演算法。該反演方法構(gòu)建的目標(biāo)函數(shù)為

φ=φd+μφm

(1)

式中，φd為數(shù)據(jù)目標(biāo)函數(shù)，φm為模型目標(biāo)函數(shù)，μ是平衡模型適用性和數(shù)據(jù)擬合誤差的一個(gè)權(quán)重參數(shù)。數(shù)據(jù)目標(biāo)函數(shù)定義為觀測數(shù)據(jù)與模型正演響應(yīng)權(quán)重差的L-2范數(shù)，即

φd=‖Wd(dfwd-dobs)‖2

(2)

式中，Wd為對角陣，由觀測數(shù)據(jù)標(biāo)準(zhǔn)誤差的倒數(shù)構(gòu)成，dfwd為正演模擬數(shù)據(jù)，dobs為觀測數(shù)據(jù)。模型目標(biāo)函數(shù)φm定義為反演密度模型和先驗(yàn)信息參考模型的權(quán)重差，即

(3)

式中，mref為先驗(yàn)信息參考模型，m為每次迭代的反演模型，ws、wx、wy和wz為權(quán)重函數(shù)。其中，空間權(quán)重函數(shù)ws控制各個(gè)網(wǎng)格單元和參考模型網(wǎng)格單元的物性特征，而權(quán)重函數(shù)wx、wy、wz控制單元模型各個(gè)方向的平滑度，wr(z)是深度或距離的權(quán)重函數(shù)，αs、αx、αy和αz為各方向的相關(guān)性系數(shù)。

為了克服重力數(shù)據(jù)三維反演垂向分辨率差及多解性問題，通常將已知物性、鉆孔和地質(zhì)信息作為約束條件融合到反演中。由式(3)可知，通過預(yù)先設(shè)置參考模型mref、邊界范圍以及權(quán)重函數(shù)三種方式可以實(shí)現(xiàn)約束反演。本研究只采用無約束三維反演。

1.2 基于三維反演的重力異常分離

為了實(shí)現(xiàn)布格重力異常分離，首先對重力數(shù)據(jù)區(qū)域進(jìn)行擴(kuò)邊處理，并利用三維反演得到地下三維密度模型，然后根據(jù)需要提取密度模型中的目標(biāo)區(qū)塊進(jìn)行正演模擬，得到局部重力異常，最后通過與布格重力異常相減得到區(qū)域重力異常。

(1)數(shù)據(jù)準(zhǔn)備。首先將實(shí)測數(shù)據(jù)通過常規(guī)數(shù)據(jù)改正和處理后得到布格重力異?！鱣;為了消除人為噪聲影響、提高三維反演的穩(wěn)定性，將布格重力異常△g向上延拓0.5～1倍點(diǎn)距。另外，三維反演效果也會(huì)受到地形起伏的影響，因此需要足夠密度的地形數(shù)據(jù)加以約束。

(2)數(shù)據(jù)擴(kuò)邊。為了削弱邊界效應(yīng)對三維反演的不利影響，獲得深部密度分布信息，通常需要對反演區(qū)域進(jìn)行擴(kuò)邊，為此研究區(qū)域的觀測布格重力數(shù)據(jù)也需進(jìn)行相應(yīng)的擴(kuò)邊處理。通常利用研究區(qū)域周邊大面積小比例尺的數(shù)據(jù)與觀測布格重力數(shù)據(jù)進(jìn)行融合拼接，或者利用重磁處理軟件(如Oasis Montaj)中的FFT算法對觀測布格重力數(shù)據(jù)進(jìn)行擴(kuò)邊得到。擴(kuò)邊范圍通常需要達(dá)到數(shù)據(jù)觀測區(qū)域的3倍以上。

(3)初始模型建立。合理的初始模型是重力三維反演的關(guān)鍵。初始模型由研究區(qū)域、擴(kuò)邊區(qū)域和延拓區(qū)域三部分組成。其中，研究區(qū)域是為了提取淺部局部異常，需要獲得較為可靠的密度模型，橫向網(wǎng)格大小通常為布格重力數(shù)據(jù)觀測點(diǎn)距或網(wǎng)格大小的0.5～2倍，垂直方向采用隨深度指數(shù)增加的網(wǎng)格，第一層網(wǎng)格厚度一般取橫向網(wǎng)格的0.5～1倍。擴(kuò)邊區(qū)域是為了提取中深部的局部異常，橫向網(wǎng)格尺寸為研究區(qū)域網(wǎng)格大小的1～3倍，而垂向網(wǎng)格尺寸與研究區(qū)域一致。延拓區(qū)域是為了減小反演的邊界效應(yīng)，提高研究區(qū)域和擴(kuò)邊區(qū)域內(nèi)反演密度模型的可靠性，初始橫向網(wǎng)格大小與擴(kuò)邊區(qū)域一致，然后呈指數(shù)遞增，垂向網(wǎng)格與研究區(qū)域一致。另外，初始模型的剩余密度一般取為零。當(dāng)研究區(qū)地形起伏較大時(shí)，初始模型需要加入地形數(shù)據(jù)。

(4)重力三維反演通常采用三維無約束反演。在布格重力異常反演數(shù)據(jù)中需給定數(shù)據(jù)誤差(數(shù)據(jù)不確定性)，通常依據(jù)數(shù)據(jù)質(zhì)量確定。當(dāng)數(shù)據(jù)質(zhì)量難于確定時(shí)，可根據(jù)多次反演試驗(yàn)的結(jié)果進(jìn)行設(shè)定。

(5)局部異常和區(qū)域異常分離。根據(jù)研究目標(biāo)需求在反演三維密度模型中將對應(yīng)區(qū)塊的密度體單獨(dú)提取，并利用三維重力正演軟件進(jìn)行正演模擬，得到局部重力異常數(shù)據(jù)，最后利用布格重力異常數(shù)據(jù)減去局部重力異常數(shù)據(jù)獲得區(qū)域重力異常數(shù)據(jù)。

2 理論模型試驗(yàn)

為了驗(yàn)證本算法的有效性，可設(shè)計(jì)一個(gè)區(qū)域重力異常和局部重力異常疊加密度模型(圖1a)。其中，區(qū)域密度模型為一個(gè)半徑3.0 km的深部球體，頂部埋深為10 km，質(zhì)心位于(0，0)，密度為1.8 g/cm3。局部密度模型為5個(gè)不同大小、位置、埋深和密度的異常體。其中，塊體1和2密度為1.5 g/cm3、塊體3和4為1 g/cm3、塊體5為1.2 g/cm3(表1)。每個(gè)塊體位置和埋深如圖1a所示。重力異常數(shù)據(jù)觀測范圍為7 km×7 km，采樣間隔0.25 km。圖1b給出重力異常疊加分布特征。由圖1可見，5個(gè)異常體所產(chǎn)生的局部重力異?；祀s在埋藏較深的大異常體所產(chǎn)生的區(qū)域重力異常之中，異常體邊界難于識(shí)別。

表1 理論模型參數(shù)Table 1 The parameter of synthetic models

對疊加重力場數(shù)據(jù)進(jìn)行三維反演獲取地下三維密度模型，然后提取研究區(qū)范圍內(nèi)3.5 km以淺密度體進(jìn)行正演模擬得到局部異常，并與疊加異常相減得到區(qū)域重力異常。同時(shí)采用二階趨勢分析和解析延拓(向上延拓3 500 m)進(jìn)行重力異常分離。圖2引用三種分離方法得到的區(qū)域重力異常與區(qū)域密度模型正演響應(yīng)的對比結(jié)果。趨勢分析法(圖2b)分離的區(qū)域異常形態(tài)與正演響應(yīng)(圖2a)基本一致，但最大異常幅值相差最大(0.595 mGal，約為7.44%)。向上延拓分離(圖2c)的區(qū)域重力異常形態(tài)與正演響應(yīng)形態(tài)有一定差別，在近北東方向有一定拉伸，呈現(xiàn)為近橢圓形，最大異常幅值相差0.567 mGal，約為7.09%。利用三維反演方法(圖2d)分離得到的區(qū)域重力異常形態(tài)與正演響應(yīng)形態(tài)基本一致，最大異常幅值相差最小，為0.416 mGal，約為5.20%。

圖3引用三種分離方法得到的局部重力異常及局部密度模型正演響應(yīng)的對比結(jié)果。三種方法提取的局部重力異常形態(tài)(圖3b，c，d)與正演響應(yīng)(圖3a)基本一致，但從細(xì)節(jié)上看，趨勢分析和向上延拓分離出的區(qū)域異常邊界不清晰。從局部重力異常最大幅值來看，趨勢分析相差最大(為0.396 mGal，約為4.44%);向上延拓相差0.377mGal，約為4.43%，而三維反演相差最小(0.294 mGal，約為3.44%)。

綜上所述，無論從形態(tài)還是從分布特征上，趨勢分解、解析延拓和三維反演方法分離的區(qū)域和局部重力異常與正演響應(yīng)均呈現(xiàn)基本一致特征，均取得較好效果。然而，從細(xì)節(jié)上利用重力三維反演分離區(qū)域和局部重力異常與理論值更接近;從數(shù)值上利用三維反演提取局部和區(qū)域重力異常與正演響應(yīng)最大幅值相差最小，具有更高準(zhǔn)確度。由此說明基于三維反演的重力異常分離算法較傳統(tǒng)的異常分離方法更有效。

3 實(shí)例分析

九瑞礦集區(qū)位于揚(yáng)子地塊中部，下?lián)P子斷裂拗陷帶的西緣，北臨大別山隆起，南靠江南褶皺造山帶，是長江中下游銅鐵金硫多金屬成礦帶的重要組成部分。區(qū)內(nèi)發(fā)育“多位一體”(斑巖型、矽卡巖型、塊狀硫化物)的城門山銅多金屬礦床，以矽卡巖型、塊狀硫化物礦床為主的武山銅礦，以及與次火山巖作用有關(guān)的洋雞山金礦等重要金屬礦床(鄧居智等，2015)，如圖4所示。該區(qū)位于九瑞—彭澤復(fù)式向斜的西段，主要褶皺自北向南有:鄧家山—通江嶺向斜、界首—大橋背斜、橫立山—黃橋向斜、大沖—丁家山背斜、烏石街—賽湖向斜、長山—城門湖背斜、新塘向斜。區(qū)內(nèi)斷裂主要為北東向、北西向、北東東向。北東向斷裂屬郯廬—贛江構(gòu)造帶中廬山穹斷束西部邊緣斷裂帶的一部分，對九瑞地區(qū)構(gòu)造格局起著重要的控制作用。該區(qū)蓋層出露的地層比較齊全，從震旦系至第四系，除缺失志留系部分地層外，其余均有出露。區(qū)內(nèi)巖漿活動(dòng)較強(qiáng)烈，小巖體星羅棋布，巖漿巖活動(dòng)有燕山期和喜山期兩個(gè)巖漿旋回。該區(qū)發(fā)現(xiàn)已知礦床近30個(gè)，受成礦地質(zhì)、斷裂構(gòu)造和巖體共同控制。

該區(qū)已開展1∶5萬比例尺重力測量。研究區(qū)位置見圖4，而圖5給出九瑞礦集區(qū)1∶5萬布格重力異常圖。布格重力異常整體呈現(xiàn)由東北向西南逐漸降低，中部出現(xiàn)明顯北東向密集梯度帶特征，與豐山洞-城門山深部斷裂相吻合。該深部斷裂帶對九瑞礦集區(qū)構(gòu)造格局起著重要的控制作用。布格重力異常無法揭示地層和構(gòu)造分布特征，需要進(jìn)行異常分離，提取局部重力異常。為此，對研究區(qū)數(shù)據(jù)通過FFT算法進(jìn)行擴(kuò)邊處理，然后實(shí)現(xiàn)無約束條件帶地形三維反演得到地下三維密度分布，并根據(jù)研究目標(biāo)提取不同深度的重力異常進(jìn)行分析和解釋。

為了進(jìn)一步說明本算法的有效性和可行性，利用三階趨勢分析、低通濾波(波長為70 km)和向上延拓(上延6 km)等傳統(tǒng)方法進(jìn)行異常分離。從圖6可見，4種方法提取的區(qū)域重力異常均表現(xiàn)為西南低，沿東北方向逐漸增大，在東北角呈現(xiàn)為高值的特征。其中，趨勢分析等值線表現(xiàn)為近北西向的弧形，數(shù)值變化激烈，三維反演、低通濾波和向上延拓重力異常等值線表現(xiàn)為近北西方向直線，幅值變化較平滑且基本接近。通過與該區(qū)區(qū)域地質(zhì)及深部構(gòu)造對比發(fā)現(xiàn)，基于三維反演的重力異常分離方法更加符合實(shí)際地質(zhì)情況。

從圖7可見，4種方法提取的區(qū)域重力異常整體呈西南向北東延伸，北西為高值，東南為低值。另外，區(qū)內(nèi)存在三條明顯高異常區(qū)，分別為中南部近東西向分布的大規(guī)模帶狀高異常區(qū)、中部近西北向延伸的高異常帶和東北部近西北向高異常帶，與該區(qū)背斜的分布特征吻合。西北部、東南部及東北部的低異常區(qū)則與向斜分布相對應(yīng)。同時(shí)，已知礦床基本出現(xiàn)在局部重力高值異常的邊界處。三維反應(yīng)提取的局部重力異常幅值變化較小，高值異常邊界清晰。

為進(jìn)一步了解不同深度重力異常分布特征，分別提取了三維密度模型中<0.5 km、0.5～1 km、1～5 km、5～10 km的密度體進(jìn)行正演模擬，得到不同深度的局部重力異常(圖8)。在淺部(圖8a)可以看出研究區(qū)存在東西向和近東西向三條重力高異常帶，分別對應(yīng)大橋-寶山復(fù)式背斜、丁家山-大沖背斜和城門山-長山背斜;四條東西向和近東西向的重力高異常帶分別對應(yīng)于江嶺-鄧家山復(fù)式向斜、黃橋復(fù)式向斜、賽城湖-烏石街復(fù)式向斜和新壙向斜;背斜和向斜交替出現(xiàn)為研究區(qū)地層的主要特征。隨著深度增加(圖8b)高值異常范圍逐漸增大，低值異常區(qū)逐漸變小。在深部(圖8c，d)表現(xiàn)為東北重力異常值高、西南低的特征。由此表明背斜在深部逐漸消失，而向斜逐漸融合為一體。結(jié)合礦床和礦點(diǎn)分布特征，發(fā)現(xiàn)向斜和背斜接觸帶(五通組和黃龍組交接處)為成礦有利區(qū)。

4 結(jié)論

利用三維反演技術(shù)實(shí)現(xiàn)重力異常分離的方法和技術(shù)流程，通過試驗(yàn)研究得出以下結(jié)論和建議:

(1)理論模型和實(shí)測數(shù)據(jù)試驗(yàn)研究表明，基于三維反演的重力異常分離方法能依據(jù)研究目標(biāo)直接提取不同深度的局部重力異常，傳統(tǒng)異常分離方法更加有效，其提取的區(qū)域和局部重力異常更加接近真實(shí)情況，同時(shí)異常邊界刻畫更加明顯，為重力異常分離提供一種新的技術(shù)手段。

(2)利用基于三維反演的重力異常分離方法提取九瑞礦集區(qū)不同深度的局部重力異常，結(jié)合地質(zhì)資料進(jìn)行分析和解釋，獲得很好的應(yīng)用效果，為九瑞礦集區(qū)深部找礦和地質(zhì)認(rèn)識(shí)提供參考和依據(jù)。