基于簡單幾何體擬合的自然電位異常反演

涂 君, 李 論, 周 軍, 李偉林, 王 成， 金 楊

(1. 成都理工大學(xué) 地球勘探與信息技術(shù)教育部重點實驗室,成都 610059;2. 電子科技大學(xué) 資源與環(huán)境學(xué)院,成都 610054)

0 引言

自然電場法(self-potential，SP)是利用巖、礦石由于電化學(xué)作用在其周圍產(chǎn)生的自然極化電場，進行找礦、填圖和解決水文地質(zhì)問題的一種被動源電法勘探方法[1]。該方法近年來被廣泛應(yīng)用于石墨礦、硫化礦等礦產(chǎn)勘探當(dāng)中。早在上世紀(jì)20、30年代，人們就開始著手自然電法的研究[2]。Yüngül[3]奠定了利用簡單形體的響應(yīng)去逼近自電異常的反演解釋方法基礎(chǔ)。對于不均勻體自電異常定量反演解釋方法經(jīng)過幾十年的發(fā)展，已形成很多分支，但一般主要分為兩大類，即復(fù)雜模型的最優(yōu)化反演與簡單形體模型的快速反演。

所謂復(fù)雜模型的最優(yōu)化反演基于嚴(yán)格二、三維正演模擬，此類方法適用于任意復(fù)雜形體的自電異常。由于此類方法中的正演計算引入了有限元或有限差分法等數(shù)值計算方法，導(dǎo)致反演計算的效率相對較低[4]。在當(dāng)前技術(shù)條件下，此類反演只能獲得異常體在地面的電流值，且無法得到模型電阻率的空間分布特征[5-6]。

所謂簡單形體模型的快速反演是基于單獨異常體的響應(yīng)分布，可用簡單形體的響應(yīng)近似代替的原理。該方法的合理性在于勘探獲得的自電異常往往是復(fù)雜的、綜合的總體異常，但該總異常是可以分解為簡單異常的疊加組合[7]，因而實際工作中的自電異?？梢宰鳛楹唵涡误w來反演[8]。

簡單形體模型的快速反演不同于復(fù)雜模型的最優(yōu)化反演，簡單形體模型的快速反演理論方法推導(dǎo)簡單、對計算機性能需要求低、計算速度快。目前一般反演技術(shù)對原始數(shù)據(jù)利用效率不高、方法不穩(wěn)定以及去噪能力不強[8-9]。針對于此，筆者提出一種基于簡單形體的自電異常半自動化反演策略。

1 方法原理

設(shè)異常體的中心埋藏深度為z，極化軸與地面的夾角為θ，圍巖的電阻率為ρ1，異常體的電阻率ρ2，異常體表面最大電位躍變綜合參數(shù)Ur，異常形體參數(shù)q，則沿X軸方向主剖面上的電位表達式為式(1)

(1)

(2)

當(dāng)異常體為球體(3D)、水平無限延伸板狀體(2D)和垂直有限延伸板狀體(3D)時，q值分別取1.5、1和0.5[7]。圖1為球體自電異常示意圖。因此簡單異常體的自電異常響應(yīng)可以根據(jù)公式(2)進行數(shù)值模擬計算求得。

圖1 球體自電異常示意圖Fig.1 A sketch showing cross-sectional view, geometries and parameters of sphere

當(dāng)異常體響應(yīng)對X軸方向進行求二階導(dǎo)數(shù)時，有：

(3)

其中：s、Vxx分別代表求導(dǎo)時的步長因子和二階導(dǎo)數(shù)。取x=0處的Vxx(0)值，可以將式(3)改寫成式(4)。

Vxx(xi,z,s,θ)=A*B

(4)

式中：

公式(4)為計算任意點處二階導(dǎo)數(shù)值計算公式。

取xi=±s處的二階導(dǎo)相加并除以xi=0處的二階導(dǎo)，化簡得到：

(5)

(6)

式中:Vxx和s的意義同上。V(x)為各點的實測數(shù)據(jù)。

至此，得到非線性方程(5)的一般表達式，方程(5)含有未知數(shù)z(異常體中心埋深)、s(二階導(dǎo)數(shù)步長因子)以及q(異常體形體參數(shù))。解方程(5)得到在不同步長因子s下的q、z值；繪制z_q曲線，出現(xiàn)交點的橫、縱坐標(biāo)即為異常體的埋深值和形體參數(shù)值q值。

利用x=0處的V(0)帶入公式(2)可以得到k的表達式：

(7)

帶回公式(2)重新整理得到電位表達式如下：

(8)

式中的zc表示由前面計算得到的異常體埋深。因為z和q值已知，所以公式(8)只含極化軸和地面的夾角θ未知。利用最小二乘法原理可以得到計算夾角θ的表達式為式(9)。

(9)

式中:V(xi)為實測值;V(xi,zc,q,θ)為當(dāng)前模型參數(shù)下理論值。令式(9)取極小值得到θ的表達式為式(10)。

θc=cot-1(C-D)

(10)

式中：

(11)

帶入zc和q值，計算公式(11)即可以得到θ的值θc。

同理將得到的zc、q、θc值設(shè)為固定值，帶入公式(2)再次利用最小二乘法得到電偶極矩k的計算表達式：

(12)

從公式(5)、公式(11)和公式(12)分別計算了異常體的中心埋深z、形體參數(shù)q、極化軸與地面的夾角θ和電偶極矩k。其基本流程見圖2。

圖2 算法流程Fig.2 Algorithm flow

2 理論模型計算

正演計算是已知模型空間求解數(shù)據(jù)空間，而反演是已知數(shù)據(jù)空間，求取模型空間；正演是基礎(chǔ)，反演是目的。筆者選擇一個典型模型按文中方法進行反演，驗證方法的正確性。

2.1 正演模擬計算

利用正演計算對理論模型自電響應(yīng)特征進行分析，為觀測數(shù)據(jù)的定性以及定量解釋奠定理論基礎(chǔ)。

表1 正演模型的埋深參數(shù)

圖3 不同埋深模型自電異常響應(yīng)Fig.3 SP anomaly aroused by different depths

模型Z/mqθ/°k/mV·mB1100.845-1000B210145-1000B3101.545-1000

正演計算中心埋深不同模型的響應(yīng)并分析，參數(shù)見表1。

圖3為模型埋深不同的正演響應(yīng)曲線，在模型其他參數(shù)相同的情況下，異常體的埋深增加使響應(yīng)極小值幅值呈非線性遞減。異常極小值對應(yīng)坐標(biāo)位置基本保持不變，極大值所對應(yīng)的坐標(biāo)位置隨著埋深增加而遠(yuǎn)離原點，并且在遠(yuǎn)離中心的位置，響應(yīng)趨于零。

正演計算模型形體參數(shù)不同的響應(yīng)特征，見表2。

圖4為模型形體參數(shù)不同下的正演響應(yīng)曲線，在模型其他參數(shù)相同的情況下，異常體的形體參數(shù)增加使響應(yīng)極小值幅值呈非線性遞減；異常極小值對應(yīng)坐標(biāo)位置基本保持不變，極大值對應(yīng)坐標(biāo)位置隨著形體參數(shù)增加而靠攏異常體中心位置。整體響應(yīng)可見球體(q=0.5)響應(yīng)的幅值最小。

圖4 不同形體參數(shù)模型自電異常響應(yīng)Fig.4 SP anomaly aroused by different shape factor

模型Z/mqθ/°k/mV·mC1101.50-1000C2101.545-1000C3101.590-1000

圖5 不同傾角模型自電異常響應(yīng)Fig.5 SP anomaly aroused by different inclination angle

正演計算模型極化軸傾角不同的響應(yīng)特征并分析(表3)。

從圖5中可以得到，隨著極化軸傾角的變化，自電異常也發(fā)生這變化[1]。具體表現(xiàn)為：

2)當(dāng)極化體為傾斜極化(0°<θ<90°)時，其電位曲線介于水平極化和垂直極化電位曲線之間。傾斜極化電位曲線的形態(tài)隨傾角的不同而不同：電位曲線的極小值已不在極化體正上方，而是向極化軸傾斜的反方向移動。θ越小，移動的距離越大；零值點情況亦是隨θ的減小向極化軸傾斜相反方向移動。零值點與原點的距離為:x0=z*tgθ。

由此可見，z一定時，θ角越大，零值點偏離原點越遠(yuǎn)；當(dāng)θ=90°時，零值點將在無窮遠(yuǎn)處。因極化軸傾斜，在傾斜一側(cè)出現(xiàn)的電位正值，是極化軸傾斜較小的標(biāo)志，且隨著傾角θ的減小，電位正值逐漸增大。在自然界中，由于水文地質(zhì)條件關(guān)系，一般極化軸近于垂直，故在金屬礦體上常觀測到負(fù)電位，只在地形切割很強的地區(qū)位于陡峭的金屬礦體上，有時能見到顯著的電位正異常。并且從公式(2)可見，電位與電偶極矩k呈正相關(guān)關(guān)系。

2.2 理論模型反演

分別設(shè)計簡單參數(shù)理論模型和一個復(fù)雜參數(shù)理論模型進行驗證。根據(jù)表4中模型參數(shù)，計算正演響應(yīng)，并作為初始實測數(shù)據(jù)，帶入反演流程，計算出的二階導(dǎo)數(shù)曲線如圖6所示。

表4 模型反演參數(shù)

圖6 自電異常二階導(dǎo)數(shù)曲線Fig.6 Second derivative of SP anomalies

得到二階導(dǎo)數(shù)值后，進行非線性方程組的求解，得到固定s下的z_q曲線(圖7)。圖7是在不同步長因子下，反演所得關(guān)于異常體埋深與形體參數(shù)的曲線圖。圖7中橫坐標(biāo)即為異常體形體參數(shù)，縱坐標(biāo)為異常體埋深參數(shù)。

圖7 模型反演z_q曲線Fig.7 The inversion z_q curves

模型Z/mqθ/°k/mV·mF14.50.955.6-1950

圖8 自電異常二階導(dǎo)數(shù)Fig.8 Second derivative of SP anomalies

從二階導(dǎo)曲線左右不對稱性，可以初步判斷異常體極化軸的傾向，二階導(dǎo)數(shù)極小值所在坐標(biāo)原點的一側(cè)為傾向方向。從圖7得到模型埋深和形體參數(shù)q分別為10 m和1.5。將得到的z、q值帶入公式(11)、公式(12)即可得到、k值分別為45°和-1000 mv.m。計算結(jié)果與設(shè)置的初始模型參數(shù)一致，初步說明該方法是可行的。

對于理論模型E1進行反演后，驗證了該方法對簡單模型能夠進行正確的反演。設(shè)計原始模型參數(shù)如表5所示，分別得到反演過程中的二階導(dǎo)數(shù)曲線圖和z_q曲線圖。

圖9 模型反演z_q曲線Fig.9 The inversion z_q curves

圖8為理論模型F1的自電異常二階導(dǎo)數(shù)，圖9為該相應(yīng)模型的反演z-q曲線。從圖9可得到zc、q值分別為4.5 m和0.9；將得到的z、q值帶入公式(11)、公式(12)得到θc、k值分別為55.6°和-1950 mV.m。計算結(jié)果與原始模型參數(shù)一致，進一步說明該方法是正確可行的。

3 探測實例

3.1 方法驗算

通過理論模型的數(shù)值模擬計算，驗證了文中所提方法的可行性。為進一步驗證該方法的實用性，選取與已知工區(qū)進行對比驗證。選取Yüngül在土耳其東南部某銅礦區(qū)采集一條自電剖面，并且Bhattacharya B B[10]對該剖面進行處理，Essa K[8]也處理了該剖面。筆者選取二階導(dǎo)步長因子s為3、5和7，對其進行處理。先計算二階導(dǎo)數(shù)，然后得到z_q曲線交匯圖，從中讀取埋深和形體參數(shù)值分別為=35.903 m和q=1。表明此異常體用水平無限板狀體(q=1)能較好擬合。將計算出來的z、q值帶入公式(11)、公式(12)，計算得到極化軸與地面的傾角為17.8201°，電偶極矩k為-12072.8 mV·m，即得到反演參數(shù)(表6)。

表6 反演結(jié)果

計算反演所得參數(shù)的正演響應(yīng)，并與實際資料和原文獻結(jié)果進行對比分析。由圖10可見，響應(yīng)結(jié)果整體吻合很好，該方法相比較于Essa K方法更能擬合原始剖面，但也存在誤差，其產(chǎn)生原因可能是地下異常體空間展布復(fù)雜引起的。綜合整條剖面，此方法整體較好的擬合了原始剖面，說明方法對實際資料也能適用。

圖10 反演響應(yīng)對比Fig.10 The comparison graph of inversion results

3.2 石墨礦區(qū)探測數(shù)據(jù)反演

測區(qū)位于楊子克拉通北緣，經(jīng)歷了結(jié)晶基底形成，褶皺基底形成，澄江湖大陸裂谷、克拉通盆地演化，內(nèi)陸盆山耦合—推覆構(gòu)造五大演化階段。區(qū)內(nèi)礦產(chǎn)主要有鐵礦、鉀長石礦、霞石鋁礦、石墨等。富含有機質(zhì)陸源碎屑及原始生物沉積形成是碳質(zhì)的主要來源[11]。

石墨礦的碳源主要為有機成因生物碳，石墨礦本身具有低阻高極化特性，而石墨礦的成礦圍巖多為電阻率較高、激化率極低的巖體，因而石墨和其他巖(礦)石相比具有明顯的低阻高極化特性，并且在成礦后具有穩(wěn)定的層位和一定的規(guī)模。石墨礦體上方自然電位值最低，勘查區(qū)巖石與礦石之間的存在電阻率和自然電位電性差異，這為物探工作的開展提供了較為理想的地球物理條件。

選取工區(qū)一條自電剖面數(shù)據(jù)進行反演解釋，該剖面含測點44個，數(shù)據(jù)采集點距40 m，對實測數(shù)據(jù)進行反演前預(yù)處理，得到反演結(jié)果如表7所示，表中角度為按李金銘[1]定義的正方向，即從右到左變化。

表7 反演結(jié)果

計算反演所得參數(shù)的正演響應(yīng)，并與實際測量資料(圖11)、長偏移距瞬變電磁法(LOTEM)反演結(jié)果和實際地質(zhì)鉆進圖進行對比分析(圖12)。

圖11中，黑色點線為經(jīng)過預(yù)處理后的測量數(shù)據(jù)，原始數(shù)據(jù)穩(wěn)定性較好，異常幅值明顯。而圖11中的紅色曲線為利用本文快速反演方法得到的反演結(jié)果，反演結(jié)果曲線更為光滑，整體趨勢與原始數(shù)據(jù)擬合較好。

圖11 反演響應(yīng)對比Fig.11 The comparison graph of inversion results

圖12分別為石墨礦體上的自電異常測量數(shù)據(jù)剖面圖與反演擬合數(shù)據(jù)，利用線源長偏移距瞬變電磁得到的石墨礦體反演響應(yīng)圖，以及利用鉆井資料得到的地質(zhì)剖面成果圖。從圖12中鉆井資料得到，石墨礦體傾角為向右大角度陡立傾斜，而且石墨礦體中心點埋深也在80 m左右，形體基本可以用垂直有限板狀體來擬合。從自電剖面上可以看出，自然電位低阻特征與礦體位置對應(yīng)較好。地面長偏移距瞬變電磁反演結(jié)果，電阻率總體呈“兩高夾一低”的特征。自電反演得到的埋深和傾角以及形體參數(shù)與瞬變電磁反演結(jié)果和鉆進資料吻合較好，所以通過圖12所示的勘探實例，進一步驗證了該方法。

在石墨礦體的平面邊界圈定方面，目前的自然電位方法發(fā)揮了重要作用，未來的石墨礦探測工作中，自然電位方法將占據(jù)重要地位。因此，有必要著力于該技術(shù)的方法機理研究，從數(shù)值模擬、實測資料的快速成像方法等方面入手[12]，進一步優(yōu)化和提升自然電位數(shù)據(jù)的解釋策略，使之可以更好的為石墨礦及其他礦產(chǎn)勘探服務(wù)。

4 結(jié)論

從理論上推導(dǎo)了該方法的可行性，設(shè)計理論模型對方法的正確性進行了驗證，應(yīng)用于實際資料反演解釋中，驗證了其可靠性。因為該方法是利用二階導(dǎo)數(shù)解構(gòu)建非線性方程，所以推導(dǎo)簡單，易于實現(xiàn)。在解傾角和電偶極距時，利用剖面全部數(shù)據(jù)，因而數(shù)據(jù)使用率高，可靠性大。相比較于復(fù)雜模型的最優(yōu)化反演，此方法方便、快速，適于初步快速反演。且在用該半自動反演方法時，加入了人為對剖面的地質(zhì)先驗信息，從而提高了數(shù)據(jù)的可靠性。

圖12 綜合對比圖Fig.12 Comprehensive comparison(a)自電反演響應(yīng)對比圖；(b)LOTEM反演剖面；(c)地質(zhì)剖面圖

致謝

感謝無人機技術(shù)項目組(項目編號：kzw027-jy)對本文的支持，感謝審稿專家對本文的審閱及其所提寶貴的修改意見，感謝物探化探計算技術(shù)期刊以及編輯部的審稿、收錄等幫助。