段成龍，閻長虹，許寶田，吳煥然，鄒明陽

南京大學(xué)地球科學(xué)與工程學(xué)院，南京，210046

內(nèi)容提要：本文介紹了在南京地鐵4號線某區(qū)間應(yīng)用跨孔地震CT技術(shù)進行現(xiàn)場數(shù)據(jù)采集的方法，并運用Matlab語言編制程序，通過最小二乘法對大型線性稀疏矩陣方程進行求解，由生成的三維網(wǎng)格數(shù)據(jù)繪制波速等值線圖，從而實現(xiàn)了孔間介質(zhì)波速層析成像的功能。根據(jù)成像圖中低速帶范圍及波形分析確定了溶洞位置和充填情況，經(jīng)checkerboard測試，并與鉆探資料和地質(zhì)雷達影像對比分析，結(jié)果表明跨孔地震CT技術(shù)用于地下巖溶探測時可靠性高，能夠?qū)崿F(xiàn)對孔間巖溶發(fā)育狀況的超前預(yù)報。

巖溶是威脅地下工程施工安全的一大隱患。在南京地鐵4號線工程勘探中發(fā)現(xiàn)，部分灰?guī)r分布地段有大規(guī)模巖溶發(fā)育，洞徑小則幾十厘米，大至4～5m，巖溶類型為覆蓋型(曾玉等，2011)，覆蓋層為第四系沉積物。為避免施工過程中發(fā)生坍塌、涌水等工程事故，必須采取有效的手段查明巖溶的具體位置和規(guī)模。傳統(tǒng)鉆探方法成本高，鉆孔間距和數(shù)量受到較大限制，因此單憑鉆探資料并不能對測區(qū)地質(zhì)異常體的分布情況做出整體判斷(李世民，2004)。若采用地面物探方法，如垂直地震剖面(VSP)法和高頻電磁波法，上覆土層可大量吸收彈性波和電磁波能量，探測深度達不到要求；直流電法因受到表面高阻層(水泥或瀝青路面)屏蔽傳導(dǎo)電流的影響，探測效果也不理想(陳易玖，1999)。相比較而言，跨孔地震CT技術(shù)在基巖中激發(fā)彈性縱波，探測深度不受影響，勘探準確性也得到了保證；此外，該技術(shù)可對孔間巖土體進行逐層掃描，并對目標探測區(qū)內(nèi)介質(zhì)的波速成像，使得孔間異常體能在圖中直觀地表現(xiàn)出來(李世民，2004)，因此地震CT技術(shù)用于探測孔間巖溶分布不失為一種較理想的方法。

對跨孔地震CT技術(shù)來說，成像分辨率是制約解釋精度的重要因素。它既受激發(fā)頻率、觀測系統(tǒng)等硬件性能的影響，也受反演軟件的影響(裴正林等，2002)。目前針對石油物探的反演軟件已得到廣泛應(yīng)用，以走時反演為例，其核心算法從分辨率極低的反投影技術(shù)(BPT)，過渡到以聯(lián)合迭代重建技術(shù)(SIRT)和正交分解(LSQR)為代表的迭代法(楊文采等，1993；楊文采，1997；羅炬等，2011)，到90年代初又提出了高分辨率級聯(lián)法(杜劍淵等，1991；楊文采等，1993)。不足的是，這些方法或針對大型反演而設(shè)計(如迭代法)，或需要較高的計算成本(如級聯(lián)法)，若盲目應(yīng)用于規(guī)模不大的工程淺層勘探很可能事倍功半。因此，本文基于廣義最小二乘法反演的思想(楊文采，1997)，在Matlab中編程反演計算介質(zhì)波速，并將反演結(jié)果在二維剖面內(nèi)成像。根據(jù)波速圖中相對低速區(qū)的范圍，對測區(qū)巖溶發(fā)育位置和大小作出判斷。

圖 1 南京地鐵4號線一期工程樺墅站—仙林東站區(qū)間探測孔相對位置及跨孔測線分布情況Fig. 1 Plan of the relative locations of boreholes and survey lines for the cross-hole seismic CT test on metro line 4, Nanjing City

1 現(xiàn)場跨孔地震CT試驗

1.1 場地工程地質(zhì)與水文地質(zhì)條件

南京地鐵4號線一期工程樺墅站—仙林東站區(qū)間(礦山及明挖段)在鉆探過程中發(fā)現(xiàn)規(guī)模較大的溶洞，溶洞為充填或半充填，充填物主要為粉質(zhì)粘土和碎石土。鉆孔深度范圍內(nèi)揭示的地層為上土下巖結(jié)構(gòu)：上覆土層以雜填土、粉質(zhì)粘土為主，厚度從0.5～20m不等；下伏基巖類型為中三疊統(tǒng)周沖組灰色中厚層灰?guī)r和泥質(zhì)灰?guī)r。研究區(qū)地下水類型隨地貌條件而異：在地勢相對較高的構(gòu)造剝蝕低山丘陵地貌單元內(nèi)，地下水類型為基巖裂隙水和巖溶水，二者存在水力聯(lián)系；在山體坡腳以下地段內(nèi)，地下水類型為孔隙潛水、基巖裂隙水和巖溶水。

表1 試驗儀器Table 1 Testing instruments

圖 2 南京地鐵4號線Q16G95—Q16Z73跨孔地震CT現(xiàn)場觀測系統(tǒng)(兩孔間距26m)Fig. 2 The field observation system of cross-hole ( in this case, boreholes are Q16G95 and Q16Z73 on metro line 4, Nanjing City) seismic CT (the distance between the two boreholes is 26m)

1.2 現(xiàn)場數(shù)據(jù)采集

測試區(qū)間內(nèi)共12個勘探孔，見圖1，分成11組進行跨孔地震CT試驗。根據(jù)勘察要求，目標是對10～35m深度處的溶洞進行探測，由于孔間距不超過30m，因此為取得盡可能完全的投影數(shù)據(jù)，勘探孔深度應(yīng)大于40m。

為具體說明跨孔地震CT的試驗步驟，這里以Q16G95-Q16Z73兩孔間的試驗為例展開討論，其余測線的處理辦法與之類似。

試驗開始前，首先應(yīng)收集鉆探資料，以便確定激發(fā)深度并檢驗成像結(jié)果的可靠性。鉆孔資料顯示，Q16G95孔基巖埋深7.4m，地下水位距地表8.8m，在9.1～12.2m段發(fā)現(xiàn)溶洞，26.7～29.8m段為破碎灰?guī)r；Q16Z73孔基巖埋深0.5m，地下水位距地表11.4m，灰?guī)r完整性較好。根據(jù)現(xiàn)場測線布置示意圖(圖2)，激發(fā)區(qū)間在Q16G95孔下8.5～20.5m。完成儀器連接后，檢查RAS-24數(shù)字地震儀(表1)與計算機間的通信是否順暢，系統(tǒng)觸發(fā)狀態(tài)是否正常，若不存在問題，即可將電火花震源(表1)置于Q16G95孔內(nèi)20.5m深度處激發(fā)縱波；為了增加數(shù)據(jù)區(qū)域覆蓋率，將12道三分量檢波器以1m的道間距排列成串(表1)，底部固定在Q16Z73孔24.5m深度處，以便同時接收同一震源點激發(fā)的縱波信號；最后通過RAS-24地震儀和計算機將這12道信號記錄下來。為保證探測分辨率，設(shè)定炮間距為1m，移動Q16G95孔中的震源至19.5m深度后再次激發(fā)地震縱波，并由地震儀和計算機記錄下縱波信號。依上述流程逐次完成激—接工作，直至地震波覆蓋目標探測區(qū)域為止。

圖 3 直射線追蹤與區(qū)域網(wǎng)格剖分示意圖Fig. 3 Schematic illustration of an example straight ray path between the source and receiver and the mesh generation of the study area

圖 4 反演計算流程圖Fig. 4 Flowchart showing the steps of the inversion calculation計算公式參考楊文采(1997)Thecalculation formulas refer to Yang Wencai (1997)

值得注意的是，試驗過程中為了使獲取的信號有較高的信噪比，激發(fā)探頭和檢波器必須處于地下水或泥漿環(huán)境中，以保持傳感器與井壁良好的耦合(楊文采，1993)。對耦合條件不佳的，可適當提高充電電壓來加大激發(fā)能量。應(yīng)當指出，加大震源能量不如改善耦合條件，這是CT采集的主要經(jīng)驗(楊文采等，1993)。此外，為保證地震儀能最大程度地記錄高頻信號，在記錄系統(tǒng)中設(shè)置前放增益值為12dB來使波形有較大的振幅動態(tài)范圍(章建宏等，2005)。

2 基于走時反演的波速層析成像

2.1 反演模型

如圖3所示，在S處發(fā)射地震波，在R處接收。目標是通過第i個地震波的初至走時ti來反演斷面內(nèi)速度分布C(r)或慢度S(r) = 1/C(r)(i= 1, 2 ,…,N，r為空間位置矢量)。設(shè)地震波的傳播路徑為Li，則

(1)

式中，dl為弧長微元。

測區(qū)地質(zhì)資料表明該區(qū)成像介質(zhì)均一性好，彈性波主要穿透灰?guī)r地層，由于不同深度的灰?guī)r彈性波速變化很小(陸基孟等，1991)，因此波的衍射和散射效應(yīng)可忽略不計。根據(jù)Fermat原理，路徑Li可由圖3中的直線唯一確定，這一方法稱為直射線追蹤法(Gu Hanming et al.，2006)。以矩形網(wǎng)格離散成像區(qū)域(圖3)，并假設(shè)網(wǎng)格波速cj為常量(下標j為網(wǎng)格號，j= 1,2,…,M)，則該網(wǎng)格的慢度sj也是常量，代入(1)式得

(2)

式中，aij是第i條射線沿路徑Lij經(jīng)過第j個網(wǎng)格時所截線段長度(圖3)。

考慮所有N條射線的情況，就得到了一個關(guān)于網(wǎng)格慢度的M元線性方程組，寫成矩陣方程的形式，即

AS=T

(3)

2.2 數(shù)據(jù)處理

依據(jù)上述思想，以南京地鐵4號線Q16G95—Q16Z73兩孔間的現(xiàn)場測試結(jié)果為例說明數(shù)據(jù)處理方法，分為以下3個步驟：

圖 5 實測信號經(jīng)40～200Hz帶通濾波后的波形記錄Fig. 5 Waveform record processed by 40～200Hz band-pass filtering to the original signal

(1) 濾波。以Q16G95 孔內(nèi)1處激發(fā)(圖2)所獲實測信號為例，進行40～200Hz的帶通濾波處理。濾波后的地震記錄如圖5所示，其中橫軸為接收通道編號，縱軸為時間(單位 s)，信號采樣間隔0.125ms，數(shù)字信號記錄格式為SEG-D 8038，模擬信號前放增益為12dB。

表 2 初至波走時拾取結(jié)果(單位ms)Table 2 First arrival picking results (unit ms)

圖 6 南京地鐵4號線Q16G95—Q16Z73斷面反演計算點Fig. 6 The inversion calculation points on the cross section of Q16G95-Q16Z73 boreholes on metro line 4, Nanjing City

(2) 拾取初至波走時。由于初至波能量很強，有明顯的波峰或波谷，因此正常情況下，初至?xí)r間取在初至波信號起跳位置。但對噪聲較大以致初至波難以辨識的情形，應(yīng)拾取第一個波峰(或波谷)對應(yīng)的時間(李世民，2004)。因初至波走時拾取往往會引入觀測誤差，從而使方程(3)的解偏離真實情況，所以要對走時數(shù)據(jù)進行調(diào)試和校正，以確保反演結(jié)果的合理性。具體方法：①在單炮記錄中，對初至走時無法確定的記錄道，利用已確定的初至波走時值，通過樣條插值重新確定，得到一連續(xù)的走時曲線；②將插值后的走時值映射到記錄道中，在此附近尋找振幅最大的點作為最終初至波走時(潘樹林等，2006)。按上述方法共拾取了156個初至波走時值，將它們合成初至波走時矢量T= (ti)156×1，如表2所示。

圖 7 南京地鐵4號線Q16G95—Q16Z73跨孔二維剖面縱波波速等值線圖(圖中所標數(shù)字為等值線對應(yīng)的波速值)Fig. 7 Contoured interval P-wave velocities calculated by seismic CT in the cross-hole plane between boreholes Q16G95 and Q16Z73 on metro line 4, Nanjing City (the number marked in the figure represents velocity of the contour line)

(3) 成像平面網(wǎng)格離散。為使溶洞在成像時有較高的分辨率，必須指定合適的網(wǎng)格尺寸。一般來說，由于炮間距和道間距均為1m，當成像區(qū)被網(wǎng)格離散后，必須保證一個網(wǎng)格至少有一條射線經(jīng)過(楊文采等，1993；張平松等，2012)。其次，成像物理分辨率的估計對網(wǎng)格尺寸的劃分也具有參考意義。以圖5為例，信號主頻取為174Hz，由于縱波在溶洞內(nèi)的速度不超過1km/s，為得到穩(wěn)定的分辨率估計，波速取上限1km/s，求得波長為5.7m。因此，沿地震波傳播方向溶洞的分辨率必大于1/4個波長，即0.7m(Fowler C M R，1990；Gu Hanming et al.，2006)。最終將網(wǎng)格尺寸取為1m×1m。

楊文采等(1993)提出，走時反演得到的圖像分辨率約為3～5個像元寬度，且像元寬度最好小于探測目標特征尺度的1/3(像元在這里等同于網(wǎng)格)。而實際需查明對工程有影響的溶洞平均特征尺度大于3m，因此上述網(wǎng)格劃分辦法是合理的。令各網(wǎng)格中心點為反演計算點，則成像平面內(nèi)312個計算點如圖6所示，圖中橫軸代表以Q16G95為原點，跨孔斷面內(nèi)任一點的水平距離，縱軸代表埋深(注：斷面內(nèi)任一點埋深已經(jīng)過地形校正，零點均設(shè)在地表處)。

表 3 由波速成像結(jié)果得到的溶洞方位及其充填情況Table 3 The position of karst caves and filled conditions by the results of velocity tomography

圖 8 南京地鐵4號線Q16G95—Q16Z73跨孔剖面checkerboard分辨率檢測結(jié)果Fig. 8 Result of the checkerboard resolution test for the cross-hole plane between boreholes Q16G95 and Q16Z73 on metro line 4, Nanjing City

2.3 反演結(jié)果和討論

確定了成像平面的尺寸、網(wǎng)格大小及初至波走時矢量后，調(diào)用事先編制好的解譯程序，得到圖7所示的二維剖面縱波波速等值線圖，其坐標軸含義同圖6，縱波波速在圖中通過不同顏色予以區(qū)分。由圖可知，波速分布于0.8～1.3km/s的區(qū)域為整個成像區(qū)的低速帶。應(yīng)當指出，成像平面經(jīng)網(wǎng)格離散后，每一網(wǎng)格所反映的是該區(qū)域的平均波速，當?shù)刭|(zhì)異常體(如巖溶發(fā)育區(qū))的尺寸小于網(wǎng)格尺寸時，其波速會受到周圍介質(zhì)的影響，因此圖7所示低速帶的邊界處波速值偏大。

進一步分析，可從圖7中獲得下述探測結(jié)果：①根據(jù)低速帶的邊界確定出溶洞位置和平均洞徑(見表3)；②低速帶波速范圍大于縱波在空氣中的速度0.43km/s(陸基孟，1991)，因此解釋為充填型溶洞。表3給出了包括Q16G95—Q16Z73跨孔在內(nèi)的5組測線上探測到的溶洞的各項具體信息。

表 4 波速成像結(jié)果與鉆探資料對比Table 4 The comparisonbetween the results of velocity tomography and drilling data

一般地球物理反演問題存在多解性，跨孔地震CT問題也不例外。為了評價反演結(jié)果的可靠性，擬綜合checkerboard測試、鉆探資料和地質(zhì)雷達探測結(jié)果來說明該問題。

圖 9 南京地鐵4號線Q16G95—Q16Z73跨孔剖面地質(zhì)雷達探測結(jié)果Fig. 9 Interpreted radar-gram showing distribution of karst caves for the cross-hole plane between boreholes Q16G95 and Q16Z73 on Metro 4, Nanjing City圖中橫軸為水平距離(單位m)；左側(cè)縱軸代表反射波的回聲時間(單位ns)；右側(cè)縱軸代表探測深度(單位m)。計算探測深度時假定電磁波傳播速度為0.1m/ns，并采用電磁波單程走時，即回聲時間的一半The horizontal axis in the figure represents horizontal distance (unit m). The left and right vertical axes represent echo time of reflected wave (unit ns) and detecting depth (unit m), respectively. Assume the velocity of electromagnetic wave is0.1m/ns and use one way travel time (i.e. half of the echo time) when calculating detecting depth

根據(jù)checkerboard(或稱檢查板、棋盤格)法檢測成像分辨率的原理(Humphreys E R et al.，1988；Zhao Dapeng et al.，1992)，在初始速度模型中加入絕對值為3%的正負相間的擾動量，對Q16G95—Q16Z73跨孔層析成像結(jié)果的分辨率進行分析。檢測結(jié)果如圖8所示，從圖中直觀地看出，成像結(jié)果中的主要特征都是可信的。

表 5 地質(zhì)雷達探測到的溶洞方位及充填情況Table 5 The position of karst caves and filled conditions byground penetrating radar (GPR)

將鉆孔周圍解譯出溶洞的4組數(shù)據(jù)與相應(yīng)的鉆探資料對比，對比結(jié)果如表4所示。從中不難看出，在容許誤差內(nèi)，解譯結(jié)果與鉆探資料有較高的一致性。

同時，還采用地質(zhì)雷達方法對CT探測結(jié)果進行了驗證。本次使用的地質(zhì)雷達是瑞典Mala Geoscience公司生產(chǎn)的第三代ProEx雷達主機及50MHz超強地面耦合天線，測點間距0.2m，雷達有效探測深度為22m(圖9)。根據(jù)反射波能量大小及同相軸特征，得到Q16G95—Q16Z73跨孔斷面內(nèi)溶洞分布如圖9所示。表5詳細列出了除該斷面外其他4組跨孔測線的雷達解譯結(jié)果，對比于由波速成像得到的解譯結(jié)果(表3)，二者基本一致。

綜上所述，利用介質(zhì)波速層析成像方法可以有效給出溶洞的大小規(guī)模和發(fā)育位置，能夠達到預(yù)期的分辨率要求，因而解譯結(jié)果可信度高。若對測區(qū)內(nèi)所有符合跨孔地震CT技術(shù)測試要求的場地進行探測，即可對測區(qū)巖溶發(fā)育情況作出準確、系統(tǒng)的評價。

3 結(jié)論

通過將跨孔地震CT技術(shù)應(yīng)用于南京地鐵4號線樺墅站～仙林東站區(qū)間巖溶探測工程，得到了以下認識：

(1)現(xiàn)場數(shù)據(jù)采集過程中，應(yīng)通過改善震源及檢波器與井壁的耦合條件、抑制振動干擾和提高波形振幅動態(tài)范圍的辦法來增強信噪比，以減小后期解譯誤差。

(2)編制Matlab程序反演網(wǎng)格波速，根據(jù)生成的三維網(wǎng)格數(shù)據(jù)得到跨孔斷面波速等值線圖。通過辨識低速帶范圍及其波速值，確定溶洞的大小規(guī)模和充填情況。

(3)波速層析成像結(jié)果與鉆探資料及地質(zhì)雷達探測結(jié)果具有較好的一致性，證明了該技術(shù)在巖溶探測方面有較高的可靠性和較好的應(yīng)用前景，能夠?qū)︺@孔間巖溶發(fā)育狀況進行超前預(yù)報。