魯凱亮, 李貅*, 戚志鵬*, 樊亞楠, 周建美,李文翰, 李賀, 張明晶, 王揚(yáng)州

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

0 引言

目前瞬變電磁反演結(jié)果為連續(xù)電阻率的空間分布，在識(shí)別電性分界面等構(gòu)造時(shí)具有一定的困難.而瞬變電磁的擬地震成像是實(shí)現(xiàn)三維成像的有效手段之一(Xue et al.，2013；戚志鵬等，2013；鐘華森等，2016).在實(shí)現(xiàn)擬地震成像之前，需要把瞬變電磁擴(kuò)散場(chǎng)轉(zhuǎn)換到虛擬波場(chǎng).只有在此基礎(chǔ)上，才能實(shí)現(xiàn)穩(wěn)定高效的瞬變電磁法擬地震解釋.

瞬變電磁擴(kuò)散場(chǎng)主要描述的是低頻電磁場(chǎng)傳播過程中的擴(kuò)散和感應(yīng)特征，具有較強(qiáng)的體積效應(yīng)，因而對(duì)電性界面的分辨能力較差(Khan et al.，2018；底青云等, 2019; Di et al.，2020；Xue et al.，2019；Li et al.，2018)；而波動(dòng)場(chǎng)能夠比較好地描述波場(chǎng)的一系列特征，所以基于波動(dòng)方程的擬地震成像可以很好地刻畫地層的界面信息(Qi et al.，2015).而直接將瞬變電磁的擴(kuò)散場(chǎng)引入至波場(chǎng)的處理流程中是不合適的，需要通過數(shù)學(xué)手段進(jìn)行變換.Lee等(1989)給出了虛擬波場(chǎng)到擴(kuò)散場(chǎng)的積分公式，實(shí)現(xiàn)了穩(wěn)定的正變換；隨后，Lee和Xie(1993)根據(jù)波場(chǎng)正變換的積分公式，嘗試進(jìn)行波場(chǎng)反變換的計(jì)算，他們指出波場(chǎng)反變換需要求解第一類的Fredholm積分方程，屬于嚴(yán)重的不適定性問題.之后，Gershenson(1997)分別使用了直接求解法、反褶積法和拉氏正弦變換法計(jì)算虛擬波場(chǎng)，并對(duì)這三種方法的優(yōu)劣進(jìn)行了比較.陳本池等(1999)使用正則化方法進(jìn)行虛擬波場(chǎng)反變換的求解，并對(duì)層狀模型進(jìn)行試算，結(jié)果表明從擴(kuò)散場(chǎng)轉(zhuǎn)換至虛擬波場(chǎng)后，在電性分界面處會(huì)出現(xiàn)“反射”和“折射”等現(xiàn)象，這一發(fā)現(xiàn)為電磁場(chǎng)的擬地震處理解釋提供了重要的理論依據(jù).近些年，李貅等(2005)，朱宏偉等(2010)，張軍等(2011)將正則化算法應(yīng)用于波場(chǎng)反變換中，取得較好的結(jié)果；戚志鵬等(2013)通過使用預(yù)條件正則化共軛梯度法，實(shí)現(xiàn)了全時(shí)域的虛擬波場(chǎng)反變換算法，通過引入地震勘探中的脈沖反褶積方法，提高了虛擬波場(chǎng)的縱向分辨率.樊亞楠等(2019)使用掃時(shí)波場(chǎng)變換算法，進(jìn)一步提高了波場(chǎng)反變換的穩(wěn)定性.雖然以上方法都取得了長(zhǎng)足的發(fā)展，但是虛擬波場(chǎng)反變換的精度和效率還有待進(jìn)一步提高.

在計(jì)算虛擬波場(chǎng)反變換時(shí)，將其離散化后的線性方程組具有高度的病態(tài)性.近年來(lái)，精細(xì)積分法作為求解高度病態(tài)線性方程組的方法，在各個(gè)領(lǐng)域得到了日益廣泛的應(yīng)用.該方法通過把求逆矩陣的過程轉(zhuǎn)化為一個(gè)求解無(wú)窮積分的過程，大大降低了求解病態(tài)線性方程組的難度.精細(xì)積分法(鐘萬(wàn)勰，1994)在計(jì)算矩陣指數(shù)函數(shù)時(shí)，具有效率高、精度高等優(yōu)點(diǎn)，已經(jīng)在各類工程問題中取得了較大進(jìn)展(Gu et al.，2001)，而目前將精細(xì)積分法引入到虛擬波場(chǎng)反變換的工作鮮有提及，本文將對(duì)這一方法的引入初做嘗試.

1 波場(chǎng)變換基本原理

1.1 基本原理

在導(dǎo)電介質(zhì)中，忽略位移電流，瞬變電磁場(chǎng)滿足擴(kuò)散方程.為了不失一般性，取f(x,y,z,t)為瞬變電磁場(chǎng)的電磁分量.

根據(jù)文獻(xiàn)(Lee et al.，1989)可知，瞬變電磁擴(kuò)散場(chǎng)轉(zhuǎn)換到虛擬波場(chǎng)的表達(dá)式為

(1)

式(1)為第一類Fredholm型積分方程，從擴(kuò)散場(chǎng)求解虛擬波場(chǎng)是一個(gè)典型的不適定問題.對(duì)式(1)進(jìn)行離散，可以得到：

(2)

其中：

(3)

離散后得到的線性方程組是高度病態(tài)的，且隨著階數(shù)的增加，矩陣的條件數(shù)急劇增大(戚志鵬等，2013)，本文采用精細(xì)積分法求解(1)式.

圖1 核函數(shù)曲線圖Fig.1 Kernel function graph

1.2 精細(xì)積分法

考慮病態(tài)線性方程組

Ax=b,

(4)

式中A為n階正定矩陣，b為n維實(shí)向量，記H=-A、r=b，則式(4)可以轉(zhuǎn)化為

Hx+r=0,

(5)

考慮如下一階微分方程

(6)

式(6)是瞬態(tài)熱傳導(dǎo)方程解的形式；穩(wěn)態(tài)的熱傳導(dǎo)方程的解具有式(5)的形式，而式(6)的解可以表示為

上式中，H是負(fù)定矩陣，r為常向量.從上式易知，當(dāng)時(shí)間t′→∞時(shí)，方程(6)趨向于方程(5)，且容易看出exp(Ht′)→0，所以式(7)中的積分項(xiàng)就逼近式(5)的解，即

(9)

所以，式(4)的解就可以寫為(富明慧，2018)

(10)

從上面的分析我們可以得出一個(gè)結(jié)論，病態(tài)線性方程組的求解可以轉(zhuǎn)化為一個(gè)求積分的穩(wěn)定過程.

(11)

1.3 迭代終止準(zhǔn)則

前文主要介紹了精細(xì)積分法，這一節(jié)給出迭代終止條件.理論上講，迭代次數(shù)越多，結(jié)果越精確，但是當(dāng)積分區(qū)間達(dá)到一定數(shù)值范圍后，隨著計(jì)算誤差的累積以及矩陣的高病態(tài)性，反而會(huì)導(dǎo)致精度隨積分區(qū)間的增加迅速下降，所以有必要在虛擬波場(chǎng)計(jì)算的過程中添加一個(gè)迭代終止條件.

目前有如下四種迭代終止條件：(1)對(duì)迭代次數(shù)設(shè)置一個(gè)迭代上限m；(2)終止條件設(shè)置為errk<ε，ε為給定的誤差容限，errk=‖xk+1-xk‖；(3)當(dāng)?shù)鷼埐钣袉握{(diào)下降和單調(diào)上升的兩個(gè)計(jì)算過程時(shí)，理想的迭代終點(diǎn)在殘差的拐點(diǎn)處，即當(dāng)errk/errk-1≥1時(shí)迭代終止；(4)在第k次計(jì)算過程中，往后逐次取n個(gè)迭代步，如果這n個(gè)迭代步的相對(duì)殘差都在增長(zhǎng)，默認(rèn)第k次的解為最優(yōu)解.

對(duì)于第一種終止條件，更多的是以經(jīng)驗(yàn)為主，要么沒有達(dá)到精度要求，要么增加不必要的計(jì)算量；第二種終止條件在應(yīng)用于病態(tài)方程組時(shí)經(jīng)常失效(富明慧和李勇息，2018)；第三種終止條件太過理想，在實(shí)際計(jì)算中的殘差曲線往往并非單調(diào)，而且經(jīng)常有波動(dòng).綜上，本文選取第四種終止條件.

設(shè)

(12)

式中n為選擇的檢測(cè)窗口，為整數(shù)，一般情況下取2，最大不超過10(富明慧和李勇息，2018).

2 高斯脈沖的波場(chǎng)反變換

2.1 算法精度

首先研究一個(gè)簡(jiǎn)單脈沖，令u(τ)=1

從圖2的結(jié)果可以看出，對(duì)于簡(jiǎn)單脈沖u(τ)=1，本文的方法與正則化共軛梯度法都能得到相對(duì)較好的結(jié)果，但是相較而言，本文結(jié)果的精度比較高.

圖2 波場(chǎng)反變換結(jié)果和相對(duì)誤差 (a) 本文方法結(jié)果； (b) 本文方法相對(duì)誤差； (c) PRCG結(jié)果； (d) PRCG相對(duì)誤差.Fig.2 Inverse wave field transformation results and relative error (a) The present method; (b) The relative error in paper； (c) The PRCG method; (d) The relative error by PRCG.

接下來(lái)本文使用單個(gè)高斯脈沖和組合高斯脈沖的方式來(lái)得到更加復(fù)雜的虛擬波場(chǎng).令u(τ)為高斯脈沖，它的表達(dá)式如下：

(13)

式中，a表示高斯脈沖的峰值大小，c表示高斯脈沖的寬度，τ0表示高斯脈沖的峰值位置，h0表示高斯脈沖在垂直方向的偏移量.通過給定不同的a、c、τ0和h0的值就可以得到不同樣式的高斯脈沖，而且還可以對(duì)上式進(jìn)行組合得到更加復(fù)雜的高斯脈沖序列.

使用精細(xì)積分法分別計(jì)算了D型、G型、Q型、A型、H型和K型模型的虛擬波場(chǎng)(圖3)，并與前人研究結(jié)果進(jìn)行了精度對(duì)比.從相對(duì)誤差圖上可以看出，本文方法的相對(duì)誤差較小，均在4%以下，而正則化共軛梯度法的相對(duì)誤差隨著模型的復(fù)雜而逐漸增大，最大相對(duì)誤差達(dá)到了將近50%，從這里可以看出本文方法具有較高的精度.

圖3 不同型模型波場(chǎng)反變換結(jié)果和相對(duì)誤差 (A) D型模型； (B) G型模型； (C) Q型模型; (D) A型模型; (E) H型模型; (F) K型模型.(a1—f1) 本文方法結(jié)果； (a2—f2) 本文方法相對(duì)誤差； (a3—f3) PRCG結(jié)果；(a4—f4) PRCG相對(duì)誤差.Fig.3 The inverse transformation results and relative errors of thedifferent model (A) D model; (B) Gmodel; (C) Q model; (D) A model; (E) H model; (F) K model; (a1—f1) The present method; (a2—f2) The relative error in paper; (a3—f3) The PRCG method; (a4—f4) The relative error by PRCG.

2.2 模型分辨率分析

令第一個(gè)脈沖的脈寬為0.02，振幅為2；第二個(gè)脈沖的脈寬為0.04，振幅為1.6；使第一個(gè)脈沖的波峰位置為τ01=0.03，改變第二個(gè)波峰τ02的位置，使其分別為0.04和0.08.結(jié)果如圖4所示：

本節(jié)研究了本文方法對(duì)Q型、A型、H型和K型模型的分辨能力(圖4).固定第一個(gè)脈沖的波峰位置為τ01=0.03，改變第二個(gè)波峰τ02的位置，使其分別為0.08和0.04，使得兩個(gè)波峰的位置越來(lái)越近.從上述結(jié)果可以看出，本文方法具有較好的分辨能力，而隨著兩個(gè)波峰的位置越來(lái)越近，正則化共軛梯度法的分辨能力越來(lái)越低.

圖4 波場(chǎng)反變換結(jié)果 (A) Q型模型； (B) A型模型; (C) H型模型; (D) K型模型; (a1—d1) 本文方法結(jié)果τ02=0.08；(a2—d2)PRCG結(jié)果τ02=0.08； (a3—d3) 本文方法結(jié)果τ02=0.04； (a4—d4) PRCG結(jié)果τ02=0.04.Fig.9 The inverse transformation results and relative errors of the different model (A) Q model； (B) A moel; (C) H model; (D) K model; (a1—d1) The present method τ02=0.08; (a2—d2) The PRCG method τ02=0.08; (a3—d3) The present method τ02=0.04; (a4—d4) The PRCG method τ02=0.04.

2.3 加噪試驗(yàn)

分別對(duì)D型、G型、Q型、A型、H型和K型模型的擴(kuò)散場(chǎng)中分別加入2%和5%的噪聲，并計(jì)算對(duì)應(yīng)的虛擬波場(chǎng).計(jì)算結(jié)果如圖5.

圖5 不同模型加噪波場(chǎng)反變換結(jié)果 (A) D型模型; (B) G型模型; (C) Q型模型; (D) A型模型; (E) H型模型; (F) K型模型; (a1—f1) 本文方法結(jié)果SNR=50∶1； (a2—f2) PRCG結(jié)果SNR=50∶1； (a3—f3) 本文方法結(jié)果SNR=20∶1； (a4—f4) PRCG結(jié)果SNR=20∶1.Fig.5 The inverse transformation results and relative errors of the different model (A) D model; (B) G model; (C) Q model; (D) A model; (E) H model; (F) K model; (a1—f1) The present method SNR=50∶1; (a2—f2) The PRCG method SNR=50∶1; (a3—f3) The present method SNR=20∶1; (a4—f4) The PRCG method SNR=20∶1.

圖5各模型的擴(kuò)散場(chǎng)信號(hào)中加入2%和5%的噪聲，測(cè)試了本文方法和正則化共軛梯度法的抗噪性.從圖中可以看出，本文的方法抗噪能力比較強(qiáng)，對(duì)于含有2%和5%噪聲的擴(kuò)散場(chǎng)，依舊能求取較為準(zhǔn)確的虛擬波場(chǎng)；對(duì)于含有2%噪聲的擴(kuò)散場(chǎng)，PRCG可以較為準(zhǔn)確地求取對(duì)應(yīng)的虛擬波場(chǎng)，對(duì)于含有5%噪聲的擴(kuò)散場(chǎng)，求得的虛擬波場(chǎng)結(jié)果較差.當(dāng)信噪比超過5%時(shí)，兩種方法的虛擬波場(chǎng)結(jié)果波動(dòng)劇烈，與真值差距較大.因此對(duì)于野外信號(hào)的采集必須控制噪聲在5%以下，并且對(duì)于局部不光滑數(shù)據(jù)不能直接求解，需要對(duì)數(shù)據(jù)進(jìn)行平滑降噪后才能計(jì)算虛擬波場(chǎng).

3 三維模型計(jì)算

為了進(jìn)一步驗(yàn)證本文使用的波場(chǎng)反變換算法的有效性，本節(jié)設(shè)計(jì)如下三維模型計(jì)算虛擬波場(chǎng)，并與正則化共軛梯度法進(jìn)行比較.圖6為均勻半空間中含有兩個(gè)低阻薄板狀異常體的模型示意圖.均勻半空間的電阻率為200 Ωm，上層異常體上表面距離地表150 m，幾何尺寸為100 m×100 m×20 m，電阻率為5 Ωm；下層異常體上表面距地表250 m，幾何尺寸為100 m×100 m×30 m，電阻率為1 Ωm.接地導(dǎo)線源的長(zhǎng)度為200 m，最短偏移距為100 m，發(fā)射電流為10A；布設(shè)13條測(cè)線，測(cè)線的線距為10 m，起始測(cè)線位置為y=-80 m，終止測(cè)線位置為y=50 m，每條測(cè)線長(zhǎng)度為500 m，點(diǎn)距為25 m，每條測(cè)線21個(gè)測(cè)點(diǎn)，接收的電磁場(chǎng)信號(hào)為Ex.本文三維正演采用擬態(tài)有限體積法(Zhou et al.，2018).

圖6 三維模型示意圖 (a) XOZ平面圖； (b) XOY平面圖.Fig.6 Schematic diagram of three-dimensional model (a) XOZ plan; (b) XOY plan.

圖7為波場(chǎng)反變換結(jié)果，本文在這里只展現(xiàn)三條測(cè)線的波場(chǎng)反變換結(jié)果，分別是測(cè)線y=-70 m，測(cè)線y=0 m，測(cè)線y=30 m.選取這三條測(cè)線的原因是測(cè)線y=-70 m能反映出深部異常信息；測(cè)線y=0 m可以反映淺部和深部的異常信息；而測(cè)線y=50 m可以反映出淺部異常信息，其他測(cè)線的結(jié)果與這三條測(cè)線的結(jié)果類似，所以在此不進(jìn)行展示.

圖7 本文方法波場(chǎng)反變換結(jié)果 (a) 測(cè)線y=-70 m； (b) 測(cè)線y=0 m； (c) 測(cè)線y=30 m.Fig.7 Results of wave field transformation in this paper (a) Line y=-70 m; (b) Line y=0 m; (c) Line y=30 m.

圖7a是測(cè)線y=-70 m的波場(chǎng)反變換結(jié)果，從圖中可以看出，除了直達(dá)波以外，還有深部異常體的上界面信息，由于瞬變電磁的高頻成分被地層迅速吸收，所以導(dǎo)致對(duì)深部異常體的分辨率較差；圖7b是測(cè)線y=0 m的波場(chǎng)反變換結(jié)果，從圖中可以看出，虛擬波場(chǎng)可以很好地反映出淺部異常體的上下界面，波形的同相軸比較明顯；對(duì)于深部異常體，圖7b也能反映異常體的上界面信息；圖7c是測(cè)線y=50 m的波場(chǎng)反變換結(jié)果，能夠反映出淺部異常體的界面信息.從三幅圖中可以看出，使用本文方法計(jì)算的結(jié)果，虛擬波場(chǎng)比較平滑，在整個(gè)時(shí)間段，無(wú)意義的波動(dòng)起伏也比較少.

圖8為使用PRCG方法得到的波場(chǎng)反變換結(jié)果.可以看出，只有測(cè)線y=-70 m的結(jié)果比較好；測(cè)線y=0 m的結(jié)果主要反映的是深部異常體上界面的部分信息；測(cè)線y=50 m對(duì)淺層異常體界面反映較好，但是由于PRCG方法的精度較低，出現(xiàn)了多余的虛擬波場(chǎng).

圖8 PRCG波場(chǎng)反變換結(jié)果 (a) 測(cè)線y=-70 m； (b) 測(cè)線y=0 m； (c) 測(cè)線y=30 m.Fig.8 Results of wave field transformation by PRCG (a) Line y=-70 m; (b) Line y=0 m; (c) Line y=30 m.

圖9是本文方法的虛擬波場(chǎng)三維切片圖，共有四張切片，對(duì)應(yīng)的位置分別為y=-50 m，y=40 m，x=120 m和x=240 m.從圖中可以看出，位置(1)和位置(2)處可以很清楚地反映異常體在xoy平面的位置，位置(3)和位置(4)處也能較好地反映出深部異常體在xoy平面的位置.通過計(jì)算，三維模型的虛擬波場(chǎng)與設(shè)計(jì)模型相符，可以進(jìn)一步說明本文方法的有效性和準(zhǔn)確性.

圖9 波場(chǎng)反變換結(jié)果三維切片圖Fig.9 Three-dimensional slice map of wavefield transformation results

4 實(shí)測(cè)數(shù)據(jù)分析

上文利用各種算例與模型對(duì)本文方法進(jìn)行了驗(yàn)證，本節(jié)對(duì)甘肅省某煤田采空區(qū)實(shí)測(cè)數(shù)據(jù)進(jìn)行處理.采區(qū)地層大部分被第四系黃土覆蓋，基本隱蔽,部分基巖裸露.地層從上到下由上白堊統(tǒng)、下白堊統(tǒng)、上侏羅統(tǒng)、中下侏羅統(tǒng)及上三疊統(tǒng)組成，中下侏羅統(tǒng)為主要煤系地層，上三疊統(tǒng)為煤系基底.采區(qū)主要含煤地層由石英細(xì)粗砂巖、砂礫巖、炭質(zhì)泥巖、炭質(zhì)粉砂巖組成.采區(qū)地表南部為丘陵，部分丘陵的白堊紀(jì)巖層出露.地表北部開闊平緩，為農(nóng)田及低矮丘陵，小沖溝發(fā)育.地表中部、北部有兩條季節(jié)性沙河.

野外工作使用V8多功能電法儀器，發(fā)射系統(tǒng)采用地面長(zhǎng)導(dǎo)線源，發(fā)射電流15 A，基頻8.3 Hz，共15條測(cè)線；每條測(cè)線400 m，點(diǎn)距10 m，共41個(gè)測(cè)點(diǎn)，沿南北向布設(shè)，接收Ex電磁信號(hào).本文選取380號(hào)測(cè)線數(shù)據(jù)進(jìn)行視電阻率成像(張瑩瑩等，2015)與波場(chǎng)反變換，其余測(cè)線不在此文贅述.

由圖10a視電阻率剖面圖可以看出，隨著深度的增加，視電阻率值總體呈現(xiàn)先減小后增大的趨勢(shì).在測(cè)線的280～320 m、340～400 m、440～500 m處，高程在800～1000 m之間，有三個(gè)低阻異常區(qū)，推斷可能為煤礦采空區(qū)塌陷，導(dǎo)致裂隙發(fā)育形成富水區(qū)域，從而形成低阻異常區(qū)，與已知的鉆孔資料較為吻合(已知鉆孔的煤線高程為1042.93 m).由圖10b可以看出，虛擬波場(chǎng)的同相軸與視電阻率的趨勢(shì)基本一致，整體呈南低北高，與煤層底板等值線圖較吻合(圖11).在測(cè)線的280～320 m、340～380 m、400～430 m，430～500 m處，由于煤層被開采，地層錯(cuò)斷較為嚴(yán)重，在虛擬波場(chǎng)圖上表現(xiàn)出橫向不連續(xù)，間斷跳躍的特征.本文通過與已知資料對(duì)比，本文的計(jì)算結(jié)果與已知資料相吻合，說明本文方法在識(shí)別電性界面方面有較好的效果.

圖10 380號(hào)測(cè)線視電阻率圖(a) 與波場(chǎng)反變換結(jié)果(b)Fig.10 Apparent resistivity profile and pseudo wave field of line 380

圖11 工區(qū)煤層底板等值線圖Fig.11 Contour map of coal seam floor in working area

5 結(jié)論

本文通過使用精細(xì)積分法，成功實(shí)現(xiàn)了從瞬變電磁擴(kuò)散場(chǎng)到虛擬波場(chǎng)全時(shí)域穩(wěn)定的反變換算法.在保證計(jì)算精度的前提下，使用指數(shù)級(jí)的積分步長(zhǎng)可以顯著提高計(jì)算效率.通過設(shè)置合理的終止迭代條件，使得本文的方法適用性更廣.

實(shí)現(xiàn)波場(chǎng)反變換后，本文對(duì)算法進(jìn)行了驗(yàn)證.通過對(duì)典型地電模型的計(jì)算，證明了本文的方法具有較好的精度(得到虛擬波場(chǎng)值與理論值的誤差小于4%)、分辨能力以及抗噪性，可以說明本文方法適用于復(fù)雜介質(zhì)模型.通過計(jì)算三維模型的虛擬波場(chǎng)值，可以看出虛擬波場(chǎng)對(duì)淺層異常體的分辨率較好，隨著深度增加，高頻成分很快被吸收，低頻成分占主導(dǎo)，波場(chǎng)幅值減小，子波寬度增加，分辨率逐漸降低.最后通過對(duì)實(shí)測(cè)數(shù)據(jù)的處理，進(jìn)一步說明了本文方法的可行性，為后續(xù)的瞬變電磁構(gòu)造成像方法打下了良好的基礎(chǔ).

從文中可以看出，當(dāng)擴(kuò)散場(chǎng)包含的噪聲越來(lái)越大時(shí)，波場(chǎng)反變換的結(jié)果與實(shí)際結(jié)果相差也就越大.因此，在后續(xù)工作中可以引入窗口掃描、相關(guān)疊加的方法進(jìn)一步提高該方法的抗噪性.

致謝感謝中山大學(xué)富明慧教授在作者研究過程中給予指導(dǎo)，同時(shí)感謝論文評(píng)審專家和編輯提出的寶貴修改意見

附錄A

(A1)

(A2)

類似地，我們可以得到：

(A3)

此時(shí)對(duì)兩邊同時(shí)乘以b，得：

(A4)

當(dāng)k→∞時(shí)，我們將得到精確解xk.但是精細(xì)積分的步長(zhǎng)ζ一般都比較小(富明慧和李勇息，2018)，如果對(duì)積分區(qū)間[0,∞)采用線性等間隔離散，計(jì)算量將大到無(wú)法接受.所以本文參考前人工作(鐘萬(wàn)勰，1994)，采用指數(shù)增長(zhǎng)的積分步長(zhǎng).

對(duì)于矩陣指數(shù)函數(shù)，鐘萬(wàn)勰提出了一種精確且高效的計(jì)算方法，在各類偏微分方程中獲得了廣泛應(yīng)用.與上文類似，設(shè)

(A5)

對(duì)上式兩邊同時(shí)乘以b，得：

x1=(I+exp(-Aζ))x0.

(A6)

由于ζ是極小量，在計(jì)算矩陣指數(shù)函數(shù)時(shí)，可使用Taylor展開的前幾項(xiàng)：

(A7)

對(duì)式(A7)在區(qū)間[0,ζ]積分，得

(A8)

綜上，x的初始值x0=F(ζ)b.

令積分步長(zhǎng)以2k-1ζ變化，則有

k=1,2,…

(A9)

且當(dāng)2kζ→Fk時(shí)，有

(A10)

(A11)

上式中Tk即為矩陣指數(shù)函數(shù)的值，這樣計(jì)算矩陣指數(shù)函數(shù)的好處在于：在計(jì)算過程中當(dāng)k很小時(shí)，能將大量和小量分離，避免小量與大量相加時(shí)出現(xiàn)數(shù)據(jù)“湮沒”，保證計(jì)算精度(富明慧和李勇息，2018).