漢克爾自適應(yīng)算法在電磁計(jì)算中的應(yīng)用

蔡 盛，柳建新，湯文武，曹創(chuàng)華，鄧小康

（1.中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司，武漢 430063；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083）

0 引言

在地球物理中，尤其是電磁法計(jì)算中，經(jīng)常會(huì)遇到貝塞爾函數(shù)的積分問(wèn)題。該積分用常規(guī)方法計(jì)算，比較復(fù)雜。十九世紀(jì)七十年代，荷蘭學(xué)者把數(shù)值線性濾波方法引入地球物理學(xué)領(lǐng)域［1］，很好地解決了這個(gè)問(wèn)題。之后，快速漢克爾濾波系數(shù)如雨后春筍，發(fā)展很快，出現(xiàn)了很多短濾波系數(shù)［2－4］。由于計(jì)算機(jī)的發(fā)展以及計(jì)算中對(duì)精度要求的不斷提高，緊接著又出現(xiàn)了很多高精度的長(zhǎng)濾波系數(shù)［5－7］。

現(xiàn)有的漢克爾長(zhǎng)濾波系數(shù)，在一個(gè)適當(dāng)?shù)姆秶畠?nèi)，精度都是很高的，相對(duì)誤差可控制在10－5以上，甚至有時(shí)精度達(dá)到了10－10以上。在地球物理中，實(shí)際應(yīng)用時(shí)需要考慮計(jì)算精度和計(jì)算需要的時(shí)間，而且這兩者同等重要。在一個(gè)可以接受的計(jì)算精度下，控制住計(jì)算時(shí)間，極其重要。對(duì)于一些精度較高的長(zhǎng)濾波器，如果能在滿(mǎn)足精度要求的前提下，應(yīng)用一種自適應(yīng)算法，得到需要的精度即結(jié)束計(jì)算過(guò)程，將節(jié)省大量的計(jì)算時(shí)間。在地球物理數(shù)值計(jì)算的過(guò)程中，經(jīng)常會(huì)遇到相關(guān)核函數(shù)的積分問(wèn)題，在數(shù)值濾波中，積分最終歸結(jié)為褶積值的計(jì)算，所以相關(guān)核函數(shù)積分時(shí)，褶積值的重復(fù)計(jì)算量很大。把褶積值存儲(chǔ)起來(lái)為之后的計(jì)算調(diào)用，結(jié)合自適應(yīng)算法在數(shù)值積分時(shí)，會(huì)節(jié)省一定的時(shí)間。

1 濾波系數(shù)的選取

在地球物理中，電磁法一維正演過(guò)程包括直流電、可控源、瞬變電磁，都含有貝塞爾函數(shù)積分式［8］，可表示為式（1）。

式中 f（λ）為積分核函數(shù)；Ji（λr）為貝塞爾函數(shù)，電磁計(jì)算中一般是零階和一階。由快速漢克爾理論，可表示為式（2）的形式：

其中 λi為采樣點(diǎn)的位置；Wi為濾波權(quán)重值。積分的精度，決定于權(quán)重值Wi。選取五種表現(xiàn)很好的快速漢克爾長(zhǎng)濾波系數(shù)如表1所示。

表1 選擇的五種高精度長(zhǎng)濾波系數(shù)Tab.1 The five high－precision long filter coefficients selected

選擇兩個(gè)有解析解的貝塞爾函數(shù)積分式，即取式（1）中積分核函數(shù)f（λ）零階時(shí)為λe－λ，用五種系數(shù)求取其各自的相對(duì)誤差，如圖1所示；一階時(shí)f（λ）為λe－2λ，用五種系數(shù)求取其各自的相對(duì)誤差，如圖2所示。用相對(duì)誤差而不用絕對(duì)誤差是因?yàn)橐粋€(gè)小的絕對(duì)誤差可能會(huì)誤導(dǎo)我們（例如如果理論計(jì)算值很小，而實(shí)際計(jì)算值有很大的誤差，也可能導(dǎo)致其絕對(duì)誤差比較?。?，因此相對(duì)誤差是更好的評(píng)價(jià)方法。

圖1、圖2中橫坐標(biāo)是求值區(qū)間，也就是公式（1）和公式（2）中的實(shí)參數(shù)r的對(duì)數(shù)，即log（r）；縱坐標(biāo)是相對(duì)誤差。由圖1和圖2可以看出，在零階積分中，801點(diǎn)系數(shù)和120點(diǎn)系數(shù)表現(xiàn)最好，而在一階積分中，801點(diǎn)系數(shù)和140點(diǎn)系數(shù)表現(xiàn)最好；801點(diǎn)系數(shù)在全求值區(qū)間表現(xiàn)比較穩(wěn)定，而120（J0）和140（J1）點(diǎn)系數(shù)受求值區(qū)間的影響較大。就計(jì)算精度而言，考慮到電磁法中進(jìn)行數(shù)值計(jì)算中的一般性，認(rèn)為801點(diǎn)系數(shù)表現(xiàn)最好，而且由于801點(diǎn)系數(shù)零階和一階具有相同的系數(shù)個(gè)數(shù)、相同的采樣坐標(biāo)，有利于統(tǒng)一進(jìn)行計(jì)算。

圖1 零階時(shí)漢克爾積分精度對(duì)比Fig.1 The relative errors of the five J0filters

2 相關(guān)核函數(shù)的自適應(yīng)算法

算法建立在已知的高精度濾波系數(shù)之上，選擇801點(diǎn)的濾波系數(shù)在一般情況下表現(xiàn)杰出，而且J0和J1濾波系數(shù)的采樣點(diǎn)的坐標(biāo)都相同，有利于算法的利用。

2.1 自適應(yīng)計(jì)算

把濾波系數(shù)按橫坐標(biāo)劃為三個(gè)部分，①固定區(qū)域；②固定區(qū)域右邊部分；③固定區(qū)域左邊部分。固定區(qū)域需要包括擾動(dòng)比較劇烈的部分［9］。在801點(diǎn)的濾波系數(shù)中，如圖3、圖4所示，規(guī)定振蕩劇烈的含有92個(gè)濾波權(quán)重值范圍為固定范圍，即從256點(diǎn)到第348個(gè)點(diǎn)。在這個(gè)區(qū)域里面，從圖3、圖4中可以看出，0階和1階的濾波器絕對(duì)值都快速趨近于“0”。自適應(yīng)計(jì)算過(guò)程的第一步，即計(jì)算固定范圍內(nèi)的褶積值；第二步，自適應(yīng)褶積從固定范圍的右邊，從左至右展開(kāi)；第三步，計(jì)算固定范圍的左邊，從右至左展開(kāi)。當(dāng)固定范圍的兩邊褶積的總和的絕對(duì)值小于或者等于截?cái)鄻?biāo)準(zhǔn)時(shí)，自適應(yīng)過(guò)程就停止了。在固定范圍之外的自適應(yīng)褶積的截?cái)鄻?biāo)準(zhǔn)，由下式確定：

復(fù)數(shù)和的實(shí)部分和虛部分，必須獨(dú)立地滿(mǎn)足這個(gè)標(biāo)準(zhǔn)，才能接受截?cái)?。截?cái)鄻?biāo)準(zhǔn)基本近似等于褶積誤差，但是不等于真正的積分誤差。

圖2 一階時(shí)漢克爾積分精度對(duì)比Fig.2 The relative errors of the five J1filters

圖3 零階801點(diǎn)濾波系數(shù)權(quán)重值截?cái)嗍疽鈭DFig.3 The truncated diagram of the 801 J0 filter coefficients

在實(shí)際中遇到的大多核函數(shù)，在固定范圍的兩側(cè)只需要少許的褶積值即可達(dá)到所需要的計(jì)算精度。然而，對(duì)于一些緩慢衰減或者振蕩的，核函數(shù)可能需要在一個(gè)較大范圍的橫坐標(biāo)之內(nèi)進(jìn)行濾波計(jì)算，這也是在一個(gè)很大范圍提供濾波權(quán)重值的原因。

2.2 相關(guān)核函數(shù)的計(jì)算

所謂相關(guān)核函數(shù)的問(wèn)題，即幾個(gè)核函數(shù)能夠用簡(jiǎn)單的代數(shù)關(guān)系聯(lián)系起來(lái)，那這樣的變換就是相關(guān)的。例如，在計(jì)算不同的循環(huán)結(jié)構(gòu)的電磁耦合率的大地二層模型時(shí)，需要計(jì)算一個(gè)零階和兩個(gè)一階漢克爾變換求值問(wèn)題，而它們的核函數(shù)幾乎是相同的。針對(duì)復(fù)變的相關(guān)核函數(shù)的0階或者1階快速高精度的漢克爾變換積分問(wèn)題，為了避免對(duì)核函數(shù)進(jìn)行重復(fù)計(jì)算，在第一次計(jì)算核函數(shù)時(shí)，應(yīng)先將計(jì)算值保存在計(jì)算機(jī)內(nèi)存中，在下一次計(jì)算相關(guān)核函數(shù)時(shí)，就可以直接調(diào)用內(nèi)存中的值，這將避免重復(fù)計(jì)算時(shí)的時(shí)間消耗。

如果一個(gè)核函數(shù)為f（λ），它的相關(guān)核函數(shù)可以表示成λi［f（λ）］j。相關(guān)核函數(shù)褶積運(yùn)算過(guò)程：首先把濾波器的權(quán)重值放在一個(gè)DATA數(shù)組里面，在直接核函數(shù)計(jì)算時(shí)，對(duì)DATA數(shù)組進(jìn)行一個(gè)初始的調(diào)用（NEW＝1），計(jì)算的核函數(shù)的值和相對(duì)應(yīng)的參數(shù)保存在計(jì)算機(jī)內(nèi)存的一個(gè)公共模塊中。當(dāng)初始調(diào)用結(jié)束之后，公共存儲(chǔ)區(qū)可以因?yàn)橄嚓P(guān)核函數(shù)之間的幾何關(guān)系進(jìn)行修改。然后，任意階的相關(guān)漢克爾變換在隨后的調(diào)用（NEW＝0）中將會(huì)獲得存儲(chǔ)區(qū)的值，不需要重復(fù)計(jì)算相關(guān)核函數(shù)的褶積值，除非當(dāng)前核函數(shù)的特性需要一些另外的計(jì)算。而這些另外的核函數(shù)褶積值，如果需要，也可以計(jì)算并且自動(dòng)添加到公共模塊區(qū)以備進(jìn)一步使用。程序流程簡(jiǎn)圖如圖5所示。

圖4 一階801點(diǎn)濾波系數(shù)權(quán)重值截?cái)嗍疽鈭DFig.4 The truncated diagram of the 801 J1 filter coefficients

圖5 相關(guān)核函數(shù)積分計(jì)算流程簡(jiǎn)圖Fig.5 The calculation diagram of the related kernels

3 實(shí)例計(jì)算與結(jié)果分析

使用一組相關(guān)核函數(shù)的漢克爾解析式，來(lái)測(cè)試相關(guān)核函數(shù)自適應(yīng)算法的精度情況。這組貝塞爾函數(shù)積分式子，都是有解析解的，可以得到算法濾波值對(duì)于精確解的相對(duì)誤差［10］。

五個(gè)貝塞爾函數(shù)積分式如下：

上面式（3）－式（7）五個(gè)積分式的核函數(shù)，分別用C1、C2、C3、C4、C5表示?？筛鶕?jù)代數(shù)關(guān)系表示成式（8）：如果用式C2（λ）表示式C3（λ），將會(huì)比式C2（λ）本身還更復(fù)雜，所以認(rèn)為方程（4）跟方程（5）是不相關(guān)的。計(jì)算結(jié)果如表2所示，其中實(shí)參數(shù)為積分式中的變量r，公差因子是在程序中設(shè)定的，但是不等于積分誤差。NEW表示相關(guān)調(diào)用情況，NEW＝1表示初始計(jì)算；NEW＝0表示相關(guān)調(diào)用。由表2可知，隨著實(shí)參數(shù)的變化計(jì)算的精度是不一樣的，這也與用濾波系數(shù)直接計(jì)算的相吻合。隨著公差因子的減小，相對(duì)誤差也是逐漸減小的，所以為了獲得較好的精度，可以設(shè)置公差因子為零。對(duì)于五個(gè)積分式的計(jì)算，精度最低的時(shí)候，相對(duì)誤差也在10－3以下，一般情況都在10－6以下，精度能夠滿(mǎn)足需求。

表2 自適應(yīng)算法計(jì)算結(jié)果表Tab.2 The calculation results of the adaptive algorithm

上面五個(gè)貝塞爾函數(shù)積分式的核函數(shù)都比較簡(jiǎn)單。當(dāng)要計(jì)算的相關(guān)核函數(shù)不是簡(jiǎn)單基本函數(shù)時(shí)，使用自適應(yīng)算法對(duì)于時(shí)間的節(jié)省將非常明顯。對(duì)于上面提到的不同循環(huán)結(jié)構(gòu)的相互耦合比的計(jì)算，對(duì)于水平M層大地模型，忽略位移電流的影響，核心計(jì)算就是下面的三個(gè)漢克爾變換積分式：

R（D，g）是一個(gè)遞歸復(fù)函數(shù)，表示M層模型參數(shù)，以式（10）表示：

其中 V1、F1由遞歸得到，其他各參數(shù)如下：

由于五種普通的循環(huán)配置，即水平共面循環(huán)、垂直循環(huán)、垂直共面循環(huán)、垂直共軸循環(huán)、水平共面循環(huán)和線元的相互耦合率，都可以用式（9）中T0、T1、T2表示出來(lái)，計(jì)算T0、T1、T2積分，用自適應(yīng)算法和用濾波系數(shù)直接計(jì)算對(duì)比，可以看出耗時(shí)情況。作者在個(gè)人計(jì)算機(jī)上用Fortran編程，取三層模型參數(shù)，即σ1＝0.01、σ2＝0.001、σ3＝0.01；h1＝100、h2＝100。用自適應(yīng)算法計(jì)算1 000個(gè)頻點(diǎn)所用時(shí)間為0.515 625s，而用801點(diǎn)濾波系數(shù)直接計(jì)算所用時(shí)間為1.343 75s。用其他模型計(jì)算的結(jié)果也基本相符，這可以從算法理論方面解釋?zhuān)撕瘮?shù)T0、T1、T2的關(guān)系是非常簡(jiǎn)單的，兩個(gè)是相同的，另一個(gè)是1／g的關(guān)系，盡管當(dāng)M＞1時(shí)R（D，g）比較復(fù)雜，但是使用了自適應(yīng)算法，在計(jì)算了第一個(gè)核函數(shù)積分之后，后面兩個(gè)可以直接調(diào)用前面的計(jì)算值，稍許修改即可，所以只需比計(jì)算單個(gè)積分多一點(diǎn)時(shí)間。

4 結(jié)論

基于801點(diǎn)快速漢克爾濾波系數(shù)的自適應(yīng)算法，對(duì)于出現(xiàn)在電磁計(jì)算里的貝塞爾函數(shù)積分問(wèn)題尤其是相關(guān)核函數(shù)積分問(wèn)題，能夠節(jié)省計(jì)算時(shí)間。這不僅是在保證精度的前提下進(jìn)行了濾波系數(shù)的截?cái)啵乙脖苊饬撕撕瘮?shù)褶積值的重復(fù)計(jì)算。由于801點(diǎn)系數(shù)對(duì)于零階和一階具有相同的橫坐標(biāo)，所以對(duì)電磁法計(jì)算中常出現(xiàn)零階和一階貝塞爾函數(shù)積分，也可以統(tǒng)一來(lái)計(jì)算。用固定長(zhǎng)度的濾波系數(shù)直接計(jì)算，要么系數(shù)太短達(dá)不到足夠的精度，要么系數(shù)太長(zhǎng)又浪費(fèi)不必要的時(shí)間，利用801點(diǎn)的自適應(yīng)算法很好地解決了這個(gè)問(wèn)題。

［1］GHOSH D P.The application of linear filter theory to the direct interpretation of geoelectrical resistivity sounding measurements［J］.Geophysical Prospecting，1971，19：192－217.

［2］KOEFOED O.A note on the linear filter method of interpreting resistivity sounding data ［J］.Geophysical Prospecting，1972，20：403－405.

［3］KOEFOED O，GHOSH D P，POLMAN G J.Computation of type curves for electromagnetic depth sounding with a horizontal transmitting coil by means of a digital linear filter［J］.Geophysical Prospecting，1972，20：406－420.

［4］DAS U C，GHOSH D P.The determination of filter coefficients for the computation of standard curves for dipole resistivity sounding over layered earth by digital linear filtering［J］.Geophysical Prospecting，1974，22：765－780.

［5］JOHANSEN H K， SORENSEN K. Fast Hankel transforms［J］.Geophysical Prospecting，1979，27：876－901.

［6］ANDERSON W L.Fast Hankel transforms using related and lagged convolutions［J］.ACM Transactions on Mathematical Software，1982，8：344－368.

［7］GUPTASARMA D，SINGH .New digital linear filters for Hankel J0and J1transforms ［J］.Geophysical Prospecting，1997，45：745－765.

［8］樸化榮.電磁測(cè)深法原理［M］.北京：地質(zhì)出版社，1990.

［9］ANDERSON W L.Computer program numerical integration of related Hankel transforms of orders 0and 1by adaptive digital filtering［J］.Geophysics，1979，44：1287－1305.

［10］張偉，王緒本，覃慶炎.漢克爾變換的數(shù)值計(jì)算與精度的對(duì)比［J］.物探與化探，2010，34（6）：753－755.