漢克爾自適應算法在電磁計算中的應用

蔡 盛，柳建新，湯文武，曹創華，鄧小康

（1.中鐵第四勘察設計院集團有限公司，武漢 430063；2.中南大學 地球科學與信息物理學院，長沙 410083）

0 引言

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

現有的漢克爾長濾波系數，在一個適當的范圍之內，精度都是很高的，相對誤差可控制在10－5以上，甚至有時精度達到了10－10以上。在地球物理中，實際應用時需要考慮計算精度和計算需要的時間，而且這兩者同等重要。在一個可以接受的計算精度下，控制住計算時間，極其重要。對于一些精度較高的長濾波器，如果能在滿足精度要求的前提下，應用一種自適應算法，得到需要的精度即結束計算過程，將節省大量的計算時間。在地球物理數值計算的過程中，經常會遇到相關核函數的積分問題，在數值濾波中，積分最終歸結為褶積值的計算，所以相關核函數積分時，褶積值的重復計算量很大。把褶積值存儲起來為之后的計算調用，結合自適應算法在數值積分時，會節省一定的時間。

1 濾波系數的選取

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

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

其中 λi為采樣點的位置；Wi為濾波權重值。積分的精度，決定于權重值Wi。選取五種表現很好的快速漢克爾長濾波系數如表1所示。

表1 選擇的五種高精度長濾波系數Tab.1 The five high－precision long filter coefficients selected

選擇兩個有解析解的貝塞爾函數積分式，即取式（1）中積分核函數f（λ）零階時為λe－λ，用五種系數求取其各自的相對誤差，如圖1所示；一階時f（λ）為λe－2λ，用五種系數求取其各自的相對誤差，如圖2所示。用相對誤差而不用絕對誤差是因為一個小的絕對誤差可能會誤導我們（例如如果理論計算值很小，而實際計算值有很大的誤差，也可能導致其絕對誤差比較?。?，因此相對誤差是更好的評價方法。

圖1、圖2中橫坐標是求值區間，也就是公式（1）和公式（2）中的實參數r的對數，即log（r）；縱坐標是相對誤差。由圖1和圖2可以看出，在零階積分中，801點系數和120點系數表現最好，而在一階積分中，801點系數和140點系數表現最好；801點系數在全求值區間表現比較穩定，而120（J0）和140（J1）點系數受求值區間的影響較大。就計算精度而言，考慮到電磁法中進行數值計算中的一般性，認為801點系數表現最好，而且由于801點系數零階和一階具有相同的系數個數、相同的采樣坐標，有利于統一進行計算。

圖1 零階時漢克爾積分精度對比Fig.1 The relative errors of the five J0filters

2 相關核函數的自適應算法

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

2.1 自適應計算

把濾波系數按橫坐標劃為三個部分，①固定區域；②固定區域右邊部分；③固定區域左邊部分。固定區域需要包括擾動比較劇烈的部分［9］。在801點的濾波系數中，如圖3、圖4所示，規定振蕩劇烈的含有92個濾波權重值范圍為固定范圍，即從256點到第348個點。在這個區域里面，從圖3、圖4中可以看出，0階和1階的濾波器絕對值都快速趨近于“0”。自適應計算過程的第一步，即計算固定范圍內的褶積值；第二步，自適應褶積從固定范圍的右邊，從左至右展開；第三步，計算固定范圍的左邊，從右至左展開。當固定范圍的兩邊褶積的總和的絕對值小于或者等于截斷標準時，自適應過程就停止了。在固定范圍之外的自適應褶積的截斷標準，由下式確定：

復數和的實部分和虛部分，必須獨立地滿足這個標準，才能接受截斷。截斷標準基本近似等于褶積誤差，但是不等于真正的積分誤差。

圖2 一階時漢克爾積分精度對比Fig.2 The relative errors of the five J1filters

圖3 零階801點濾波系數權重值截斷示意圖Fig.3 The truncated diagram of the 801 J0 filter coefficients

在實際中遇到的大多核函數，在固定范圍的兩側只需要少許的褶積值即可達到所需要的計算精度。然而，對于一些緩慢衰減或者振蕩的，核函數可能需要在一個較大范圍的橫坐標之內進行濾波計算，這也是在一個很大范圍提供濾波權重值的原因。

2.2 相關核函數的計算

所謂相關核函數的問題，即幾個核函數能夠用簡單的代數關系聯系起來，那這樣的變換就是相關的。例如，在計算不同的循環結構的電磁耦合率的大地二層模型時，需要計算一個零階和兩個一階漢克爾變換求值問題，而它們的核函數幾乎是相同的。針對復變的相關核函數的0階或者1階快速高精度的漢克爾變換積分問題，為了避免對核函數進行重復計算，在第一次計算核函數時，應先將計算值保存在計算機內存中，在下一次計算相關核函數時，就可以直接調用內存中的值，這將避免重復計算時的時間消耗。

如果一個核函數為f（λ），它的相關核函數可以表示成λi［f（λ）］j。相關核函數褶積運算過程：首先把濾波器的權重值放在一個DATA數組里面，在直接核函數計算時，對DATA數組進行一個初始的調用（NEW＝1），計算的核函數的值和相對應的參數保存在計算機內存的一個公共模塊中。當初始調用結束之后，公共存儲區可以因為相關核函數之間的幾何關系進行修改。然后，任意階的相關漢克爾變換在隨后的調用（NEW＝0）中將會獲得存儲區的值，不需要重復計算相關核函數的褶積值，除非當前核函數的特性需要一些另外的計算。而這些另外的核函數褶積值，如果需要，也可以計算并且自動添加到公共模塊區以備進一步使用。程序流程簡圖如圖5所示。

圖4 一階801點濾波系數權重值截斷示意圖Fig.4 The truncated diagram of the 801 J1 filter coefficients

圖5 相關核函數積分計算流程簡圖Fig.5 The calculation diagram of the related kernels

3 實例計算與結果分析

使用一組相關核函數的漢克爾解析式，來測試相關核函數自適應算法的精度情況。這組貝塞爾函數積分式子，都是有解析解的，可以得到算法濾波值對于精確解的相對誤差［10］。

五個貝塞爾函數積分式如下：

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

表2 自適應算法計算結果表Tab.2 The calculation results of the adaptive algorithm

上面五個貝塞爾函數積分式的核函數都比較簡單。當要計算的相關核函數不是簡單基本函數時，使用自適應算法對于時間的節省將非常明顯。對于上面提到的不同循環結構的相互耦合比的計算，對于水平M層大地模型，忽略位移電流的影響，核心計算就是下面的三個漢克爾變換積分式：

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

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

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

4 結論

基于801點快速漢克爾濾波系數的自適應算法，對于出現在電磁計算里的貝塞爾函數積分問題尤其是相關核函數積分問題，能夠節省計算時間。這不僅是在保證精度的前提下進行了濾波系數的截斷，而且也避免了核函數褶積值的重復計算。由于801點系數對于零階和一階具有相同的橫坐標，所以對電磁法計算中常出現零階和一階貝塞爾函數積分，也可以統一來計算。用固定長度的濾波系數直接計算，要么系數太短達不到足夠的精度，要么系數太長又浪費不必要的時間，利用801點的自適應算法很好地解決了這個問題。

［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］樸化榮.電磁測深法原理［M］.北京：地質出版社，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］張偉，王緒本，覃慶炎.漢克爾變換的數值計算與精度的對比［J］.物探與化探，2010，34（6）：753－755.