多層巖土環境的溶洞智能識別方法研究

巴澤群 * 周玉鋒 周其勛 ** 周立成 * 劉逸平 * 劉澤佳 * 湯立群 *

*(華南理工大學土木與交通學院，亞熱帶建筑與城市科學全國重點實驗室，廣州 510640)

?(廣州廣佛肇高速公路有限公司，廣州 510289)

**(廣州市高速公路有限公司，廣州 510030)

喀斯特地貌廣泛分布在我國西南部和中部地區。這種地貌是由地下水對可溶性巖石進行溶蝕、沉積等作用形成的，因此容易產生地下溶洞等地質現象。在實際工程中，溶洞會影響樁基施工安全，嚴重時可能導致地面坍塌。因此，在施工前了解地下環境并排除溶洞等安全隱患，有助于確保施工安全并減少經濟損失。

目前，傳統溶洞探測方法包括鉆孔取樣法、彈性波CT 法[1]、地震反射波法[2-4]等機械類方法以及電磁波跨孔CT 法[5]、探地雷達法[1]、高密度電法[1]、電阻率層析法[6-8]、瞬態電磁場法[1,9]等電磁類方法。這些方法技術成熟且應用廣泛，但仍存在著一些缺點，如主要以二維截面探測為主，數據反演和解釋需要依賴專業知識和專家經驗等。亦有方法可進行三維溶洞形態反演[10]，但該方法建立在已對溶洞進行定位的前提下，且數據反演較復雜。因此，傳統溶洞探測方法在探測效率、經濟成本、時間成本和智能化方面亟待改進。

在地質勘探領域，已有部分學者利用智能方法開展了一些探索性的研究工作。例如，Zhang 等[11]使用U-NET 網絡實現地球物理勘探信息的智能反演；Di 等[12]基于支持向量機/多層感知機等機器學習方法實現地震數據的智能解釋和反演；Zheng 等[13]用主成分分析、K 聚類和線性判別式分析方法，使用聲速信號、電阻率、密度、孔徑數據等測井相信息進行地層類型的預測和識別；Xiong 等[14]使用深度自編碼網絡推斷礦物儲備情況；Huang 等[15]基于交叉小波變換處理振動時序信號獲得振動相位譜和功率譜，并通過卷積神經網絡來進行爆破位置的預測；Li 等[16]使用地震屬性數據結合隨機森林、遺傳算法–支持向量機等機器學習方法，實現了煤層地下結構的預測；Yang等[17]基于卷積神經網絡和支持向量機等算法，通過振動信號的時頻特征對地震源深度進行了識別。這些研究使用了人工智能方法代替人工數據反演與解釋，表明了人工智能方法在地質勘探領域的可行性和應用潛力。

本文作者團隊在前期對高效、智能的溶洞探測技術進行了初步研究，提出了一種適用于單一土層環境的基于機器學習算法的球形溶洞定量探測方法[18]。在此基礎上，本文將研究場景推廣到更符合實際工程的多層巖土環境，并設計了基于(bidirectional long short-term memory, BiLSTM)神經網絡的深度學習模型來實現多層巖土環境下的溶洞三維定量智能探測。本文技術路線如圖1所示，主要步驟如下：（1）借鑒淺層地震反射波法，建立含有單個無填充球形孔洞的多層巖土結構模型，并計算樁錘激振下地表的加速度響應信號；（2）基于應力波傳播的基本原理對響應信號進行分析以提取特征數據；（3）針對不同溶洞工況進行大量建模，構建用于深度學習的數據集；（4）設計并優化深度學習模型，實現多層巖土環境下溶洞的三維位置和大小的定量智能識別。本文提出的深度學習模型可不依賴于專業知識和專家經驗，為高效、智能化的溶洞探測方法的工程應用奠定基礎。

圖1 技術路線Fig.1 Technical roadmap

1 含溶洞的多層巖土結構有限元建模

在實際工程中，大部分地下巖土體為分層的多層結構，具有多樣的巖土層種類和復雜的巖土層分布，這導致多層巖土環境下對溶洞大小和位置的預測變得非常困難。由于獲得實測數據需要高昂的經濟成本和時間成本，本文采用有限元方法對多層巖土環境下的大量溶洞工況進行建模，從而構建不同溶洞工況下的響應數據庫，為基于深度學習的溶洞定量識別方法的探索奠定基礎。

1.1 有限元模型

本文的建模過程參考了廣東省廣州市廣佛肇高速公路SG03 段的地質勘探數據，根據其地勘數據對巖土層種類、巖土層組合以及各層巖土的厚度變化范圍進行建模。本文建立的有限元模型如圖2(a)所示。其中，整個巖土體模型為一個邊長30 m 的立方體，分為8 層，按照從上往下的順序進行編號，1～6 層為土層，7～8 層為巖層。為簡化建模過程，土層和巖層均設置為各向同性、均勻密實的材料，且土層采用Mohr–Columb 本構模型，巖層采用Drucker–Prager 本構模型，詳細參數見表1。同時，為簡化模型，假設土體內含一個無填充、真空的球形孔洞。激振源使用一個高度0.5 m，直徑1.5 m，重約6.8 t 的圓柱形樁，其材料采用線彈性本構，密度為7800 kg/m3，彈性模量為200 GPa，泊松比為0.3。激振形式采用樁錘自由下落的方式，撞擊地表時的速度設為4.42 m/s。如圖2(b)所示，為了獲取地表的加速度響應，以激振點為圓心，在距離其10 m 的采樣圓上均勻布置6 個加速度傳感器，按照順時針方向依次編號為1～6 號，且傳感器的采樣頻率為1000 Hz，采樣時間為0.6 s。值得注意的是，為了防止應力波在模型邊界上的反射，土體側面和底面設置為無限元邊界。

表1 巖土層分布示例Table 1 Example of rock-soil structure

圖2 有限元模型示意圖Fig.2 Diagram of finite element model

1.2 工況設置

如表2 所示，球形孔洞參數包括洞頂與激振線之間的水平偏移距離、洞頂深度和孔洞直徑。其中，洞頂水平偏移距離的取值范圍為1～10 m，取值間隔為1 m；洞頂深度取值范圍為5～21 m，取值間隔為3 m；孔洞直徑取值范圍為1～10 m，取值間隔為1 m。上述3 個取值變量依照取值要求隨機組合，生成大量涵蓋不同孔洞工況的模型。同時，如表3 所示，設置9 種典型的多層巖土結構，與表2 的孔洞模型隨機組合，即可構建9 種不同巖土結構下含有不同孔洞工況的海量模型。

表2 溶洞工況設置Table 2 Settings of karst cave parameters

表3 巖土層結構設置Table 3 Settings of rock–soil structures

通過統計，每個有限元算例在搭載了Intel Core i5-10400F 型CPU 的計算機上的計算時間平均為5 h 以上。因此，構建數據集需要大量的算力和時間開銷。針對有限元算例計算時間長的問題，本文將有限元算例進行等效工況擴充。具體地，如圖2(c)所示，虛線圓框代表溶洞，藍色圓和綠色圓分別代表兩組加速度傳感器，藍色對應x – y坐標系，綠色對應x′-y′坐標系?？梢钥闯?，以不同坐標系描述的同一個溶洞具有不同的空間位置。因此，對于一個有限元算例，本文在其加速度采樣的圓圈上均布了288 個節點，即每隔1.25°方位角布置一個加速度傳感器，進而通過旋轉變換采樣不同的傳感器組合（6 個為一組）的方式來擴充溶洞工況數據，即每一個有限元算例可以等效擴充為288 個工況。最終，每種巖土結構對應的溶洞工況數大致相同，平均約包含2880 條數據工況，總數據量為25 344 條。

2 特征工程

為了提高深度學習模型的識別效果，本文基于應力波反射的基本原理，對地表加速度響應信號進行特征工程，力求挖掘與預測目標相關性強的特征并作為模型的輸入，從而提高模型的預測精度。

地下多層巖土結構中存在許多波阻抗界面，如巖土層分界面、溶洞表面等，當應力波遇到這些界面時，會發生反射現象。因此，可通過計算特征界面的應力波反射信號的行程與時間，捕捉受溶洞和地下巖土體結構影響的加速度信號，進而進行特征工程。

首先，對于不含水的各向同性致密巖土體，其縱波波速的計算公式為

式中，vp為縱波波速，ν為泊松比，E為彈性模量，ρ為密度?；谑剑?），結合巖土體結構中各層的厚度，可以估算不同界面應力波反射信號的行程和走時。應力波傳播路徑示意圖如圖3所示，其中圖3(a)展示了應力波被激發后傳播至巖土層交界面被反射回地表并被傳感器接收的過程；圖3(b)展示了應力波被激發后傳播至溶洞頂點被反射回地表并被傳感器接收的過程?？疾鞄r土層界面反射情況時，取傳感器與激振點連線的中位線與各巖土層交界點的位置作為反射參考點，以計算應力波的行程?？疾烊芏幢砻娣瓷淝闆r時，近似取溶洞洞頂作為反射參考點進行計算。以一列應力波的傳播過程分析作為示例，對于任意傳感器m，其自激振作用起到接收到巖土層界面反射信號的時間可按式（2）進行計算

圖3 應力波路徑示意Fig.3 Diagram of stress wave propagation

對于任意傳感器m，其自激振作用起到接收到溶洞洞頂反射應力波的時間可按公式(3)進行估算

其中，x，y，z代表溶洞洞頂反射點的空間坐標；xr，yr和zr代表巖土層反射點r 的空間坐標；xm，ym和zm代表傳感器m 的空間坐標。hi代表應力波經過巖土層i的垂直距離，vpi代表巖土層i中縱波的波速?；谑剑?）和式（3），本文計算了各個巖土層交界面和溶洞洞頂的應力波反射信號的接收時間，如表4 所示。其中包含了洞頂深度為2 m，6 m 和10 m 的溶洞模型的計算結果，所含溶洞的半徑均為2 m。圖4(a)給出了不同溶洞工況下的多層巖土體結構加速度響應，其中黑色曲線為無溶洞情況下的多層巖土結構的加速度響應時程曲線，紅色、藍色和綠色曲線為溶洞深度遞增的多層巖土結構的響應時程曲線，且4 條圖線所對應的巖土體分層情況一致。

表4 反射信號走時Table 4 Travel time of reflected stress wave

圖4 加速度響應信號對比Fig.4 Contrast of acceleration response signals

綜合表4 和圖4(a)的結果可知，巖土層界面反射信號和溶洞洞頂反射信號的接收時間與加速度時序信號的波峰波谷時間接近，雖然無法準確說明每個波峰波谷包含哪些反射波成分，但是可以推斷在加速度曲線的極值數據中蘊含著與巖土層界面位置和溶洞深度等相關的信息。同時，若溶洞直徑增大，其反射面增大，響應信號波峰波谷對應的加速度大小以及所對應的時間會相應發生變化。因此，可以提取加速度信號的極值數據作為輸入數據的特征。具體地，將每個工況下的6 個傳感器的時序信號的極值點全部提取出來，并按照出現的先后順序進行排列，形成極值序列數據。通過統計，平均每個加速度時序信號包含19 個極值點，故極值特征的數據格式為[19, 6]。由以上的分析可知，極值特征與溶洞的深度和大小有較強的相關性。

圖4(b)展示了對應于一個含有溶洞的多層巖土結構的一組傳感器的時序信號，其中溶洞洞頂的水平偏移距離為5 m，溶洞半徑為2 m。觀察不同傳感器的時序信號可知，由于傳感器空間方位的不同，傳感器接收的響應信號也有一定的差異，這種差異來源于傳感器與溶洞的不同距離?？梢酝茢?，不同空間方位的傳感器采集的加速度時序信號蘊含著溶洞的位置特征信息，即可以提取一組6 個不同空間位置傳感器的時序信號，作為輸入數據的時空特征。具體地，首先提取每個傳感器采樣信號的0.15～0.60 s 時間段的加速度數據，然后按照傳感器編號順序即空間位置順序排列時序信號，得到格式為[6, 451]的時序數據作為后續深度學習模型的輸入。

3 多層巖土環境的溶洞智能識別模型

深度學習是人工智能研究領域的重要方法，其基本原理是通過搭建人工神經元組成的神經網絡，建立輸入與輸出的映射關系，解決各種模式識別問題。目前，深度學習方法廣泛應用于文本翻譯、圖像識別、時間序列預測等領域。其中，BiLSTM 作為循環神經網絡（recurrent neural network, RNN）的一種改進模型，有著強大的時間序列數據學習能力?？紤]到本文提取了以時間序列數據為主的輸入特征，因此引入BiLSTM 網絡，以BiLSTM 網絡為核心模塊設計了雙數據通道的分離架構網絡模型，對多層巖土環境下的不同溶洞工況進行定量識別。

3.1 BiLSTM 的原理

BiLSTM 由兩條方向相反的長短期記憶（long short-term memory, LSTM）神經網絡構成，其結構如圖5(a)所示，可以從正反兩個方向學習時間序列的特征。LSTM 網絡的基本單元為LSTM 神經元，其基本原理如圖5(b)所示。其中，LSTM 神經元由遺忘門ft，輸入門it，輸出門ot和細胞臨時狀態Ct～組成。LSTM 神經元將上一時刻的輸出值ht?1和本時刻的輸入值xt進行矩陣拼接，并將拼接矩陣激活后得到遺忘門、輸入門、細胞臨時狀態和輸出門。遺忘門控制上一時刻細胞狀態Ct?1的信息流入；輸入門決定了本時刻的細胞臨時狀態Ct～對本時刻細胞狀態Ct的影響程度；輸出門決定了LSTM 神經元的輸出值。將LSTM 神經元連接起來，通過細胞狀態的傳播，即可實現序列數據的長時記憶，從而全面、完整地提取時間序列的特征。

圖5 BiLSTM 示意圖Fig.5 Diagram of BiLSTM

3.2 網絡架構

本文采用BiLSTM 網絡作為多層巖土環境下溶洞識別深度學習模型的核心模塊，基于該網絡設計了如圖6(a)所示的雙數據通道融合架構、如圖6(b)所示的雙數據通道分離架構和如圖6(c)所示的單數據通道架構。如圖6(a)所示，兩支通道的輸入數據分別為時序信號序列和極值特征，其分別通過BiLSTM 網絡層提取特征，隨之與多層全連接網絡連接，在相鄰網絡層間設置隨機失活層（dropout）以避免網絡過擬合，最后進行兩個通道特征矩陣的拼接融合，并通過全連接層得到最后的預測結果。預測結果為4 個標簽X,Y,Z和D，分別對應了溶洞洞頂的三維坐標和溶洞直徑。其網絡各層神經元個數及激活函數設置如表5 所示。如圖6(b)所示，雙通道分離架構的輸入數據分為兩支，分別為時序信號序列和極值特征，兩支數據通道的訓練同步進行。且兩支數據通道的網絡層神經元權重、偏置等網絡參數的更新優化互不影響，其中以時序信號序列為輸入的分支對應的輸出為溶洞洞頂水平坐標X,Y；以極值特征為輸入的分支對應的輸出為溶洞洞頂深度Z和溶洞的直徑D。該網絡架構的各個網絡層神經元個數與激活函數設置如表6 所示。圖6(c)上部展示了以時序信號序列為輸入的單數據通道架構，其通過BiLSTM 層提取特征后，隨之連接多層全連接層，同時在相鄰全連接層間設置隨機失活層（dropout），最終輸出溶洞洞頂的空間坐標X，Y，Z和溶洞的直徑D；圖6(c)下部為以極值特征為輸入的單數據通道架構，其網絡設置與以時序信號序列為輸入的單通道架構相同。單通道架構的各層神經元個數與激活函數設置如表7所示。

表5 雙通道融合架構參數設置Table 5 Settings of dual-channel fusion architecture

表6 雙通道分離架構參數設置Table 6 Settings of dual-channel separated architecture

表7 單通道架構參數設置Table 7 Settings of single-channel architecture

圖6 不同網絡架構Fig.6 Different model architectures

3.3 模型訓練

基于9 種巖土結構下共25 344 個工況的數據集，以8:2 的數據量比例劃分訓練集與測試集。其中，訓練集數據量為20 275 條，測試集數據量為5069 條。從訓練集中劃分出10%的數據作為驗證集，并在訓練集上訓練數據，在驗證集上觀察模型訓練效果，在測試集上檢驗預測效果。在訓練過程中，選擇平均平方誤差（meansquare error, MSE）作為損失函數，其表達式為

其中，ypi表示第i個數據的預測值，yti表示第i個數據的真實值，n代表總數據量。網絡優化器選擇Adam 優化器，其算法由Kingma 等[19]提出，能夠實現自適應學習率調節，并且初始學習率設置為0.001。網絡訓練輪次上限設置為500 次，收斂條件為：若連續20 次訓練迭代后，驗證集的MSE 減小值小于0.1，則模型達到收斂，訓練終止。平均每個網絡模型訓練10 次，每次訓練結果相互獨立。

圖7(a)展示了雙通道融合架構訓練過程中損失函數的變化情況，可以看出損失函數穩定下降，收斂程度良好；圖7(b)展示了雙通道分離架構的情況，其損失函數下降速度先快后慢，最終收斂；圖7(c)和圖7(d)分別展示了單通道架構（時序信號序列）和單通道架構（極值特征）的訓練情況，其中單通道架構（極值特征）的收斂狀態下的損失函數值約為3，相比其他網絡訓練效果較差。由圖7 可知，4 種網絡架構模型都能較好收斂。雙通道架構的收斂效果更優，以雙通道分離架構網絡模型的收斂程度為最好。

圖7 各網絡訓練效果可視化Fig.7 Visualization of training process

3.4 評價標準

本文提出的多層巖土環境下溶洞識別模型的識別目標為溶洞的三維位置和大小，因此評價網絡性能的標準是溶洞參數的預測準確率，對溶洞預測準確與否的標準所示為

其中，r表示容許誤差，下標pred 和true 分別表示預測值和真實值。當預測結果同時滿足式（5）～式（8）時，可判定模型完成了溶洞工況的準確預測；當預測結果滿足上述任一公式條件時，可判定該模型對于溶洞的對應參數預測準確。模型的預測準確率計算式為

其中，Ntrue表示預測準確的樣本數，Ntotal表示參與預測的總樣本數。

4 結果與討論

本節對基于BiLSTM 的雙數據通道分離架構、雙數據通道融合架構和單數據通道架構模型開展對比分析，各模型在3 m 容差下的預測結果如表8 所示。其中，雙通道分離架構網絡在3 m 容差下的綜合預測準確率最高，達到了98.23%，高于雙通道融合架構網絡的95.95%綜合預測準確率、以時序信號作為輸入的單通道架構網絡的95.51%綜合預測準確率和以極值特征為輸入的單通道架構的17.99%綜合預測準確率。對比雙通道分離架構和雙通道融合架構的預測結果，可知雙通道分離架構充分利用了時序信號與溶洞坐標X和Y的相關性以及極值特征與溶洞深度Z和溶洞大小D的相關性，且雙通道分離架構無需特征融合，規避了雙通道融合架構特征融合困難、依賴調參的缺點。觀察雙通道分離架構網絡模型與以極值特征為輸入的單通道架構網絡模型的預測結果，可知極值特征與溶洞深度Z和溶洞大小D的相關性較強，與溶洞水平位置X和Y的相關性較弱。觀察基于時序信號的單通道架構網絡綜合預測結果，其預測準確率相較雙通道融合架構網絡僅有不到0.5%的差距；而與雙通道分離架構網絡的綜合預測準確率相比，存在約2.7%的差距?？梢酝茢鄷r序信號包含了足夠的信息，具有不同水平方位、深度和大小的溶洞對溶洞表面反射應力波的行程和走時會產生影響，因此通過對不同空間位置加速度傳感器信號進行完整采樣，能夠從中挖掘溶洞的三維位置和大小特征，從而實現對溶洞的三維位置和大小的較高精度的識別，而加入極值特征有助于提高識別準確率，說明了本文特征工程的有效性。

表8 預測結果Table 8 Results of prediction

綜上所述，采用時序信號和極值特征的雙通道輸入有助于提高模型預測準確率，其中，雙通道分離架構模型能充分利用輸入特征與相應識別目標的強相關性，因此其能達到最優的預測準確率。

5 結論

本文提出了一種基于深度學習的多層巖土環境的溶洞智能三維定量探測方法，所得到的主要結論如下：（1）地表加速度響應時序信號包含足夠的信息，以時序信號為輸入的單通道架構模型在3 m 容差下對溶洞位置和大小的識別準確率達到了95.51%，加入極值特征可提高模型識別準確率；（2）基于BiLSTM 的雙通道分離架構網絡模型在3 m 容差范圍內對多層巖土環境下溶洞的三維位置坐標和直徑大小的綜合預測準確率達到了98.23%，實現了多層巖土環境下溶洞的定量智能預測。本文研究的不足之處在于：（1）巖土結構和溶洞的組合工況可能性極多，本文受到算力的限制，只能在一些有限的組合工況下進行方法的驗證；（2）僅進行了基于有限元數據的可行性研究，未在實測數據上對方法進行驗證。