電偶極子切分算法研究

胡瑞華，林 君，孫彩堂, 劉長(zhǎng)勝, 周逢道

(吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院, 長(zhǎng)春 130061)

0 引言

當(dāng)前電磁法勘探已成為解決地質(zhì)問(wèn)題的重要手段[1-2]，學(xué)者們通過(guò)正反演理論和技術(shù)來(lái)研究地質(zhì)模型的電磁響應(yīng)和地下結(jié)構(gòu)的復(fù)現(xiàn)，考慮到探測(cè)系統(tǒng)的發(fā)射功率和分辯能力，野外作業(yè)時(shí)常采用人工源激發(fā)達(dá)到獲取高信噪比的電磁數(shù)據(jù)的要求。磁性源和電性源是常被采用的兩種人工源，其中電性源使用一定長(zhǎng)度的導(dǎo)線向地下供一定強(qiáng)度的交變電流來(lái)實(shí)現(xiàn)。在以電性源為場(chǎng)源的電磁探測(cè)方法中，其正演常需要研究基于電偶極子激發(fā)的場(chǎng)響應(yīng)，并根據(jù)響應(yīng)規(guī)律來(lái)解釋野外實(shí)測(cè)數(shù)據(jù)(即反演)[3～12]。這種正反演通常都是假定發(fā)射端與接收點(diǎn)距離遠(yuǎn)大于供電兩端的長(zhǎng)度，即要求導(dǎo)線長(zhǎng)度滿足偶極子的條件。文獻(xiàn)[4]指出，當(dāng)收發(fā)距大于等于10倍的導(dǎo)線長(zhǎng)度時(shí)，導(dǎo)線可被看作是電偶極子。然而在實(shí)際應(yīng)用中由于受地形地貌、發(fā)射功率等多種因素的制約，導(dǎo)致長(zhǎng)導(dǎo)線并不滿足電偶極子的條件，而使基于電偶極子的場(chǎng)計(jì)算公式不再適用，在這種情況下若仍然按照電偶極子激勵(lì)理論開展數(shù)據(jù)解釋，將造成較大偏差甚至是錯(cuò)誤的地質(zhì)結(jié)構(gòu)的復(fù)現(xiàn)。為了解決這個(gè)問(wèn)題，在正反演中可以先將長(zhǎng)導(dǎo)線源切分成若干個(gè)電偶極子，然后對(duì)切分的電偶極子的響應(yīng)進(jìn)行求和，以達(dá)到計(jì)算長(zhǎng)導(dǎo)線源場(chǎng)響應(yīng)的目的。

針對(duì)如何將長(zhǎng)導(dǎo)線源切分成多個(gè)電偶極子的方法，本次研究并實(shí)現(xiàn)了兩種切分算法，①借鑒二分法思想設(shè)計(jì)的遞歸切分算法；②是借鑒窮舉法思想設(shè)計(jì)的窮舉切分算法。通過(guò)對(duì)兩種切分算法的研究表明：遞歸切分算法切分出的電偶極子數(shù)目相對(duì)較多；窮舉切分算法在窮舉步長(zhǎng)較小(如0.1 m)時(shí)能較準(zhǔn)確地切分出真實(shí)電偶極子的位置。將兩種算法應(yīng)用有于CSAMT電場(chǎng)的計(jì)算中，計(jì)算結(jié)果表明：兩種切分算法都可應(yīng)用于實(shí)際場(chǎng)的計(jì)算。

1 切分算法原理

1.1 觀測(cè)點(diǎn)與長(zhǎng)導(dǎo)線源的拓?fù)浣Y(jié)構(gòu)

人工源野外測(cè)量時(shí)，觀測(cè)點(diǎn)離長(zhǎng)導(dǎo)線源有一定的距離，接收和發(fā)射之間的是三維拓?fù)浣Y(jié)構(gòu)[9,11]，常見(jiàn)的拓?fù)浣Y(jié)構(gòu)有平地形的地面發(fā)射地面接收、起伏地形的地面發(fā)射地面接收、地面發(fā)射井中接收等(圖1)。

圖1 觀測(cè)點(diǎn)與長(zhǎng)導(dǎo)線源的拓?fù)浣Y(jié)構(gòu)Fig．1 Topology structure of survey point and long lead source

圖1中長(zhǎng)方體代表大地，AB表示長(zhǎng)度為AB的長(zhǎng)導(dǎo)線發(fā)射源，A和B分別代表發(fā)射源的兩極，P1、P2、P3分別表示與AB共面、共線和地下的觀測(cè)點(diǎn)。P1、P2、P3與AB的三種位置關(guān)系與野外測(cè)量時(shí)的施工布局一致。偶極子的位置常用發(fā)射源兩極(導(dǎo)線兩端)的中點(diǎn)代表。若空間坐標(biāo)系中A極、B極、測(cè)點(diǎn)P三點(diǎn)坐標(biāo)分別取為A(xA，yA，zA)、B(xB，yB，zB)、P(xP，yP，zP)，其中x、y、z表示點(diǎn)的坐標(biāo)分量，若用C表示AB的中點(diǎn)，則接收點(diǎn)與發(fā)射源的距離(發(fā)收距)表示為P點(diǎn)到C點(diǎn)的距離PC(常用r字符表示)。根據(jù)A、B坐標(biāo)及代數(shù)關(guān)系很容易計(jì)算出C點(diǎn)的坐標(biāo)、發(fā)收距PC和發(fā)射極距AB。另外，偶極子的條件是收發(fā)距大于等于N(通常最小值可取為10)倍的導(dǎo)線長(zhǎng)度，即PC≥N*AB。實(shí)際應(yīng)用中A、B、P三點(diǎn)坐標(biāo)已知，而N參數(shù)的確定可在綜合考慮工區(qū)、收發(fā)系統(tǒng)參數(shù)以及電磁數(shù)據(jù)質(zhì)量的情況下選定。

1.2 二分法遞歸切分算法

二分法算法是一種快速搜尋算法，算法的核心是從有序數(shù)據(jù)列中按照折半查找的方式尋找滿足條件的結(jié)果[12-14]。借鑒二分法思想研究切分長(zhǎng)導(dǎo)線源為多個(gè)電偶極子的算法如下：

1) 根據(jù)長(zhǎng)導(dǎo)線兩端A點(diǎn)、B點(diǎn)的坐標(biāo)計(jì)算其中點(diǎn)C的坐標(biāo)和A點(diǎn)到B點(diǎn)的距離AB；根據(jù)P點(diǎn)、C點(diǎn)的坐標(biāo)計(jì)算P點(diǎn)到C點(diǎn)的距離PC，轉(zhuǎn)到步驟2)。

2) 如果PC≥N*AB則記錄下C點(diǎn)的位置，算法結(jié)束；否則轉(zhuǎn)到步驟3)。

3)用C點(diǎn)坐標(biāo)替換B點(diǎn)的坐標(biāo)，轉(zhuǎn)到步驟1)。

4)用C點(diǎn)坐標(biāo)替換A點(diǎn)的坐標(biāo)，轉(zhuǎn)到步驟1)。

此算法是遞歸算法，共分四步。第一步計(jì)算發(fā)射距AB和收發(fā)距PC，為算法的下一步作準(zhǔn)備；第二步判斷AB是否為電偶極子，若是電偶極子，則AB的中點(diǎn)就是電偶極子的位置，算法結(jié)束，它是遞歸算法的出口；第三步當(dāng)AB不是電偶極子時(shí)，搜尋CA段導(dǎo)線電偶極子的位置；第四步當(dāng)搜尋完CA段導(dǎo)線的電偶極子時(shí)，再搜尋CB段導(dǎo)線的電偶極子。每遞歸一次就能找到一個(gè)電偶極子，當(dāng)整個(gè)遞歸算法結(jié)束時(shí)，便找到了所有電偶極子的位置點(diǎn)(圖2)。

圖2 二分法遞歸切分算法流程圖Fig．2 Flow of binary segmentation method

1.3 窮舉切分算法

窮舉算法是一種試算遍歷算法，算法的核心是從取值區(qū)間對(duì)每個(gè)取到的值(或組合值)進(jìn)行試算，通過(guò)判斷試算結(jié)果是否滿足要求，從而找到問(wèn)題的解[12,15]，這種算法常被用在方程少而未知數(shù)多的問(wèn)題中。將長(zhǎng)導(dǎo)線源切分為多個(gè)電偶極子的問(wèn)題，本質(zhì)上是查找AB上哪些點(diǎn)是偶極子的位置，但我們僅知道偶極子的條件PCnew≥N*AnewBnew(式中AnewBnew表示AB上的一段導(dǎo)線，Cnew為其中點(diǎn))，這樣必須通過(guò)1個(gè)方程(PCnew=N*AnewBnew)求解AnewBnew的中點(diǎn)Cnew的位置(即要確定Cnew的3個(gè)坐標(biāo)分量)，顯然方程少而未知數(shù)多，無(wú)法直接求解出方程的解，即無(wú)法定位偶極子的位置，因此這個(gè)問(wèn)題屬于窮舉算法能解決的范疇。借鑒窮舉算法的思想，設(shè)計(jì)將長(zhǎng)導(dǎo)線源切分為多個(gè)電偶極子的算法如下：

1)根據(jù)切分步長(zhǎng)以A或B為參考點(diǎn)(算法描述中以A為參考點(diǎn))對(duì)AB進(jìn)行切分，得到切分點(diǎn)列表C=(C1,C2,…，Cn);設(shè)置兩個(gè)計(jì)數(shù)器變量i、j，并使i= 0，j= 0，轉(zhuǎn)到步驟2)。

2)改變i使i=i+1；計(jì)算點(diǎn)P到點(diǎn)Ci的距離PCi，點(diǎn)Ci到點(diǎn)B的距離CiB，點(diǎn)Ci到點(diǎn)A的距離CiA，轉(zhuǎn)到步驟3)。

3)如果PCi

4)改變i使i=i-1，如果CiB≥CiA，則記錄下Ci；用Ci+j的坐標(biāo)替換A點(diǎn)的坐標(biāo)；i=i+j，j=0，轉(zhuǎn)到步驟6)；否則轉(zhuǎn)到步驟5)。

5)記錄下AB的中點(diǎn)，結(jié)束算法。

6)改變j使j=j+1，轉(zhuǎn)到步驟2)。

算法共由六步組成。第一步是參考窮舉算法的原理根據(jù)切分步長(zhǎng)對(duì)AB長(zhǎng)導(dǎo)線源進(jìn)行切分，切分的小段中除最后一小段外，其他的切分小段長(zhǎng)度是等長(zhǎng)的(若AB能被步長(zhǎng)整除，則所有切分出的小段長(zhǎng)度都是等長(zhǎng)的)，并記下所有的切分點(diǎn)Ci，偶極子的位置就存在于這些切分點(diǎn)中和最后一個(gè)偶極子的中點(diǎn)，為了能找到偶極子的位置，用i變量跟蹤切分點(diǎn)的位置，用j變量跟蹤上一個(gè)偶極子B端到Ci的切分點(diǎn)個(gè)數(shù)；第二步計(jì)算PCi、CiB、CiA；第三步和第四步判斷PCi是否同時(shí)滿足偶極子的條件和CiB是否為原AB的最后一段，在兩個(gè)條件都滿足情況下，則記下Ci，然后計(jì)算出下一個(gè)待求偶極的A端位置，然后執(zhí)行第六步進(jìn)入查找下一次偶極子的循環(huán)中；另外，在第四步中若CiB是原AB切分段的最后一段，則偶極子的位置直接取CiB的中點(diǎn)，算法結(jié)束(圖3)。

圖3 窮舉切分算法流程圖Fig．3 Flow of exhaustion segmentation method

2 算法應(yīng)用

使用兩種算法對(duì)地表長(zhǎng)導(dǎo)線源AB，其中A極坐標(biāo)為(0，0，0)，B極坐標(biāo)為(500，0，0)，觀測(cè)點(diǎn)P坐標(biāo)為(50，500，0)進(jìn)行偶極子的切分，偶極子的條件取為收發(fā)距是10倍的發(fā)射距。切分的偶極子位置對(duì)比如圖4所示。

圖4 兩種算法切分結(jié)果對(duì)比Fig．4 Segmentation results comparison to two algorithms

窮舉切分算法切分步長(zhǎng)分別選擇0.1 m、1 m和5 m，用這三種步長(zhǎng)切分的偶極子數(shù)目都是9個(gè)。從圖4可以看出，用步長(zhǎng)0.1 m和步長(zhǎng)1 m切分出的偶極子位置差別很小，與理論上的準(zhǔn)確偶極子位置基本一致；而用步長(zhǎng)5m切分出的偶極子位置與用步長(zhǎng)0.1 m切分出的偶極子位置差別明顯。如果將0.1 m步長(zhǎng)切分出的偶極子看作是真實(shí)(理論上)的偶極子，那么切分步長(zhǎng)越大切分出的偶極子的位置與真實(shí)偶極子的位置誤差就越大。采用遞歸算法切分出的偶極子數(shù)目為15個(gè)，多于0.1 m步長(zhǎng)的窮舉法切分的偶極子數(shù)目，表明遞歸算法切分的偶極子長(zhǎng)度要小于用窮舉法切分的偶極子的長(zhǎng)度。

為了進(jìn)一步了解場(chǎng)計(jì)算中需要將長(zhǎng)導(dǎo)線源切分為多個(gè)電偶源的必要性及算法的有效性，本次研究以可控源音頻大在電磁法(CSAMT)一維電場(chǎng)計(jì)算為例對(duì)其說(shuō)明。與切分算法相關(guān)的參數(shù)：發(fā)射源極距AB長(zhǎng) 2 000 m，A點(diǎn)坐標(biāo)為(0，0，0)，B點(diǎn)坐標(biāo)為(2 000，0，0)，觀測(cè)點(diǎn)P位于離AB垂直距離 7 000 m 處，坐標(biāo)為(0，7 000，0)，偶極子的條件取收發(fā)距是10 倍的發(fā)射距。使用文獻(xiàn)[1]中場(chǎng)的計(jì)算公式求解一個(gè)三層H型地電模型(電阻率：300 Ω·m，50 Ω·m，1 000 Ω·m；厚度：300 m，200 m)x方向的電場(chǎng)，并分別使用文中的兩種切分算法和不切分三種方式計(jì)算出的電場(chǎng)結(jié)果如圖5所示。

圖5是頻率-電場(chǎng)曲線圖，從圖5中可以看出，用兩種算法切分的偶極子數(shù)目都為 4 個(gè)，且用遞歸算法切分的偶極子計(jì)算出的電場(chǎng)曲線與用窮舉算法在步長(zhǎng)分別取 0.1 m、1 m 切分的偶極子計(jì)算出的電場(chǎng)曲線基本重合，這表明了遞歸切分算法可用于實(shí)際應(yīng)用中；而步長(zhǎng)取 50 m 切分的偶極子計(jì)算出的電場(chǎng)曲線與步長(zhǎng)取 0.1 m 切分的偶極子計(jì)算出的電場(chǎng)曲線明顯分離，表明用較大步長(zhǎng)切分偶極子計(jì)算出的電場(chǎng)會(huì)帶來(lái)較大的誤差；不做任何切分計(jì)算得到的場(chǎng)曲線與用0.1 m步長(zhǎng)切分后計(jì)算得到的場(chǎng)曲線不重合，這表明在不滿足偶子條件的情況下，若仍然按照偶極子激發(fā)公式計(jì)算場(chǎng)將會(huì)得到不正確的結(jié)果。

3 結(jié)論

采用長(zhǎng)導(dǎo)線源激勵(lì)時(shí)，為了避免收發(fā)距較小時(shí)源長(zhǎng)度對(duì)場(chǎng)計(jì)算的影響，需要考慮到將長(zhǎng)導(dǎo)線源切分為多個(gè)偶極子源，這樣長(zhǎng)導(dǎo)線源場(chǎng)的計(jì)算就轉(zhuǎn)化為多個(gè)偶極子場(chǎng)的疊加。遞歸切分算法和窮舉切分算法都可以實(shí)現(xiàn)對(duì)長(zhǎng)導(dǎo)線的切分。在實(shí)際物探應(yīng)用中，遞歸切分算法可直接使用，窮舉算法的切分步可取 0.1 m或 1 m。

圖5 兩種切分算法和不切分電場(chǎng)計(jì)算結(jié)果對(duì)比Fig．5 Electric field calculating comparison to two segmentation algorithms and no segmentation

參考文獻(xiàn)：

[1] 湯井田, 何繼善.可控源音頻大地電磁法及其應(yīng)用[M]．長(zhǎng)沙: 中南大學(xué)出版社,2005．

[2] 何繼善．可控源音頻大地電磁法[M]．長(zhǎng)沙:中南大學(xué)出版社,1990.

[3] 王剛,王書民, 雷達(dá)，等.CSAMT場(chǎng)源效應(yīng)試驗(yàn)研究[J].物探化探計(jì)算技術(shù),2011,33(5):527-530.

[4] 趙廣茂,李志華,朱旭東，等.長(zhǎng)導(dǎo)線源CSAMT一維正演研究[J].鐵道工程學(xué)報(bào),2010,143(8):21-24.

[5] 殷長(zhǎng)春. 可控源音頻磁大地電流法一維正演及精度評(píng)價(jià)[J].長(zhǎng)春地質(zhì)學(xué)院學(xué)報(bào),1994( 4):438- 453.

[6] 孫婭,何展翔,柳建新，等. 長(zhǎng)導(dǎo)線源頻率域電磁測(cè)深場(chǎng)源靜態(tài)位移的模擬研究[J].石油地球物理勘探,2011,46( 1):149-154.

[7] 李勇,林品榮,徐寶利.電離層影響層狀介質(zhì)長(zhǎng)導(dǎo)線源的電磁場(chǎng)[J].物探化探計(jì)算技術(shù),2009,31( 3):183-188.

[8] 劉昱,李云哲,黃子瑩.水平長(zhǎng)導(dǎo)線源瞬變電磁三維積分方程模擬[J].現(xiàn)代地質(zhì),2010,24(2):397-402.

[9] 屈有恒,張貴賓,晉風(fēng)明.傾斜線源的三維電場(chǎng)數(shù)值模擬研究[J].物探化探計(jì)算技術(shù),2007,29(5):431-435.

[10] 許振奎,張勝業(yè).層狀介質(zhì)中磁性源頻率測(cè)深的一種視電阻率算法[J].勘探地球物理進(jìn)展,2004,27(5):333-336.

[11] 徐建華,朱德懷,胡文寶，等.多層介質(zhì)中偶極子場(chǎng)的系數(shù)遞推關(guān)系[J].石油地球物理勘探,1994,29(1):69-74.

[12] ALSUWAIYEL M H. Algorithms Design Techniques and Analysis[M].Beijing: Publishing House of Electronics Industry,2003.

[13] 王海濤,朱洪.改進(jìn)的二分法查找[J].計(jì)算機(jī)工程,2006,32(10):60-62.

[14] 陳寶平.遞歸算法的設(shè)計(jì)模式與調(diào)試[J].電子科技,2011,24(9):28-33.

[15] 孫義欣,馮娜.窮舉法在程序設(shè)計(jì)中的應(yīng)用[J].計(jì)算機(jī)時(shí)代,2012,8:50-52.