朱智清，王道坤，汪偉民，吳殿偉，施衛(wèi)力，張 賽

(1.機械工業(yè)勘察設計研究院有限公司，陜西西安 710043；2.安徽省皖北煤電集團有限責任公司，安徽宿州 234000；3.中鐵上海工程局集團有限公司，上海 201900；4.安徽省建設工程測試研究院有限公司，安徽合肥 230051)

0 引言

隨著我國經(jīng)濟的高速發(fā)展，各種類型的城市建筑物數(shù)量逐年增長，城市土地資源緊缺，城市地下空間利用與開發(fā)是緩解這一狀況的重要措施(趙鐠等，2017；楊洋等，2019；閻浩等，2020)。然而城市淺層地質(zhì)條件復雜，地下孤石等地質(zhì)災害體較為常見且危害巨大，為降低城市地下空間開發(fā)過程中的安全隱患，必須詳細查明地下孤石的分布情況(施有志等，2019)，目前孤石探測仍處于發(fā)展階段(吳成平和胡祥云，2007；黨如姣，2012；徐佩芬等，2012)。通常情況下孤石分布在特定的地層中(馮濤，2007；謝壯，2010；楊雄飛和袁鴻，2011；楊書江，2014)，在構造條件簡單的情況下可以借助鉆探查明一部分建構筑物基礎區(qū)域的地質(zhì)情況，但大多數(shù)情況下大量采用鉆孔勘察會大大提高成本和工期，而充分利用工程物探手段，往往能取得較好的效果(李萬倫等，2018)。

在工程物探領域中，地震層析成像技術有著不可替代的重要作用(王妙月，1997；趙永貴等，2000；孫宇等，2003；肖寬懷等，2003；郭慧麗和徐佩芬，2011；張凱等，2013)。這種方法借鑒了醫(yī)學上X射線斷面掃描診斷的原理，利用地震波在不同方向投射的波場信息，得出被測區(qū)域內(nèi)巖體波速分布規(guī)律，從而對地下介質(zhì)內(nèi)部結構進行精細成像(劉暢等，2020)。與常規(guī)工程物探方法相比，地震層析具有很高的分辨率，有助于全面細致地對巖體進行質(zhì)量評價，查明地質(zhì)異常體(譚顯江等，2013；舒森和胡金星，2016；秦晶晶等，2018；張雨飛，2019；何耀京，2019)。

目前，地震層析成像技術已在孤石探測問題上引起重視并得到了研究(李紅立等，2016)，但針對地震波在含孤石介質(zhì)中的傳播特征及探測效果的研究較少。本文以正演模擬為基礎，構建含孤石的地質(zhì)模型，分析地震波在含孤石介質(zhì)中的傳播特征，根據(jù)地震波的傳播特征利用層析成像技術對孤石進行速度成像，并通過廈門某軌道工程實例進行了驗證。

1 基本原理

地震層析成像通常使用兩個或多個鉆孔，經(jīng)過多點激發(fā)、多點接收，在被測區(qū)域內(nèi)形成致密的射線交叉網(wǎng)絡(成谷等，2002；胡明順等，2009；呂子強，2009；孫嘉妤，2019)。這些射線的傳播時間包含了被測區(qū)域內(nèi)介質(zhì)的波速信息。為了求解巖體的波速分布規(guī)律，首先需要將被測區(qū)域劃分成許多規(guī)則的矩形網(wǎng)格，網(wǎng)格化方法及觀測方式如圖1所示。

圖1 地震層析成像測區(qū)網(wǎng)格化方法

不同接收關系下地震波的到達時間用數(shù)學方程可表示為：

(1)

式(1)中，i為射線條數(shù)，i=1，2，…，n，n為射線總條數(shù)；j為劃分成規(guī)則矩形單元個數(shù)，j=1，2，…，m，m為矩形網(wǎng)格總個數(shù)；Ti為第i條射線接收的時間；vj(x,y)為坐標x、y的第j個矩形單元的速度值(待求)，Sj(x,y)為坐標x、y的第j個矩形單元的慢度值；lij為第i條射線通過第j個單元矩形的長度。

當矩形網(wǎng)格劃分得足夠小時，可以認為每個網(wǎng)格內(nèi)的介質(zhì)是均勻的，即波速單一。因此公式(1)可寫成離散化形式：

(2)

公式(2)實際上為線性方程組，可表示為：

(3)

其矩陣表達式為：

Tn×1=Ln×mSm×1

(4)

實際計算中，網(wǎng)格很小，穿過每個網(wǎng)格的射線條數(shù)很多，但每條射線穿過的總網(wǎng)格數(shù)都小于2 m，所以公式(4)中的L是一個大型的稀疏矩陣。為了求解這個稀疏矩陣，首先將被測區(qū)域中的v或s賦予一個初值，然后將所得到的投影值殘差沿射線方向均勻地反投影過去，并不斷對網(wǎng)格中的數(shù)值進行校正，一直到滿足精度條件為止。上述過程用數(shù)學公式可表示為：

(5)

公式(5)可計算出一系列向量S(1)，S(2)，…，當k足夠大時，這個向量收斂于所求的值。這個值即為被測區(qū)域內(nèi)波速分布，常用等值線圖進行表示。

2 數(shù)值模擬

地震波場數(shù)值模擬是通過計算機模擬地震波在給定模型中的傳播規(guī)律，具有方便、靈活成本低等特點。隨著計算機技術的發(fā)展，地震波場的數(shù)值模擬技術已經(jīng)成為地球物理勘探的重要手段(陳建，2015)。

2.1 單個孤石體的數(shù)值模擬

由于地震層析成像的研究對象為直達波走時，為了簡化模型提高計算效率，采用聲波方程模擬含孤石均勻地層中的地震波場特征。三維正演模型如圖2a所示，模型在X、Y、Z方向上的大小分別是50 m、50 m、50 m。其中XY方向為水平地面方向，Z方向為垂直方向。孤石為直徑4 m的球形體，球心位于X、Y方向上25 m，Z方向上15 m處。模型在X、Y和Z方向上空間網(wǎng)格大小均為0.5 m。模型的彈性參數(shù)如表1。震源采用頻率為120 Hz的零相位雷克子波爆炸震源，模型周圍添加PML吸收邊界，減小模型邊界對有效波場的影響。觀測系統(tǒng)布置如圖2b，采用單孔激發(fā)，單孔接收觀測系統(tǒng)。震源位置從(5，5，0)到(5，5，30)，炮間距為1 m。檢波器位置從(45，45，0)到(45，45，30)，道間距為1 m。

表1 模型介質(zhì)參數(shù)表Table 1 Parameter of model media

圖2 三維模型與觀測系統(tǒng)

對地震波場三維空間傳播特征進行分析。圖3為第1炮不同時刻波場快照。其中，左邊為三維顯示，右邊為對應的剖面圖。通過分析波場快照可以得出，在20 ms時刻，震源激發(fā)產(chǎn)生的直達波在均勻介質(zhì)中呈球形擴散；在40 ms時刻，直達縱波在球形擴散過程中遇到孤石異常，在接觸點發(fā)生繞射現(xiàn)象，繞射波相位與原波場相反，且能量減弱。

圖3 不同時刻三維地震波場快照與切片

圖4為前6炮地震記錄。其中，同相軸1為直達波；同相軸2為孤石產(chǎn)生的繞射波。二者速度均為800 m/s。通過對比信號的能量可知，直達波的能量強，孤石繞射波能量弱于直達波，但頻率略高于直達波，這與波場快照的分析吻合。前12道檢波器可接收到分離的直達波與繞射波，后18道接收到的是直達波與繞射波疊加后的信號。

圖4 地震記錄

建立好觀測系統(tǒng)后對地震信號進行初至時間拾取。初至時間是地震層析成像反演工作中最重要的基礎數(shù)據(jù)，其精度直接影響層析的精度，因此初至拾取工作需要反復多次的進行。圖5為第1、5、10、15、25炮每一道信號繪制對應的初至時間曲線。

圖5 不同炮地震記錄初至時間曲線

將探測區(qū)域劃分成規(guī)則正方形網(wǎng)格，網(wǎng)格邊長為0.5 m，對到時數(shù)據(jù)按公式(5)進行迭代成像，結果如圖6。成像結果中，黑色虛線表示數(shù)值模擬設置的孤石位置。通過層析成像結果可知，探測區(qū)域內(nèi)存在明顯的高速異常，異常范圍較大。異常位置在X方向23～27 m，Z方向-13～-16 m處，異常速度為2300 m/s以上，與建立的模型基本吻合，但在高速異常周圍存在較大的速度漸變。說明層析成像能很好地反映孤石的尺寸與位置，但成像精度與反演網(wǎng)格大小有關。

圖6 單個異常成像結果

2.2 兩個孤石體的數(shù)值模擬

為研究分析反演網(wǎng)格大小對地震層析成像的精度影響和孤石模型可以區(qū)分的縱向距離，設置垂直方向組合的孤石模型。

在單個孤石模型的基礎上添加另外一個孤石模型，其它參數(shù)保持不變。三維正演模型如圖7所示的。孤石1為直徑4 m的球形體，球心位于X、Y方向上25 m，Z方向上15 m處。孤石2為直徑4 m的球形體，球心位于X、Y方向上25 m，Z方向上21 m處。兩個孤石中心在垂直方向上間距為6 m。

圖7 三維模型與觀測系統(tǒng)

將探測區(qū)域劃分成規(guī)則正方形網(wǎng)格，網(wǎng)格邊長為0.2 m，對采集的數(shù)據(jù)按上述方法進行迭代成像，結果如圖8。成像結果中，黑色虛線表示數(shù)值模擬設置的孤石位置。通過層析成像結果可知，探測區(qū)域內(nèi)存在兩個明顯的高速異常，異常區(qū)域相對于單個孤石模型的反演結果明顯收斂。異常1在X方向23～27 m，Z方向-13～-17 m范圍內(nèi)，異常速度為2300 m/s以上；異常2在X方向23～27 m，Z方向-20～-23 m范圍內(nèi)，異常速度為1700 m/s左右。由于異常1和異常2距離近，波場相互干涉，導致反演結果與實際模型在速度特征與實際大小上存在一定誤差。

圖8 兩個異常成像結果

數(shù)值模擬結果表明，地震層析成像能很好地反映孤石異常的位置與速度特征，通過減小反演網(wǎng)格的大小可以明顯提升反演結果的精度；當存在相鄰孤石異常時，孤石異常的縱向距離大于孤石半徑時可以有效區(qū)分。

3 工程實例

本次選取實際工程實例來進一步說明地震層析成像的效果。廈門地質(zhì)條件特殊，某軌道交通施工區(qū)間存在較多的孤石，對工程建設造成了較大的困難。根據(jù)前期的巖土勘測資料，淺層的巖土主要為碎石填土和花崗巖。碎石素填土主要由黏性土、碎石角礫組成，也有部分塊石和殘積黏性土，成分不均，縱波速度均值為850 m/s?；◢弾r以中等風化花崗巖為主，土質(zhì)較堅硬，縱波速度均值為2400 m/s；速度的差異作為此次地球物理探查的物性基礎。震源采用電火花震源，16道孔中檢波器作為信號接收儀器。圖9為實測地震記錄與拾取的初至。

圖9 地震記錄初至拾取

將震源與檢波器所圍成的區(qū)域采用0.2 m×0.2 m的網(wǎng)格進行劃分。對提取到的時間初至，采用公式(5)中的算法進行迭代，當?shù)螖?shù)足夠多，反演結果基本不再變化時，即為地震層析成像得到的最終結果。圖10為聯(lián)合所有鉆孔不同深度的平面反演結果圖。對各異常區(qū)域解釋如下：

圖10 異常解釋

異常區(qū)域1在水平方向上4～16 m處，豎直方向上-8～-14 m處，異常區(qū)范圍較大。包含2～3個高速異常中心，高速異常中心呈橢球形，其地震波波速2100 m/s以上，解釋為多個孤石組成的孤石群；異常區(qū)域2在水平方向上24～32 m處，豎直方向上-8～-12 m處，異常區(qū)范圍較大。主要為1個高速異常中心，高速異常中心呈橢球型，其地震波波速2100m/s以上，解釋為單個孤石異常；異常區(qū)域3在水平方向上38～43 m處，豎直方向上-10～-14 m處，異常區(qū)范圍較大。主要為2個高速異常中心，高速異常中心呈橢球形，其地震波波速2100 m/s以上，解釋為兩個孤石異常。

經(jīng)現(xiàn)場鉆孔施工驗證，在劃分的異常區(qū)域內(nèi)均有孤石存在，且孤石的大小與規(guī)模與反演結果基本吻合，說明地震層析成像法對被測區(qū)域內(nèi)速度差異反映敏感，能有效用于城市地下工程建設的孤石探測，防范孤石地質(zhì)災害帶來的生命財產(chǎn)安全損失。

4 結論

(1)文章通過數(shù)值模擬方法得到了不同孤石異常模型的地震波響應與層析成像結果，說明地震層析成像方法能很好地反映孤石異常的位置與速度特征，通過減小反演網(wǎng)格的大小可以明顯提升反演結果的精度。

(2)由于實際孤石異常復雜，當存在相鄰孤石異常時，孤石異常的縱向距離大于孤石半徑時可以有效區(qū)分相鄰孤石。

(3)現(xiàn)場工程實踐有效查明了廈門某軌道交通施工區(qū)間內(nèi)孤石的分布與規(guī)模，說明地震層析成像是一種經(jīng)濟有效的地質(zhì)異?？辈旆椒?。