本發(fā)明涉及地球物理勘探，尤其是涉及一種基于DBIM的瞬變電磁電導(dǎo)率反演方法。

背景技術(shù)：

瞬變電磁法(Transient Electromagnetic Methods)又稱時間域電磁方法(Time Domain Electromagnetic Methods)，簡稱TDEM或TEM。航空電磁(AEM,Airborne Electromagnetics)法是在地面電磁法基礎(chǔ)上發(fā)展起來的一種空中測量的電磁法，被廣泛用來進(jìn)行資源勘探，并且其非常適合探測一些復(fù)雜地形和人員難以達(dá)到的地區(qū)。其中半航空瞬變電磁(Semi-Transient Electromagnetic Methods)是為了克服全航空瞬變電磁的探測深度有限所發(fā)展出的新型航空電磁法，其是在地面鋪設(shè)幾公里長度的電性源，然后利用無人機(jī)或者直升機(jī)在空中接收磁場，通過分析接收到磁場的信息，來反演出地下介質(zhì)的情況。

航空電磁法是基于巖石電性及磁性差異，利用電磁感應(yīng)原理，以固定翼飛機(jī)或直升機(jī)等飛行器作為運(yùn)載工具，實時地球物理探測的勘探方法，其方法具有高效、經(jīng)濟(jì)、適應(yīng)性強(qiáng)等特點(diǎn)，能夠廣泛應(yīng)用于地面大面積的礦產(chǎn)普查、水工環(huán)境普查、詳查和精細(xì)測量等工作。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的在于提供一種在瞬變電磁系統(tǒng)中能準(zhǔn)確反演地層電導(dǎo)率，在頻率域中直接進(jìn)行反演，不需要進(jìn)行頻時轉(zhuǎn)換，在數(shù)據(jù)處理過程比較簡單，只需要準(zhǔn)確提取相應(yīng)頻率的頻譜即可，不僅適用于半航空瞬變電磁系統(tǒng)，同樣適用于全航空瞬變電磁系統(tǒng)的基于DBIM的瞬變電磁電導(dǎo)率反演方法。

本發(fā)明包括以下步驟：

1)讀取觀測數(shù)據(jù)，具體步驟如下：

將實際接收機(jī)接收數(shù)據(jù)作為輸入數(shù)據(jù)。

在步驟1)中，所述輸入數(shù)據(jù)包括電壓數(shù)據(jù)或磁場數(shù)據(jù)。

2)建立初始模型，具體步驟如下：

根據(jù)已有的信息建立初始模型，初步定義下列初始參數(shù)：

ε＝(ε₁,ε₂...ε_n)；σ＝(σ₁,σ₂...σ_n)；h＝(h₁,h₂...h_n)

其中，ε：為地下分層模型的介電常數(shù)矩陣，ε₁,ε₂...ε_n：為地下第一層、第二層…第n層模型的介電常數(shù)；σ：為地下分層模型的電導(dǎo)率矩陣，σ₁,σ₂...σ_n：為地下第一層、第二層…第n層模型的電導(dǎo)率；h：為地下分層模型的電導(dǎo)率矩陣，h₁,h₂...h_n：為地下第一層、第二層…第n層模型的厚度；

在步驟2)中，所述初始模型包括分層信息、每層的厚度、每層介質(zhì)電導(dǎo)率、介電常數(shù)等信息。

3)更新模型參數(shù)，具體步驟如下：

為迭代過程中更新模型的參數(shù)過程，每一次迭代的前一次都會得到一個模型的更新量，這里是δσ＝(δσ₁,δσ₂...δσ_n)，其中δσ為模型電導(dǎo)率更新矩陣，δσ₁,δσ₂...δσ_n為地下第一層、第二層…第n層的電導(dǎo)率更新量；得到電導(dǎo)率更新量之后，有：

σ^k+1＝σ^k+δσ^k+v

其中，σ^k+1為第k+1次迭代所需要的電導(dǎo)率參數(shù)矩陣；σ^k為第k次迭代過程中所需要的迭代矩陣；δσ^k+1為第k次迭代過程中所得到的模型電導(dǎo)率更新矩陣，其中第一次迭代，模型的更新量為0；

4)計算模型場值，具體步驟如下：

得到模型后，對模型的場值進(jìn)行計算，采用頻率域的計算方法：

E＝<G^EJ；J>

H＝<G^HJ；J>

其中，E和H表示電場和磁場；G^EJ和G^HJ為電流源電場并矢Green函數(shù)和磁場Green函數(shù)，<；>為點(diǎn)乘積分和；

5)計算誤差，具體步驟如下：

得到模型的場值后，對實際測量的場值與模型計算的場值做差，得到兩者之間的誤差；

6)計算Frechet導(dǎo)數(shù)，具體步驟如下：

Frechet導(dǎo)數(shù)矩陣由以下方程求出：

其中為對應(yīng)于地下介質(zhì)所產(chǎn)生的等效源，為由源J'所產(chǎn)生的場值的變化量，得到去譜域的解后對其進(jìn)行Fourier逆變換得到其頻域的解。如下為在y方向源的作用下所產(chǎn)生的磁場的x分量：

δH_x＝δH_x1+δH_x2+δH_x3

δH_x為磁場的x方向分量，δH_x1、δH_x2、δH_x3分別為其分量的第一項、第二項與第三項，具體表達(dá)式如下：

上式中J_y表示y方向的源的強(qiáng)度；y_-l與y_l表示沿y方向源的起點(diǎn)與終點(diǎn)；y為接收點(diǎn)的y坐標(biāo)；x'表示源的x坐標(biāo)，x、y表示接收點(diǎn)的x和y坐標(biāo)；表示為沿著第n層的上表面積分到n層的下表面，即第n層厚度的積分；i_rip(r)V_rip(r)表示位于z'的1V串聯(lián)脈沖電流源(或1A并聯(lián)脈沖電壓源)在深度z處的電流和電壓，下標(biāo)p＝e,h分別表示橫磁波與橫電波；J₀、J₁表示第0類Bessel函數(shù)和第1類Bessel函數(shù)；其中，k_x和k_y分別為x方向和y方向的波數(shù)；δσ_n為地下第n層的電導(dǎo)率參數(shù)變化量。

對于上述的場值變化δH_x＝F·δσ，其中F即為地層參數(shù)的Frechet導(dǎo)數(shù)。

7)計算更新量，具體步驟如下：

對于上述地層模型，建立如下成本函數(shù)，

其中，C表示成本函數(shù)(Cost Function)；δf表示場值的誤差；F為Frechet導(dǎo)數(shù)矩陣；δσ^k+1為第k次迭代所得到的k+1次的電導(dǎo)率更新量；f^obs為測量到的場值；σ^k為第k次迭代所得到的電導(dǎo)率矩陣；γ²表示正則化系數(shù)；||||²表示矩陣的二范數(shù)。

為了使成本函數(shù)取得最小值，等效于解如下方程：

其中，為F矩陣的轉(zhuǎn)置，若其為復(fù)數(shù)矩陣，則表示為共軛轉(zhuǎn)置。

由上式可以得到每一次迭代過程中電導(dǎo)率的更新量δσ^k+1。

8)判斷收斂條件，具體步驟如下：

若未達(dá)到收斂條件，則返回步驟3)更新數(shù)據(jù)的參數(shù)；若達(dá)到收斂條件，則結(jié)束反演過程。

本發(fā)明應(yīng)用于瞬變電磁快速反演，是利用瞬變電磁系統(tǒng)接收到的信號的頻譜信息在頻率域基于變形波恩迭代方法(DBIM Distorted Born iterative methods)的一種快速反演方法。

在瞬變電磁系統(tǒng)中，發(fā)射信號為雙極性方波或者是半正弦信號，同時利用接收線圈接收一次場或者是二次場的時域信息，本發(fā)明可以使用任意發(fā)射信號，只要準(zhǔn)確記錄發(fā)射信號和接收信號，就可以準(zhǔn)確的求得對應(yīng)發(fā)射信號下的理論接收波形；對于時域的接收信號對其進(jìn)行簡單的去直流，去噪等操作后做FFT運(yùn)算，得到其對應(yīng)頻率的頻譜信息，就可以用來進(jìn)行反演運(yùn)算。

本發(fā)明不僅適用于半航空瞬變電磁系統(tǒng)，而且適用于全航空瞬變電磁系統(tǒng)。本發(fā)明從瞬變電磁的頻譜信息出發(fā)，只要提取接收信號的準(zhǔn)確頻譜信息，就可以用本發(fā)明進(jìn)行反演。本發(fā)明基于DBIM方法建立迭代反演過程，最后反演的結(jié)果能很好與實際數(shù)據(jù)相吻合，從而可以大大提高瞬變電磁系統(tǒng)的計算速度和反演精度。

附圖說明

圖1為本發(fā)明的系統(tǒng)反演流程圖。

圖2為瞬變電磁地下模型分層圖。

圖3為本發(fā)明理論模型5層的反演效果圖。

圖4為本發(fā)明理論模型10層的反演效果圖。

圖5為本發(fā)明實際數(shù)據(jù)效果圖。

具體實施方式

以下將結(jié)合附圖對本發(fā)明的技術(shù)方案及其突出效果作進(jìn)一步說明。

參見圖1，本發(fā)明的具體實施步驟如下：

(1)讀取觀測數(shù)據(jù)

將實際接收機(jī)接收數(shù)據(jù)作為輸入數(shù)據(jù)，這些數(shù)據(jù)可以是電壓數(shù)據(jù)，也可以是接收到的磁場數(shù)據(jù)。

(2)建立初始模型

可根據(jù)已有的信息建立一個簡單的初始模型，初始模型需要包括分層信息，每層的厚度，每層介質(zhì)電導(dǎo)率、介電常數(shù)等信息，參見圖2，在圖2中，Z₁、Z₂、…Z_n-1為層界面的高度，線源可采用電流源。

初步定義下列初始參數(shù)：

ε＝(ε₁,ε₂...ε_n)；σ＝(σ₁,σ₂...σ_n)；h＝(h₁,h₂...h_n)

ε：為地下分層模型的介電常數(shù)矩陣，ε₁,ε₂...ε_n：為地下第一層、第二層…第n層模型的介電常數(shù)；σ：為地下分層模型的電導(dǎo)率矩陣，σ₁,σ₂...σ_n：為地下第一層、第二層…第n層模型的電導(dǎo)率；h：為地下分層模型的電導(dǎo)率矩陣，h₁,h₂...h_n：為地下第一層、第二層…第n層模型的厚度。

(3)更新模型參數(shù)

為迭代過程中更新模型的參數(shù)過程，每一次迭代的前一次都會得到一個模型的更新量，這里是δσ＝(δσ₁,δσ₂...δσ_n)，其中δσ為模型電導(dǎo)率更新矩陣，δσ₁,δσ₂...δσ_n為地下第一層、第二層…第n層的電導(dǎo)率更新量。得到電導(dǎo)率更新量之后，有：

σ^k+1＝σ^k+δσ^k+1

其中σ^k+1為第k+1次迭代所需要的電導(dǎo)率參數(shù)矩陣；σ^k為第k次迭代過程中所需要的迭代矩陣；δσ^k+1為第k次迭代過程中所得到的模型電導(dǎo)率更新矩陣。其中第一次迭代，模型的更新量為0。

(4)模型場值計算

得到模型后，就要對模型的場值進(jìn)行計算，這里采用頻率域的計算方法：

E＝<G^EJ；J>

H＝<G^HJ；J>

這里E和H表示電場和磁場；G^EJ和G^HJ為電流源電場并矢Green函數(shù)和磁場Green函數(shù)，<；>為點(diǎn)乘積分和。

(5)誤差計算

得到模型的場值后，需要對實際測量的場值與模型計算的場值做差，得到兩者之間的誤差。

(6)Frechet導(dǎo)數(shù)計算

其中Frechet導(dǎo)數(shù)矩陣可由以下方程求出：

其中為對應(yīng)于地下介質(zhì)所產(chǎn)生的等效源，為由源J'所產(chǎn)生的場值的變化量，得到去譜域的解后對其進(jìn)行Fourier逆變換得到其頻域的解。如下為在y方向源的作用下所產(chǎn)生的磁場的x分量：

δH_x＝δH_x1+δH_x2+δH_x3

δH_x為磁場的x方向分量，δH_x1、δH_x2、δH_x3分別為其分量的第一項、第二項與第三項，具體表達(dá)式如下：

上式中J_y表示y方向的源的強(qiáng)度；y_-l與y_l表示沿y方向源的起點(diǎn)與終點(diǎn)；y為接收點(diǎn)的y坐標(biāo)；x'表示源的x坐標(biāo)，x、y表示接收點(diǎn)的x和y坐標(biāo)；表示為沿著第n層的上表面積分到n層的下表面，即第n層厚度的積分；i_rip(r)V_rip(r)表示位于z'的1V串聯(lián)脈沖電流源(或1A并聯(lián)脈沖電壓源)在深度z處的電流和電壓，下標(biāo)p＝e,h分別表示橫磁波與橫電波；J₀、J₁表示第0類Bessel函數(shù)和第1類Bessel函數(shù)；其中，k_x和k_y分別為x方向和y方向的波數(shù)；δσ_n為地下第n層的電導(dǎo)率參數(shù)變化量。

對于上述的場值變化δH_x＝F·δσ，其中F即為地層參數(shù)的Frechet導(dǎo)數(shù)。

(7)更新量計算

對于上述地層模型，建立如下成本函數(shù)，

其中C表示成本函數(shù)(Cost Function)；δf表示場值的誤差；F為Frechet導(dǎo)數(shù)矩陣；δσ^k+1為第k次迭代所得到的k+1次的電導(dǎo)率更新量；f^obs為測量到的場值；σ^k為第k次迭代所得到的電導(dǎo)率矩陣；γ²表示正則化系數(shù)；||||²表示矩陣的二范數(shù)。

為了使成本函數(shù)取得最小值，等效于解如下的方程

其中：為F矩陣的轉(zhuǎn)置，若其為復(fù)數(shù)矩陣，則表示為共軛轉(zhuǎn)置。

由上式可以得到每一次迭代過程中電導(dǎo)率的更新量δσ^k+1。

(8)收斂條件判斷

判斷收斂條件，若未達(dá)到收斂條件，則返回步驟(3)更新數(shù)據(jù)的參數(shù)；若達(dá)到收斂條件，則結(jié)束反演過程。

本發(fā)明在理論上具有較高的反演精度。圖3為反演5層結(jié)果對比圖，從圖3中可以看出，反演結(jié)果可以很準(zhǔn)確反演出地下電導(dǎo)率參數(shù)，并且反演的初值可以任意設(shè)定。這里僅使用1個觀察點(diǎn)和20個頻率點(diǎn)對分層情況進(jìn)行反演。圖4為反演10層的結(jié)果對比圖，假設(shè)地下為5層的分層模型，而實際分10層來反演，結(jié)果對比圖如圖3，可以看出在10層結(jié)果也是可以很好與地下分層情況相符合的。圖5為山東昌邑地區(qū)的實測數(shù)據(jù)的效果圖，為山東昌邑某地區(qū)的鐵礦采集區(qū)的一條半航空瞬變電磁數(shù)據(jù)反演結(jié)果圖。