低頻可控震源激發(fā)信號特征及高精度動態(tài)監(jiān)測應用

王力偉 李麗 楊微 王寶善 黃柳芳

摘要：通過實驗測試了國產KZ28LF低頻可控震源在線性升頻、偽隨機編碼和單頻3種激發(fā)模式下的信號特征及在地下介質高精度動態(tài)監(jiān)測中的應用效果。結果表明：①線性升頻和偽隨機編碼激發(fā)模式可穩(wěn)定激發(fā)1.5～12Hz的低頻信號，近場信號互相關系數大于0.95的分別占81%和63%，信號能量在9～12Hz較強、在1.5～9Hz較弱，低頻外傳能量較弱，而單頻激發(fā)時高次諧波能量較強，信號重復性較差。②相干法和反褶積法均能壓縮震源子波，恢復介質真實響應，相干法的結果重復性和信噪比最高，反褶積法保幅和保相效果最好。③疊加后探測距離可達20km，波速變化探測精度達10-。

關鍵詞：低頻可控震源;聚類分析;反褶積;加權相位疊加;波速變化;高精度動態(tài)監(jiān)測

中圖分類號：P315.3文獻標識碼：A文章編號：1000-0666（2021）04-0622-13

0引言

利用地震波高精度動態(tài)監(jiān)測地下介質物理性質及其變化是人們認識地球物理過程的重要手段。實現地下一定目標深度的高精度動態(tài)監(jiān)測，需要能穩(wěn)定激發(fā)、高度重復、寬頻帶、高信噪比和高時間服務精度信號的地震震源。天然震源如背景噪聲和重復地震已廣泛應用于與火山噴發(fā)（Bren—guier et al，2014）、斷裂帶破壞和愈合（Peng， -Zion，2006;Chaves et al，2020）、流體運移（Niu et al，2003）、大地震（Liu et al，2014;Li et al， 2017）和地殼介質受外力加卸載（Hillers et al，2015;Mao et，2019）等相關波速變化研究中。但背景噪聲易受噪聲源變化及其能量分布影響（Zhan et al，2013），重復地震會受到定位精度和時空分布限制（Liet al，2017）。與天然震源相比，人工震源激發(fā)位置、深度和激發(fā)時間都精確已知，且可與流動臺陣組成高質量觀測系統(tǒng)，是地下介質高精度動態(tài)監(jiān)測的首選震源（Chen et al，2017）。一些學者使用壓電陶瓷、電動落錘、氣槍和精密常時可控震源等不同類型的人工震源，以10-5～10-的精度觀測到與地震相關的短期波速變化（Niu et al，2008;張元生等，2017;Wang et al，2020），以及與大氣壓（Silver et al，2007）、溫度（Wang et，2020）、水文（Yang et al， 2018）、固體潮（mamura et al，2003）等相關的不同尺度的波速變化。

可控震源是另一種廣泛應用于地球物理勘探的人工主動震源（Sallas，2010），根據驅動方式的不同，可分為液壓式、電磁式以及精密可控震源3種，與精密常時可控震源相似，其激發(fā)信號精密可控，可激發(fā)出具有一定持續(xù)時間、瞬時能量密度低、高度重復的地震波（陶知非，1995）。其中，液壓式可控震源還具有掃描信號可精確設計，激發(fā)信號低頻能量強、頻帶寬、分辨力高和作業(yè)方便等優(yōu)點，已是全球陸上使用最廣泛的機械震源（Krohn et al，2010）?？煽卣鹪匆驯粐L試用于局部尺度的地下介質時變特征研究，如Clymer和McEvilly（1981）利用可控震源觀測到可能與地震相關的0.2～0.5ms的走時變化;Korneev等（2000）利用橫波可控震源通過連續(xù)監(jiān)測圣安德烈斯斷裂內部介質波速變化來研究地震成核過程，發(fā)現斷裂帶引起的波速變化達6%;McCartney 和Cox（2013）利用可控震源作為重復加載源，研究土層的液化響應;Pevzner等（2011）利用可控震源連續(xù)監(jiān)測CO2的驅油過程。

常規(guī)可控震源由于機械和液壓系統(tǒng)的限制，激發(fā)信號頻帶范圍通常為8～300Hz，很難激發(fā)低于8Hz的低頻信號（Sallas，2010;Wei，Phillips，2013）。低頻可控震源是指工作頻率低頻極限達6Hz以下的震源，因低頻信號傳播過程中不易耗散和衰減、穿透力強、有效傳播距離遠，在區(qū)域尺度探測中具有較大優(yōu)勢。中石油東方地球物理公司在KZ28可控震源車基礎上，生產的KZ28LF新型低頻震源成功實現了1.5Hz的低頻信號激發(fā)，是當時全球唯一針對3Hz以下低頻設計的可控震源（陶知非，徐小剛，2017）。與其它震源相比，KZ28LF具有更寬的激發(fā)頻帶、更低的信號畸變和更高的信號重復性，穩(wěn)定的低頻信號使區(qū)域尺度結構成像和高精度動態(tài)監(jiān)測成為可能。

為測試KZ28LF激發(fā)信號特征及其在高精度動態(tài)監(jiān)測中的應用效果，筆者于2014年11月在河北徐水開展了連續(xù)激發(fā)和接收實驗。本文分析了其激發(fā)信號特征;使用波形互相關和聚類分析方法評估了激發(fā)信號的重復性和不良激發(fā)班次的剔除方法;測試了互相關、移動窗口互相關、相干法和反褶積4種信號壓縮方法對線性升頻、偽隨機編碼和單頻3種激發(fā)模式的使用效果;測試了不同信號疊加方式下低頻信號的傳播距離;最后使用互相關時延檢測方法分析了其在地下介質高精度動態(tài)監(jiān)測中的應用效果。

1工作原理和數據

液壓式可控震源集電子、機械和液壓系統(tǒng)為一體，主要包括電控箱體和振動器（圖1）兩部分（Sallas，2010;Wei et al ，2010;陶知非，徐小剛，2017）。其中電控箱體產生線性或非線性參考信號，經電液伺服閥生成液壓油流信號，高壓液壓油交替進入振動器重錘內的上腔室和下腔室，并驅動重錘做上下往復振動，而活塞桿與上板、立柱和振動平板等剛性連接，因此，重錘振動信號便通過與大地緊密耦合的振動平板傳入大地。電控箱體通過實時監(jiān)控設置在振動器重錘和平板上的加速度傳感器等反饋信號，按照Sallas模型精確控制振動器的工作狀態(tài)（Sallas，2010），確保振動器的實際振動信號與輸入的參考信號相符。

液壓式可控震源的工作狀態(tài)一般使用震源特征信號來描述，主要包括輸入的理論參考信號，重錘加速度計記錄的重錘信號，平板加速度計記錄的平板信號以及根據重錘加速度和平板加速度矢量合成的地面力信號（Sallas，1984）。地面力信號代表被平板下的近地表彈性介質轉換為的彈性地震波。因此，數據處理中一般使用地面力信號作為參考源信號進行子波壓縮處理（Brittle etal，2001），也可在震源車附近架設臺站，用記錄的近場信號作為參考信號進行子波壓縮處理（Poletto et al ，2016）。

為測試低頻可控震源激發(fā)信號特征及其在區(qū)域地下結構和斷裂帶應力變化的監(jiān)測能力，2014年11月，筆者聯合中國石油集團東方地球物理公司，在河北省保定市徐水區(qū)北奧特車試車場內開展低頻可控震源連續(xù)激發(fā)實驗（圖2）。實驗采用1臺載荷為26t的KZ28LF低頻可控震源，該震源最大靜載壓重達281kN、液壓峰值力達276kN、重錘質量為4445kg、平板質量為1724kg、激發(fā)頻帶1.5～86Hz、全流量最低頻率低于3Hz（陶知非等，2010）。實驗采取獨立掃描作業(yè)方式（表1），使用3種震源掃描信號進行連續(xù)激發(fā)（Cunningham ，1979;Dean，2014）。為增加振動平板與地表土層的耦合效果，實驗采取空震的方式壓實淺表土層。接收系統(tǒng)采用流動地震臺，跨安新斷裂和徐水斷裂（楊承先，1983），分別布設在偏移距2、5、10和20km位置（圖2），其中D20km臺站為基巖臺，其它為土層臺。為記錄震源車近場震源時間函數，在距離震源車10m處布設了參考地震臺（D0km）。觀測地震計為GURALP -40T（10m、20km距離處）和GURALP -6TD（2、5和10km距離處），儀器頻帶范圍分別為2s～100Hz，20s～100Hz，數據采集器為rEFTEK -130B，采樣率為200Hz，GPS采取連續(xù)授時模式。

2研究方法

液壓式可控震源得益于掃描信號精確可控、機械系統(tǒng)和耦合條件穩(wěn)定不變的優(yōu)點，能持續(xù)穩(wěn)定的激發(fā)出具有一定持續(xù)時間的、帶限的、瞬時能量密度低的高度重復信號。但因機械、電氣和平板耦合等多種原因，外傳信號會出現出力不均、畸變大和基波能量弱等現象（Nagarajappa，wilkin-son，2010）。為分析其激發(fā)信號特征，本文采用短時傅里葉變換時頻分析方法，分別分析了理論參考信號、實際出力信號和近場參考信號在線性升頻、偽隨機編碼和單頻3種信號掃描模式下的信號特征，然后利用D0km臺記錄的垂直分量波形，使用波形互相關和聚類分析方法研究了激發(fā)信號的重復性和不良班次的剔除方法。

2.1激發(fā)信號特征及數據篩選方法

聚類分析是根據一個類中對象間的相似程度，將類中的對象分成更多子類的過程，其目的是使類內部對象間的距離最小，不同類之間的距離最大（Sneath et al ，1973），目前已廣泛應用于相似地震識別（王偉濤，王寶善，2012）、震相拾?。≧owe et al ，2002）、精定位（Son et al ，2018）和地震成核過程研究（Zaliapin，Ben-Zion，2013）。

聚類分析方法使用基于樹狀結構的凝聚式層次聚類方法（Rowe et al ，2002）。首先，使用移動窗口互相關方法，計算了任意兩個事件之間的歸一化最大互相關系數，構建了相似性矩陣c，同時計算了每一次事件與將全部事件線性疊加構建的參考事件之間的最大互相關系數，互相關計算窗口為0～20s，移動步長為0.005s，事件之間振幅的差異被忽略，因為波形相位的一致性對于疊加和互相關時延估計更重要（Wang et al ，2008）。其次，計算了距離矩陣dn=1-c，這個矩陣包含所有即將聚類事件之間的相似距離，并再次計算事件之間、類之間、事件和類之間新的距離，反復迭代，直至所有事件均歸為一類為止。兩個事件之間距離由距離矩陣中對應元素確定，但是兩類之間的距離使用3種定義模式（Rowe et al ，2002），最嚴格的是最大距離法，即類間距離是兩個類之間所有事件的最大距離：

最寬松的是最短距離法，即類間距離使用兩個類間所有事件的最小距離：

式中：C、C為不同的類;x1和x1為類中不同對象;d（x，x）為對象間的距離;D為類間距離。當所有事件歸為一類后，通過設定相似度距離閾值a，便可得到此標準下的類別數和分屬各類的事件（Lance ，williams，1967;Rowe et al ，2002）。使用聚類分析方法避免了互相關計算時參考波形的選取，能更合理地得到震源車激發(fā)信號的穩(wěn)定性分布特征。

聚類分析前以近場信號的起跳時刻為參考零點，使用移動窗口互相關方法提取了各次事件的精確激發(fā)時刻，相關系數計算窗長20s，提取閾值0.955，濾波范圍單頻激發(fā)信號為1～2Hz，其它為1～13Hz，并裁剪了激發(fā)零時刻前5s、后25s總長度為30s的事件波形用于后續(xù)分析。通過波形互相關方法，共提取出線性升頻激發(fā)事件431個、偽隨機編碼激發(fā)事件121個和單頻激發(fā)事件400個，經與震源車激發(fā)班報進行對比，發(fā)現通過波形互相關方法不但可以100%的提取出激發(fā)班報給出的事件目錄，還能剔除一些發(fā)生嚴重畸變的激發(fā)班次。

2.2震源子波壓縮方法

可控震源激發(fā)信號為長時間的連續(xù)調頻信號，因此，為恢復介質的真實響應，必須對接收到的信號進行子波壓縮處理（Wapenaar，Fokke-ma，2006;楊微等，2013）。本文選用D0km臺記錄的近場信號作為子波壓縮的參考信號，分別測試了互相關、短時互相關、相干和反褶積4種方法在線性升頻、偽隨機編碼和單頻3種激發(fā)模式中的應用效果，4種信號壓縮方法的計算基本原理為：

互相關法壓縮子波是基于可控震源掃描信號的自相關函數可近似為零相位Klauder子波的假設（Klauder，1960），利用互相關法將長時間記錄的掃描信號壓縮為一個類似脈沖信號：

式中：U（r，t）是儀器記錄項;S（s，t）是近場震源項;G（r，s，t）是介質經驗格林函數，其包含直達體波、面波及多次波等介質的全部響應信息。互相關計算相當于經過掃描信號頻譜的2次方濾波，互相關結果可以壓制相關噪聲，但其結果與掃描信號的形狀有關，會改變各頻率的能量（Brittle et al， 2001; Wapenaar et al， 2010）

短時互相關法與互相關法計算過程類似，但每次僅計算第i個窗口w（n）內，儀器記錄信號U（r，t）和近場震源記錄信號S1（s，t）對應部分的互相關C（r，s，t），并以一定的步長滑動互相關計算窗口，遍歷整個記錄道，然后將每個窗口內的互相關函數疊加得到平均互相關結果（Bensen et al ，2007;Tal et al ，2011）：

相干法是根據地下介質信號傳播的褶積模型和互相關原理，將儀器記錄信號U（r，t）和近場震源記錄信號S（s，t）在頻率域表示為：

式中：分子項與頻率域互相關相同;分母是臺站記錄信號和震源車激發(fā)近場信號的振幅譜。因此，通過頻域振幅譜相除，相干法去除了信號的振幅信息，僅保留了相位信息，能顯著壓制噪聲和道間振幅不一致的影響（Cao，2016）。

反褶積法同樣基于地下信號傳播的褶積模型，在頻域將震源車至臺站間的介質格林函數G（r，s，ω）表示為：

為避免除法數值運算的不穩(wěn)定引起的信號畸變，采用頻率域水準反褶積方法，引入水準比例因子c（Mehta et al ，2007）：

式中：max為求最大值。通過水準比例因子c將低于該閾值的振幅提升到閾值水平，以保證數值運算的穩(wěn)定性。最后再通過傅里葉反變換到時間域，得到震源車到接收臺站之間的格林函數G（r，s，t）與相干法類似，反褶積法也能顯著消除震源子波不一致和噪聲的影響，反褶積結果只與掃描信號的頻帶有關，與振幅譜的形狀無關，且保幅和拓頻效果較好，廣泛應用于地震資料的高保真數據處理（Brittle et al ，2001）。

2.3信號增強方法

低頻可控震源單次激發(fā)信號能量較弱，為增強有效信號、壓制噪聲，使用線性疊加法和加權相位疊加法來提高信噪比，以測試其最遠探測距離。研究表明，線性疊加有用信號隨疊加次數N成N倍增長，而隨機噪聲信號隨疊加次數N成√N倍增長，因此線性疊加法可使結果信噪比增加N倍（蔣生淼等，2017）。直接線性疊加法給記錄的每個采樣點賦予了相同的疊加權重，加權相位疊加法可根據信號和噪聲的相關程度大小，分別賦予不同的權重，可更有效地壓制噪聲，增強遠距離臺站微弱信號的檢測效果（Ken-nett，2000）。

2.4走時變化測量方法

為測試低頻可控震源在地下介質高精度動態(tài)監(jiān)測中的應用效果，使用互相關時延檢測方法分別測試了3種激發(fā)模式下走時變化的檢測精度。該方法采用滑動窗口互相關技術（Wang et al ，2008），通過測量不同窗口內比較波形和參考波形的互相關系數，最大互相關系數對應的走時延遲即為走時變化d，進而得到一定窗口范圍內的平均dt及其誤差估計。采用互相關方法進行波形時延檢測的誤差理論下限可由Cramer -Rao公式估算（Niu et al， 2008; Wang et al， 2008）

式中：f為震源的中心頻率;B為頻寬比;SNR為信噪比;p為波形的互相關系數;T為計算互相關的窗長。從式（8）可見，提高激發(fā)信號的主頻、頻帶寬度、信噪比和互相關系數均可大大降低走時測量的誤差。

因反褶積法可以有效抑制源致諧波干擾、拓展頻帶寬度，恢復的振幅和相位響應更接近介質真實的格林函數，且可消除由于機電系統(tǒng)、平板耦合和近地表介質物理參數變化等震源處的變化而造成的震相到時的微弱變化（Ikuta et al ，2002;Lebedev ，Beresnev，2004），因此，測量信號使用信噪比較高的D2km臺與震源車附近參考臺的垂直分量反褶積處理結果。為進一步降低震源附近近場信號變化對走時變化測量的影響，對使用聚類分析提取的優(yōu)良激發(fā)班次：①采用波形互相關技術將參考臺的波形對齊。②將裁剪出的D2km臺的事件波形和參考臺的事件波形進行反褶積，線性升頻和偽隨機編碼濾波范圍為1～13Hz，單頻為1～2Hz。③通過疊加構建走時計算的參考波形和比較波形，參考波形是將全部波形線性進行疊加得到，比較波形是將相鄰震次進行線性疊加得到，且僅使用比較波形和參考波形的互相關系數大于0.91的事件進行走時變化測量。④選擇信噪比和互相關系數均較高的4.1～4.5s窗口范圍計算3種激發(fā)模式下的平均dt、平均互相關系數和標準差。線性升頻和偽隨機編碼模式dt計算窗長為0.4s，單頻模式dt計算窗長為0.9s，窗口移動步長0.005s，為獲取更高精度的走時變化，對互相關計算結果進行了5倍采樣率的余弦函數內插（Wang et al ，2008）。

3研究結果

3.1激發(fā)信號特征

線性開頻、偽隨機編碼和單頻模式下激發(fā)信號頻譜特征如圖3所示，從圖3a可見，線性升頻激發(fā)模式下，實際出力信號和近場信號激發(fā)頻率隨時間變化線性增大，信號主頻在1.5～12Hz，除能量較強的基波外，還存在2階、3階等高階諧波，諧波能量較弱。近場信號能量集中在9～12Hz，在1.5～9Hz低頻段，近場信號能量不如出力信號強。雖然震源車在前0～15s加在平板上的力較大，且75%振動時間集中在1.5～9Hz，而在9～12Hz僅振動5s，占25%，低頻信號外傳能力仍需改進。從圖3b可見，偽隨機編碼激發(fā)模式與線性升頻激發(fā)模式類似，實際出力信號主頻集中在1.5～12Hz，在4Hz和11Hz存在2個能量主峰，近場信號能量集中在9～12Hz，能量峰值隨時間變化特征與出力信號符合較好，無明顯諧波干擾，抗干擾性強，編碼效果較好，而1.5～9Hz低頻信號外傳能量相對較弱。從圖3c可見，單頻激發(fā)模式下，理論掃描信號能量都積累到1.5Hz，頻率域呈現單一峰值，主頻集中在1.5Hz，但由于震源是機械震源，除1.5Hz基頻外，出力和近場信號還存在能量較強的高階諧波，諧波頻率是基波頻率的整數倍，高能量的諧波耗散了大部分激發(fā)能量，降低了基波出力的信噪比。

綜上，從實際出力信號看，偽隨機編碼激發(fā)模式抗干擾能力最強，線性升頻模式次之，單頻模式最差，其諧波干擾最嚴重。這是由機械震源的Sallas力學模型決定的，如液壓系統(tǒng)的非線性、機械系統(tǒng)的窄帶限制、輸出信號振幅和相位畸變、平板的低剛度、近地表耦合條件的非線性、激發(fā)參數等，均可導致可控震源在振動出力時包含許多諧波畸變，尤其是在低于20Hz的低頻段（Weietal，2010，2013）。從近場記錄信號看，因基波和諧波能量一起輻射進入大地，線性升頻和偽隨機編碼模式下近場信號能量均集中在9～12Hz，外傳效果較好。1.5～9Hz低頻信號以及單頻激發(fā)模式1.5Hz的低頻信號受能量較大的源致諧波干擾影響較大，外傳能量相對較弱。

3.2重復性分析及數據篩選

圖4a給出了3種激發(fā)模式下，每次激發(fā)的近場信號相對通過疊加構建的參考信號的互相關系數分布。從圖4a-1可見，線性升頻激發(fā)時重復性最好，相關系數在0.95以上的占81%，0.90以上的占88%，信號畸變發(fā)生在第48到100次激發(fā)之間，畸變發(fā)生時段較集中，且因實際出力信號畸變造成，調整激發(fā)參數后，激發(fā)波形非常穩(wěn)定。偽隨機編碼激發(fā)重復性次之（圖4a-2），相關系數在0.95以上的占63%，0.90以上的占68%，信號畸變發(fā)生在第80次激發(fā)之后。單頻激發(fā)重復性最差（圖4a-3），17%的相關系數在0.85以上、42%在0.8以上、57.7%在0.75以上，相關系數在整個激發(fā)時段較離散。

相似度距離閾值α的不同選取會影響類和類內事件的數目，為評價低頻震源車激發(fā)信號的重復性，筆者從0.8至0.995，每隔0.005計算了不同相似度閾值α下的聚類結果，并統(tǒng)計了類內對象數目最多的前4個類以及剩余事件數目占全部事件的百分比，根據類的數目和對象數的分布情況選擇最優(yōu)的閾值a（圖4b），其中線性升頻和偽隨機編碼的類間距離采用嚴格的最大距離法，單頻由于相似性較差，采用寬松的最短距離法。聚類結果顯示，線性升頻激發(fā)模式下a取0.925，87.7%的事件能聚為一類（圖4b-1），其它類的數目雖較多，但每一類的事件數卻較少，可見采用線性升頻激發(fā)模式時，激發(fā)信號高度重復。偽隨機編碼激發(fā)模式下a取0.945，66.1%的事件能聚為一類（圖4b-2），且大多是前84次激發(fā)的信號，后37次激發(fā)信號不穩(wěn)定，可能因振動系統(tǒng)或平板耦合問題造成信號畸變較大。單頻激發(fā)模式下a取0.885，51.3%的事件能聚為一類（圖4b-3），但其它事件也大多自成一類，說明單頻模式雖能穩(wěn)定激發(fā)1.5Hz的低頻信號，但重復性較差。圖4c顯示了該閾值下，數量最多的前5個類的波形特征，其中，第1類是高質量的激發(fā)班次，其它多為背景噪聲、無出力、畸變大等不良激發(fā)班次。

可控震源的失效模式可分為瞬間崩潰式失效和漸變失效（陶知非，徐小剛，2017），可控震源激發(fā)時通過多種實時質控方式保證激發(fā)信號質量（王秋成等，2019）。從聚類分析結果看，KZ28LF低頻可控震源通過實時質控及時調整激發(fā)參數，可以穩(wěn)定激發(fā)高度重復的線性升頻和偽隨機編碼信號，使利用疊加技術增加信噪比進行區(qū)域尺度探測和介質時變特征研究成為可能。本文使用聚類篩選的378個線性升頻激發(fā)事件、80個偽隨機編碼激發(fā)事件和205個單頻激發(fā)事件用于后續(xù)研究。

3.3震源子波壓縮方法

圖5給出采用不同信號壓縮方式D2km臺垂直分量波形的處理結果，數據處理包括去均值、去傾斜、帶通濾波（線性升頻和偽隨機編碼1～13Hz帶通濾波，單頻1～2Hz帶通濾波），然后將D2km臺記錄的事件波形和D0km臺相應的事件波形進行子波壓縮處理，并將50次計算結果進行線性疊加。從圖5a可見，對于線性升頻激發(fā)模式，4種信號壓縮方法均能正確壓縮震源子波，得到介質響應信息，互相關和短時互相關法結果相似，恢復信號的體波振幅較大，面波振幅較小，相干法和反褶積法計算結果相似，但相干法得到的體波振幅較大，面波振幅相當，反褶積法計算結果的保幅和保相效果最好，體波和不同階次的面波信號幅度恢復均較好（圖5a-1）。從圖5b可見，對于偽隨機編碼激發(fā)模式，4種方法均能正確壓縮震源子波，互相關和短時互相關法結果相似，能正?；謴统鲶w波信號，但恢復的面波信噪比較低，而相干法和反褶積法均能較好地恢復信號，但相干法得到的體波振幅較大，面波振幅相當，反褶積法計算結果的保幅和保相效果最好。從圖5c可見，對于單頻激發(fā)模式，相干法和反褶積法均能正確壓縮震源子波，但相干法得到的信噪比更高，互相關和短時互相關法壓縮結果較差?？傊?種方法中，相干法和反褶積方法使用范圍更廣，通過計算恢復后信號的波形互相關系數和信噪比，發(fā)現相干法得到結果的重復性和信噪比最高，反褶積法保幅和保相效果最好。

3.4信號傳播距離

在線性升頻激發(fā)模式下，低頻可控震源經反褶積處理后的信號經378次線性疊加和加權相位疊加后的記錄剖面圖如圖6所示。從圖中可見，采用加權相位疊加后，信噪比相對直接線性疊加得到了較大提高，經過378次疊加后，20km范圍內的臺站能記錄到有效信號。顧廟元等（2016）在地學長江計劃安徽實驗中，使用EV56低頻可控震源和更長的儀器觀測測線，發(fā)現通過時頻率相位加權疊加，200多次疊加后可在50km范圍內提取到有效信號。

3.5走時變化測量

本文使用D2km臺垂直分量反褶積恢復的4.0～5.1s窗口內的面波信號進行走時變化測量（圖6紅色方框），從垂向-徑向的質點運動圖可見（圖6c），其為基階瑞雷面波（Godfrey et al ，2017）。為提高比較波形的信噪比和相關系數，對相鄰震次的信號進行了線性疊加，圖7a-1、b-1和c-1展示了4.0～5.1s窗口內不同疊加次數的比較波形和參考波形的互相關系數分布圖。由圖可見，當線性升頻模式疊加15次，偽隨機編碼和單頻模式疊加20次可顯著降低相關系數小于0.9的事件比例。線性升頻模式平均激發(fā)間隔為44.7s，15次疊加相當于11.2min的平均;偽隨機編碼模式平均激發(fā)間隔29.4s，20次疊加相當于9.8min的平均;單頻激發(fā)模式平均激發(fā)間隔82.3s，20次疊加相當于27.4min的平均。最后使用信噪比和互相關系數均較高的4.1～4.5s窗口范圍計算了3種激發(fā)模式下的平均dt、平均互相關系數和標準差。

由圖7a-2、b-2和c-2可見，3種激發(fā)模式d集中在±10ms之間，平均走時約為4.3s，相對波速變化和相對走時變化關系為：8v/v=-8t/t，因此，相對走時變化和相對波速變化范圍約為0.23%。從dt互相關系數來看，掃頻和編碼模式互相關系數相對較高，大多在0.95以上，單頻互相關系數較低，大多在0.84以上。

從dt測量誤差來看（圖7a-3、b-3和c-3），線性升頻和單頻模式誤差為0～9ms，偽隨機編碼模式的誤差為0～7ms，線性升頻模式誤差小于其均值2.1ms的占比為72.3%，小于1ms的占比為42.9%，偽隨機編碼模式誤差小于其均值2.9ms的占比為63.1%，小于1ms的占比為17.5%，單頻模式誤差小于其均值4.5ms的占比為56.7%，小于1ms的占比為1.6%。因測量窗口中心走時是4.3s，因此1ms、2ms、3ms和4.5ms的走時測量誤差分別相當于2.3×10-4、4.6×10-4、6.9×10-4和1.0×10-3的相對走時變化和相對波速變化測量誤差。從測量誤差來看，線性升頻模式因波形相似性高，測量精度最高，偽隨機編碼模式其次，單頻模式因波形相似性低且頻帶窄，測量誤差最大。

本文線性升頻和偽隨機編碼模式f0≈5Hz，B≈8，SNR≈50，p≈0.97，T≈0.4s，代入式（8）計算得到時延估計誤差的理論值約為5.6×10-4s，而觀測使用的RefTek -130B數采鐘差小于10μs（Wang et al ，2020），因測量窗口中心走時是4.3s，因此相對走時變化的理論誤差為1.3×104，可見實測誤差和理論誤差相近。從圖7a-2、b-2和c-2可見，不同激發(fā)模式觀測到的dt隨時間變化存在約25min的周期，單頻激發(fā)模式存在約1h的周期，將線性升頻模式相鄰班次的疊加次數增加至35次，約相當于26min的平均，dt計算結果也存在約1h的周期成分，不同激發(fā)模式觀測到的類似現象也印證了觀測到的信號的真實性。

4討論

本文實驗表明，低頻可控震源能以線性升頻和偽隨機編碼模式穩(wěn)定激發(fā)頻率在1.5～12Hz的低頻信號，激發(fā)數據經過嚴格篩選后，波形互相關系數分別優(yōu)于0.925和0.945，疊加后探測距離可達20km，能以（2～7）×10-4的測量精度探測到2.3×10-3的相對波速變化。低頻可控震源經過精細處理后，能達到其它人工震源相似的重復性和波速變化探測精度，不僅如此，低頻可控震源掃描信號可精確設計，作業(yè)方便，對場地要求低，且能以低能量密度穩(wěn)定激發(fā)低頻信號，特別適合區(qū)域尺度的斷裂帶精細結構及斷裂帶內部大震成核過程的應力變化研究（Korneev et al ，2000）。同時，由于低頻可控震源瞬時能量密度低，場地破壞小，也可適用于滑坡觀測，謝凡等（2020）通過背景噪聲互相關干涉測量發(fā)現，滑坡失穩(wěn)前存在大于1%波速相對變化。低頻可控震源發(fā)展迅速，是當今勘探裝備的研究熱點，國內低頻可控震源在KZ28LF之后又有了發(fā)展，LFV3激發(fā)頻率范圍為1.5～120Hz，EV56激發(fā)頻率范圍為1.5～190Hz，輸出力平均畸變降低了10%，不用特殊設計的掃描信號就可以穩(wěn)定激發(fā)1.5Hz的低頻信號（陶知非等，2018）。因此，隨著技術的進步，低頻可控震源在高精度動態(tài)監(jiān)測方面的應用將越來越廣泛。

5結論

本文使用國產低頻可控震源KZ28LF，通過野外實驗測試了其激發(fā)信號特征及在地下介質高精度動態(tài)監(jiān)測中應用效果，得到的結論和認識如下：

（1）低頻可控震源能持續(xù)穩(wěn)定地激發(fā)頻率在1.5～12Hz的線性升頻信號和偽隨機編碼信號;出力信號和近場信號除有能量較強的基波外，還存在高階諧波，偽隨機編碼諧波能量較弱。近場信號9～12Hz能量較強，1.5～9Hz相對較弱，低頻信號外傳能量較弱。

（2）線性升頻激發(fā)模式的波形重復性最好，81%的互相關系數大于0.95，偽隨機編碼激發(fā)模式其次，63%的互相關系數大于0.95，單頻激發(fā)模式波形重復性最差，17%大于0.85。線性升頻和偽隨機編碼激發(fā)模式波形畸變時段集中，可通過實時質控提高激發(fā)波形的重復性。

（3）相干法和反褶積方法可用來壓縮線性升頻、偽隨機編碼和單頻3種類型的激發(fā)信號，相干法得到結果的重復性和信噪比最高，反褶積法保幅和保相效果最好。

（4）經過378次疊加后，有效信號傳播距離達20km。

（5）剔除不良激發(fā)班次后，線性升頻和偽隨機編碼2種模式能以（2～7）×10-4的測量精度探測到2.3×10-3的相對波速變化。

在實驗過程中，中國石油集團東方地球物理公司提供了低頻可控震源，中國地震局地震科學探測臺陣技術中心提供了流動觀測設備，東方地球物理公司陶知非總工和祝楊高級工程師指導并參加了野外實驗。審稿專家為本文提出了寶貴修改意見，在此一并表示衷心感謝！

參考文獻：

顧廟元，姚佳琪，張偉，等。2016。地學長江計劃安徽實驗中低頻可控震源地震波信號提取方法評估[J].中國地震，32（2）：356-378.

蔣生淼，王寶善，張云鵬，等。2017。噪聲對氣槍信號疊加效果的影響及自動數據篩選方法[J].地震研究，40（4）：534-542。

陶知非，蘇振華，趙永林，等。2010?？煽卣鹪吹皖l信號激發(fā)技術的最新進展[J].物探裝備，20（1）：1-5。

陶知非，徐小剛，駱飛，等。2018。LFV3低頻可控震源技術與應用效果[J].物探裝備，28（1）：1-3。

陶知非，徐小剛。2017。再論可控震源輸出信號畸變的問題[J].物探裝備，27（6）：351-355。

陶知非。1995?？煽卣鹪磁c高分辨率地震勘探[J].石油物探裝備，5（4）：1-9.

王秋成，王梅生，孫哲，等。2019?？煽卣鹪磩討B(tài)掃描高效采集實時質控要點分析[J].物探裝備，29（3）：141-146。

王偉濤，王寶善。2012?；诰垲惙治龅亩喑叨认嗨频卣鹂焖僮R別方法及其在汶川地震東北端余震序列分析中的應用[J].地球物理學報，55（6）：1952-1962。

謝凡，夏開文，黃會寶，等。2020。基于多重散射波波速變化的滑坡實時監(jiān)測方法與應用研究[J].巖石力學與工程學報，39（11）2274-2282.

楊承先。1983。河北徐水一文安地區(qū)第三紀伸展斷裂的地質特征及其活動性分析[J].地震，（5）：34-36，45。

楊微，王寶善，葛洪魁，等。2013。精密控制機械震源特征及信號檢測方法[J]中國石油大學學報（自然科學版），37（1）：50-69。

張元生，王寶善，陳颙，等。2017。兩次強震發(fā)生前后主動源觀測走時數據的變化[J].地球物理學報，60（10）：3815-3822。

Bensen G D.Ritzwoller MH.BarminMp. et a.2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dis-persion measurements [ J].Geophysieal Journal International， 169（3）：1239-1260.

Brenguier F， Campillo M， Takeda T， et al.2014. Mapping pressurized vol-canic fluids from induced crustal seismic velocity drops[ J]. Sci-ence，345（6192）：80-82.

Brittle KF， Lines L R， Dey A K.2001.Vibroseis deconvolution： a com-parison of cross-correlation and frequeney-domain sweep deconvolution]. Geophysical Prospecting， 49（6）： 675-686.

CaoHT.2016. Comparisons of seismic interferometry by cross correla-tion， deconvolution and cross coherence [ D]. Houghton： Michigan Technological University

Chaves EJ， Schwartz SY， Abercrombie R E.2020. Repeating earthquakes record fault weakening and healing in areas of megathrust postseismicslipLJ]. Science Advances， 6（32）： 1-8.

Chen Y， Wang B S， Yao H J.2017. Seismic airgun exploration of conti-nental crust structures[ J]. Science China Earth Seiences， 60（10）1739-1751.

Clymer R W， McEvilly TV.1981. Travel-time monitoring with vibroseis LJ. Bulletin of the Seismological Society of America， 71（6）： 1903-1927.

Cunningham A B.1979. Some altemate vibrator signalsJ]. Geophysics 44：1901-1921

Dean T.2014. The use of pseudorandom sweeps for vibroseissurveys[J Geophysical Prospecting， 62（1）： 50-74.

Godfrey HJ， Fry B， Savage M K.2017. Shear-wave veloeity structure of the Tongariro Voleanie Centre， New Zealand Fast Rayleigh and slow Love waves indicate strong shallow anisotropylJ]. Journal of volean-ology and Geothermal Research 336：33-50

Hillers G， Ben-Zion Y， Campillo M， et al.2015. Seasonal variations of seismie velocities in the San Jacinto fault area observed with ambient seismic noise [ J. Geophysical Journal International， 202（2）： 920-932.

Ikuta R. Yamaoka K， Miyakawa K， et al.2002. Continuous monitoring of propagation veloeity of seismie wave using ACROSS[ J]. Geophysical Research Letters 29（13）： 1627.

Kennett BLN.2000. RESEARCH NOTE： Stacking three-component seis-mograms[ J. Geophysical Journal International， 141（1）： 263-269

Klauder JR.1960. The design of radar signals having both high range res-olution and high velocity resolution[ J]. Bell System Technieal Jour-nal，39（4）：809-820.

Komeev VA， MeEvilly TV， Karageorgi E D.2000. Seismological studies at Parkfield VIll： Modeling the observed travel-time changes[J]Bulletin of the Seismological Society of America， 90（3）： 702-708

Krohn C， Johnson M， Ho R.2010.Vibroseis productivity： shake and go [ J. Geophysical Prospecting， 58（1）： 101-122

Lance GN， Williams WT.1967. A general theory of classificatory sorting strategies：1. hierarchical systems[ J]. The Computer Journal， 9（4 373-380

Lebedev av.beresney IA.2004. Nonlinear distortion of signals radiated by vibroseissources[ J].Geophysies， 69（4）： 877-1103.

Li L， Niu F L， Chen oF. et a.2017. Post seismic veloeity changes a-long the 2008 M7.9 Wenchuan earthquake rupture zone revealed by S coda of repeating events [ J]. Geophysical Journal International 208（2）：1237-1249.

Liu ZK. Huang JL， Peng zG， et al.2014. Seismic velocity changes in the epicentral region of the 2008 Wenchuan earthquake measured from three-component ambient noise correlation techniques [ J Geophysical Research Letters， 41（1）： 37-42.

MaoSJ， Campillo M， Van der Hilst R D， et al.2019. High temporal reso-lution monitoring of small variations in erustal strain by dense seis-mic Arrays[ J]. Geophysical Research Letters， 46（1）： 128-137.

MeCartney JS， Cox B R.2013. Role of strain magnitude on the deforma-tion response of geosynthetic - reinforced soil layers[ J].Geosynthel-ics International， 20（3）： 174-190

Mehta K.Snieder r.graizer V.2007. Downhole receiver funetion： a case study [ j]. Bulletin of the Seismological Society of America， 97（5）1396-1403.

Nagarajappa N， Wilkinson D.2010. Source measurement effect on high fidelity vibratory seismic separation[ J]. Geophysical Prospecting， 58（1）：55-67.

Niu F L， Silver PG， Daley TM， et al.2008.Preseismic velocity changes observed from active source monitoring at the Parkfield saFod drill site[J].Nature ，454（7201）：204-208.

Niu F L. silver pg.nadeau r m. et al.2003. Migration of seismic scat-terers associated with the 1993 Park field aseismie transient event [J].Nature ，426（6966）：544-548.

Peng ZG， Ben - Zion Y.2006.Temporal changes of shallow seismic velocity around the KaradereDiizce Branch of the North Anatolian Fault and strong ground motion [ J]. Pure and Applied Geophysics 163（2）：567-600.

Pevzner R， Shulakova V， Kepie A， et al.2011. Repeatability analysis of land time-lapse seismic data： CO2 CRC Otway pilot project case study [ J]. Geophysical Prospecting， 59（1）： 66-77.

Poletto F， Schleifer A， Zgaue F， et al.2016. Acquisition and deconvolution of seismic signals by different methods to perform direct ground force measurements [ J]. Journal of Applied Geophysies， 135（SI）191-203.

Rowe CA， Aster R C， Borchers B， et al.2002. An automatic， adaptive algorithm for refining phase picks in large seismie data sets[].Bul-letin of the Seismological Society of America， 92（5）： 1660-1674.

Sallas J.1984. Seismic vibrator control and the down going P-wave[Geophysics， 49：732-740

Sallas J.2010. How do hydraulic vibrators work？ A look inside the black box[ J]. Geophysical Prospecting， 58（1）： 3-17.

Silver PG， Daley TM， Niu F L， et a.2007. Active source monitoring of cross-well seismic travel time for stress - induced changes [ J.Bul-letin of the Seismological Society of America， 97（1）： 281-293

Sneath PHA， Sokal RR. Freeman WH.1973. Numerical taxonomy. The principles and practice of numerical classification[ M]. United states of America： W H freeman and company

Son M， Cho C S， Shin JS， et al.2018. Spatiotemporal distribution of events during the first three months of the 2016 Gyeongju， Korea Earthquake Sequence[ J. Bulletin of the Seismological Society of A-merica ，108（1）：210-217.

Tal B. Bencze A. Zoletnik s. et al.2011.Cross - correlation based time delav estimation for turbulent fow velocity measurements. Statistical considerations[ J]. Physics of Plasmas， 18（12）： 122304.

Wang Bs， Yang w， Wang w T， et al.2020.Diumal and semidiumnal P-and S-Wave velocity changes measured using an airgun source[J]Journal of Geophysical Research： Solid Earth， 125（1）：e2019 JB018218.

Wang b S， Zhu P， Chen Y， et al.2008. Continuous subsurface veloeity measurement with coda wave interferometry [ J. Journal of Geophysi-cal research： Solid Earth， 113（ B12）： B12313.

Wapenaar K， Draganov D， Snieder R. et al.2010.Tutorial on seismic interferometry： Part 1- Basic principles and applications [ J]. Geo-physics ，75（5）：195-209.

Wapenaar K， Fokkema J.2006. Green's function representations for seis-mic interferometry [ J]. Geophysics， 71（4）： S133-SI46

WeiZH， Phillips TF， Hall MA.2010. Fundamental discussions on seis-mic vibrators[ J]. Geophysics， 75（6）： W13-W25

WeiZH， Phillips TF.2013. On the generation of low frequeneies with modern seismic vibrators[ J. Geophysics， 78（2）： WA91-WA97.

Yamamura K， Sano O， Utada H， et al.2003.Long-term observation of in situ seismie velocity and attenuation [ J]. Journal of Geophysica Research： Solid Earth， 108（ B6）： 2317.

Yang w， Wang Bs， Yuan SY， et al.2018. Temporal variation of seismic-wave velocity associated with groundwater level observed by a down-hole airgun near the Xiaojiang Fault Zone [ J]. Seismological Research Letters， 89（3）： 1014-1022

Zaliapin I， Ben-Zion Y.2013. Earthquake clusters in southemn California Il： Classification and relation to physical properties of the crust[J]Journal of Geophysical Research： Solid Earth， 118（6）： 2865-2877.

Zhan Zw， Tsai VC， Clayton .2013. Spurious velocity changes caused by temporal variations in ambient noise frequency content J. Geophysical Joumnal International， 194（3）： 1574-1581

Signal Characteristics of Low Frequency Vibroseis and Its Application to High Precision Dynamic Monitoring

WANG Liwei.2.3. LI Li， YANG Wei， WANG Baoshan， HUANG Liufang

（1. Institute of Geophysics， China Earthquake Administration， Beiing 100081， China）

（2. CEA Key Laboratory of earthquake Monitoring and Disaster Mitigation Technology，Guangdong Earthquake Agency， Guangzhou 510070， Guangdong， China）

（3. Guangdong Provincial Key Laboratory of earthquake Early Warning and Safety Diagnosis ofMajor Projects， Guangdong Earthquake Agency， Guangzhou 510070， Guangdong， China）

（4. School of Earth and Space Sciences， University of Science and Technology of China， Hefei 230026， Anhui， China）

Abstract

The low frequency vibroseis can generate broad band， low energy density and high repeatable signals with fre-quencies as low as 1.5 Hz， which is particularly suitable for high resolution exploration and regional seismological research. In order to evaluate its signal characteristics and feasibility in high precision dynamic monitoring， we con-duct a field experiment on KZ28 LF low frequency vibroseis. We found that it can stably generate low frequency sig-nals in 1.5～12 Hz in upsweeping and pseudorandom sweeping mode， the near field energy is strong in 9.12 Hz and weak in 1.5-9 Hz， and the waveform correlation coefficient is greater than 0.95， accounting for 81o and 63o respectively， but it is still difficult to put sufficient energy into the ground in low frequency band. The signal repeatability is low and the harmonic energy is strong in single frequency excitation mode. Both the coherent and de-convolution method can successfully compress the source wavelet and restore the interpretable record， but the re-peatability and signal-to - noise ratio of coherent method is higher， and the deconvolution result can keep a high-er fidelity in phase and amplitude. After phase - weighted stacking， the max recording offsets can reach more than 20 km， and the monitoring precision of seismic velocity changes can reach up to .10. The result can provide a reference for high resolution monitoring of active fault zones， earthquake nucleation process and land slip.

Keywords： low frequency vibroseis; cluster analysis; deconvolution; phase weighted stack; velocity changes; high precision dynamic monitoring