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

（1.中國煤炭地質總局湖北煤炭地質局，湖北 武漢 430070；2.湖北煤炭地質物探測量隊，湖北 武漢 430200）

1 引言

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

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

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

2 電磁波CT成像基本原理

根據電磁波理論，電磁波傳輸過程中，耗散介質的初始場強E0與接收場強E關系為［7］：

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

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

假設R是n條電磁波射線在m個網格單元中的射線元組成的m×n維矩陣；X是電磁波吸收系數(shù)組成的矩陣，X＝（α1，α2，α3，α4，…，αn）T；D為接收場強的觀測值，D＝（d1，d2，d3，d4，…，dn）T，則所有觀測點的發(fā)射和接收之間得到的線性方程組可以寫成：RX＝D，通過合適的迭代算法求解上述線性方程組，就能夠得到兩個鉆孔間介質的吸收系數(shù)，進一步得到吸收系數(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）等，不同方法的成像結果各有區(qū)別［13］。

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

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

3.1 反演方法的影響

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

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

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

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

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

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

3.2 初始場強E0的影響

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

4 相鄰孔對的拼接

4.1 拼接思路

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

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

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

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

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

4.2 反距離加權的局部多項式插值法

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

4.2.1 局部多項式插值方程

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

文中采用二次插值，由n個采樣點的值可以得到方程組：

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

求解多項式系數(shù)矩陣，表達式如下：

通過最小二乘求解多項式系數(shù)：

式（9）中，P為每個采樣點通過反距離加權法確定的權重，由此求得局部多項式插值方程，進而得到內插點的數(shù)值。

4.2.1 反距離加權法

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

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

4.3 相鄰孔對的拼接

以上述貴州巖溶勘查項目為例，電磁波CT鉆孔孔對如圖4所示，根據鉆孔在平面上的位置關系，將92對鉆孔劃分為縱橫交錯的23條剖面，以下選擇兩種不同異常下的剖面拼接效果進行分析，一是團塊狀及不規(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對鉆孔原始剖面數(shù)據拼接成圖。由圖5可以看出，整個剖面數(shù)據的整體水平基本一致，但是在k7、k8、k9、k10鉆孔位置存在比較明顯的吸收系數(shù)突變現(xiàn)象，圖上等值線圖連續(xù)，這種現(xiàn)象無法從地質上進行解釋。圖5（b）為經過反距離加權的局部多項式插值法處理后成圖效果，通過對比可以看出，經處理后拼接的等值線變得圓滑，電磁波吸收系數(shù)異常的位置及大小均未發(fā)生明顯變化，但各孔對拼接位置等值線不連續(xù)的現(xiàn)象消失，各孔對之間等值線過渡自然，且與實際鉆孔情況吻合。

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

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

圖8 驗證地質剖面示意圖Fig.8 Schematic diagram of verified geological profile

綜上所述，采用反距離加權的局部多項式插值法對實測井間電磁波CT數(shù)據拼接是可行且有效的。

5 結論及討論

5.1 結論

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

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

3）單孔對成像在綜合成果解釋時對近直立異常，特別是鉆孔附近的直立異常無法識別，通過反距離加權的局部多項式插值法進行拼接處理，可以準確識別單孔對無法識別的近直立異常，同時多個剖面綜合解釋還可以準確識別跨越多個孔對，有一定跨度的水平異常，對提高電磁波CT成果解釋具有實際應用價值。

5.2 討論

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