球對稱有限水體氣泡振蕩模型及其對陸地水體氣槍研究的意義

徐逸鶴，王寶善，王偉濤

1 中國地震局地球物理研究所，北京 100081 2 英國劍橋大學地球科學系，劍橋 CB3 0EZ 3 愛爾蘭都柏林高等研究院宇宙物理學學院地球物理學分部，都柏林 D02Y006 4 中國科學技術大學地球和空間科學學院，合肥 230026

0 引言

利用人工震源向地下發(fā)射地震波，是開展地殼及淺部成像的重要方法.在地震活動性較弱的地區(qū)，人工源探測可增加震源點位，提升成像所需的路徑覆蓋，服務相應的成像要求(陳颙和朱日祥，2005; Li et al., 2006; Etgen et al., 2009; Zhang et al., 2011).地震波的振幅會因為幾何擴散和黏彈性衰減而隨著傳播距離減小，所以人工震源需要具有足夠的能量來保證清晰的遠距離或深部成像.炸藥震源是人工震源的重要代表，其能量密度高，使用數(shù)百到數(shù)千公斤級的炸藥可以獲得近千公里地殼結構的成像(高銳等，2006；Zhang et al., 2011).但是由于炸藥震源對環(huán)境的負面影響及本身的安全性問題，尋找替代震源成為了主動源探測的一個重要的研究方向(李麗等，2004；Alekseev et al., 2005; 羅桂純等，2006；陳颙等，2007a；李穩(wěn)等，2020).其中，之前用于海洋勘探的氣槍震源具有能量大、重復性高、環(huán)境友好等特點，具有深部探測的潛力(Dragoset, 2000；丘學林等，2007；趙明輝等，2008；王偉濤等，2017；吳志強等，2019).

氣槍震源通過在水體中瞬間釋放預先充在氣槍槍體中的高壓(可達一百多個大氣壓)氣體產生了類似于炸藥水下爆炸所產生的高壓氣泡，進而引起水體振蕩，產生地震波.氣槍震源單次激發(fā)能量高，相當于0.5～1級的地震(張尉等，2009)，并且氣槍震源具有很好的重復性，可以通過疊加多次激發(fā)的結果進一步提高探測距離(林建民等，2008；唐杰等，2009).目前，氣槍震源已在不同陸地水體中有成功性的應用，包括在水庫中的固定激發(fā)(陳颙等，2007b；Wang et al., 2012)以及在河中船拖拽的流動激發(fā)(田曉峰等，2016；張云鵬等，2016；徐逸鶴等，2016)，并獲得了體波、面波走時層析成像的結果(田曉峰等，2016；張云鵬等，2016；She et al., 2018；Zhang et al., 2020).此外，得益于氣槍震源的高重復性，氣槍還可以進行精細的時變監(jiān)測 (Wang et al., 2020).

對于氣槍震源激發(fā)過程和氣泡振蕩的研究，是深化陸內水體氣槍震源分析的重要方面.氣槍震源的激發(fā)過程可以大致分為兩個階段：(1)槍體中高壓氣體通過氣閥釋放到水中，形成初始的高壓氣泡；(2)高壓氣泡在水中的自由振蕩.第一個階段會產生一個尖銳、高振幅的信號，通常稱為主壓力脈沖或主脈沖，而第二個階段會產生一個振幅逐漸減小的周期振蕩信號，稱為氣泡脈沖.同時，主脈沖和氣泡脈沖在水面的反射會產生一個時間略微滯后的負脈沖信號.這三個脈沖信號構成了氣槍震源源區(qū)信號的主要特征.

經典的氣槍震源激發(fā)過程研究通常采用解析方法，優(yōu)勢是計算快速，缺點是需要事先評估氣槍激發(fā)物理過程的影響因素的主次，然后進行合理的簡化(Ziolkowski, 1970; Schulze-Gattermann, 1972; Ziolkowski et al., 1982; Laws et al., 1990; 王立明，2010).通常只考慮高壓氣泡的自由振蕩階段，而將高壓氣泡的形成階段簡化為氣泡壁運動的初始條件，氣泡上浮、變形等次要因素也予以忽略.模型核心是求解球形氣泡在無限可壓縮無黏流體中的振蕩過程 (Lamb, 1923; Gilmore, 1952; Keller and Kolodner, 1956; Ziolkowski, 1970; Plesset and Prosperetti, 1977; Brennen, 2013; de Graaf et al., 2014).

經典氣槍模型和實測海洋氣槍信號有很好的一致性，說明相關簡化符合海洋激發(fā)壞境中的實際物理過程.但是在陸內水體氣槍研究中，觀測到的陸內氣槍源區(qū)信號與經典氣槍模型的預測存在較大偏差(唐杰等，2009；唐杰，2010).該差別很可能是由于陸內水體較小而導致的，但是對于陸內水體如何影響氣槍激發(fā)過程還缺乏系統(tǒng)的研究.通過數(shù)值模擬對不同水體形狀、水深、氣槍位置的模擬可以解釋部分趨勢，但是由于使用的是預設的子波，更多地關注壓力波在水體中傳播及其轉換為地震波過程，不能研究有限水體對氣泡振蕩本身的影響(孫楠等，2017；胡久鵬等，2017)，因此無法解決臨界水體大小的問題.

研究有限水體對氣槍震源初始激發(fā)信號的影響是目前亟待解決的關鍵問題之一.在經典氣槍理論的球對稱氣水雙層模型基礎上增加一層固體，即中心為球形氣泡，中間層為球形水體(中空)，外層為無限大的固體，可以對有限水體中氣泡的振蕩進行研究.該模型也經過了適當簡化，但這一簡單模型的優(yōu)點是可以避免其他因素的干擾，更準確地分析水體大小這一因素的獨立影響.該有限水體氣槍模型在水體半徑足夠大時，氣泡振蕩過程接近無限水體模型，而球腔半徑逐漸減少時，固液界面的影響會逐漸增強.當氣泡振蕩過程開始與無限水體模型出現(xiàn)明顯偏差時，可以將對應水體大小視為臨界水體.當水體大于該臨界大小時，可以使用經典氣槍模型(無限水體)來模擬氣槍信號的震源時間函數(shù)；當水體小于臨界大小時，需要考慮有限水體對氣泡振蕩過程本身的改變.

利用解析方法求解該有限水體模型存在較大的困難，原因是固液邊界的反射波使氣泡振蕩方程中出現(xiàn)了時延項.Doinikov等(2018)得到的有限水體模型的氣泡振蕩方程包括11項，比無限水體的氣泡振蕩方程增加了近2倍(4項)，且每一項的形式也更為復雜.此時使用理論解的難度已經高于數(shù)值解，并且理論解具有很弱的擴展性.

使用常規(guī)數(shù)值方法對該問題進行精確求解也同樣存在一定的難度，包括：(1)氣體、水體、固體滿足不同的控制方程，且水體模擬多用歐拉描述，即同一位置在不同時刻的運動狀態(tài)，而固體模擬多用拉格朗日描述，即同一質點在不同時刻的運動狀態(tài)；(2)氣水界面的高密度差、高壓力差容易產生數(shù)值振蕩；(3)精確地描述和追蹤氣水界面需要特殊的方法(如Level-set，Volume of Fluid等)；(4)氣泡尺度(1 m)與模擬區(qū)域尺度(百米)相比相差較大；(5)氣泡振蕩過程的準確模擬通常需要以很小的時間步長(約10-5s)計算1 s (King, 2016)；(6)常規(guī)模擬方法多采用不可壓縮流體，對高壓氣泡的振蕩過程模擬可能不準確(董明榮等，2016).因此使用常規(guī)數(shù)值模擬方法不僅較難實現(xiàn)，計算量大，而且還可能會犧牲結果的精度.

采用歐拉有限元通過組合拉格朗日和歐拉描述，可以對水下爆炸產生的高壓氣泡進行較好地模擬(Liu et al., 2018b). 該方法的第一步是采用拉格朗日描述的有限元方法對流體力學方程在軸坐標下進行求解.本研究利用有限水體氣槍模型的球對稱屬性，在球坐標下重新推導了有限元數(shù)值格式.由于水體在有限的空間內振蕩，用拉格朗日有限元替代歐拉有限元，進一步簡化了問題的求解.由于氣泡也為球對稱，可以將模擬區(qū)域的邊界設在氣液邊界和固液邊界上，既將三種介質的模擬簡化為單一介質的模擬，又可以精確、顯式地追蹤氣泡(即氣液邊界)的變化.氣泡和固體的作用以邊界條件的形式施加在水體上.

本文將介紹水體控制方程及其球對稱條件下拉格朗日有限元的推導過程，并以此分析氣泡半徑的振蕩曲線、氣泡振蕩周期和振幅與水體大小的關系以及氣槍振蕩過程受影響的機制.通過數(shù)值模擬方法，對氣槍在有限水體中的激發(fā)過程，不同激發(fā)參數(shù)對水體臨界大小的影響進行了討論分析.

1 有限水體內的氣泡振蕩模型及其有限元模擬

1.1 有限水體氣泡振蕩模型及控制方程

氣槍在陸內水體中的激發(fā)問題可以近似為一個球形高壓氣泡在充滿水的球形空腔中的振蕩問題，空腔外是無限大的固體介質(圖1).由于整個問題的球對稱性，氣泡和水體只存在徑向運動.將復雜的氣槍震源過程近似為一個球形高壓氣泡在水中的自由振蕩問題是氣槍和水下爆炸研究中的常規(guī)處理，與實際氣槍實驗數(shù)據具有很好的一致性(Ziolkowski, 1970, 1998; Johnson, 1994)，具體的近似包括(Ziolkowski, 1970)：(1)忽略重力，因為高壓下氣泡運動的加速度遠大于重力加速度；(2)忽略水的黏性，因為黏性項的貢獻遠小于速度項的貢獻(約百萬分之一)；(3)忽略氣泡的表面張力，因為其遠小于高壓氣泡產生的壓力.

圖1 有限水體內氣槍激發(fā)模型R和RS分別是氣泡半徑和水體半徑，其中R=R(t)是時間的函數(shù)，RS是常數(shù).Fig.1 Air gun model in the finite water bodyR and RS are the radii of the air bubble and the water body, respectively. Of which R=R(t) is a function of time and RS is a constant.

水體運動的控制方程由動量守恒、質量守恒和介質狀態(tài)方程組成.在拉格朗日描述下，無黏可壓縮球對稱水體的動量守恒方程為

(1)

其中ρ,u,p分別為水體的密度、速度和壓力，t和r為時間和徑向距離，D/Dt為物質導數(shù).在拉格朗日描述下，坐標跟隨物質點移動，因此任意指定區(qū)域內的質量不變(區(qū)域可能會變形)，

(2)

以下壓力和密度的經驗公式使這一關于ρ,u,p三個變量的方程組達到閉合 (Ziolkowski, 1970)，

(3)

其中p∞和ρ∞是無擾動狀態(tài)下水體的壓力和密度.n和B是經驗參數(shù)，并沿用Ziolkowski(1970)選擇：n=8,B=2500大氣壓.

氣泡和固體的作用通過邊界條件的形式施加于水體.氣液邊界處壓力連續(xù)(忽略氣泡的表面張力)，所以

(4)

其中R=R(t)是氣泡半徑，p0和R0是初始時刻的氣泡壓力和半徑，γ=1.13為經驗常數(shù) (Ziolkowski, 1970).剛性固體邊界的速度為零，因此，u|r=RS=0.

(5)

初始的氣泡體積等于氣槍容積，壓力等于氣槍槍壓，水體無擾動狀態(tài)的壓力為氣槍所在深度的靜水壓(大氣壓+ρg×氣槍深度).

1.2 基于拉格朗日描述的顯式有限元的數(shù)值實現(xiàn)

有限元方法首先使用形函數(shù)對速度等變量進行近似，然后求解控制方程對應的弱解.研究區(qū)域(從水體內邊界到外邊界)等分為N個單元，其中第i個單元對應從ri到ri+1的球殼，而ri=ri(t) 是第i個節(jié)點的位置.記ui,ai為節(jié)點i上的速度、加速度，其他位置上的值可由插值得到，

(6)

其中Ni(x)是第i個節(jié)點的權重函數(shù)，或稱為形函數(shù)，取如下的線性格式：

(7)

壓力和密度則定義在單元上，在每個單元內為常數(shù) (Benson, 1992).

弱解并不要求處處滿足動量方程，而是滿足動量方程的某種加權積分，采用形函數(shù)作為權函數(shù)的方法稱為伽遼金方法.使用伽遼金方法，動量方程(式(1))轉化為

(8)

(9)

節(jié)點位置更新后，可以得到單元的新體積，進而根據單位質量守恒得到單元的密度.單元的壓力根據狀態(tài)方程(3)由計算得到，然后可以進行下一個時間步的計算.

1.3 數(shù)值模擬的參數(shù)選擇

時間步長的選擇由空間網格間隔和最快速度決定，過大的時間步長會導致數(shù)值模擬的不穩(wěn)定.水中速度有兩個：聲波速度(1500 m·s-1)和激波速度.根據前人的研究結果，高壓氣泡產生的激波速度也約為1500 m·s-1(Liu et al., 2018b)，因此本文中采用min(Δr)/(4×1500)為時間步長，其中min(Δr)代表最小的空間網格，經驗因子4是為了保證數(shù)值模擬的穩(wěn)定.因為空間網格在模擬過程中會隨著水體變形，所以時間步長在每個時間步后也會動態(tài)更新.

空間間隔通常由壓力波場的最短波長決定，以保證對波場的準確描述.但激波的存在導致壓力場存在間斷，不能預先確定合適的最短波長.空間間隔為0.02 m到2 m的9個不同的實驗結果表明Δr≤0.05 m的氣泡振蕩曲線沒有明顯區(qū)別(圖2)，因此空間間隔選擇為Δr=0.05 m.

圖2 空間間隔測試黑色實線為氣泡半徑隨時間的變化曲線，2條黑色垂直虛線標注的是氣泡收縮的時間(由Δr=0.02 m的氣泡半徑曲線確定)，空間間隔標注在右側.該測試中使用的水體半徑為300 m.Fig.2 Test of the node intervalBlack solid lines show the bubble radius as a function of time. The two black dashed lines mark the time of the bubble contraction, which is determined by the bubble radius obtained with Δr=0.02 m. The node interval is labelled on the right side. The radius of the water body used in the test are 300 m.

2 有限水體內高壓氣泡振蕩的模擬結果

為了驗證數(shù)值模擬結果的正確性，本文將大水體的有限水體氣泡振蕩曲線與無限水體的結果進行對比.采用的參數(shù)如下：氣槍體積為1500 in3(0.02458 m3，1 in3=1.63871×10-5m3)，槍壓為15 MPa，靜水壓以沉放深度為10 m計算，水體半徑為300 m，模擬時長為0.5 s.按反射波的速度為1500 m·s-1計，兩個模型在約0.4 s之前均為自由振蕩，振蕩過程應保持一致.無限水體的氣泡振蕩曲線使用Johnson(1994)的簡化公式計算得到.兩者對比表明，有限水體結果與無限水體結果的振蕩周期和氣泡大小在前兩個振蕩周期內具有較好的一致性(圖3).兩者之間的少許不同可能是由于Johnson模型的簡化導致.

圖3 與無限水體氣槍模型的對比黑線為無限水體條件下氣槍模型計算得到的氣泡振蕩曲線，紅線為本文有限水體模型的計算結果.有限水體的水體半徑為300 m，空間間隔為0.05 m.Fig.3 Comparison with the air gun model in infinite water bodiesThe black line is the bubble radius computed by the air gun model in an infinite water body. The red line is the bubble radius computed using our air gun model in a finite water body. The radius of the finite water body is 300 m. The node interval is set to 0.05 m.

為了研究有限水體效應作用于氣泡振蕩過程的影響，本文計算了14個算例，其中水體半徑從10 m變化到500 m.氣泡振蕩曲線的對比表明，當水體半徑大于300 m時，整個0.5 s長的氣泡振蕩曲線趨于一致，且與無限水體的結果一致(圖4).水體半徑為200 m時，氣泡振蕩曲線在第二個周期開始偏離大水體的結果.當水體半徑小于100 m時，振蕩曲線的第一個周期也開始偏離大水體的結果，并且振蕩周期隨水體半徑減小而縮短.該結果證實了大水體中的氣泡振蕩曲線趨于無限水體的結果，而在小水體中出現(xiàn)明顯偏離.

圖4 不同大小水體中的氣泡振蕩曲線垂直虛線是RS=500 m對應的氣泡收縮時間.Fig.4 Bubble radius in water bodies of different sizeVertical dashed lines are the time of bubble contraction corresponding to RS=500 m.

圖5 氣泡振蕩周期和水體半徑的關系黑色虛線為水體趨于無窮大時氣泡振蕩周期的極限.Fig.5 Period of bubble oscillation as a function of the radius of the water bodyThe black dashed line marks the limit of the period of bubble oscillation when the size of the water body approaches infinity.

使用兩個氣泡振蕩的關鍵參數(shù)，氣泡振蕩周期(第一個振蕩周期)和最大氣泡半徑，可以定量地確定有限水體效應的臨界水體半徑.振蕩周期的結果表明，水體半徑小于100 m時，振蕩周期隨著水體變大逐步增加，而半徑大于100 m時，振蕩周期穩(wěn)定于0.17 s(圖5).最大氣泡半徑呈現(xiàn)相似的趨勢，水體半徑小于60 m，氣泡最大半徑從0.78 m (RS=10 m)逐漸增加到1.17 m (RS=60 m)；水體半徑大于60 m，氣泡最大半徑穩(wěn)定在1.17 m(圖6). 雖然兩個關鍵參數(shù)得到的臨界水體半徑有些許不同(100 m和60 m)，不過均再次支持了小水體對氣泡振蕩過程存在明顯的影響.

圖6 氣泡最大半徑和水體半徑的關系Fig.6 Maximum bubble radius as a function of the radius of the water body

圖7 有限水體氣泡振蕩曲線(紅線)與“無限”水體(黑線)的對比“無限”水體的模擬結果由水體半徑為500 m的有限水體模型得到，因為最快的反射波(0.667 s)在最大模擬時長(0.5 s)之后才會到達氣泡壁.紅色豎線標記了有限水體中的氣泡振蕩最早偏離“無限”水體曲線的時間，對應的速度為1338 m·s-1.Fig.7 Comparison of bubble oscillations in finite water bodies (red lines) to those in the ‘infinite’ water body (black lines)We used the finite water body simulation with the radius of 500 m as the result of ‘infinite’ water body.The fastest reflected waves arrive at the bubble at 0.667 s, later than our longest simulation time (0.5 s). The vertical red lines mark the earliest time when the bubble oscillations in finite water bodies deviate from the ‘infinite’ result, corresponding to the velocity of 1338 m·s-1.

氣泡振蕩曲線偏離自由振蕩的時間可以提供理解有限水體效應的關鍵證據.以500 m水體的模擬結果作為參考(在前0.5 s處于自由振蕩階段)，10 m到400 m水體的結果的最初偏離時間與水體半徑存在線性的關系，對應于速度為1338 m·s-1的應力波從氣泡出發(fā)、經固液界面反射后返回到氣泡的時間(圖7).

水中壓力場快照可以用來查看壓力波的傳播速度，進一步找到速度為1338 m·s-1的壓力波的來源.半徑300 m的水體的波場快照展示一個傳播速度為1338 m·s-1的壓力波(圖8).這一壓力波是氣槍槍壓和靜水壓之間的強壓力間斷產生的激波.雖然激波后續(xù)的波形存在明顯的數(shù)值振蕩，但其激波的前緣與1338 m·s-1的速度參考線吻合很好.

圖8 壓力波場快照.為了更好地顯示激波，我們選擇聚焦于10 m以內的壓力場虛線是速度為1338 m·s-1的參考線.Fig.8 Snapshots of pressure wavefield. For better illustration of shock waves, we zoom in on the pressure field within the distance of 10 mThe dashed line is the reference line of velocity of 1338 m·s-1.

3 討論

球對稱有限水體氣槍模型的有限元模擬表明，小水體和大水體的氣泡振蕩曲線存在系統(tǒng)性區(qū)別.不同尺度的水體中氣泡振蕩曲線在振蕩初期具有很好的一致性，而后出現(xiàn)偏離，且偏離發(fā)生的時間與水體的半徑存在線性關系.該線性關系對應的反射波速度與實驗中激波速度吻合.氣泡振蕩周期和氣泡最大半徑表明，當水體半徑小于60～100 m時，氣泡振蕩過程會與無限水體氣槍模型存在較大的差別.當水體半徑大于100 m時，氣泡振蕩過程與無限水體模型比較接近，這時使用無限水體模型計算震源子波不會有明顯的誤差.

3.1 氣槍在陸內水體的激發(fā)過程及其在海洋中激發(fā)的區(qū)別

氣槍在海洋和陸內水體中的過程具有一定區(qū)別，對其進行對比利于理解陸內水體的氣槍激發(fā)模型.本質而言，氣槍在海洋和陸內水體的激發(fā)過程具有很高的相似性：先在氣槍內充入高壓空氣，然后在水中釋放.釋放時氣槍內外有幾十到幾百大氣壓的壓力差，高壓力差會激發(fā)在水中傳播的壓力間斷——激波，產生主壓力脈沖(Liu et al., 2018a).氣泡兩側的壓力差還會將水體推離氣槍所在位置，導致高壓空氣泡的膨脹及其壓力降低.當氣泡內外壓力相等時，水體會由于慣性的原因繼續(xù)往外移動，但負壓力差會降低氣泡擴張的速度，直到氣泡擴張速度降為零，此時氣泡膨脹到最大.隨后，氣泡由于負壓力差開始收縮，經過氣泡內外壓力相等的位置時，會同樣由于慣性的原因繼續(xù)收縮，但收縮速度會開始減慢.當氣泡收縮到最小時，氣泡壓力差為正，然后開始新一輪的膨脹過程.這種周期性的氣泡膨脹收縮過程會產生壓力波.這些壓力波會穿過水底和岸邊轉換為地震波，同時也有部分壓力波的能量會在水面和固液邊界之間發(fā)生反射或多次反射后再轉換為地震波.

海洋勘探中，往往利用的是氣槍高壓氣體釋放時產生的高頻信號，用來提高勘探的分辨率.而在陸內水體激發(fā)時，為了實現(xiàn)更遠距離的探測，所分析的信號為氣泡振蕩產生的相對低頻的信號.就氣泡振蕩而言，海洋和陸內水體氣槍的主要區(qū)別在于氣槍離固液邊界的距離.海水深度為幾百米到幾公里，而陸內水體的深度通常為十幾到幾十米.由于氣槍信號的振幅衰減與傳播距離呈正相關，小水體意味著反射波的振幅更強，更可能對氣泡振蕩本身產生影響.

反射波影響氣泡振蕩的原理是反射波到達氣泡壁時會改變氣泡內外的壓力差，進而改變氣泡振蕩的過程.如圖9所示，當反射激波(壓力間斷)到達氣泡壁時，會迅速增加氣泡壁上的壓力，在短時間內(10 ms)形成負加速度(圖9b)，降低氣泡膨脹的速度，使得氣泡振蕩過程偏離了自由振蕩(圖9a).雖然之后的二次反射(到時約0.014 s)又短暫提高了速度，但已無法將氣泡膨脹的過程校正回原本的程度.

圖9 反射激波到達氣泡壁時氣泡半徑、加速度和氣泡外壁上的壓力(a) 氣泡半徑； (b) 氣泡加速度，內嵌圖展示的是水體半徑為5 m和300 m的氣泡壁上的壓力差.作為參考的水體半徑為300 m的結果在前0.015 s仍處于自由振蕩階段.黃色區(qū)域顯示的是反射激波到達氣泡壁的階段.Fig.9 Bubble radius, acceleration and pressure on the bubble wall when shock wave reflected to the bubble(a) Bubble radius; (b) Bubble acceleration. The inset map shows the difference in pressure on the bubble wall between the radius of the water body of 5 m and 300 m. The result with RS=300 m is still in undisturbed bubble oscillation within the first 0.015 s. The yellow region highlights the period when the reflected shock wave arrives the bubble.

3.2 氣槍激發(fā)參數(shù)的影響

氣槍的槍容、槍壓和沉放深度是氣槍激發(fā)的重要參數(shù)，三者決定了無限水體內氣泡振蕩過程，也在有限水體模型中起著重要作用.

改變槍容、槍壓和沉放深度對氣泡振蕩周期產生了明顯的影響(圖10).當水體較大時，單獨將槍容從1500 in3增加到8000 in3會導致氣泡振蕩的周期從0.174 s增加到了0.3035 s，降低槍壓(15 MPa到7 MPa)會減小振蕩周期到0.1355 s，增加沉放深度會減小振蕩周期到0.1445 s，因此氣泡振蕩周期的改變與氣槍的槍容、槍壓呈正相關，與沉放深度呈負相關，換言之，如果想要降低氣槍信號的氣泡振蕩階段的主頻(低頻信號衰減更慢，探測距離更遠)，可以采用增加槍容(或者采用多槍同時激發(fā))、增加槍壓、或者減小沉放深度的方法.這三種方法的一致作用是導致氣泡最大半徑的增加，因此可以簡單地理解為大氣泡的振蕩周期長，小氣泡的振蕩周期短.

氣槍釋放能量的公式為(Ronen, 2002; 陳颙等, 2007b)

(10)

其中P是槍壓，V是槍容，P0是靜水壓力(主要由沉放深度決定).該公式表明，導致氣泡最大半徑增加的參數(shù)改變也會增加氣槍能量，因此水體半徑較大時，氣泡振蕩周期、氣泡最大半徑和氣槍釋放能量三者呈正相關.

當水體半徑很小時(<20 m)，改變這三個參數(shù)對于氣泡振蕩周期的影響很小(圖10).可能意味著在這個水體尺度下，有限水體效應已經占據了主導作用，改變氣槍本身的激發(fā)參數(shù)已較難控制氣槍激發(fā)頻率.

其次，仔細觀測曲線的變化可以發(fā)現(xiàn)，大氣泡的臨界水體半徑比小氣泡的更大，意味著在同等尺度的水體中，大氣泡受到的有限水體效應更明顯.

3.3 理論模型與陸內水體氣槍實驗的對比分析

在已經開展的陸內水體氣槍實驗中，在不同大小的水體中均開展了氣槍激發(fā)，并在近場記錄了震源信號.將有限水體模型模擬結果與實測結果進行對比，可以驗證相關結果的可靠性.

在云南賓川大銀甸人工水井中的實驗是典型的小水體：水井半徑2.5 m，井深20 m.氣槍沉放深度為12 m，距離水體邊界的垂向和橫向距離分別是8 m和2.5 m，平均距離約為5 m(胡久鵬等，2017).根據圖5中展示的振蕩周期曲線的趨勢，對應的周期約在0.2～0.3 s, 即33～50 Hz，與實測氣槍信號的氣泡振蕩主頻(約35 Hz)有很好的一致性(胡久鵬等，2017).這說明在小尺度的陸內水體中，需要考慮有限水體效應，且本文的有限水體模型可以很好地估計氣泡振蕩的頻率.

另一個小水體的例子是北京昌平開展的馬刨泉實驗(Wang et al., 2010).馬刨泉是一個水面直徑為25～37 m的人工水池，氣槍沉放深度為11 m，沉放處的水深為14.5 m.氣槍與水體邊界的垂向(到水底)和橫向(到岸邊)距離約為3.5 m和15 m，取其平均約為10 m.根據圖5中展示的振蕩周期的曲線，水體半徑為10 m的球對稱有限水體氣槍信號的振蕩周期約為0.05 s，即20 Hz，與實測的主頻(20～25 Hz)具有很好的一致性(Wang et al., 2010).

新疆呼圖壁人工水池是一個軸對稱的碗形水體，水面半徑為48.5 m，深15 m(魏斌等，2016).氣槍位于水體中心，距岸邊約48.5 m，距水底5 m，平均距離約25 m，對應的振蕩頻率為6.6 Hz.實測氣槍信號由高頻(40～50 Hz)的主脈沖和低頻(約3 Hz)的氣泡振蕩組成(胡久鵬等，2017)，理論預測的頻率低了一半.可能的解釋是呼圖壁實驗使用了6條氣槍，氣泡能量更大，導致振蕩周期降低.

賓川大銀甸水庫水面直徑更大(約1 km)，但由于水面形狀不規(guī)則，且氣槍不在水庫中心(離岸邊30 m處)，氣槍距離水體邊界最近為3 m(水底)，離岸邊最近為 30 m，最遠距離為1200 m，情況更為復雜.實測的氣槍信號中有多個頻率，分別為5，10，20，40 Hz，與有限水體模型預測的振蕩周期對比(圖10a)，分別對應于水體半徑為40，20，10，5 m.可能的解釋是：(1)氣槍離水體邊界的距離變化產生多個振蕩周期；(2)最低的振蕩周期與水體自由振蕩有關，而更高頻的信號由主脈沖或諧振產生.進一步檢驗具體原因可能需要對非規(guī)則的水體中的高壓氣泡振蕩問題進行數(shù)值模擬.

3.4 未來可能的改進

本文描述的數(shù)值模擬方法雖然可以很好地研究有限水體對氣槍振蕩過程的影響，但激波后存在明顯的數(shù)值振蕩(圖8)，可以通過增加數(shù)值黏性來解決 (von Neumann and Richtmyer, 1950; Brode, 1955; Persson and Peraire, 2006).此外，將固體介質從剛體改為彈性體，應該可以更準確地模擬氣泡的振幅，并給出準確的震源子波.這兩個改進可以在球對稱的框架下進行.

非規(guī)則水體對氣泡振蕩的影響不能繼續(xù)使用球對稱框架，需要在二維或三維條件下求解水體的控制方程.該情況下將網格依附在氣泡邊界上的做法可能較為復雜，更可行的方法是將氣泡、水體、固體都包括在計算域中，利用Level set或Volume of Fluid method (VOF)等方法來描述邊界(King, 2016; Liu et al., 2018b).

目前氣槍模型中主脈沖被簡化為高壓氣泡的初始壓力、體積和速度.為了更好地描述主脈沖，可以采用數(shù)值模擬方法對氣槍槍體和氣閥建模，模擬高壓氣體噴射出槍體形成高壓氣泡(或者泡沫流)的過程.這一研究可能需要比氣泡振蕩模擬更密的空間網格，可以單獨模擬后作為初始條件加入到有限水體氣槍模型中.

4 結論

充液球腔的球對稱有限水體模型可以較好地解決陸內水體氣槍的關鍵參數(shù)：無限水體模型失效的臨界水體大小.使用球坐標系下的一維顯式有限元方法的數(shù)值模擬結果表明，臨界水體大小約為100 m.氣槍在小于這個尺度的水體中激發(fā)時，會受到固液邊界反射激波的影響而偏離自由振蕩過程.臨界水體大小受到初始槍壓、槍容和沉放深度等因素的影響，與槍壓、槍容呈正相關，與沉放深度呈負相關，且三者可統(tǒng)一為與氣槍能量呈正相關.有限水體氣槍模型可以幫助理解陸地小水體對氣槍震源激發(fā)過程的影響，定量地衡量有限水體效應.未來可以在該模型的基礎上考慮更復雜的情況，比如氣泡在非規(guī)則水體中的振蕩過程.

致謝謹此祝賀陳颙先生從事地球物理教學科研工作60周年．