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

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

0 引言

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

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

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

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

1 方法原理

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

(1)

(2)

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

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

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

(3)

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

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

(4)

式中：

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

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

(5)

(6)

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

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

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

(7)

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

(8)

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

(9)

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

θc=cot-1(C-D)

(10)

式中：

(11)

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

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

(12)

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

圖2 算法流程Fig.2 Algorithm flow

2 理論模型計(jì)算

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

2.1 正演模擬計(jì)算

利用正演計(jì)算對(duì)理論模型自電響應(yīng)特征進(jìn)行分析，為觀測(cè)數(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

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

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

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

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

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

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

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

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

2.2 理論模型反演

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

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

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

得到二階導(dǎo)數(shù)值后，進(jìn)行非線性方程組的求解，得到固定s下的z_q曲線(圖7)。圖7是在不同步長(zhǎng)因子下，反演所得關(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)曲線左右不對(duì)稱性，可以初步判斷異常體極化軸的傾向，二階導(dǎo)數(shù)極小值所在坐標(biāo)原點(diǎn)的一側(cè)為傾向方向。從圖7得到模型埋深和形體參數(shù)q分別為10 m和1.5。將得到的z、q值帶入公式(11)、公式(12)即可得到、k值分別為45°和-1000 mv.m。計(jì)算結(jié)果與設(shè)置的初始模型參數(shù)一致，初步說(shuō)明該方法是可行的。

對(duì)于理論模型E1進(jìn)行反演后，驗(yàn)證了該方法對(duì)簡(jiǎn)單模型能夠進(jìn)行正確的反演。設(shè)計(jì)原始模型參數(shù)如表5所示，分別得到反演過(guò)程中的二階導(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。計(jì)算結(jié)果與原始模型參數(shù)一致，進(jìn)一步說(shuō)明該方法是正確可行的。

3 探測(cè)實(shí)例

3.1 方法驗(yàn)算

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

表6 反演結(jié)果

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

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

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

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

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

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

表7 反演結(jié)果

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

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

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

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

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

4 結(jié)論

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

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

致謝

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