多孔對(duì)電磁波CT剖面拼接的影響因素及拼接方法

江振寅，姜 杰，鄭 軍，田必林

（1.中國(guó)煤炭地質(zhì)總局湖北煤炭地質(zhì)局，湖北 武漢 430070；2.湖北煤炭地質(zhì)物探測(cè)量隊(duì)，湖北 武漢 430200）

1 引言

井間電磁波CT成像是高效探測(cè)地下介質(zhì)電磁特性的地球物理方法［1］，其主要原理就是利用一定頻率的電磁波分別在兩個(gè)鉆孔發(fā)射和接收，根據(jù)不同位置接收到場(chǎng)強(qiáng)的大小，對(duì)兩孔之間電磁波吸收系數(shù)進(jìn)行處理來(lái)確定地下不同介質(zhì)分布，以此來(lái)精確地界定地下目標(biāo)體的屬性［2］。井間電磁波CT法以其抗干擾能力強(qiáng)，不受地形限制等優(yōu)勢(shì)［3］，在工程、煤炭、油氣和金屬礦產(chǎn)等資源勘探開(kāi)發(fā)領(lǐng)域以及水文工程方面都得到了廣泛的應(yīng)用，并取得了較好的效果［4］。

通過(guò)收集的資料來(lái)看，近年來(lái)井間電磁波CT研究工作多集中于單對(duì)孔反演理論研究，在實(shí)際工程勘察工作中，常常需要在同一個(gè)工作區(qū)完成幾對(duì)鉆孔乃至幾十對(duì)鉆孔的電磁波CT成像測(cè)量，對(duì)巖溶、破碎、構(gòu)造等地質(zhì)結(jié)構(gòu)進(jìn)行勘查，目前最常用的做法就是對(duì)每一個(gè)孔對(duì)的剖面進(jìn)行單獨(dú)成像單獨(dú)解釋［5］。但單孔對(duì)成像往往不能滿足綜合解釋的需要，有時(shí)需要將相鄰的孔對(duì)相互連接成為一個(gè)完整剖面，而在連接孔處的成像結(jié)果卻常常不吻合，相鄰孔對(duì)之間的背景值差異過(guò)大，無(wú)法連接，以至于連接剖面的成像結(jié)果遭到質(zhì)疑［5］。針對(duì)這類(lèi)問(wèn)題，相關(guān)資料顯示，2010年雷旭友［6］提出多源多孔對(duì)電磁波CT聯(lián)合反演成像的方法，將采集的原始數(shù)據(jù)進(jìn)行簡(jiǎn)單的組合，對(duì)因迭代次數(shù)、隨機(jī)誤差等方面因素所引起的相鄰孔對(duì)成像結(jié)果不吻合現(xiàn)象具有較好的消除作用，反演的成果可在數(shù)值上綜合反映地下地質(zhì)體共同的特性，對(duì)研究地質(zhì)體延伸、橫向變化具有實(shí)用價(jià)值，但該方法最后并未得到推廣應(yīng)用。

本文從單孔對(duì)入手，從反演方法和初始場(chǎng)強(qiáng)E0兩個(gè)方面對(duì)可能導(dǎo)致相鄰鉆孔拼接不吻合的因素進(jìn)行分析，探索尋求實(shí)現(xiàn)多孔對(duì)電磁波CT剖面拼接的方法，并結(jié)合工程勘察中既有數(shù)據(jù)，分析拼接效果，為今后多孔對(duì)電磁波CT成果解釋和判斷提供參考。

2 電磁波CT成像基本原理

根據(jù)電磁波理論，電磁波傳輸過(guò)程中，耗散介質(zhì)的初始場(chǎng)強(qiáng)E0與接收?qǐng)鰪?qiáng)E關(guān)系為［7］：

其中，E是觀測(cè)點(diǎn)接收到的場(chǎng)強(qiáng)（dB）；E0為初始場(chǎng)強(qiáng)（dB）；a為介質(zhì)對(duì)電磁波吸收系數(shù)；r為收發(fā)點(diǎn)間距（m）；f（θ）為收發(fā)天線方向因子；θ為接收點(diǎn)處天線與電場(chǎng)方向夾角［8］。一般情況下E0值是未知的，且介質(zhì)電磁參數(shù)變化較大時(shí)，E0值隨位置而變化［9］。

將式（1）兩邊經(jīng)過(guò)變換后取對(duì)數(shù)可以得到ar＝，寫(xiě)成離散化形式，得到第i個(gè)觀測(cè)點(diǎn)的方程＝di，其中di＝，aj是第j（j＝1，2，3，……）個(gè)網(wǎng)格電磁波的吸收系數(shù)，rij是第i條射線在第j個(gè)網(wǎng)格中的射線長(zhǎng)度，di是第i（i＝1，2，3，……）個(gè)觀測(cè)點(diǎn)的觀測(cè)值［10］。

假設(shè)R是n條電磁波射線在m個(gè)網(wǎng)格單元中的射線元組成的m×n維矩陣；X是電磁波吸收系數(shù)組成的矩陣，X＝（α1，α2，α3，α4，…，αn）T；D為接收?qǐng)鰪?qiáng)的觀測(cè)值，D＝（d1，d2，d3，d4，…，dn）T，則所有觀測(cè)點(diǎn)的發(fā)射和接收之間得到的線性方程組可以寫(xiě)成：RX＝D，通過(guò)合適的迭代算法求解上述線性方程組，就能夠得到兩個(gè)鉆孔間介質(zhì)的吸收系數(shù)，進(jìn)一步得到吸收系數(shù)的等值線圖［11］。

解上述線性方程組的方法有很多［12］，比如共軛梯度法CG（Conjugate Gradient）、代數(shù)重建法ART（Algebraic Reconstruction Technique）、奇異值分解法SVD（Singular Value Decomposition）、聯(lián)合迭代法SIRT（Simultaneous Iterative Reconstruction Technique）、最小二乘QR分解法LSQR（Least Squares QR－decomposition）等，不同方法的成像結(jié)果各有區(qū)別［13］。

3 影響相鄰孔對(duì)數(shù)據(jù)不吻合的因素

從上述電磁波基本原理可以看出，影響電磁波CT成像誤差的主要因素包括反演方法的選擇和初始場(chǎng)強(qiáng)E0的處理。

3.1 反演方法的影響

反演方法包括算法本身和反演參數(shù)的設(shè)置。2019年趙威對(duì)電磁波CT幾種常用的成像方法進(jìn)行了對(duì)比分析［14］，確定SIRT是最合適的方法。實(shí)際工作中，同一個(gè)工作區(qū)相鄰孔對(duì)的反演成像通常采用相同的反演算法，但反演參數(shù)的設(shè)置可能略有差別，因此下文中僅對(duì)反演參數(shù)的設(shè)置對(duì)相鄰孔對(duì)成像誤差的影響進(jìn)行說(shuō)明。

電磁波CT反演參數(shù)主要包括迭代次數(shù)、射線角度、初始模型等。參數(shù)較多，以迭代次數(shù)為例，建立一個(gè)簡(jiǎn)單吸收系數(shù)模型，如圖1所示，模型中ZK1－ZK2、ZK2－ZK3為兩對(duì)相鄰孔對(duì)，ZK2為共用鉆孔，模型背景值均設(shè)置為0.20Np／m，設(shè)置高吸收異常體和低吸收異常體各一處，大小為6m×4m，吸收系數(shù)分別為0.3Np／m、0.1Np／m，異常體中心深度分別為5m、15m。

圖1 吸收系數(shù)模型Fig.1 Model of the absorption coefficient

對(duì)上述模型中兩對(duì)孔ZK1－ZK2、ZK2－ZK3分別采用SIRT法成像，除迭代次數(shù)外，其他反演參數(shù)均保持一致，提取成像結(jié)果中共用鉆孔ZK2處反演數(shù)據(jù)，進(jìn)行對(duì)比，計(jì)算相對(duì)誤差。表1為該模型不同迭代次數(shù)成像結(jié)果在共用孔ZK2處的數(shù)據(jù)誤差。

表1 共用孔處不同迭代次數(shù)相對(duì)誤差（單位：%）Table 1 Relative error of different stacking times at shared holes（unit：%）

從表1中可以看出，理想情況下，反演參數(shù)完全相同時(shí)共用孔處數(shù)據(jù)的誤差最小達(dá)到了0；隨著迭代次數(shù)的增大，相對(duì)誤差逐漸增大。實(shí)際工作中存在一些干擾，可能達(dá)不到理想的最小誤差0，但總的來(lái)說(shuō)，反演過(guò)程中統(tǒng)一反演方法和反演參數(shù)可以使相鄰孔對(duì)在共用鉆孔位置的誤差最小。

3.2 初始場(chǎng)強(qiáng)E0的影響

在射線條件下，影響井間電磁波CT成像的因素中，E0值的確定是比較關(guān)鍵的問(wèn)題，不僅因?yàn)镋0是未知的，而且當(dāng)介質(zhì)電磁參數(shù)變化較大時(shí)，E0的值是隨位置而不斷變化的。張輝等［16］詳細(xì)研究了初始場(chǎng)強(qiáng)對(duì)成像結(jié)果的重大影響，并進(jìn)一步探究了該影響特征所呈現(xiàn)的規(guī)律性，得出不同場(chǎng)強(qiáng)會(huì)使成像結(jié)果形態(tài)不同。實(shí)際資料處理時(shí)，可以通過(guò)觀察初始場(chǎng)強(qiáng)與反演結(jié)果之間的關(guān)系，從理論和經(jīng)驗(yàn)兩方面共同確定合理的初始場(chǎng)強(qiáng)值來(lái)部分消除對(duì)反演結(jié)果的影響，但如何準(zhǔn)確獲取各點(diǎn)的初始場(chǎng)強(qiáng)以得到更準(zhǔn)確的反演結(jié)果，目前還沒(méi)有一個(gè)比較成熟的方法。

4 相鄰孔對(duì)的拼接

4.1 拼接思路

根據(jù)以上分析，相鄰孔對(duì)共用鉆孔位置的誤差并不能從根本上完全消除，但是可以通過(guò)統(tǒng)一反演算法和反演參數(shù)、確定合理的初始場(chǎng)強(qiáng)值來(lái)最大限度地降低誤差。

以貴州某巖溶勘查項(xiàng)目中電磁波CT既有數(shù)據(jù)為例，圖2為兩個(gè)相鄰孔對(duì)各項(xiàng)反演參數(shù)進(jìn)行統(tǒng)一、多次嘗試確定初始場(chǎng)強(qiáng)的基礎(chǔ)之上反演得到的成果圖。從圖2中可以看出，對(duì)于單對(duì)孔剖面而言，兩個(gè)孔對(duì)的反演成果與實(shí)際鉆孔均能吻合，但當(dāng)兩個(gè)孔對(duì)的反演數(shù)據(jù)直接進(jìn)行拼接后，成圖如圖3（a）所示，在共用鉆孔位置并不能較好地吻合，出現(xiàn)明顯的等值線不連續(xù)現(xiàn)象。

圖2 相鄰孔對(duì)單孔等值線Fig.2 Contour map of adjacent holes versus single holes

圖3 相鄰孔對(duì)拼接效果Fig.3 Picture of adjacent hole pairs splicing effect

經(jīng)過(guò)統(tǒng)計(jì)，勘查區(qū)內(nèi)92對(duì)孔均與圖3（a）中情況相似，對(duì)于單孔對(duì)而言，反演成果均能夠跟實(shí)際鉆孔吻合，但相互拼接時(shí)，卻在共用孔位置出現(xiàn)明顯的等值線不連續(xù)現(xiàn)象。按照上述情況，如果將兩個(gè)相鄰孔對(duì)在共用孔處的數(shù)據(jù)進(jìn)行圓滑，而剖面其他位置的數(shù)據(jù)盡量保持不變，則可以實(shí)現(xiàn)相鄰孔對(duì)的拼接，亦能夠保證各剖面中異常不變。通過(guò)多次試驗(yàn)發(fā)現(xiàn)，采用反距離加權(quán)的局部多項(xiàng)式插值方式能夠較好地實(shí)現(xiàn)相鄰孔對(duì)的拼接，如圖3（b）所示。

4.2 反距離加權(quán)的局部多項(xiàng)式插值法

局部多項(xiàng)式插值法是一種局部加權(quán)最小二乘擬合，根據(jù)有限個(gè)已知采樣點(diǎn)數(shù)據(jù)，采用多個(gè)多項(xiàng)式來(lái)進(jìn)行曲面的擬合，每一個(gè)插值點(diǎn)的預(yù)測(cè)值都對(duì)應(yīng)一個(gè)多項(xiàng)式，利用局部多項(xiàng)式插值法得到的曲面，能夠較好地描述短程變異；在一定區(qū)域內(nèi)，不同采樣點(diǎn)對(duì)內(nèi)插點(diǎn)的影響不同，因此在計(jì)算多項(xiàng)式時(shí)，各點(diǎn)的權(quán)重也存在差異。通常情況下，臨近區(qū)域內(nèi)采樣點(diǎn)的權(quán)重會(huì)隨內(nèi)差點(diǎn)與采樣點(diǎn)之間的距離的增加而減小［17］，反距離加權(quán)法是目前最廣泛采用的方法之一。

4.2.1 局部多項(xiàng)式插值方程

假設(shè)采樣點(diǎn)的坐標(biāo)為（xi，yi），其中（i＝1，2，…，n），內(nèi)插點(diǎn)為N，坐標(biāo)（XN，YN）為以N為原點(diǎn)建立局部坐標(biāo)系，則（xi，yi）表示為［18］：

文中采用二次插值，由n個(gè)采樣點(diǎn)的值可以得到方程組：

其中i＝1，2，…，n。寫(xiě)成矩陣形式為：

求解多項(xiàng)式系數(shù)矩陣，表達(dá)式如下：

通過(guò)最小二乘求解多項(xiàng)式系數(shù)：

式（9）中，P為每個(gè)采樣點(diǎn)通過(guò)反距離加權(quán)法確定的權(quán)重，由此求得局部多項(xiàng)式插值方程，進(jìn)而得到內(nèi)插點(diǎn)的數(shù)值。

4.2.1 反距離加權(quán)法

反距離加權(quán)法用來(lái)反映不同采樣點(diǎn)與內(nèi)插點(diǎn)的相關(guān)程度，確定的原則與該數(shù)據(jù)點(diǎn)到內(nèi)插點(diǎn)的距離di（i＝1，2，…，n）有關(guān)，di越大，數(shù)據(jù)點(diǎn)對(duì)內(nèi)插點(diǎn)的影響越小，權(quán)重也越小，反之權(quán)重越大［19，20］。本文以距離平方的反比為權(quán)，權(quán)函數(shù)的表達(dá)式為：

其中，di為數(shù)據(jù)點(diǎn)到插值點(diǎn)的距離。

4.3 相鄰孔對(duì)的拼接

以上述貴州巖溶勘查項(xiàng)目為例，電磁波CT鉆孔孔對(duì)如圖4所示，根據(jù)鉆孔在平面上的位置關(guān)系，將92對(duì)鉆孔劃分為縱橫交錯(cuò)的23條剖面，以下選擇兩種不同異常下的剖面拼接效果進(jìn)行分析，一是團(tuán)塊狀及不規(guī)格形狀異常剖面的拼接（圖5），二是近直立異常的剖面拼接（圖6）。

圖4 井間電磁波CT拼接剖面示意Fig.4 Schematic diagram of cross well electromagnetic wave CT splicing profile

圖5 P2剖面處理前后拼接等值線Fig.5 Splicing contour map before and after P2profile processing

圖6 P22剖面處理前后拼接等值線Fig.6 Splicing contour map before and after P22profile processing

圖5（a）為P2剖面5對(duì)鉆孔原始剖面數(shù)據(jù)拼接成圖。由圖5可以看出，整個(gè)剖面數(shù)據(jù)的整體水平基本一致，但是在k7、k8、k9、k10鉆孔位置存在比較明顯的吸收系數(shù)突變現(xiàn)象，圖上等值線圖連續(xù)，這種現(xiàn)象無(wú)法從地質(zhì)上進(jìn)行解釋。圖5（b）為經(jīng)過(guò)反距離加權(quán)的局部多項(xiàng)式插值法處理后成圖效果，通過(guò)對(duì)比可以看出，經(jīng)處理后拼接的等值線變得圓滑，電磁波吸收系數(shù)異常的位置及大小均未發(fā)生明顯變化，但各孔對(duì)拼接位置等值線不連續(xù)的現(xiàn)象消失，各孔對(duì)之間等值線過(guò)渡自然，且與實(shí)際鉆孔情況吻合。

圖6 為P 22剖面處理前后的等值線圖。其中，圖6（a）為P22剖面5對(duì)鉆孔原始剖面數(shù)據(jù)拼接成圖，圖中K37－K35、K32－K28、K28－K24、K24－K22四對(duì)孔剖面數(shù)據(jù)的整體水平基本一致，但K35－K32孔對(duì)的數(shù)據(jù)較兩側(cè)都有很大的差異，該孔對(duì)在整個(gè)剖面中呈直立的異常存在；經(jīng)過(guò)曲線擬合后數(shù)據(jù)成圖如圖6（b）所示，等值線變得圓滑，各孔對(duì)拼接位置等值線錯(cuò)斷現(xiàn)象消失，但K35－K32孔對(duì)仍然保持近似直立的異常存在。

基于單孔來(lái)看，K35－K32剖面等值線形態(tài)與其他孔對(duì)的剖面形態(tài)區(qū)別不大，整體呈上高下低的趨勢(shì)，只是在多孔對(duì)拼接在一個(gè)剖面中才會(huì)發(fā)現(xiàn)該孔較相鄰兩側(cè)的異常。為了分析此類(lèi)現(xiàn)象是否僅由剖面拼接而造成的虛假異常，筆者收集了全區(qū)存在類(lèi)似異常的剖面，共5條，投影到平面位置上有明顯的對(duì)應(yīng)關(guān)系，詳見(jiàn)圖7。后期經(jīng)過(guò)驗(yàn)證為地下暗河，邊界與電磁波CT近直立異常一致（圖8）。

圖8 驗(yàn)證地質(zhì)剖面示意圖Fig.8 Schematic diagram of verified geological profile

綜上所述，采用反距離加權(quán)的局部多項(xiàng)式插值法對(duì)實(shí)測(cè)井間電磁波CT數(shù)據(jù)拼接是可行且有效的。

5 結(jié)論及討論

5.1 結(jié)論

1）通過(guò)對(duì)反演方法、初始場(chǎng)強(qiáng)兩個(gè)主要影響因素的分析認(rèn)為，電磁波CT相鄰孔對(duì)拼接誤差受反演方法及初始場(chǎng)強(qiáng)的影響，并不能從根本上完全消除，但是可以通過(guò)統(tǒng)一反演算法、反演參數(shù)、多次試驗(yàn)確定合理的初始場(chǎng)強(qiáng)值來(lái)最大限度地降低誤差。

2）目前常用的單孔對(duì)成像解釋能夠識(shí)別兩個(gè)鉆孔之間的大部分形態(tài)異常，通過(guò)反距離加權(quán)的局部多項(xiàng)式插值法進(jìn)行拼接后的多孔對(duì)電磁波CT剖面，其結(jié)果與原單孔剖面的擬合度較高，在拼接圖上，異常均能夠完全遵循原始剖面的特征，共用孔處等值線連續(xù)、過(guò)渡自然，成像結(jié)果可信度高。

3）單孔對(duì)成像在綜合成果解釋時(shí)對(duì)近直立異常，特別是鉆孔附近的直立異常無(wú)法識(shí)別，通過(guò)反距離加權(quán)的局部多項(xiàng)式插值法進(jìn)行拼接處理，可以準(zhǔn)確識(shí)別單孔對(duì)無(wú)法識(shí)別的近直立異常，同時(shí)多個(gè)剖面綜合解釋還可以準(zhǔn)確識(shí)別跨越多個(gè)孔對(duì)，有一定跨度的水平異常，對(duì)提高電磁波CT成果解釋具有實(shí)際應(yīng)用價(jià)值。

5.2 討論

影響多孔對(duì)電磁波CT剖面拼接的因素是復(fù)雜的，目前尚沒(méi)有成熟的解決方法，這些難題都需要在以后的工作中進(jìn)一步深入研究。本文多孔對(duì)電磁波CT剖面拼接是基于既有單孔對(duì)反演的成果數(shù)據(jù)進(jìn)行的，各單孔對(duì)經(jīng)反演參數(shù)統(tǒng)一和合理初始場(chǎng)強(qiáng)反演后，整體背景相差不大，且單孔對(duì)反演成果均能夠跟實(shí)際鉆孔吻合，采用反距離加權(quán)的局部多項(xiàng)式插值法拼接能夠達(dá)到理想的拼接效果，在此基礎(chǔ)上亦可以嘗試將拼接后的數(shù)據(jù)賦予三維空間位置坐標(biāo)，然后采用三維成像軟件進(jìn)行成像，結(jié)合三維切片初步實(shí)現(xiàn)三維展示。