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

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

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

0 引言

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

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

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

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

1 方法原理

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

(1)

(2)

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

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

當異常體響應對X軸方向進行求二階導數時，有：

(3)

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

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

(4)

式中：

公式(4)為計算任意點處二階導數值計算公式。

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

(5)

(6)

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

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

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

(7)

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

(8)

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

(9)

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

θc=cot-1(C-D)

(10)

式中：

(11)

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

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

(12)

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

圖2 算法流程Fig.2 Algorithm flow

2 理論模型計算

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

2.1 正演模擬計算

利用正演計算對理論模型自電響應特征進行分析，為觀測數據的定性以及定量解釋奠定理論基礎。

表1 正演模型的埋深參數

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

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

正演計算中心埋深不同模型的響應并分析，參數見表1。

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

正演計算模型形體參數不同的響應特征，見表2。

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

圖4 不同形體參數模型自電異常響應Fig.4 SP anomaly aroused by different shape factor

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

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

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

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

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

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

2.2 理論模型反演

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

表4 模型反演參數

圖6 自電異常二階導數曲線Fig.6 Second derivative of SP anomalies

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

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

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

圖8 自電異常二階導數Fig.8 Second derivative of SP anomalies

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

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

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

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

3 探測實例

3.1 方法驗算

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

表6 反演結果

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

圖10 反演響應對比Fig.10 The comparison graph of inversion results

3.2 石墨礦區探測數據反演

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

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

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

表7 反演結果

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

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

圖11 反演響應對比Fig.11 The comparison graph of inversion results

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

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

4 結論

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

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

致謝

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