周建美，劉文韜，齊彥福，李貅*，戚志鵬，魯凱亮,劉正銀,邢濤

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安 710054 2 長安大學(xué)地球物理場多參數(shù)綜合模擬實驗室(中國地球物理學(xué)會重點實驗室)，西安 710054 3 山東省交通規(guī)劃設(shè)計院，濟南 250031 4 北京探創(chuàng)資源科技有限公司，北京 100071

0 引言

瞬變電磁法作為一種重要的電磁勘探方法(底青云等，2019;Xue et al.,2020)，廣泛應(yīng)用于金屬礦產(chǎn)勘探(Di et al.,2018;底青云等，2020)、煤田災(zāi)害探測(王若等,2018;Jiang et al.，2019;Chang et al.，2019)、水文地球物理研究(Danielsen et al.，2003)和隧道超前地質(zhì)預(yù)報(Xue et al.,2018)等領(lǐng)域.理論上瞬變電磁發(fā)射波形設(shè)計為一個階躍波，接收不含一次場的地下感應(yīng)渦流產(chǎn)生的純二次場信號(Nabighian and Macnae，1991).由于儀器制造上的限制(Davis and Macnae，2008)或探測目標(biāo)的需要(殷長春等，2015)，實際中采用多種發(fā)射波形，包括半正弦波、三角波、梯形波(Liu，1998)和其他復(fù)雜波形(Davis and Macnae，2008;Chen et al.，2015).

發(fā)射波形對瞬變電磁響應(yīng)的影響得到了廣泛的研究.Liu(1998)、陳曙東等(2012)利用自由空間中回線作為目標(biāo)體，推導(dǎo)了不同波形的解析解，并研究了不同發(fā)射波形對機載瞬變電磁響應(yīng)的影響.Chen等(2015)開發(fā)了一個多波形(MULTIPULSE)系統(tǒng)，該系統(tǒng)在半周期內(nèi)間斷發(fā)射一個高功率半正弦波和一個低功率梯形波，保證探測深度的同時提供較高的近表面分辨率.殷長春等(2015)以半正弦波和梯形波為例，系統(tǒng)研究了幾種典型地下目標(biāo)情況下機載TEM方法的測量能力，并提供了參數(shù)組合以實現(xiàn)最佳勘探效果.不同發(fā)射波形的瞬變電磁響應(yīng)不同，因此在開發(fā)瞬變電磁三維正演算法時需要考慮發(fā)射波形的影響(Zeng et al.，2019).

計算全波形瞬變電磁正演響應(yīng)的方法主要有以下幾種：(1)褶積算法(殷長春等，2013)，對階躍響應(yīng)與發(fā)射波形的時間導(dǎo)數(shù)進行褶積，可以得到任意發(fā)射波形的全時瞬變電磁響應(yīng).該方法可以直接利用高效的階躍響應(yīng)正演算法，但該算法需要獲得精確的電流二階導(dǎo)數(shù)，實際中發(fā)射電流數(shù)據(jù)采樣率有限，難以獲得電流二階導(dǎo)數(shù)的精確值，嚴重影響了褶積結(jié)果的計算精度(齊彥福等，2017)；(2)時域有限差分法(Commer and Newman,2004；孫懷鳳等，2013)，直接對連續(xù)發(fā)射電流波形進行采樣并將其包含到控制方程式中，從而得到全時響應(yīng).該方法不需要求解大規(guī)模線性方程組，內(nèi)存需求相對較低，但為了保持其數(shù)值穩(wěn)定性，只能采用非常小的時間步長，從而需要較長的計算時間(孫懷鳳等，2013)；(3)后退Euler法(Um et al.，2010；齊彥福等，2017；Zeng et al.，2019)，對發(fā)射波形的處理方法與時域有限差分法類似，雖然該方法是無條件穩(wěn)定的，但仍需要謹慎選取時間步長，以在計算效率和準(zhǔn)確性之間取得有效的平衡(Li et al.,2018).

對于三維階躍響應(yīng)計算，基于Krylov子空間投影的模型降階算法(Zaslavsky et al.，2011;B?rner et al.，2015;Zhou et al.，2018；Qiu et al.，2019)比時域有限差分法和后退Euler法具有更高的計算效率(Zaslavsky et al.，2011;B?rner et al.，2015).結(jié)合直接法求解器，一次時間域正演的模型降階算法的計算量可以減少為一次系數(shù)矩陣分解和上百次矩陣回代(Zhou et al.,2018).最近Zhou等(2020)將模型降階算法擴展到了全波形瞬變電磁正演響應(yīng)的計算，通過指數(shù)求積法則(Hochbruck and Ostermann,2010)，將on-time響應(yīng)和off-time響應(yīng)分別表示為關(guān)于空間離散系數(shù)矩陣的phi函數(shù)和指數(shù)函數(shù)，采用基于位移逆Krylov子空間投影的模型降階算法(Zhou et al.,2018)實現(xiàn)了高精度正演計算.但該算法在計算on-time時間內(nèi)每個時刻的正演響應(yīng)時都需要構(gòu)建一個新的位移逆Krylov子空間；而為了保證on-time響應(yīng)的正演精度，一般需要在on-time范圍內(nèi)計算上百個時刻的響應(yīng)，導(dǎo)致該算法在on-time時間內(nèi)的正演響應(yīng)計算非常耗時.

由于儀器設(shè)備的限制，目前實際工作中一般只對off-time數(shù)據(jù)進行處理解釋，因此對于實際需求而言，只需要精確模擬任意發(fā)射波形的off-time響應(yīng)即可.Smith和Neil(2013)指出，on-time發(fā)射波形非精確采樣，也可以得到精確的off-time響應(yīng).Zhou等(2020)同樣發(fā)現(xiàn)，在一定范圍內(nèi)改變on-time發(fā)射波形采樣時刻數(shù)量，對off-time響應(yīng)的精度沒有顯著影響.由于模型降階算法(Zhou et al.,2020)on-time響應(yīng)的計算時間與on-time范圍內(nèi)的采樣時刻數(shù)量近似成正比，因此減少on-time范圍內(nèi)的采樣時刻數(shù)量，就能夠減少正演計算時間.基于頻譜分析理論，時間域信號是頻率域信號的疊加，而頻率域信號的幅值隨時間快速衰減，因此off-time某一時刻的響應(yīng)只和特定頻率范圍內(nèi)的信號有關(guān).基于這一原則，本文提出了一種定量化評價on-time采樣波形精度的方法，能夠在保證off-time響應(yīng)精度的同時，顯著的減少on-time采樣時刻，最大程度的減少正演計算時間.本文首先簡述全波形瞬變電磁模型降階算法；然后給出定量化評價on-time采樣波形精度的方法；最后通過典型發(fā)射波形的正演計算驗證本文算法的精度和計算速度.

1 正演算法

1.1 控制方程離散

全波形瞬變電磁法對應(yīng)的時間域Maxwell方程，忽略位移電流：

(1)

式中，E是電場強度矢量，B是磁感應(yīng)強度矢量，t是時間，σ是電導(dǎo)率，μ是磁導(dǎo)率，s是源項，考慮如圖1所示的任意發(fā)射波形的電流源，在on-time時間段，0

圖1 任意發(fā)射波形示意圖Fig.1 Schematic diagram of the arbitrary transmitting waveform

采用簡單的自然邊界條件(Haber,2014)，利用擬態(tài)有限體積方法(周建美等，2018)對方程(1)進行空間離散，消去電場，得到關(guān)于磁場的空間離散控制方程：

(2)

(3)

式中，tk=tk-1+hk，k=1,2,…,n,hk是時間步長,t0=0,b0=0.

對于on-time時間段，公式(3)中g(shù)不為零.利用指數(shù)求積法則(Hochbruck and Ostermann,2010)對公式(5)中的積分項進行離散化，并引入函數(shù)φj(z),j=0,1,2,…，滿足φj(0)=1/j!且具有遞歸關(guān)系φ0(z)=exp(z)，φj+1(z)=[φj(z)-φj(0)]/z，則on-time磁場響應(yīng)可以表示為(Zhou et al.,2020)：

bk=bk-1+hk(-Abk-1+gk-1)+hkφ2(-hkA)

×[-hkA(-Abk-1+gk-1)+gk-gk-1]，

0

(4)

對于off-time時間段，公式(3)中g(shù)為零.為了避免在求解每個時刻的響應(yīng)時都重新構(gòu)建子空間，統(tǒng)一以b(toff)為起點(B?rner et al.,2015;Zhou et al.,2018)，則off-time時間段任意時刻磁場的求解不再依賴于上一個時刻的磁場，簡化后的公式(3)可以轉(zhuǎn)化為：

bk=e-(tk-toff)Ab(toff),tk≥toff，

(5)

綜合公式(4)和(5)，全時域磁場響應(yīng)可表示為：

(6)

其中r=φ2(-hkA)[-hkA(-Abk-1+gk-1)+gk-gk-1]，q=e-(tk-toff)Ab(toff).可以看出，公式(6)中大部分的計算量集中于r和q的求解.它們可以寫成統(tǒng)一的矩陣與向量乘積的函數(shù)形式

y=f(A)v，

(7)

對于on-time 時段，

f(A)=φ2(-hkA),

v=-hkA(-Abk-1+gk-1)+gk-gk-1,

0

(8)

對于off-time時段，f(A)=e(toff-tk)A,v=b(toff),tk≥toff.

(9)

1.2 位移逆Krylov子空間算法

采用基于位移逆Krylov子空間投影的模型降階算法 (van den Eshof and Hochbruck,2006;Botchev,2016;Zhou et al.,2018)求解公式(7)的近似解，主要計算量是構(gòu)建合適的Krylov子空間.對于給定的位移量γ和子空間階數(shù)m，采用Arnoldi 方法(Zhou et al.,2018)構(gòu)建Krylov子空間Km((I+γA)-1,v).根據(jù)子空間的正交基Vm和上Hessenberg矩陣Hm+1,m得到公式(11)的模型降階解為(Hochbruck and Ostermann,2010)：

(10)

其中，Hm,m是上Hessenberg矩陣Hm+1,m中除去第m+1列，e1=[1,0,…,0]T.公式(10)中矩陣的維度m遠遠小于原系統(tǒng)中系數(shù)矩陣A的維度，可以快速求解.注意到公式(7)中向量v在off-time時段是常數(shù)，而在on-time時段隨著發(fā)射時間變化，因此on-time時段的計算比off-time時段的計算復(fù)雜.為了得到整體高效的算法，對這兩個時段分別采用不同的子空間近似策略進行求解.

在off-time時段，把公式(10)代入公式(9)中，得到off-time時段的磁場響應(yīng)：

(11)

在on-time時段，把公式(10)代入公式(8)中，而后代入公式(6)，可得on-time時段的磁場響應(yīng)：

bk≈bk-1+hk(-Abk-1+gk-1)+hkVmkφ2

(12)

1.3 頻譜分析

在on-time時段進行密集時刻采樣，可以得到精確的全波形正演響應(yīng)(Zhou et al.,2020).但在on-time時段每個計算時刻tk需要重新構(gòu)建子空間Kmk((I+γkA)-1,vk)，導(dǎo)致on-time響應(yīng)的計算效率低下.如果只關(guān)注off-time響應(yīng)，則on-time采樣時刻可以顯著減少，從而減小總的正演計算時間.

由傅里葉變換理論，任意一個時間域信號，都可以分解為無窮多個不同頻率的正弦信號的疊加(胡廣書,2003)，因此可以從頻率域信號出發(fā)，分析接收的時間域信號的特點.如圖1所示，瞬變電磁發(fā)射波形是一個連續(xù)非周期信號.實際工作中對發(fā)射波形進行離散化采樣，得到時間域中一個離散非周期信號.對該信號進行快速傅里葉變換,即可得到該發(fā)射波形對應(yīng)的頻譜分布(胡廣書,2003).下面從頻譜分析角度，根據(jù)一維介質(zhì)中電磁場的分布特征，定量化的給出on-time時段的波形采樣規(guī)則.

(1)根據(jù)電磁場的傳播理論，忽略位移電流，一維介質(zhì)中頻率域磁場的通解為(湯井田和何繼善，2005)：

B=B0e-z/δeiφ，

(13)

Vp=ωδ.

(14)

對于瞬變電磁法，在發(fā)射源關(guān)斷后t時刻接收到的時間域信號中，頻率為f的場的相位傳播距離為

z1=Vpt=2πfδt.

(15)

(2)常見發(fā)射波形的頻譜中，高頻信號的幅值相比低頻信號的幅值小(牛之璉，2007)，對于測量的時間域信號的貢獻有限，可以進一步將比較波形的頻譜范圍限定為大于最大幅值A(chǔ)max的1%的頻率范圍.本文令A(yù)c=0.01Amax為截止振幅.

基于上述兩條原則，可以在保證off-time響應(yīng)計算精度的同時，極大的減少on-time采樣時間的數(shù)量.

2 數(shù)值實例

通過對比基于位移逆Krylov子空間模型降階的全波形算法(Zhou et al.,2020)(簡稱為SAI)和本文結(jié)合頻譜分析的位移逆Krylov子空間快速算法(簡稱為SAI-SA)計算典型的半正弦波、梯形波和VTEM波的正演響應(yīng)，測試本文算法在計算off-time響應(yīng)的精度和速度.SAI算法對on-time波形進行密集時刻采樣(Zhou et al.,2020)，SAI-SA算法對on-time波形進行基于頻譜分析的稀疏時刻采樣.本文計算設(shè)備為32 G內(nèi)存、四核主頻3.6 GHz的Intel i7 CPU的臺式電腦.

2.1 半正弦波

首先采用均勻半空間模型，通過與一維解析解結(jié)果(殷長春等，2013)對比，驗證本文算法的計算精度.半空間模型的電阻率為10 Ωm，空氣層電阻率設(shè)置為10-6Ωm，采用中心回線裝置的航空裝置，發(fā)射線圈是邊長6.66 m的正八邊形，飛行高度為30 m.發(fā)射波形設(shè)置為單位電流強度、脈沖寬度為4 ms的半正弦波(圖2a中黑色曲線所示)，接收off-time時間段[0.01,10]ms的dbz/dt響應(yīng).

圖2 半正弦波形和采樣數(shù)nton=6的采樣波形對比(a)半正弦波形(黑色)和采樣波形(灰色)；(b)半正弦波形(黑粗線)和采樣波形(灰粗線)頻譜分布，虛線對應(yīng)截止頻率，黑色細實線對應(yīng)截止振幅；(c)采樣波形的頻譜分布與半正弦波形的頻譜分布的比值.Fig.2 Comparison of half-sine waveform and sampling waveform with nton=6(a)Half-sine waveform (black)and sampling waveform (gray);(b)Spectrum of half-sine waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and half-sine waveform.

圖3 半正弦波形和采樣數(shù)nton=7的采樣波形對比(a)半正弦波形(黑色)和采樣波形(灰色)；(b)半正弦波形(黑粗線)和采樣波形(灰粗線)頻譜分布，虛線對應(yīng)截止頻率，黑色細實線對應(yīng)截止振幅；(c)采樣波形的頻譜分布與半正弦波形的頻譜分布的比值.Fig.3 Comparison of half-sine waveform and sampling waveform with nton=7(a)Half-sine waveform (black)and sampling waveform (gray);(b)Spectrum of half-sine waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and half-sine waveform.

圖4 半正弦波形和采樣數(shù)nton=8的采樣波形對比(a)半正弦波形(黑色)和采樣波形(灰色)；(b)半正弦波形(黑粗線)和采樣波形(灰粗線)頻譜分布，虛線對應(yīng)截止頻率，黑色細實線對應(yīng)截止振幅；(c)采樣波形的頻譜分布與半正弦波形的頻譜分布的比值.Fig.4 Comparison of half-sine waveform and sampling waveform with nton=8(a)Half-sine waveform (black)and sampling waveform (gray);(b)Spectrum of half-sine waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and half-sine waveform.

接下來采用采樣數(shù)nton=7的SAI-SA算法計算半正弦波的正演響應(yīng)，同時采用密集采樣數(shù)nton=100的SAI算法計算的結(jié)果(Zhou et al.,2020)作為三維高精度結(jié)果進行比較.結(jié)果如圖5所示，圖5a為三維正演結(jié)果，圖5b為三維正演結(jié)果相對一維解析解的誤差，其中實線為一維解析解，圓圈為采樣數(shù)nton=100的SAI算法的計算結(jié)果，叉點為采樣數(shù)nton=7的SAI-SA算法的計算結(jié)果.由圖5可知，密集采樣的SAI算法能夠得到on-time和off-time時間段內(nèi)的高精度正演結(jié)果；而稀疏采樣的SAI-SA算法的正演結(jié)果雖然在on-time時間段的誤差很大，但在off-time時間段內(nèi)與密集采樣的SAI正演結(jié)果基本一致，相對一維解析解的相對誤差均小于5%.

圖5 半正弦波模型的dbz/dt響應(yīng)(a)及誤差對比(b)實線為一維解析解，圓圈為采樣nton=100的SAI算法的正演結(jié)果，叉點為采樣nton=7的SAI-SA算法的正演結(jié)果.Fig.5 dbz/dt response (a)and error comparison (b)of the half-sine waveform modelThe line is 1D analytical solution,circle is result of SAI algorithm with nton=100 and cross is result of SAI-SA algorithm with nton=7.

表1給出了密集采樣的SAI算法和稀疏采樣的SAI-SA算法的計算時間對比.由表1可知，采用密集采樣nton=100的SAI算法，由于每次計算on-time響應(yīng)時需要構(gòu)建一次新的Krylov子空間，導(dǎo)致總的系數(shù)矩陣回代次數(shù)高達1622，總的正演時間為217.1 s；而采用稀疏采樣nton=7的SAI-SA算法，計算on-time響應(yīng)只需要構(gòu)建7次新的Krylov子空間，從而總的系數(shù)矩陣回代次數(shù)減少為171，總的正演時間減少為39.5 s，計算速度提高為密集采樣算法的5.5倍.

表1 SAI-SA算法和SAI算法計算半正弦波模型的時間對比Table 1 Solution time of SAI-SA algorithm and SAI algorithm for a half-sine waveform model

2.2 梯形波

接下來考察梯形波的正演.采用與半正弦發(fā)射波形相同的半空間模型和發(fā)射接收裝置系統(tǒng).發(fā)射波形設(shè)置為單位電流強度，脈沖寬度為4 ms、上升沿和下降沿均為0.2 ms的梯形波(圖6a中黑色曲線)，接收off-time時間段[0.01,10]ms的dbz/dt響應(yīng).利用SAI-SA算法計算off-time響應(yīng).采用與半正弦發(fā)射波形相同的空間離散網(wǎng)格.先確定on-time時間段的采樣波形.對發(fā)射波形采用線性等間隔采樣，顯然當(dāng)采樣數(shù)nton=20(不包括初始時刻t=0)時，采樣波形與原梯形波完全重合(圖6a中灰色圓圈)，頻譜分布也完全一致(圖6b).注意到梯形波的規(guī)則形狀，如果采用線性不等間隔采樣，采樣數(shù)nton=3(不包括初始時刻t=0)，既可以保證采樣波形與原梯形波完全重合(圖7a)，頻譜分布也完全一致(圖7b).

圖6 梯形波形和采樣數(shù)nton=20的采樣波形對比(a)梯形波形(黑色)和采樣波形(灰色)；(b)梯形波形(黑粗線)和采樣波形(灰粗線)頻譜分布，虛線對應(yīng)截止頻率，黑色細實線對應(yīng)截止振幅；(c)采樣波形的頻譜分布與梯形波形的頻譜分布的比值.Fig.6 Comparison of trapezoid waveform and sampling waveform with nton=20(a)Trapezoid waveform (black)and sampling waveform (gray);(b)Spectrum of trapezoid waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and trapezoid waveform.

圖7 梯形波形和采樣數(shù)nton=3的采樣波形對比(a)梯形波形(黑色)和采樣波形(灰色)；(b)梯形波形(黑粗線)和采樣波形(灰粗線)頻譜分布，虛線對應(yīng)截止頻率，黑色細實線對應(yīng)截止振幅；(c)采樣波形的頻譜分布與梯形波形的頻譜分布的比值.Fig.7 Comparison of trapezoid waveform and sampling waveform with nton=3 (a)Trapezoid waveform (black)and sampling waveform (gray);(b)Spectrum of trapezoid waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and trapezoid waveform.

但在on-time時間段采用線性不等間隔采樣時，需要重新考慮最優(yōu)化位移量γon的選取.當(dāng)on-time時間段采用線性等間隔采樣，SAI算法的最優(yōu)化位移量γon與采樣間隔成正比(Zhou et al.,2020).當(dāng)在on-time時間段采用線性不等間隔采樣時，每個采樣間隔對應(yīng)一個最優(yōu)化位移量，但采用多個不同位移量會增加系數(shù)矩陣的LU分解次數(shù)，從而增大計算時間.由于位移量對采樣間隔不靈敏(van den Eshof and Hochbruck,2006；Zhou et al.，2018)，位移量γon選取可以采用兩種簡單的方法：(1)選取某個采樣間隔對應(yīng)的最優(yōu)化位移量作為總的最優(yōu)化位移量，比如以第一個時間間隔對應(yīng)的最優(yōu)化位移量作為總的最優(yōu)化位移量，對于采樣數(shù)nton=3，采用SAI-SA算法，求解on-time響應(yīng)需要構(gòu)建的3次Krylov子空間，子空間階數(shù)分別為m1=19，m2=55，m3=19；(2)選取off-time時間段對應(yīng)的最優(yōu)化位移量γoff作為總的最優(yōu)化位移量，對于采樣數(shù)nton=3，采用SAI-SA算法，求解on-time響應(yīng)構(gòu)建的3次Krylov子空間的階數(shù)分別為m1=16，m2=30，m3=16.兩種方法需要構(gòu)建的Krylov子空間階數(shù)相差不大，對應(yīng)計算量也相差不大，而且都表現(xiàn)為大的采樣間隔對應(yīng)大的子空間階數(shù).但選取γoff作為總的最優(yōu)化位移量時，求解on-time時間段正演響應(yīng)需要LU分解的系數(shù)矩陣與求解off-time時間段正演響應(yīng)需要LU分解的系數(shù)矩陣相同，因此可以在求解過程中減少1次系數(shù)矩陣的LU分解，表現(xiàn)出顯著的計算優(yōu)勢.因此本文選取γon=γoff來計算線性不等間距稀疏采樣波形的正演響應(yīng).

接收off-time時間段[0.01,10]ms仍然離散為31個對數(shù)等間隔時刻.分別采用密集采樣數(shù)nton=100的SAI算法(Zhou et al.,2020)、采樣數(shù)nton=20的線性等間隔采樣的SAI-SA算法、采樣數(shù)nton=3的線性不等間隔采樣的SAI-SA算法和一維解析算法(殷長春等，2013)計算梯形波模型的正演響應(yīng).結(jié)果如圖8所示，圖8a為正演計算結(jié)果，圖8b為三維正演結(jié)果相對一維解析解的誤差，其中實線為一維解析解，圓圈為nton=100的SAI算法的計算結(jié)果，三角形為nton=20的SAI-SA算法的計算結(jié)果，叉點為nton=3的SAI-SA算法的計算結(jié)果.由圖8可知，nton=100的密集采樣的SAI算法能夠得到on-time和off-time時間段內(nèi)的高精度正演結(jié)果；而nton=20和nton=3的稀疏采樣的SAI-SA算法雖然在on-time時間段的誤差很大，但由于nton=20和nton=3的稀疏采樣的頻譜分布與nton=100的密集采樣的頻譜分布是相同的，都與原梯形波完全重合，因此三種采樣得到的off-time時間段內(nèi)的正演響應(yīng)是完全相同的，相對一維解析解的相對誤差均小于5%.

圖8 梯形波模型的dbz/dt響應(yīng)(a)及誤差對比(b)實線為一維解析解，圓圈為采樣nton=100的SAI算法的正演結(jié)果，叉點為采樣nton=3的SAI-SA算法的正演結(jié)果，三角形為采樣nton=20的SAI-SA算法的正演結(jié)果.Fig.8 dbz/dt response (a)and error comparison (b)of the trapezoid waveform modelThe line is 1D analytical solution,circle is result of SAI algorithm with nton=100,cross is result of SAI-SA algorithm with nton=3,and triangle is result of SAI-SA algorithm with nton=20.

表2給出了密集采樣的SAI算法和稀疏采樣的SAI-SA算法的計算時間對比.由表2可知，采用密集采樣nton=100的SAI算法，由于總的系數(shù)矩陣回代次數(shù)高達1546，總的正演時間為213 s；采用線性等間隔稀疏采樣nton=20的SAI-SA算法，總的系數(shù)矩陣回代次數(shù)減少為303，總的正演時間減少為56.5 s，計算速度提高為密集采樣算法的3.8倍；采用線性不等間隔稀疏采樣nton=3的SAI-SA算法，不僅總的系數(shù)矩陣回代次數(shù)減少為171，而且減少了一次LU分解，因此總的正演時間減少為36.5 s，計算速度提高為密集采樣算法的5.8倍.

表2 SAI-SA算法和SAI算法計算梯形波模型的時間對比Table 2 Solution time of SAI-SA algorithm and SAI algorithm for a trapezoid waveform model

2.3 VTEM波

進一步采用復(fù)雜VTEM波的三維模型(齊彥福等，2017)的正演響應(yīng)驗證本文算法的有效性.三維模型如圖9所示，電阻率為100 Ωm的均勻半空間中存在一個電阻率為1 Ωm的低阻異常體，異常體上頂面距離地表50 m，異常體尺寸為100 m×100 m×200 m，空氣層電阻率設(shè)置為106Ωm，在三維目標(biāo)體的正上方采用與半正弦發(fā)射波形相同的發(fā)射接收裝置系統(tǒng).

圖9 三維模型Fig.9 3-D model

圖10 VTEM波形和采樣數(shù)nton=26的采樣波形對比(a)VTEM波形(黑色)和采樣波形(灰色)；(b)VTEM波形(黑粗線)和采樣波形(灰粗線)頻譜分布，虛線對應(yīng)截止頻率，黑色細實線對應(yīng)截止振幅；(c)采樣波形的頻譜分布與VTEM波形的頻譜分布的比值.Fig.10 Comparison of VTEM waveform and sampling waveform with nton=26 (a)VTEM waveform (black)and sampling waveform (gray);(b)Spectrum of VTEM waveform (black thick line)and sampling waveform (gray thick line),dotted line is cut-off frequency,black thin line is cut-off amplitude;(c)Ratio of spectrum amplitude of sampling waveform and VTEM waveform.

接下來采用采樣數(shù)nton=26的SAI-SA算法計算VTEM發(fā)射波形的正演響應(yīng)，同時采用密集采樣數(shù)nton=200的SAI算法計算的結(jié)果(Zhou et al.,2020)作為三維高精度結(jié)果進行比較.結(jié)果如圖11所示，圖11a為三維正演結(jié)果，其中實線為三維隱式時間步迭代算法的計算結(jié)果(齊彥福等，2017)，圓圈為采樣數(shù)nton=200的SAI算法的計算結(jié)果，叉點為采樣數(shù)nton=26的SAI-SA算法的計算結(jié)果.圖11b為nton=26的SAI-SA算法的計算結(jié)果相對nton=200的SAI算法的計算結(jié)果的誤差.由圖11可知，與半正弦波模型正演和梯形波模型正演結(jié)果一致，密集采樣的SAI算法能夠得到on-time和off-time時間段內(nèi)的高精度正演結(jié)果；而稀疏采樣的SAI-SA算法的正演結(jié)果雖然在on-time時間段的誤差很大，但在off-time時間段內(nèi)與密集采樣的SAI正演結(jié)果基本一致，相對密集采樣的正演結(jié)果的相對誤差小于5%.

圖11 VTEM波模型的dbz/dt響應(yīng)(a)實線為三維BDF2算法的正演結(jié)果，圓圈為采樣nton=200的SAI算法的正演結(jié)果，叉點為采樣nton=26的SAI-SA算法的正演結(jié)果;(b)為nton=26的SAI-SA算法的正演響應(yīng)相對nton=200的SAI算法的正演響應(yīng)的誤差.Fig.11 dbz/dt response of the VTEM waveform model(a)The line is result of 3D BDF2 algorithm,circle is result of SAI algorithm with nton=200,cross is result of SAI-SA algorithm with nton=26;(b)Relative error of response of SAI-SA algorithm with nton=26 to response of SAI algorithm with nton=200.

表3 SAI-SA算法和SAI算法計算VTEM波模型的時間對比
Table 3 Solution time of SAI-SA algorithm and SAI algorithm for VTEM waveform model

表3給出了密集采樣的SAI算法和稀疏采樣的SAI-SA算法的計算時間對比.由表3可知，采用密集采樣nton=200的SAI算法，由于總的系數(shù)矩陣回代次數(shù)高達3114，總的正演時間為903.9 s；采用線性不等間隔稀疏采樣nton=26的SAI-SA算法，不僅總的系數(shù)矩陣回代次數(shù)減少為436，而且減少了一次LU分解，因此總的正演時間減少為154.4 s，計算速度提高為密集采樣算法的5.9倍，在保證off-time響應(yīng)計算精度的同時，表現(xiàn)出非常好的加速效果.

2.4 不同發(fā)射波形的探測能力分析

為了比較不同發(fā)射波形的勘探能力，參考Commer和Newman(2004)，設(shè)計如圖12的三維垂直接觸帶模型.空氣電阻率為106Ωm，地表下方是厚度為50 m、電阻率為100 Ωm的覆蓋層，覆蓋層下方由兩部分不同電阻率的垂直接觸帶構(gòu)成，電阻率分別為500 Ωm和800 Ωm.垂直接觸帶中間存在一個電阻率為2 Ωm的三維復(fù)雜形狀的低阻體，沿x方向的長度為400 m，沿y方向的寬度為100 m，沿z方向的厚度為500 m，其具體位置信息如圖12所示.采用中心回線的航空裝置，發(fā)射線圈是邊長6.66 m的正八邊形，飛行高度為30 m.采用直立六面體網(wǎng)格對計算區(qū)域進行剖分，最小網(wǎng)格為10 m，網(wǎng)格擴大因子為1.36，三個方向的網(wǎng)格數(shù)為66×76×60.

SAI算法(nton=200)SAI-SA算法(nton=26)求解次數(shù)求解時間(s)求解次數(shù)求解時間(s)矩陣分解228.11 13.9矩陣回代3114833.6436116.0搜索位移量γ7.57.5求解方程(7)和(8)23110.9574.8其他23.812.2總計903.9154.4

圖12 三維垂直接觸帶模型Fig.12 3D model of a conductor at a vertical contact

為了便于比較，采用具有相同的發(fā)射能量(發(fā)射電流與時間軸圍成的面積相同)的半正弦波和梯形波形，其中半正弦波的持續(xù)時間為4 ms，發(fā)射電流為1 A；梯形波的上升沿和下降沿均為0.2 ms，持續(xù)時間為4 ms，發(fā)射電流為0.58 A.圖13為半正弦波和梯形波的off-time時段[0.01,50]ms時間范圍內(nèi)?bz/?t響應(yīng)的多測道圖.從圖13可以看出：半正弦波和梯形波的off-time時段?bz/?t響應(yīng)均能對地下低阻異常體有明顯的反映，在x=-25 m區(qū)域，對應(yīng)模型圖12中低阻異常體最淺部位置，多測道曲線呈下凹特征，顯示明顯的低阻異常.

圖13 不同發(fā)射波形off-time階段dbz/dt響應(yīng)的多測道圖(a)半正弦波；(b)梯形波.Fig.13 Profiles of off-time dbz/dt responses for different transmitting waveforms(a)Half-sine wave;(b)Trapezoid wave.

圖14為x=-25 m處正上方半正弦波與梯形波在off-time階段的?bz/?t響應(yīng)曲線，其中虛線為電阻率100 Ωm的半空間的正演響應(yīng)曲線，實線為圖12所示三維模型的正演響應(yīng)曲線，灰色為半正弦波的響應(yīng)曲線，黑色為梯形波的響應(yīng)曲線.由圖14可知，半正弦波和梯形波的響應(yīng)曲線對異常體均有明顯的反映：早期響應(yīng)相比半空間響應(yīng)的數(shù)值大，顯示為低阻異常；晚期響應(yīng)相比半空間響應(yīng)的數(shù)值小，則表現(xiàn)為高阻地層的影響.對比圖14中的曲線幅值可知，在激發(fā)能量相同時，早期的梯形波響應(yīng)相比半正弦波響應(yīng)的幅值大，因此在相同噪聲條件下，梯形波可以得到信噪比更好的早期響應(yīng)；但晚期的梯形波響應(yīng)曲線和半正弦波響應(yīng)曲線基本重合，說明發(fā)射波形對晚期響應(yīng)影響小.

圖14 x=-25 m處不同發(fā)射波形off-time階段的dbz/dt響應(yīng)對比灰色為半正弦波響應(yīng)，黑色為梯形波響應(yīng).Fig.14 Comparison of off-time dbz/dt responses at x=-25 m for different transmitting waveformsGray for half-sine wave and black for trapezoid wave.

3 結(jié)論

本文采用位移逆Krylov子空間模型降階算法結(jié)合頻譜分析理論實現(xiàn)了考慮發(fā)射波形的瞬變電磁off-time響應(yīng)的三維快速正演.利用指數(shù)求積法則，將on-time響應(yīng)和off-time響應(yīng)表示為統(tǒng)一的矩陣與向量乘積的函數(shù)形式，采用位移逆Krylov子空間模型降階算法實現(xiàn)有效求解.為了減少計算時間，從頻譜分析角度，定量化的給出了on-time時段的波形采樣規(guī)則，減少了on-time波形采樣數(shù)，從而在保證off-time響應(yīng)計算精度的同時，極大的減少了總的正演計算時間.

半正弦波模型、梯形波模型以及復(fù)雜的VTEM波模型的正演計算結(jié)果表明，對于規(guī)則變化的發(fā)射波形，比如半正弦波，采用線性等間隔采樣，能夠?qū)崿F(xiàn)好的加速效果；對于具有典型拐點的發(fā)射波形，比如梯形波和VTEM波，采用線性不等間隔采樣，能夠?qū)崿F(xiàn)更好的加速效果.三種發(fā)射波形的模型正演結(jié)果，相比密集采樣的常規(guī)正演結(jié)果(Zhou et al.,2020)，在保證off-time響應(yīng)計算精度的同時，都實現(xiàn)了5倍以上的加速，表現(xiàn)出非常穩(wěn)定的加速效果.

線性不等間隔采樣能夠最大化的實現(xiàn)加速效果，但對于復(fù)雜的發(fā)射波形，比如VTEM，需要進行手動選取采樣點，很難快速直接的得到最優(yōu)化的線性不等間隔采樣波形.未來進一步工作考慮引入有效的自適應(yīng)選取線性不等間隔采樣波形的方法.