基于DBIM的瞬變電磁電導率反演方法與流程

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

背景技術：

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

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

技術實現要素：

本發明的目的在于提供一種在瞬變電磁系統中能準確反演地層電導率，在頻率域中直接進行反演，不需要進行頻時轉換，在數據處理過程比較簡單，只需要準確提取相應頻率的頻譜即可，不僅適用于半航空瞬變電磁系統，同樣適用于全航空瞬變電磁系統的基于DBIM的瞬變電磁電導率反演方法。

本發明包括以下步驟：

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

將實際接收機接收數據作為輸入數據。

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

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

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

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

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

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

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

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

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

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

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

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

E＝<G^EJ；J>

H＝<G^HJ；J>

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

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

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

6)計算Frechet導數，具體步驟如下：

Frechet導數矩陣由以下方程求出：

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

δH_x＝δH_x1+δH_x2+δH_x3

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

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

對于上述的場值變化δH_x＝F·δσ，其中F即為地層參數的Frechet導數。

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

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

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

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

其中，為F矩陣的轉置，若其為復數矩陣，則表示為共軛轉置。

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

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

若未達到收斂條件，則返回步驟3)更新數據的參數；若達到收斂條件，則結束反演過程。

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

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

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

附圖說明

圖1為本發明的系統反演流程圖。

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

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

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

圖5為本發明實際數據效果圖。

具體實施方式

以下將結合附圖對本發明的技術方案及其突出效果作進一步說明。

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

(1)讀取觀測數據

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

(2)建立初始模型

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

初步定義下列初始參數：

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

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

(3)更新模型參數

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

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

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

(4)模型場值計算

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

E＝<G^EJ；J>

H＝<G^HJ；J>

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

(5)誤差計算

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

(6)Frechet導數計算

其中Frechet導數矩陣可由以下方程求出：

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

δH_x＝δH_x1+δH_x2+δH_x3

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

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

對于上述的場值變化δH_x＝F·δσ，其中F即為地層參數的Frechet導數。

(7)更新量計算

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

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

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

其中：為F矩陣的轉置，若其為復數矩陣，則表示為共軛轉置。

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

(8)收斂條件判斷

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

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