管道穿越巖溶段綜合探測方法研究及應用

DOI：10.3969/j.issn.1001-2206.2024.06.007

摘" " 要：管道定向鉆穿越工程需提前查明巖溶裂隙地層發育狀況，為穿越工程設計提供依據?？缈讖椥圆▽游龀上窈碗姶挪▽游龀上穹椒ㄔ趲r溶裂隙調查中具有明顯優勢，但存在單一反演多解性和局限性的問題。提出了一種新的彈性波CT和電磁波CT聯合反演方法，實施過程簡單，反演收斂穩定，克服了單一反演多解性和局限性的問題。該方法的技術路徑是將電磁波觀測場強數據等效轉換為彈性波旅行時，然后同彈性波數據進行聯合迭代反演。數值模擬分析表明，能夠對直徑為孔距1/10的溶洞實現準確成像，比單獨反演精度更高。實際應用中，對16對彈性波數據和13對孔電磁波數據實施了三維聯合CT反演，對巖溶裂隙發育的空間展布進行了精細刻畫，取得了較好的效果。

關鍵詞：管道穿越；巖溶探測；電磁波CT；彈性波CT；聯合反演

Abstract：Early identification of the development status of karst fissures is necessary for the pipeline crossing project， providing a reference for the design of crossing engineering. Though cross-hole elastic wave tomography and electromagnetic wave tomography methods have obvious advantages in investigating karst fissures， they are limited by the simplification and multiple solutions of the single inversion. This article proposes a new joint inversion method involving elastic wave CT and electromagnetic wave CT， which overcomes the above problems with an easy operation process and stable convergence.The technical path of this method includes equivalently converting electromagnetic wave observation field strength data into elastic wave travel time and performing joint iterative inversion with elastic wave data. Numerical simulations prove its ability to more accurately image karst caves having a diameter of only one-tenth of the hole spacing than single inversion. In practical applications， the three-dimensional joint CT inversion is conducted based on 16 pairs of elastic wave data and 13 pairs of hole electromagnetic wave data， accurately depicting the space distribution of karst fissures development.

Keywords：pipeline crossing; karst detection; electromagnetic wave CT; elastic wave CT; joint inversion

定向鉆穿越施工一旦遇到巖溶裂隙地層，會有大量泥漿漏失，致使施工過程出現鉆具卡鉆、鉆桿斷裂、回拖卡管、管道變形等極端風險。為保障巖溶裂隙地層穿越施工的安全性，提高穿越工程的成功率，盡量規避管道穿越施工的風險，保證管道穿越一次性成功回拖，查明溶洞裂隙發育的大小、形狀、位置及充填屬性十分重要。

當前用于巖溶裂隙探測的方法主要有鉆探法和物探法。鉆探法只能進行點樣勘查，無法做到快速大面積調查。物探法包括地面物探和鉆孔物探。地面物探法主要有：直流電法中的高密度電阻率法，電磁法中的瞬變電磁法、探地雷達法，地震波法中的折射波法、反射波法、面波法等[1-2]。

由于巖溶裂隙不良地質體賦存屬性（形態、位置、充填）非常復雜，常規地面物探方法難以取得較好效果，目前主要采用高分辨率跨孔層析成像（Computed Tomography，CT）探測技術，包括電磁波CT、電阻率CT、彈性波CT等[3-5]。各種方法的層析反演原理基本相同，但利用的物理場信息不同，因此各種跨孔CT方法具有不同特征的異常表征能力。其中，跨孔彈性波CT和跨孔電磁波CT的適用場景更為廣泛且分辨率更高。根據彈性波反演的速度場和電磁波反演的吸收系數分布，結合地質資料進行綜合分析，可以有效解釋斷裂破碎帶的邊界與產狀、巖土分界面，以及溶洞的邊界、產狀、發育與分布情況[6-10]。

當前，跨孔彈性波CT和跨孔電磁波CT反演方法已比較成熟。當采用彈性波CT和電磁波CT聯合勘探時，主要技術路徑是采用獨立反演，但容易出現因結果矛盾而導致難以解釋的問題。對反演結果進行數據融合是一種傳統的改進方法，能夠在一定程度上提高解釋準確度。2019年陳湘華等[11]對于電磁波和彈性波層析成像探測的聯合分析方法進行了研究。

此外，近年來有學者對跨孔彈性波和跨孔電磁波聯合CT反演方法開展了研究，試圖改善單一反演方法多解性和局限性的問題。20世紀70年代，Vozoff等[12]首次提出地球物理聯合反演理論，對直流電阻率法和大地電磁法展開了聯合反演工作。Gallardo等[13]首次利用交叉梯度函數來構建不同物性參數之間的結構耦合關系，令交叉梯度函數為零，來達到不同模型空間結構相似性一致的目的。李桐林等[14]進行了部分區域約束下的交叉梯度多重地球物理數據聯合反演研究。張镕哲等[15]開展了大地電磁、重力、磁法和地震初至波走時的交叉梯度聯合反演研究。2022年師學明等[16]開展了跨孔電磁波與地震波CT交叉梯度聯合反演算法研究與應用，取得了一定的成效。

然而，基于交叉梯度理論的聯合反演仍然存在諸多難題，如收斂速度慢、抗干擾能力不高等問題。為此，本文提出一種新的彈性波和電磁波聯合CT反演方法。其基本思路是將電磁波觀測場強數據等效轉換為彈性波旅行時。在每次迭代反演過程中，遍歷所有數據，達到相互約束，共同修正模型速度參數，實現聯合反演的目的。首先，本文詳細推導了該方法原理基本公式；然后，開展了有效性檢驗的數值模擬分析；最后，在平度輸油管道定向鉆穿越巖溶段的探測中進行了應用。

1" " CT成像方法原理

按照波的傳播理論劃分，層析反演分為射線層析和波動方程層析，后者又稱為散射層析（全波形反演）。射線層析的理論基礎穩定、方法比較簡單、干擾因素較少、應用效果好；盡管信息量比較少，如果采用誤差較小的反演算法，并且能夠充分利用可觀測空間和介質的先驗信息，仍能夠獲得理想的效果。波動方程層析理論能夠根據地震波旅行時、振幅、相位、頻率等因素獲得地層介質更多的信息量，提高了分辨率；但是，實際應用中卻存在很多問題，比如散射數據的提取較為困難，而且不易消除影響波形的干擾因素等，所以應用上不如射線層析反演廣泛[17]。

為便于闡述本文提出的聯合反演方法，下文首先簡要介紹基于射線理論的彈性波走時CT成像和電磁波吸收衰減CT成像原理，在此基礎上導出聯合反演數據處理方法流程。

1.1" " "彈性波走時CT成像原理

彈性波走時CT成像是根據接收到的初至旅行時數據來反演該剖面的速度分布v（x， y）或慢度s（x， y）=1/v（x， y）。在高頻近似下，彈性波傳播路徑近似為射線，假設第i條彈性波傳播路徑為Li，其旅行時為Ti，如圖1所示，則有如下彈性波路徑和走時關系式：

式（1）是一曲線積分，ds為射線弧長微元，v（x， y）和Li均為未知，Ti為已知。這實際上是一個非線性問題。在速度場變化不大的情況下，可以把射線路徑近似看作是直線，即Li為直線，實際上地下地質情況是復雜的，射線路徑也往往是曲線。現將反演區域離散化，如圖1所示，假如離散化后的單元數目為N。第1至第N單元的慢度依次記為s1、s2、…、sN。這樣，第i條射線的旅行時Ti表示為：

式中：rij是第i條射線穿過第j個網格的射線長度。 當有大量射線（如M條射線）穿過反演區域時，根據式（2）就可以得到關于未知量sj（j =1， 2，…， N）的M個方程（i =1， 2，…， M），M個方程組合成一線性方程組，即：

寫成矩陣形式如下：

式中：R=（rij）M×N稱作距離矩陣；T=（Ti）M×1為旅行時向量，即檢波器接收得到的初至旅行時；S=（si）N×1為慢度列向量。

通過求解方程組（4）就可以得到離散慢度分布，從而實現跨孔區域的速度場反演成像。

1.2" " 電磁波吸收衰減CT成像原理

跨孔電磁波CT由發射端、接收端和主機3部分組成。發射端有效電磁波初始輻射場強為E0，與發射機輻射功率、天線周圍介質有關；接收端觀測點場強為Er，觀測沿孔軸方向的電場分量，二者關系如下：

式中：β為介質對電磁波吸收系數，r為收發點距離，f（θ）為發射天線方向性因子，θ為發射天線（激發孔軸）與電磁波射線路徑的夾角，φ為接收點處天線（接收孔軸）與電場方向的夾角。

跨孔電磁波測量采用的是半波對稱偶極天線，因此方向性因子通常按照均勻介質中的公式f（θ）=cos（π/2cosθ）/sinθ計算。一般地，激發和接收天線相互平行，則有θ和φ互為余角。則式（5）簡化為：

式中：[f] （θ）=cos（π/2 cosθ）。

為更好與實際生產接軌，對式（6）進行一定處理，把公式兩端都取10為底的對數，然后乘以20，則式（6）改為：

將式（7）左端定義為觀測點場強讀數，用Di表示；右端第一項定義為初始場強分貝數，用D0表示。實際場強大小用測點電壓來衡量，因此場強E0和Er計量單位為V，換算為分貝后，Di和D0計量單位為dBV，吸收系數β的單位為dBV/m。如圖1所示，根據式（7）得到第i條射線場強觀測值公式為：

式中：ri為第i條射線的長度。其中，D0是未知的，是影響反演精度的重要參數，在固定頻率測量過程中通常是變化不大的，因此可將其視作常數。實踐中，可先求出初始場強D0，再進行層析反演。

線性擬合法是求取初始場強廣泛應用的有效方法之一。在均勻介質假設條件下，式（8）可以寫成變量為射線長度ri，截距為初始場強D0的線性關系式：

式中：Mi=Di-20lg[[f]（θ） ri-1]，A=-20（ln10）-1[β]。利用式（9），選擇無地下異常體影響且只反映背景值的數據進行線性擬合，即可求出初始場強D0，也可獲得背景場吸收系數[β]。

按圖1所示非均勻介質條件離散化，把式（8）重新整理可得：

式中：rij是第i條射線穿過第j個網格的長度，βj（j=1， 2， …，N）是第j個網格的吸收系數。

式（10）右端定義為第i條射線場強衰減量，用Ui表示，即：

當有大量射線（如M條射線）穿過反演區域時，根據式（11）就可以得到關于未知量βj（j=1， 2， …，N）的M個方程（i=1， 2， …，M），M個方程組合成一線性方程組為：

寫成矩陣形式如下：

式中：R=（rij）M×N稱作距離矩陣；U=（Ui）M×1為衰減量向量，即接收天線測量校正后得到的場強衰減量；β=（βi）N×1為吸收系數列向量。

通過求解方程組（13）就可以得到離散系數分布，從而實現跨孔區域電磁吸收系數場反演成像。

1.3" " 彈性波和電磁波聯合CT反演

在傳統彈性波和電磁波獨立層析成像原理的基礎上，提出了一種聯合層析成像的算法，具體步驟如下。

在均勻介質和高頻近似條件下，式（2）和式（11）可寫為：

式中：[s]為成像區域的平均慢度，[β]為成像區域的平均吸收系數。

由上述兩個公式可得第i條彈性波射線走時表達式，即：

利用式（16）可以將跨孔電磁波的衰減量轉換為彈性波走時，轉換后再進行層析反演，可獲得慢度量綱的場分布。式中的平均吸收系數[β]可以根據式（9）線性擬合求得，平均慢度[s]可按下式計算。

按式（16）將電磁波衰減量進行轉換后，聯合彈性波走時數據，可進行聯合反演，達到相互約束，提高成像分辨率的目的，聯合反演的結果為慢度量綱。

值得注意的是，在層析成像過程中距離矩陣R往往為大型無規則的稀疏矩陣（R中每行都有N個元素，而射線只通過所有N個像元中一小部分），而且常是病態的。實際應用中要反復求解方程組，來得到重建區域的慢度場，本文采用廣泛使用的聯合迭代重建技術（SIRT）進行反演，其迭代初始值使用反投影重建技術（BPT）求取[18]。

2" " 彈性波和電磁波聯合CT反演數值模擬

2.1" " 模型參數

為了驗證本文提出的聯合CT反演方法有效，設計了如圖2所示的跨孔數值模型，模型跨孔距離為30 m，孔深為50 m，模型中含有兩個直徑均為3 m的圓形異常體，其圓心坐標分別為（x1=15 m，z1=15 m）和（x2=15 m，z2=35 m）。淺部為充氣溶洞，深部為充水溶洞，用于檢測不同方法的探測分辨能力。模型詳細物性參數如表1所示。

2.2" " 波場模擬

在所建模型基礎上，彈性波和電磁波傳播過程均采用有限差分法模擬，激發間距和接收間距均為1 m，總激發點數51炮，總接收點數51道。彈性波模擬中，采用主頻800 Hz的雷克子波作為震源子波，采樣間隔0.2 ms，記錄長度40 ms。電磁波模擬中，采用主頻6 MHz的雷克子波為脈沖源，激發和接收均為Y方向電場（Ey），采樣間隔0.25 ns，記錄長度700 ns。

由于模型中的2個溶洞深度分別為15 m和35 m，因此圖3所示記錄所對應的激發點深度分別為15 m和35 m，正好與溶洞中心深度相同。

在彈性波記錄中，拾取所有地震道的縱波初至時間，如圖4所示。對電磁波記錄進行包絡計算，取包絡最大幅值作為觀測場強，所有激發接收點的場強大小如圖5所示。可以看出，不同充填類型的溶洞對彈性波和電磁波的影響是不一樣的，這是CT反演識別異常的基礎。電磁波模擬中，激發和接收采用的都是Ey分量，因此在估算初始場強時，取方向性因子f?（θ）=1。然后根據式（16），將電磁波場強轉換為圖6紅線所示的等效初至時間。

2.3" " 聯合CT反演

對彈性波初至時間和電磁波電場強度分別按走時層析和相對衰減層析的反演結果如圖7（a）和圖7（b）所示。結果顯示，彈性波的CT分辨率比電磁波分辨率高；彈性波對充氣溶洞更敏感，而電磁波對充水溶洞更敏感，這與彈性波場和電磁波場的傳播特性是一致的。此外，盡管空氣的電磁波吸收系數很小，但由于電磁波在充氣溶洞壁上產生散射，導致透射能量減小，使得電磁波反演結果顯示充氣溶洞為高吸收系數，這與模型正好相反?？梢?，電磁波CT無法區分溶洞是充氣還是充水。采用圖7（c）所示聯合CT反演方法，對充水和充氣溶洞均得到較好成像，這對提高巖溶地層探測分辨率具有重要意義。

3" " 管道穿越巖溶段探測應用

為了查明某輸油管道工程平度市黃同河穿越區巖溶及裂隙發育情況，為后期設計、施工提供基礎資料，采用了跨孔彈性波CT和跨孔電磁波CT綜合勘探的方法。

根據鉆孔資料揭示，場地在勘察深度范圍內的地層主要為第四系全新統-上更新統陸相沉積砂土及粉質黏土，下部為下元古代荊山群（Pt1j），上覆一定厚度的人工填土。包含如下地層類型：①素填土，②粉質黏土，③粗砂，④層全風化大理巖，⑤層強風化大理巖，⑥層中風化大理巖，⑦層中風化片麻巖。部分鉆孔揭露存在巖溶裂隙發育情況。

本次探測在穿越段設計了鉆孔8個，根據孔位實施了16對跨孔彈性波CT、13對跨孔電磁波CT的數據采集。鉆孔及CT剖面布設如圖8所示。

在上述二維數值模擬的基礎上，為進一步檢驗聯合CT反演策略的有效性，選取ZK1-ZK2剖面進行不同方法CT反演對比，結果見圖9。

由圖9可知，彈性波和電磁波CT數據采集時，均將ZK1作為激發孔，ZK2作為接收孔。相比圖9（b）所示的電磁波CT結果，圖9（a）所示的彈性波CT分辨率更高。鉆孔揭露巖溶分布情況與彈性波CT基本吻合，與電磁波CT結果不完全一致。圖9（c）所示的聯合CT反演結果，與鉆孔揭露情況一致性更好。

數據處理過程：首先將彈性波和電磁波分別進行獨立的二維CT反演和聯合CT反演，然后結合鉆孔資料檢查反演結果的正確性。根據鉆孔資料對所有鉆孔反演結果進行分析，經統計得到如下認識：素填土、粉質黏土的縱波波速Vp1000～1 500 m/s，電磁波相對吸收系數0.5～0.6 dB/m；全風化大理巖縱波波速Vp 2 000～3 000 m/s，電磁波吸收系數0.3～0.5 dB/m；中風化片麻巖、中風化硅質大理巖、中風化大理巖的縱波波速Vp 3 500～5 000 m/s，電磁波吸收系數0.2～0.4 dB/m；各風化層的波速和電磁波吸收系數存在極為明顯的差異。巖溶裂隙帶的縱波波速，因充填物和充填程度的不同在1 500～2 500 m/s之間，電磁波吸收系數在0.4～0.6 dB/m之間。

本文提出的聯合反演策略不僅適用于二維反演，也可用于三維反演中。為了獲得更加全面準確的巖溶空間分布信息，利用本文提出的聯合反演方法，對本應用采集獲得的所有跨孔數據開展了多孔三維聯合CT反演。圖10為三維聯合反演結果在高程80、85、90、95 m的水平切片圖。

反演結果揭示巖溶主要在東部ZK1、ZK2、ZK3和ZK4較發育，并在高程95 m水平形成連通。可見設計管道線路穿越溶洞發育區，因此建議此處實施定向穿越施工前，進行局部注漿處理。

4" " 結論

1）本文提出了一種新的彈性波CT和電磁波CT聯合反演方法，實施過程簡單，反演收斂穩定，克服了單一反演多解性和局限性的問題。數值模擬分析表明，能夠對直徑為孔距1/10的溶洞實現準確成像，比單獨反演精度更高。

2）在平度穿越工程巖溶探測中，對16對彈性波數據和13對孔電磁波數據實施了三維聯合CT反演，對巖溶裂隙發育的空間展布進行了精細刻畫，取得了較好的應用效果。

3）基于本文提出的聯合CT反演策略，今后可以拓展更多類型的CT源數據進行聯合反演，有望進一步提高探測精度。

參考文獻

[1]" 張健，馮旭亮，岳想平. 綜合物探方法在隱伏巖溶探測中的應用[J]. 物探與化探，2022，46（6）：1 403-1 410.

[2]" 張學亮，謝濤，周煒，等. 等值反磁通瞬變電磁和微動勘探在淺部巖溶探測中的應用[J]. 煤田地質與勘探，2023，51（12）：157-166.

[3]" 趙武陽. 跨孔地震波層析成像在巖溶探測中的應用研究[D]. 桂林：桂林理工大學，2021.

[4]" 趙威.電磁波CT幾種常用成像方法應用效果對比[J].工程地球物理學報，2019，16（5）：749-754.

[5]" 李陽陽. 基于測井約束反演的跨孔電阻率CT在城市巖溶探測中的應用[D]. 濟南：山東大學，2020.

[6]" 王運生，王家映，顧漢明. 彈性波CT關鍵技術與應用實例[J]. 工程勘察，2005，33（3）：66-68.

[7]" 羅術，金俊俊，甄大勇，等. 基于數值模擬分析的彈性波CT巖溶探測能力研究與應用[J]. 工程地球物理學報，2023，20（3）：330-336.

[8]" 陳春飛，沈曉武，張秉政. 基于電磁波層析成像技術的巖溶探測正演模擬及應用研究[J]. 工程地球物理學報，2021，18（1）：98-106.

[9]" 朱鑫磊，楊磊，馮光福，等. 地磁波CT和微動技術在盾構穿越巖溶地層中的綜合應用研究[J]. 工程地球物理學報，2022，19（5）：619-629.

[10] 王薇，鄧小虎，金聰，等. 電磁波CT揭露重大工程巖溶發育特征——以某地鐵巖溶勘察為例[J]. 科學技術與工程，2020，20（34）：13977-13982.

[11] 陳湘華，王啟明.基于電磁波和彈性波層析成像探測的聯合分析方法[J].科學技術與工程，2019，19（16）：304-312.

[12] VOZOFF K，JUPP D L B.Joint inversion of geophysicaldata[J].Geophysical Journal Royal Astronomical Society，1975，42（3）：977-991.

[13] GALLARDO L A，MEJU M A. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints[J]. Joural of Geophyssical Research：Solid Earth，2004，109（B3）：3311-3315.

[14] 李桐林，張镕哲，樸英哲. 部分區域約束下的交叉梯度多重地球物理數據聯合反演[J].地球物理學報，2016，59（8）：2979-2988.

[15] 張镕哲，李桐林，鄧海，等.大地電磁、重力、磁法和地震初至波走時的交叉梯度二維聯合反演研究[J].地球物理學報，2019，62（6）：2139-2149.

[16] 師學明，商祥，柳思龍.跨孔電磁波與地震波CT交叉梯度聯合反演算法研究及應用[C]//2022年中國地球科學聯合學術年會論文集. 北京：中國地球物理學會，2022：19-22.

[17] 何云川，黃金強. 基于初至波層析的全波形反演[C]//2022年中國地球科學聯合學術年會論文集. 北京：中國地球物理學會，2022：83-86.

[18] 楊艷，秦克偉，張東，等.一種改進的近地表層析成像SIRT算法[J].武漢大學學報（理學版），2009，55（2）：201-205.

基金項目：中石化石油工程設計有限公司科研課題“巖溶裂隙地段定向鉆勘察及施工技術應用研究”（KY2023）。

作者簡介：牟曉東（1972—），男，山東濰坊人，高級工程師，2004年畢業于中國海洋大學環境工程專業，碩士，現主要從事巖土工程勘察、工程物探、地質災害勘察與治理，以及巖溶發育區定向鉆穿越場地適宜性評價。Email：mouxd1972@163.com

收稿日期：2024-08-19