火箭發(fā)射的次聲信號分析?

戴翊靖 滕鵬曉 呂 君 姬培鋒 楊 軍 程 巍

(1 中國科學院聲學研究所噪聲與振動重點實驗室 北京 100190)

(2 中國科學院大學 北京 100049)

0 引言

次聲是頻率低于20 Hz 的聲波，在大氣中傳播損耗小，傳播距離遠。因此在一些自然或人為的突發(fā)事件的遠程監(jiān)測中，利用次聲是一種重要的手段。地震、閃電、火山噴發(fā)等自然現(xiàn)象以及火箭發(fā)射、爆炸等人類活動都伴隨著次聲的產(chǎn)生［1?4]。因此對次聲信號的研究有助于進一步認識上述諸多現(xiàn)象的物理特性。火箭發(fā)射是一類具有重要研究價值的次聲源，由于發(fā)射時間和地點是確定的，對研究次聲在大氣中的傳播路徑與大氣層結構有重要意義。

火箭的飛行任務由不同的階段組成，例如點火、助推飛行和殘骸載入等，而火箭在不同階段會產(chǎn)生不同頻率和不同強度的次聲波。因此，每次火箭發(fā)射都可以看到沿著軌跡，不同階段所輻射的次聲波。根據(jù)目前的文獻記錄，次聲監(jiān)測陣最遠可采集到相距4500 km的大型火箭發(fā)射信號［5]。這些遠程信號的頻率通常處于0.01～10 Hz，持續(xù)時間為幾分鐘。受平流層水平風在沿傳播方向的分量的影響，距離約1000 km 處次聲接收陣采集的信號的振幅從幾十mPa到大于2 Pa不等［6?8]。Blom等［9]使用多個次聲接收陣監(jiān)測一次火箭發(fā)射，采用傳播模型法很好地將火箭點火預計到達時間與觀測相匹配。Cotton等［10]使用次聲陣采集到了在軌道高度大型火箭經(jīng)過的信號，并結合氣體動力學對幅度進行計算與分析。在短距離范圍內(nèi)的靜態(tài)火箭發(fā)動機測試中，發(fā)現(xiàn)次聲信號持續(xù)時間與燃燒時間、發(fā)動機推力幅度相關，而頻率成分與排氣速度、推力和火箭發(fā)動機類型等因素相關［9,11]。

2019年12月27 ，我國在中國文昌航天發(fā)射場發(fā)射了長征五號遙三運載火箭。本文將著重研究此次火箭發(fā)射產(chǎn)生的次聲信號，并結合已有的研究進行分析比較。本文從采集信號中識別出遠程火箭發(fā)射次聲監(jiān)測中難以檢測的聲爆事件。通過經(jīng)驗模態(tài)分解(Empirical mode decomposition,EMD)預處理與Fisher 檢測算法，對火箭發(fā)射初始階段運行狀態(tài)進行分析，并估計了點火與聲爆信號的預計達到時間與方位，驗證了識別的結果。根據(jù)采集信號聲譜圖，發(fā)現(xiàn)火箭點火發(fā)射和火箭飛行達到超聲速時聲爆激發(fā)的次聲信號的頻帶較寬，兩事件間的信號低頻段(0.4～5 Hz)能量較弱的現(xiàn)象。結合火箭加速飛行過程，給出了低頻段能量較弱現(xiàn)象的解釋。

1 Fisher檢測理論

Fisher 檢測是一種常用的次聲檢測方法，其本質是方差比率檢測。根據(jù)對多通道信號能量進行分析，可以較快地找出聲源信號。Fisher 檢測法是通過構造F統(tǒng)計量，結合陣元的位置信息，計算Fisher比值與預先設定的閾值進行比較與分析。當用構建的方差比計算某個時間窗的Fisher 比值大于設定的閾值時，認為此時間窗內(nèi)的數(shù)據(jù)與常規(guī)隨機噪聲有顯著性差別，將它判定為信號［12]。

假設次聲監(jiān)測陣列有N個傳感器，記為X=[x1(t),x2(t),···,xN(t)]T。對于長度為T的時間窗，N個通道的待估計信號的能量估計可表示為

在時刻t，所有通道數(shù)據(jù)的隨機變化反映出了噪聲的存在?？捎媒M內(nèi)方差表示時刻t的這種變化：

則對于長度為T的時間窗，N個通道的平均預測噪聲能量可表示為

而Fisher 比為組間方差和組內(nèi)方差的比值，記作F，表示為

由此也可以得到信噪比和Fisher 比值之間的關系：

頻率波數(shù)(F-K)法是一種陣列信號處理方法，可以用來估計信號的方位角與視速度。將Fisher比值結合慢度網(wǎng)格和陣元的時延信息，可以較好地提升F-K 法估計方位角的分辨能力［13?14]。對于單頻的F-K功率譜密度可表示為

對于未知的信號，計算某一時間窗內(nèi)的頻譜按照給定慢度網(wǎng)格中的慢度矢量的Fisher比值。之后選擇Fisher 比值最大的慢度網(wǎng)格，這樣最有可能包含信號的來波方向就被選定，表示如下：

其中，Ps表示待測信號的功率，而噪聲功率由傳感器所記錄的總體功率Pt減去Ps計算而得。

2 火箭發(fā)射及次聲信號特征

2.1 次聲監(jiān)測陣部署

利用非線性漸進波動方程(Nonlinear progressive wave equation,NPE)仿真了聲波在預期聲源位置向茂名方向在大氣中傳播的相對能量強度分布，如圖1所示，確定影區(qū)范圍。根據(jù)影區(qū)、地形、天氣狀況以及預期射向等因素，在距離發(fā)射點約20 km 和230 km 處分別搭建了一個機動式次聲監(jiān)測陣以采集期望的信號。次聲監(jiān)測陣與發(fā)射點的相對位置如圖2所示。

圖1 聲波在大氣中的傳播損失分布Fig.1 The distribution of sound wave propagation loss in atmosphere

圖2 火箭發(fā)射點與次聲監(jiān)測陣相對位置Fig.2 The relative position of the rocket launch point and the infrasound arrays

對于監(jiān)測陣A，由于地形的限制，此點采用3元陣，陣列孔徑為170 m，發(fā)射點位于A 陣北偏東174.7°，距離約20 km。對于觀測陣B，此點采用4元陣，陣列孔徑為218 m，發(fā)射點位于B陣北偏東108°，距離約230 km。A、B 兩陣陣元的布放的相對位置如圖3(a)、圖3(b)所示。每個陣元包括一個次聲傳感器和一套降風噪設備，所用數(shù)據(jù)采集器的采樣率為100 Hz，通過無限傳輸將數(shù)據(jù)傳送到處理中心。設備的布設如圖3(c)所示。

圖3 次聲采集陣列信息Fig.3 Infrasound receiving array information

2.2 次聲信號特征

本次實驗分析選取采集信號中長度為1.5×105點的一段數(shù)據(jù)，即火箭點火后1500 s 的數(shù)據(jù)。A、B兩陣采集到的火箭發(fā)射所產(chǎn)生的次聲信號的時域波形如圖4所示?？梢悦黠@地看出，A陣采集到信號在火箭發(fā)射的初始時期，由于距離聲源較近，采集到的信號的幅度較大，信噪比較大，隨著時間的推移，采集到信號的振幅逐漸減小，反映出了火箭向著遠離A 陣的方向運動；B 陣采集的原始信號中火箭發(fā)射的幅度較小，在1200 s 左右信號振幅有一個明顯的跳變，根據(jù)實驗記錄與陣元對應的到達時間差分析，判斷為爆破干擾聲，與火箭發(fā)射次聲信號無關。

圖5(a)、圖5(b)分別為A、B兩陣采集的第一路(即圖4(a)、圖4(b)中的第一個)信號的聲譜圖。A陣距離聲源較近，信號能量主要集中在火箭點火后50～700 s的時間段內(nèi)；B 陣距離聲源較遠，信號能量主要集中在750～1200 s 之間。聲譜圖上顯示火箭發(fā)射初期，能量分布的頻帶較寬且值較大，隨著時間的推移，能量主要分布的頻帶逐漸減小，表明隨著傳播距離的增加，信號逐漸衰減，且衰減的速率與頻率的大小呈正相關。

圖4 A、B 兩陣采集的信號Fig.4 Received signals from Array A and Array B

圖5 A、B 兩陣采集信號的聲譜圖Fig.5 The spectrograms of received signals from Arrays A and Arrays B

由圖5(a)、圖5(b)中框1所示，在次聲接收陣采集到信號的初始階段，除了本底噪聲，均存在一個頻帶較寬、能量較大的信號。A 陣距離聲源較近，信號約在火箭點火后56 s 到達，由A 陣與發(fā)射點的距離可求得信號傳播的平均速度約為344 m/s，判斷接收的信號為直達波。B 陣距離聲源較遠，信號約在火箭點火后795 s到達，由B陣與發(fā)射點的距離可求得信號傳播的平均速度約為289 m/s，根據(jù)經(jīng)驗判斷該信號是通過平流層波導傳播到達的信號。此階段火箭點火開始爬升，距離地面較近，發(fā)動機排氣尾焰引起劇烈的空氣擾動，產(chǎn)生一個較強的次聲信號。發(fā)射初期火箭從靜止狀態(tài)開始做加速運動，速度較慢，需要較長的時間遠離地面，故此階段引起空氣擾動的現(xiàn)象將持續(xù)一段時間，據(jù)觀察本次火箭發(fā)射此過程持續(xù)約12 s。

圖5(a)中，在框1 之后的一段時間內(nèi)，0.4～5 Hz 的低頻頻帶能量較低。而在框2 處，在這個低頻頻帶上有明顯的能量分布。在火箭點火后60～70 s 之間，火箭飛行速度加速至超聲速，會產(chǎn)生聲爆現(xiàn)象，根據(jù)表1 中的數(shù)據(jù)與陣列的位置，可計算兩者距離約為21 km，傳播的時間大約為62 s。則聲爆信號會在點火后122～132 s 之間傳遞至A 陣，與框2 所在的時間相符。因此可以判斷，框1 為火箭點火沖擊波激發(fā)的低頻聲信號，框2 為聲爆激發(fā)的次聲信號。在框1 和框2 之間，信號主能量分布在7 Hz 以上，此時火箭逐漸升高，產(chǎn)生的信號主要由排氣尾焰與大氣作用產(chǎn)生，其所含頻率成分較高，在大氣傳播中衰減較大，因此只有A 陣采集到此信號。在156 s 之后的時間段，信號能量主要集中在5 Hz 以下，且持續(xù)時間較長。圖5(b)中，與圖5(a)在0.4～5 Hz 的低頻上有相似的變化趨勢，但在框1和框2之間，沒有明顯的信號能量分布，該現(xiàn)象由次聲遠距離傳播衰減所致。

表1 火箭超聲速時的參數(shù)Table 1 Parameters of rocket at transonic and supersonic speeds

3 次聲數(shù)據(jù)分析

經(jīng)過對信號時域波形與聲譜圖的分析，采用EMD 做預處理，保留8 階固有模態(tài)函數(shù)(Intrinsic mode functions,IMF)分量，EMD 計算量小，且在無參考信號的情況下能較好地提高F-K Fisher 檢測算法的分辨能力。對處理后數(shù)據(jù)做時域Fisher檢測，得到的結果如圖6所示。

由圖6所示，檢測出A 陣火箭信號持續(xù)時間約為560 s，B陣信號持續(xù)時間約為400 s。A、B兩陣信號計算的Fisher 比在開始的一段時間內(nèi)均有一個下降的階段，結合時域信號波形與聲譜圖中低頻頻帶能量較低的區(qū)域，判斷此段信號中包含火箭運動狀態(tài)的存在差異。信號幅度與能量衰減歸因于漸遠聲源的傳播損失。

圖6 Fisher 檢測結果Fig.6 The results of Fisher detection

3.1 火箭點火起飛

從聲譜圖與時域Fisher 檢測的結果中，均能在目標信號的初始階段分辨出與背景噪聲明確的界限。圖7 是F-K Fisher 檢測計算出此段時間聲源方位角在慢度網(wǎng)格上的結果。圖7(a)中包含56 s 的時間窗估計的聲源方位為107.8°(正北方向為坐標軸0°方向，順時針計算)，由于發(fā)射點位于觀測點北偏東108°，與計算所得的方位角能夠很好的對應；圖7(b)中包含795 s 的時間窗估計的聲源方位為174.8°，與發(fā)射點位置(北偏東174.7°)相符，兩者均指向發(fā)射點。

圖7 F-K Fisher 檢測計算的火箭點火時間段的慢度網(wǎng)格Fig.7 Slowness grids of time period of motor ignition calculated by F-K Fisher detection method

在火箭發(fā)射起始，聲源距地面較近，對于相同能量的聲源，采集信號的能量由于地面反射等現(xiàn)象會高于高空中無反射的情況，圖6(a)中包含56 s 的時間窗計算的Fisher 值較大，判斷是火箭對地面噴射的過程。隨后的一段時間為火箭遠離地面的升空過程，此時火箭排氣尾焰逐漸遠離地面，兩者的作用減小，產(chǎn)生次聲信號的較弱，計算的Fisher 比值較小。

3.2 連續(xù)運動次聲源分析

圖8 是F-K Fisher 檢測計算的方位角與視速度，圖中檢測到信號后50 s 的時間窗計算的方位角沒有明顯的變化，且指向發(fā)射場，判斷此段時間內(nèi)火箭處于點火爬升的階段，其水平位移較小。

圖8(a)視速度在檢測信號開始階段為367 m/s，假設當?shù)氐慕?jīng)驗聲速為340 m/s，可計算出采集信號的俯仰角為22.11°，之后視速度有一個先變大后變小的過程，對應的俯仰角也是先變大后變小，火箭水平方向加速度逐漸變大，且快于豎直方向加速度的變化。在50～130 s 內(nèi)，視速度緩慢增加，聲源高度逐漸上升，而此時A、B 兩陣檢測方位角變化較小，結合火箭加速過程與Fisher 比的變化，判斷此時對應火箭加速至超聲速，為火箭助推飛行的階段，進一步驗證了126 s的信號為聲爆事件。

圖8 方位角與視速度Fig.8 Azimuth and apparent velocity

在火箭信號持續(xù)的時間里，檢測出的方位角存在著連續(xù)的變化趨勢。A 陣所測得的聲源的方位角由初始方向逐漸向著90°的方向偏移，B 陣所測得的方位角由174.8°不斷減小，兩者結合，反映了火箭向東飛行的運動狀態(tài)。B 陣視速度在一段時間內(nèi)緩慢上升也能反映出火箭升高的過程。

3.3 聲爆前后的信號

從信號的聲譜圖與Fisher 檢測的結果可以得出，在此次監(jiān)測實驗采集的連續(xù)次聲信號中觀察到了聲爆信號。圖5 中，點火起飛信號與聲爆信號之間有一段頻率成分較高且能量較低的信號，判斷為火箭排氣尾焰與大氣作用產(chǎn)生的信號。此時火箭飛行速度低于聲速，雖然在此高度環(huán)境氣壓較大，氣體密度較大，但是殼體與大氣的作用較弱，產(chǎn)生次聲的能量較弱，而未能被監(jiān)測陣采集到。

聲爆信號之后，如圖9所示，上方的紅色曲線隨著時間而降低，反映出隨著距離增加，高頻信號衰減更快，符合次聲傳播規(guī)律；下方的箭頭指聲爆信號預計到達A 陣的時間，此時在低頻頻帶內(nèi)有較大的能量，之后信號成分中有了穩(wěn)定且持續(xù)的低頻能量。此階段火箭持續(xù)加速飛行，殼體與大氣的作用逐漸增強，根據(jù)頻譜特征可以判斷此過程產(chǎn)生信號的低頻成分能量較大。

圖9 A 陣經(jīng)過EMD 處理的信號的聲譜圖Fig.9 The spectrogram of a EMD processed signal from Arrays A

3.4 定位結果

圖10 顯示了兩個次聲接收陣交叉定位的結果，此圖采用通用橫軸墨卡托投影(Universal transverse Mercator,UTM)將大地坐標系轉換為笛卡爾坐標系。結合聲速分析結果，次聲接收陣采集的信號為300 km以內(nèi)的火箭發(fā)射信號，因此在分析時忽略了地球弧度的影響。圖中定位的3 個點分別對應著點火、聲爆信號和采集信號結束時刻，黑線是用已公開資料中此次火箭發(fā)射的相關數(shù)據(jù)做出的實際投影軌跡。圖10 中兩陣對火箭發(fā)射初始階段3個關鍵節(jié)點的定位與實測軌跡基本相符，但由于在方位角的計算過程中受到慢度網(wǎng)格、處理窗長以及次聲在大氣中傳播特性等因素的影響，存在一定的誤差。本文采用的是320×320 的慢度網(wǎng)格做F-K 分析，根據(jù)Fisher檢測法計算的視速度的大小，以及對應的方位角，求出兩個次聲接收陣估計方位角的誤差。A 陣的計算誤差為±0.49°，B 陣的計算誤差為±0.40°。

圖10 定位結果Fig.10 Positioning results

4 結論

本文結合Fisher 檢測算法對火箭發(fā)射產(chǎn)生的次聲信號進行監(jiān)測與分析，從采集信號的聲譜圖中，觀察到火箭點火和聲爆信號，并發(fā)現(xiàn)兩事件間的信號分布在低頻頻帶能量較弱的現(xiàn)象。

根據(jù)單個次聲接收陣所采集的信號的特征，結合檢測算法，得到運行中的火箭的方位角與視速度信息，驗證了火箭點火與聲爆信號的預計到達時間與采集到的時間相匹配，并且判斷了火箭爬升過程以及后續(xù)的飛行方向。結合火箭速度變化的特征，給出了低頻頻帶能量較弱的現(xiàn)象的解釋，火箭超聲速飛行時，殼體與大氣的作用產(chǎn)生次聲波信號的低頻成分能量較大。也可以根據(jù)此現(xiàn)象估計聲爆事件出現(xiàn)的時間段，有助于進一步了解火箭發(fā)射次聲信號與估計火箭的運行狀態(tài)。使用兩個監(jiān)測陣的方位角信息解出火箭發(fā)射初始階段3個關鍵節(jié)點的位置，與公布的火箭飛行軌跡基本相符。

對此種次聲信號的監(jiān)測有助于研究火箭發(fā)射階段的關鍵事件及其信號特征，對某些敏感地區(qū)發(fā)射場的次聲全天候偵察有實際意義。