本發(fā)明具體涉及一種電性源瞬變電磁電場響應(yīng)成像方法，屬于瞬變電磁地球物理勘探方法技術(shù)領(lǐng)域。

背景技術(shù)：

瞬變電磁法(Transient electromagnetic method，TEM)被廣泛應(yīng)用于煤田水患、礦產(chǎn)資源和地下資源探測中，是一種重要的地球物理勘探方法。電性源瞬變電磁法屬于瞬變電磁法的一個(gè)分支，其采用電性源向地下注入電磁場信號，在源關(guān)斷期間采集來自地下地質(zhì)體的感應(yīng)信號以獲得其分布信息。

瞬變電磁場信號的處理和解釋是瞬變電磁法的關(guān)鍵技術(shù)之一。通過地面采集的電磁場信號，分析地下介質(zhì)的電性分布是一個(gè)典型的地球物理反問題。通常，瞬變電磁場響應(yīng)通過視電阻率、“煙圈”等效反演和S反演等快速成像方法或者基于最小二乘的OCCAM、Marquadt-Levenberg反演獲得地下介質(zhì)的電阻率分布。快速成像方法的成像精度有限，而基于最小二乘的反演算法受到瞬變電磁場響應(yīng)計(jì)算的復(fù)雜性和計(jì)算量大的限制，發(fā)展精度較高的高維反演難度較大。為此，探索新的瞬變電磁解釋技術(shù)具有重要的現(xiàn)實(shí)意義。

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

因此，本發(fā)明目的是目的在于提出一種對電性源瞬變電磁法電場響應(yīng)進(jìn)行成像的新方法。

具體的，本發(fā)明的方法包括以下步驟：

步驟A電性源瞬變電磁響應(yīng)虛擬子波解析解的推導(dǎo)；

步驟A具體包括：

步驟A1將場值進(jìn)行拉普拉斯變換；

步驟A2進(jìn)行非線性變換s＝p²，得到波動(dòng)場與擴(kuò)散場之間的等價(jià)表達(dá)式；

步驟A3反拉普拉斯變換；

步驟B電性源瞬變電磁響應(yīng)虛擬波場數(shù)值求解；

首先對公式一所示的積分方程進(jìn)行離散；

其中t為時(shí)間，(x,y,z)為空間坐標(biāo)，q被稱為虛擬時(shí)間，其量綱為s^1/2；

得到線性方程組：d＝Gm， 公式二

采用正則化的方法求解公式二的線性方程組，引入的模型向量的二范數(shù)作為正則化項(xiàng)，并采用L曲線法求取最佳正則化因子；

步驟C基于虛擬子波峰值速度的電性源瞬變電磁數(shù)據(jù)成像；

在步驟A和步驟B的基礎(chǔ)上，對電性源瞬變電磁場電場響應(yīng)進(jìn)行成像，分析虛擬子波峰值隨虛擬時(shí)間的傳播速度，估算地下介質(zhì)的電阻率，將基于所提取的虛擬子波的峰值虛擬時(shí)刻的傳播速度進(jìn)行成像。

進(jìn)一步的，所述步驟A的推導(dǎo)過程具體為：

在均勻半空間下，單位偶極源激發(fā)所產(chǎn)生的軸向電場響應(yīng)為：

其中，c²＝ρ/μ，μ＝4π10^-7H/m，erf為誤差函數(shù)。令

則公式三可以表示為：

對F₁(t)關(guān)于時(shí)間t進(jìn)行拉普拉斯變換，得到

與F₁(t)所對應(yīng)的虛擬波場的拉普拉斯變換為：

根據(jù)拉普拉斯變換恒等式

其中H表示Heavi side函數(shù)，令a＝r/c，對公式九進(jìn)行拉普拉斯反變換，得到

對公式六中的F₂(t)進(jìn)行拉普拉斯變換，我們得到

與F₂(t)所對應(yīng)的虛擬波場的拉普拉斯變換為：

對式公式十三進(jìn)行拉普拉斯反變換，可以得到：

將公式十一和公式十二進(jìn)行結(jié)合，得到與公式三所表示的均勻半空間下軸向電場解析式所對應(yīng)的虛擬波場的解析表達(dá)式：

化簡，得到：

進(jìn)一步的，所述步驟B中對公式一所示的積分方程進(jìn)行離散具體為：

積分區(qū)間表達(dá)為[a,b]，采用中點(diǎn)公式，將其離散為n個(gè)積分子區(qū)間，每個(gè)子區(qū)間的中點(diǎn)為q₁,q₂,...,q_n，則

其中，

假設(shè)

公式一可以離散為：

其中，

G_i,j＝G(t_i,q_j)Δq 公式二十一

m_j＝m(q_j) 公式二十二

則我們可以得到線性方程組：

d＝Gm 公式二。

進(jìn)一步的，所述步驟B中，采用奇異值分解法、Tikhonov零階正則化、一階正則化、截?cái)郤VD或阻尼SVD進(jìn)行正則化的求解。

進(jìn)一步的，所述步驟C中，電阻率和傳播速度的獲取方法具體為：

首先考慮均勻半空間情形，分別計(jì)算電阻率為200Ω·m,500Ω·m的均勻半空間下，單位偶極源激發(fā)所產(chǎn)生的下降沿階躍響應(yīng)經(jīng)波場變換得到的虛擬子波；

根據(jù)虛擬子波所滿足的波動(dòng)方程，可知虛擬子波的傳播速度為虛擬子波的峰值時(shí)刻等于；

通過計(jì)算峰值時(shí)刻隨偏移距的變換關(guān)系，得到虛擬子波的傳播速度，進(jìn)而得到地下介質(zhì)的電阻率，通過計(jì)算虛擬子波的峰值時(shí)刻對時(shí)間的偏導(dǎo)數(shù)，得到虛擬子波峰值時(shí)刻的傳播速度。

本發(fā)明的有益效果在于：本專利提出一種對電性源瞬變電磁法電場響應(yīng)進(jìn)行成像的新方法，解決了傳統(tǒng)的瞬變電磁快速成像方法和基于最小二乘的反演方法，在進(jìn)一步提高瞬變電磁響應(yīng)的解釋精度方面遇到瓶頸的問題。該方法基于擴(kuò)散場與波動(dòng)場之間的數(shù)學(xué)積分變換，計(jì)算量小、計(jì)算速度快，為瞬變電磁場響應(yīng)的擬波場精細(xì)解釋奠定了基礎(chǔ)。

附圖說明

圖1a至圖1d為采用模型二范數(shù)作為正則化項(xiàng)的求解結(jié)果，其中圖1a為零階Tikhonov正則化解；圖1b為TSVD解；圖1c為DSVD解；圖1d為不同截?cái)帱c(diǎn)的TSVD解；

圖2a至圖2d為采用一階模型粗糙度作為正則化項(xiàng)的求解結(jié)果,其中圖2a為一階Tikhonov正則化解；圖2b為TSVD解；圖2c為DSVD解；圖2d為不同截?cái)帱c(diǎn)的TSVD解；

圖3a至圖3d為采用二階模型粗糙度作為正則化項(xiàng)的求解結(jié)果,其中圖3a為二階Tikhonov正則化解；圖3b為TSVD解；圖3c DSVD解；圖3d為不同截?cái)帱c(diǎn)的TSVD解；

圖4a、圖4b為不同電阻率均勻半空間的虛擬子波解析解，圖4a中200 Ohm-m；圖4b中500Ohm-m；

圖5a、圖5b為不同電阻率均勻半空間的虛擬子波數(shù)值解，圖5a中200Ohm-m；圖5b中500Ohm-m；

圖6a、圖6b為Q型模型虛擬子波數(shù)值解，圖6a中ρ₁＝50Ohm-m,ρ₂＝200Ohm-m；圖6b中ρ₁＝200Ohm-m,ρ₂＝100Ohm-m；

圖7a、圖7b為三層模型虛擬子波數(shù)值解，圖7a中ρ₁＝200Ohm-m,ρ₂＝50Ohm-m,ρ₃＝200Ohm-m；圖7b中ρ₁＝50Ohm-m,ρ₂＝500Ohm-m,ρ₃＝50Ohm-m。

具體實(shí)施方式

下面結(jié)合附圖對本發(fā)明的具體實(shí)施方式進(jìn)行說明：

瞬變電磁場數(shù)據(jù)的擬波場解釋是瞬變電磁新解釋技術(shù)的探索方向之一。瞬變電磁場在地下介質(zhì)中的傳播滿足麥克斯韋方程組。受瞬變電磁信號的頻帶范圍的限制，其以擴(kuò)散方式傳播，傳播過程可以亥姆霍茨方程表征。Bragg(1968)和Filippi(1969)在數(shù)學(xué)上，從擴(kuò)散方程和波動(dòng)方程出發(fā)，推導(dǎo)了滿足擴(kuò)散方程的數(shù)學(xué)量與滿足波動(dòng)方程的數(shù)學(xué)量之間的數(shù)學(xué)積分關(guān)系式，若將該式應(yīng)用于滿足擴(kuò)散場的電磁場E與滿足波動(dòng)方程的虛擬子波U的轉(zhuǎn)換，可以得到

其中t為時(shí)間，(x,y,z)為空間坐標(biāo)，q被稱為虛擬時(shí)間，其量綱為s^1/2。Lee(1993)采用射線追蹤的方法，利用式(1)對低頻電磁場響應(yīng)進(jìn)行了成像，Lee(2006)利用頻率域電磁場數(shù)據(jù)，獲得虛擬子波U實(shí)現(xiàn)了頻率域電磁場響應(yīng)的成像處理。在國內(nèi)，李貅等(2006)將回線源瞬變電磁場響應(yīng)轉(zhuǎn)換為虛擬子波，并采用克西霍夫偏移方法實(shí)現(xiàn)了回線源瞬變電磁響應(yīng)的成像。

本發(fā)明提出一種新的電性源瞬變電磁響應(yīng)擬波場成像方法。該方法基于式(1)所推導(dǎo)的擴(kuò)散場與波場之間的積分轉(zhuǎn)換，能夠?qū)Φ叵陆橘|(zhì)的電阻率信息進(jìn)行精確成像。具體包括：(1)電性源瞬變電磁響應(yīng)虛擬子波解析解的推導(dǎo)；(2) 電性源瞬變電磁響應(yīng)虛擬波場數(shù)值求解方法；(3)基于虛擬子波峰值速度的電性源瞬變電磁數(shù)據(jù)成像。

(模塊1)：電性源瞬變電磁響應(yīng)虛擬子波解析解的推導(dǎo)

式(1)所示積分變換為第一類Fredholm積分方程，無法直接對其進(jìn)行反變換求取均勻半空間下與電磁場響應(yīng)對應(yīng)的虛擬子波的解析解。我們通過以下推導(dǎo)步驟對電性源瞬變電磁場響應(yīng)的虛擬子波解析解進(jìn)行推導(dǎo)：

1)將場值進(jìn)行拉普拉斯變換；

2)進(jìn)行非線性變換s＝p²，得到波動(dòng)場與擴(kuò)散場之間的等價(jià)表達(dá)式；

3)反拉普拉斯變換。

下面為具體推導(dǎo)過程：

在電性源瞬變電磁法測量中，若觀源軸向方向上的電場分量。在均勻半空間下，單位偶極源激發(fā)所產(chǎn)生的軸向電場響應(yīng)為：

其中，c²＝ρ/μ，μ＝4π10^-7H/m，erf為誤差函數(shù)。令

則式(2)可以表示為：

對F₁(t)關(guān)于時(shí)間t進(jìn)行拉普拉斯變換，得到

相應(yīng)地，我們得到與F₁(t)所對應(yīng)的虛擬波場的拉普拉斯變換為：

根據(jù)拉普拉斯變換恒等式

其中H表示Heavi side函數(shù)，令a＝r/c，對式(8)進(jìn)行拉普拉斯反變換，可以得到

同樣地，對式(5)中的F₂(t)進(jìn)行拉普拉斯變換，我們得到

相應(yīng)地，我們得到與F₂(t)所對應(yīng)的虛擬波場的拉普拉斯變換為：

對式(12)進(jìn)行拉普拉斯反變換，可以得到：

將式(10)和式(11)進(jìn)行結(jié)合，即可得到與式(2)所表示的均勻半空間下軸向電場解析式所對應(yīng)的虛擬波場的解析表達(dá)式：

化簡，得到：

(模塊2)：電性源瞬變電磁響應(yīng)虛擬子波數(shù)值求解方法

本模塊面向電性源瞬變電磁法虛擬子波提取。模塊1推導(dǎo)了均勻半空間虛擬子波的解析解。然而在層狀模型或者二、三維模型下，無法求得虛擬子波的解析解，需要采用數(shù)值方法進(jìn)行求解。本模塊給出了電性源瞬變電磁響應(yīng)虛擬子波數(shù)值求解的最優(yōu)方法。首先對式(1)所示的積分方程進(jìn)行離散。在實(shí)際算例中，積分區(qū)間是有限的，將其表達(dá)為[a,b]。采用中點(diǎn)公式，將其離散為n個(gè)積分子區(qū)間，每個(gè)子區(qū)間的中點(diǎn)為q₁,q₂,...,q_n，則

其中，

若假設(shè)

此時(shí)，式(1)可以離散為：

其中，

G_i,j＝G(t_i,q_j)Δq (20)

m_j＝m(q_j) (21)

則我們可以得到線性方程組：

d＝Gm (22)

因此，虛擬子波的求解問題轉(zhuǎn)換為了式(23)所示線性方程組進(jìn)行求解的過程。積分表達(dá)式(1)的病態(tài)性使得線性方程組(23)的求解過程具有病態(tài)性，需要采用正則化的方法求解。為此，引入的模型向量的二范數(shù)作為正則化項(xiàng)，并采用L曲線法求取最佳正則化因子。正則化的求解過程采用奇異值分解法(Singular value decomposition,SVD)進(jìn)行。考慮了Tikhonov零階正則化、一階正則化、截?cái)郤VD(truncated SVD,TSVD)和阻尼SVD(Damped SVD)等方法進(jìn)行求解。下面為求解結(jié)果：

圖1a至圖1c分別給出了引入模型二范數(shù)作為正則化項(xiàng)，分別采用Tikhonov正則化，TSVD和DSVD的求解結(jié)果，圖中縱坐標(biāo)表示振幅，曲線表示真實(shí)虛擬子波，數(shù)據(jù)點(diǎn)表示數(shù)據(jù)求解結(jié)果。對于Tikhonov正則化和DSVD，正則化因子由L曲線獲得；對于TSVD，其截?cái)帱c(diǎn)的位置由L曲線獲得，在圖中均以λ標(biāo)識為正則化因子。圖1a中零階Tikhonov正則化求解結(jié)果的晚期穩(wěn)定，但是早期偏離真實(shí)值較大。圖1c中DSVD求解結(jié)果在早期和晚期均不穩(wěn)定。圖1b中的TSVD求解結(jié)果的穩(wěn)定性要優(yōu)于另外兩種方法，由于奇異值的截?cái)啵缙诘慕獾牟环€(wěn)定性得到了壓制，但是解的晚期出現(xiàn)了輕微的振蕩。為此，計(jì)算了由L曲線獲得的最佳正則化參數(shù)附近的TSVD解，將其繪制于圖1b中。隨著截?cái)帱c(diǎn)后延，考慮的奇異值增多，所得到的TSVD解早期的穩(wěn)定性降低。但是，所獲得的虛擬子波的幅值更接近與真實(shí)值。相反地，隨著截?cái)帱c(diǎn)提前，考慮的奇異值減小，得到的TSVD解總體而言穩(wěn)定性增加?？傮w而言，采用模型二范數(shù)作為正則化項(xiàng)時(shí)，所得到的正則化解早期的穩(wěn)定性較差。

圖2和圖3分別一階模型粗糙度和二階模型粗糙度作為正則化項(xiàng)得到的虛擬子波提取結(jié)果。當(dāng)采用高階模型粗糙度時(shí)，需要采用廣義SVD(Generalized SVD)方法進(jìn)行求解。由圖可知，當(dāng)采用模型粗糙度作為正則化項(xiàng)時(shí)，Tikhonov方法和TGSVD方法均能得到的穩(wěn)定的解。但是，采用阻尼奇異值分解得到的解的效果仍然不理想。通過對比Tikhonov方法和TGSVD方法所獲得的解，認(rèn)為TGSVD對虛擬子波的峰值和峰值時(shí)刻的提取效果更好。另外，通過對比圖2d和圖3d的TGSVD的提取結(jié)果可知，采用二階模型粗糙度作為正則化項(xiàng)時(shí)，所提取的虛擬子波的峰值和峰值時(shí)刻更為準(zhǔn)確。當(dāng)采用一階模型粗糙度作為正則化項(xiàng)時(shí)，所提取的虛擬子波的晚期振蕩較弱。

通過上述分析可知，波場變換的數(shù)值解并不能準(zhǔn)確的恢復(fù)出虛擬子波，而僅是虛擬子波的一個(gè)近似。求解病態(tài)反問題的難點(diǎn)在于反問題本身，而非所采用的正則化算法。因此，對病態(tài)反問題的理解和對解的認(rèn)識往往比單純求得反問題的穩(wěn)定解更為重要。本模塊通過采用均勻半空間下的瞬變電磁響應(yīng)及其相對應(yīng)的虛擬子波作為參照，給出了不同正則化方法的解的特征。在后續(xù)的成像 中，所關(guān)心的是所提取的虛擬子波的峰值時(shí)刻，因此在這些正則化方法中，采用二階模型粗糙度并采用TGSVD所獲得的解的特性更符合成像需要。

(模塊3)：電性源瞬變電電磁法虛擬波場成像

本模塊在模塊1和模塊2的基礎(chǔ)上，對電性源瞬變電磁場電場響應(yīng)進(jìn)行成像。在瞬變電磁場數(shù)據(jù)解釋中，成像的目的是獲得地下介質(zhì)的電阻率特性。均勻半空間下，虛擬子波隨虛擬時(shí)間的傳播速度與地下介質(zhì)的電阻率相關(guān)。通過分析虛擬子波峰值隨虛擬時(shí)間的傳播速度可以估算地下介質(zhì)的電阻率。為了借鑒地震數(shù)據(jù)處理中的成像技術(shù)，本模塊將基于所提取的虛擬子波的峰值虛擬時(shí)刻的傳播速度進(jìn)行成像。

首先考慮均勻半空間情形，分別計(jì)算電阻率為200Ω·m,500Ω·m的均勻半空間下，單位偶極源激發(fā)所產(chǎn)生的下降沿階躍響應(yīng)經(jīng)波場變換得到的虛擬子波。圖4為偏移距1000m至7000m范圍內(nèi)虛擬子波的解析解，圖5為相對應(yīng)的數(shù)值解，圖4中線段表示所拾取的虛擬子波的峰值時(shí)刻隨偏移距的變化，圖5中虛線對應(yīng)圖4中的線段。由于此處考慮的是虛擬子波的峰值時(shí)刻，因此對圖4和圖5的虛擬子波的幅值進(jìn)行了歸一化處理。根據(jù)虛擬子波所滿足的波動(dòng)方程，可知虛擬子波的傳播速度為另外，均勻半空間下，虛擬子波的解析表達(dá)式表明，虛擬子波的峰值時(shí)刻等于

由此可知，通過計(jì)算峰值時(shí)刻隨偏移距的變換關(guān)系，即可得到虛擬子波的傳播速度，進(jìn)而得到地下介質(zhì)的電阻率。通過計(jì)算線段所示虛擬子波的峰值時(shí)刻對時(shí)間的偏導(dǎo)數(shù)，即可得到虛擬子波峰值時(shí)刻的傳播速度。表1給出了不同電阻率的均勻半空間下虛擬子波的傳播速度及相對應(yīng)的電阻率。

表1

由表1可知，在不同的電阻率情形下，通過計(jì)算虛擬子波的峰值時(shí)刻的傳 播速度，可以得到較為準(zhǔn)確的地下介質(zhì)的電阻率。在低阻介質(zhì)中，虛擬子波的峰值時(shí)刻的傳播速度較小，在高阻介質(zhì)中虛擬子波峰值時(shí)刻的傳播速度較大，這與電磁場在地下介質(zhì)的擴(kuò)散速度的基本規(guī)律是一致的。

數(shù)值算例

以數(shù)值算例來說明本專利所提出的電性源瞬變電磁電場響應(yīng)成像方法及其應(yīng)用效果。采用了兩層模型和三層模型驗(yàn)證了本專利提出的成像方法的有效性。

兩層模型算例。圖6給出了Q型和D型模型虛擬子波的數(shù)值解。圖6a所示虛擬子波由參數(shù)為ρ₁＝50Ω·m，h₁＝800m，ρ₂＝200Ω·m的Q型模型的瞬變電磁場響應(yīng)生成，圖6b所示虛擬子波由參數(shù)為ρ₁＝200Ω·m，h₁＝800m，ρ₂＝100Ω·m的D型模型的瞬變電磁場響應(yīng)生成。圖中線段表示所拾取虛擬子波的峰值時(shí)刻隨偏移距的變化。與均勻半空間不同，Q型模型下所拾取的線段分為兩段。由圖6a計(jì)算兩段虛擬子波的峰值時(shí)刻的傳播速度分別為和計(jì)算得到相應(yīng)的電阻率為53.89Ω·m和208.24Ω·m。由圖6a計(jì)算兩段虛擬子波的峰值時(shí)刻的傳播速度為和相對應(yīng)的兩層模型的電阻率為200.52Ω·m和107.79Ω·m。計(jì)算得到的電阻率能夠很好的與模型的真實(shí)電阻率吻合。

三層模型算例。圖7a和圖7b分別給出了H型模型和K型模型的虛擬子波的數(shù)值解。H型模型的參數(shù)為ρ₁＝200Ω·m，h₁＝800m，ρ₂＝50Ω·m，h₂＝200m,ρ₃＝200Ω·m，K型模型的參數(shù)為ρ₁＝50Ω·m，h₁＝800m，ρ₂＝500Ω·m，h₂＝200m,ρ₃＝50Ω·m。由圖7可知，所拾取的虛擬子波的峰值時(shí)刻隨偏移距的變換可以分為三條線段，H型模型相對應(yīng)的虛擬子波的峰值時(shí)刻的傳播速度分別為和計(jì)算得到的響應(yīng)的電阻率為204.86Ω·m，53.72Ω·m和195.52Ω·m。K型模型相對應(yīng)的虛擬子波的峰值時(shí)刻的傳播速度分別為和計(jì)算得到的電阻率分別為53.74Ω·m，488.67Ω·m和190.52Ω·m。H型模型與K型模型各層的電阻率均得到了較好的恢復(fù)。由此可知，本發(fā)明所提出的成像方法能夠?qū)Φ叵陆橘|(zhì)的電阻率進(jìn)行較精確成像。

以上所述是本發(fā)明的優(yōu)選實(shí)施方式，應(yīng)當(dāng)指出，對于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明所述原理的前提下，還可以作出若干改進(jìn)和潤飾，這些改進(jìn)和潤飾也應(yīng)視為本發(fā)明的保護(hù)范圍。