趙 罡 ,孫乃葳 ,申 珅 ,楊益新
(1.西北工業(yè)大學 航海學院,陜西 西安,710072;2.中國船舶集團有限公司 第705 研究所,陜西 西安,710077)
迄今為止聲脈沖仍是水下通信及目標探測的主要手段。在水下傳播過程中,受信道多徑及目標閃爍等多種因素的作用,聲脈沖信號會發(fā)生衰落和起伏等變化,這一變化對工作在淺海水域的魚/水雷等主戰(zhàn)武器,以及無人水下航行器(unmanned undersea vehicle,UUV)、微型潛艇、蛙人探測聲吶等水下特種裝備的作用尤其顯著,對水下通信質量及目標檢測有效性帶來較大影響[1-3]。為消除這一影響,人們在針對水聲衰落信號的自適應處理方面開展了較多工作[4-5],但對于提供驗證環(huán)境相關的仿真方法和模型研究卻所做不多。實際水聲信道衰落情況復雜多變[6],外場試驗本身又存在成本高昂難以開展的問題,在包括魚/水雷、UUV及特種作戰(zhàn)裝備在內的聲探測武器研發(fā)過程中,使用外場試驗的方式對信號處理方法進行改進和驗證顯然是不現(xiàn)實的??紤]到信噪比(signal noise ratio,SNR)在目標檢測中的重要性[7],建立與實際情況有較高符合度的水聲雙程信道目標回波時域起伏仿真模型,并提供一種可有效利用實航數據的模型參數反演方法,是水下目標精細回波仿真所必須考慮的問題,在水下聲探測武器的研制過程中具有重要的意義。
首先以經典水聲單程信道模型為基礎進行雙程信道模型的建立: 考慮到一定時空不變條件下信道的對稱性將去程信道與返程信道特性等同,然后在直達和多徑信號強度2 個方向,通過積分計算二者聯(lián)合的全概率表達式,最后建立收發(fā)同置應用背景下的水聲往返雙程信道模型。
水聲信道的信號傳輸主要受界面、介質和障礙物的反射、衍射等作用,其信道物理過程復雜,難以解析描述,常利用經驗數據擬合概率分布函數來表示信號幅度變化等重要特征,針對水聲環(huán)境特點,可通過直達信號和多徑信號的加性合成來表示,即
式中:z(t)為直達信號;d(t)為多徑信號。
魚/水雷、UUV 及特種作戰(zhàn)裝備常用高頻短脈沖作為目標聲源探測信號,其頻率高且傳播距離較近,因此開展信道物理過程復雜性建模時可做如下合理假設: 1)信道的直射分量受水體中障礙物遮擋的影響遠小于界面反射的多徑效應;2)傳播衰減造成的短脈沖內慢起伏可忽略。此時可認為z(t)保持包絡不變、而d(t)受界面及溫度梯度的影響變化明顯,則r(t)符合Ricean 分布且滿足條件概率分布
式中: 忽略時間參量t;表示d(t)平均散射功率;I0(·)表示第1 類零階修正貝塞爾函數。
雙程信道模型的建立需聯(lián)合考慮信號去返2 次行程,即分別設定r1(t)、r2(t)為去程信號和返程信號,則二者對應的直達和多徑成分滿足
r2(t)是由水下目標體對去程信號r1(t)反射產生的,故其直達成分z2(t)的包絡與目標反射特性相似呈對數正態(tài)分布,可用r1(t)的條件概率表達,即
其中,表示去程信號的平均散射功率。
文中所用實航數據以聲周期方式保存,相對于聲在水中傳播時1 個聲周期尋的短脈沖時間片段內可認為目標回波過程是確定的,其強度并不發(fā)生閃爍,此時雙程水聲信道模型可簡化為
水聲雙程信道參數的反演問題中,非啟發(fā)式的非線性反演方法如遺傳算法(genetic algorithm,GA)、粒子群算法(particle swarm optimization,PSO)等受算法本身設計限制只能以概率密度函數(probability density function,PDF)差異作為代價函數開展反演[10],其對實航數據時域真值信息利用不足。為有效利用實航數據時域采樣先驗結果[11],流程中以馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法為主要方法實現(xiàn)水聲雙程信道模型抽樣,并進一步結合Bayes 方法的最大后驗估計(maximum a posteriori,MAP)完成參數反演過程,而反演參數的有效性判斷主要通過全局誤差評價模型對比分析海試數據、反演結果分布特性差異來進行。水聲雙程信道反演算法基本流程如圖1所示。
圖1 水聲雙程信道反演算法基本流程Fig.1 Basic process of underwater acoustic go-back channel inversion algorithm
一般而言,Bayes-MCMC 反演的基本思路是通過采樣構造一條在給定真值樣本條件下穩(wěn)態(tài)分布服從所需模型參數向量后驗概率密度的馬爾科夫鏈,當數據鏈收斂時使用采樣數據計算所求模型參數集的各矩特征量,數據鏈同時作為以反演結果為參數的正演模型數據樣本[12]。MCMC 方法的核心在于構造轉移核接收函數,其設計原則與采樣器直接相關,最為通用的采樣器包括Metropolis-Hastings(M-H)、吉布斯采樣(Gibbs sampling,GS)方法等。GS 方法作為M-H 方法的特例是對參數集逐一完成單維度正交采樣,其效率雖然較高但要求準確的參數條件建議分布,而水聲反演問題往往準確的模型參數向量先驗條件概率密度難于確定;M-H 方法使用任意建議分布在整個參數集空間整體采樣以生成數據樣本,且建議分布取對稱分布時可進一步簡化轉移核接受函數形式,具有較好的實現(xiàn)性,因此文中采用M-H 采樣方法[13]。
由于雙程信道實航數據分布特性已知的情況與通常MCMC 采樣先驗概率未知的假設條件不同[14]: 當先驗概率未知時,MCMC 采樣一般是獨立并發(fā)抽取2 組隨機序列,并以2 組序列收斂情況一致性作為停止采樣的判據[15],或者通過對數據采樣過程的設計達到降低計算量提高反演效率的目的[16-17];而雙程信道實航數據分布特性已知時,可據此對每次抽取結果的有效性進行判斷,并使用抽樣分布與實航數據分布的絕對差異或相似程度作為采樣狀態(tài)改變的判定條件。結合圖1 反演算法特征和M-H 采樣存在參數多維耦合條件下命中效率低的問題,為保證數據的高效抽樣,文中基于MCMC 進行了優(yōu)化,提出一種M-H 自適應單維度串行采樣算法,即通過M-H 方法進行精度自適應調節(jié)的單維度串行采樣過程。
首先給出如下定義: 定義模型參數向量m={m1,m2,···,mn}的取值空間為 ?;定義m中mi的補集為?mmi={m1,m2,···,mi-1,mi+1,···,mn}。
具體反演算法步驟如下:
1)設定目標反演偏差 εT、反演偏差增量 Δε、最大反演次數Rm,初始化反演偏差ε0=∞、反演次數Rc=0、反演參數序號i,j=1;
2)在 ?上隨機抽取或采用全局尋優(yōu)方法生成模型參數向量m={m1,m2,···,mn};
4)對每一個被接受的模型參數ms,使用計算模型正演PDF;
5)對模型正演PDF 及真值PDF 進行同步離散采樣并計算其面積誤差ε,當不滿足誤差收斂條件即ε >ε0時回到步驟3);否則接受,同時令i=i+1 及ε0=ε;
6)判斷 εT與 ε0關系: 若ε0<εT,即目標反演偏差已經達到,則令m=并停止反演;若εT<ε0,則令m=,同時判斷若i<n,即尚有模型參數分量未執(zhí)行M-H 采樣,則回到步驟3);否則反演次數Rc加1,同時令i=1;
7)當Rc>Rm時即達到單輪最大反演次數時,為避免陷入局部優(yōu)化令εT=εT+Δε進行反演精度調整,并將反演次數計數置0 后回到步驟3);否則直接回到步驟3)。
算法步驟1)中的各初始設定參數的選擇可根據反演時預期的結果精度設定,但精度要求越高、反演次數越多,雖然會得到更為精確的反演結果,也會帶來更大的計算壓力,對算法執(zhí)行效率產生影響;算法步驟3)中的對稱分布僅僅是為滿足簡化M-H 采樣接受條件的需求,常用如均勻分布、高斯分布等都可使用??偠灾?基于MCMC 方法的M-H 自適應單維度串行采樣算法可綜合考慮精度和效率因素,通過初始參數設置靈活控制參數采樣過程,且在算法終止條件上通過對全局優(yōu)化代價函數結果的判斷以達到有效收斂的目的,可以有效避免因初始參數設置不當而引起的反演結果偏差過大的問題。
為驗證反演算法的有效性,首先對經典分布模型參數進行反演驗證。Ricean 分布是描述海洋信道時常用的統(tǒng)計模型,當使用直達信號與多徑信號功率比即K因子進行描述時其表達式為
此時反演模型流程步驟中參數MAP 估計的主要任務是完成參數和K的估計。記真值樣本為Si,構建Ricean 模型參數反演轉移核接受函數為
式中,Ri(m1)及Ri(m2)為符合式(7)在設定模型參數條件下的離散化采樣結果,其中m1和m2分別代表迭代采樣過程中前后2 次由σ和K構成的參數向量。
由于各采樣獨立,v為各維數據誤差協(xié)方差矩陣Cd特征值積乘,通??墒褂肎A、模擬退火(simulated annealing,SA)、PSO 等全局尋優(yōu)方法獲取極大似然估計。考慮到文中方法Cd矩陣是由真值和采樣數據PDF 函數值殘差構成,同時方法本身存在參數收斂性判斷,且v值與數據矩陣相似性反向相關在反演過程中逐漸趨小,因此可在收斂性判斷同時進行v的迭代更新,而其初值只需在0~1 范圍任意取值即可,不失一般性可將v設定為1。
常規(guī)MCMC 反演一般采用針對2 組獨立并行采樣序列的非參數檢驗方法,如Kolmogorov-Smirnov(K-S)檢驗等的結果作為收斂判斷和算法終止條件,但在文中方法中MCMC 作為參數采樣過程主要用于提高采樣效率,而對反演過程的終止條件判斷可針對整體反演結果進行評價?;诖?設計針對模型參數反演結果的全局誤差模型如下
式中:p、q分別代表反演及真值樣本數據PDF,(p||q)代表2 個PDF 全局差異程度,當其小于εT時,認為反演目標達成,終止采樣,進行上述判斷的參數區(qū)間取Z2即二西格瑪水平(2-Sigma Level)。
全局誤差模型可對仿真采樣結果進行定量評價,但對于多組采樣及反演結果還需要有相似性的定性評價。在機器學習中,KL(Kullback-leibler)散度常用來對2 個概率分布的相似性進行評價,其定義為
KL 散度越小,代表2 個分布越相似,2 個相同分布的KL 散度值為0。文中將直接使用KL 散度指標對多組反演全局誤差模型計算結果進行定性評價。
3.3.1 信道模型真值樣本生成及評價
在水深一定的條件下,對于幾公里內較近的傳輸距離,由于直達信道的存在以及多途情況不那么嚴重,以Ricean 描述的水聲信道表現(xiàn)為較小的σ值和較大的K值。取一組典型值采用MCMC 采樣方法生成符合由該組參數定義的Ricean 分布仿真樣本如圖2(a)所示,縱坐標r為采樣信號幅值;其采樣結果及其與真值PDF 和累積分布函數(cumulative distribution function,CDF)對比如圖2(b)和圖2(c)所示,其中縱坐標p為某幅值信號出現(xiàn)概率。
圖2 MCMC 采樣結果及其與真值分布對比Fig.2 MCMC Sampling results and comparison with their true distributions
仿真結果分析表明,與真值模型設定參數相比,各組采樣均可較快收斂。σ2參數仿真結果誤差0.4%,K參數仿真結果誤差0.5%,PDF 全局誤差1.2%。仿真采樣生成的Ricean 分布樣本與設定模型具有較高一致性,可以作為參數反演真值。
3.3.2 信道模型參數反演
使用所設計參數抽樣及結果評價模型,以上節(jié)所生成數據樣本作為真值對Ricean 分布模型參數開展反演,為比對反演結果將上節(jié)給定參數PDF作為擬合結果。反演設定條件為: 目標反演偏差εT=1%、反演偏差增量Δε=1%、最大反演次數Rm=5;每趟次參數采樣數5 000 個。
參數迭代反演過程以及結果對應的Ricean 分布PDF/CDF 與真值分布對比如圖3 所示。
圖3 模型參數采樣過程及反演結果與真值分布對比Fig.3 Sampling process of model parameter and comparison of inversion results with true distributions
反演結果分析表明,與正演模型設定參數相比,參數 σ2反演結果誤差為0.8%,參數K反演結果誤差為1.5%,PDF 全局誤差為2.6%。各參數均有較好的收斂一致性,單維度串行采樣模型具有較高的反演精度和效率,可用于水聲雙程信道模型參數反演。
不同于單程信道測量信號,水下自導武器聲尋的脈沖通常為矩形或梯形包絡的脈沖信號,其長度一般在10 ms 量級,而回波脈沖受傳播、吸收衰減及信道衰落等作用,在接收端信號幅值往往較低并常采用時變或自動增益技術以提升其SNR,這些因素的影響在信道模型參數反演時都需考慮。一般來說,為了得到原始信號數據,需要基于補償曲線消除增益作用影響,同時為了得到較為準確的脈沖起始和終止時刻,副本相關等基本信號處理過程也是必須的。
在一段實際回波信號中選擇出反演所需的脈沖回波后,還需要對信號進行幅值歸一化處理,同時考慮到采樣及計算均在離散點上進行,所選擇PDF 值域及積分計算值域范圍均會對結果帶來一定的影響。即使經過了預處理過程,受水文及界面環(huán)境的影響,弱回波脈沖信號被本底噪聲部分淹沒,仍會對反演結果帶來一定的影響,下面的反演就是針對SNR 條件不同的幾組回波展開的。以某湖上試驗數據為例,采樣時信道傳輸距離為500~2 000 m,2 組不同SNR 條件真值采樣如圖4 所示。
圖4 不同SNR 下脈沖回波時域信號Fig.4 Pulse echo time domain signals under different SNRs
反演過程的M-H 抽樣模型轉移核接受函數形式及參數定義同式(8),不同之處在于符合設定模型參數的離散化采樣結果生成時,采取式(6)所定義不考慮目標閃爍條件的水聲往返雙程信道時域起伏模型。
反演結果評價同樣采用式(9)和式(10)模型開展,但是考慮到真實信號的模糊性,當信號樣本屬于較低SNR 條件(小于3 dB)時反演設定條件及相關評價準則放寬至目標反演偏差εT=0.03,其他不變;同時為了避免模型陷入單個參數采樣的局部優(yōu)化,在采樣算法的步驟5)設定單參數最大無效反演5 輪次后,以新的隨機參數向量重新開始反演。
基于實際數據樣本對無目標閃爍條件的水聲往返雙程信道時域起伏模型參數開展反演。設定σ2隨機采樣區(qū)間[0.1 1],K隨機采樣區(qū)間[0.5 3],反演最終結果如下。
4.3.1 不同SNR 信號反演
表1 不同SNR 信號參數反演結果對比Table 1 Comparison of parameter inversion results under different SNR signals
2 組反演結果對應的Ricean 分布PDF/CDF 與真值分布對比如圖5 所示。
圖5 不同SNR 信號反演結果與真值分布對比Fig.5 Comparison of inversion results with true distributions under different SNR signals
反演結果表明,SNR 水平對于反演結果影響較大: 在相同發(fā)射接收水平下,SNR 越高反映信道質量越好,代表直達信號強度的K因子越大,而表征多途信號功率的 σ2參數越小。此外,在較高SNR 條件下信號起伏分布規(guī)律與正演模型相似度更好,表現(xiàn)為用于評價PDF 整體差異的結果更小,表征分布相似性的KL 散度指標也更小。同時方差分析表明不同信號條件下,各評價參數均有較好的收斂一致性。最后,反演結果與真值的總體差異也表明所建立的水聲雙程信道仿真模型是有效的,尤其是在較高SNR 條件下。
4.3.2 多樣本反演
文中最后選擇了來自于某海域實航試驗的9 組在時空上有一定聯(lián)系的實際信號進行了模型參數反演。其中,前8 組數據信道傳輸距離自1 000~2 000 m 非等間隔分布,SNR 為1~2 dB,第9 組數據作為對照SNR 條件約3 dB,信號時域形態(tài)如圖6 所示。
圖6 多樣本脈沖回波時域信號Fig.6 Multiple groups of pulse time domain echo
反演初始設定及過程參數與4.3.1 節(jié)相同,均進行100 趟反演。1~8 組反演結果對應的Ricean分布PDF/CDF 與真值分布對比如圖7 所示。
1~8 組反演結果如表2 所示。
各參數的變化趨勢如圖8 所示。
作為對照,對第5、9 組數據按同樣條件進行了100 趟反演,參數反演過程結果對比見圖9。
第5、9 組信號參數反演結果的對比如表3所示。
表3 第5、9 組信號參數反演結果對比Table 3 Comparison of inversion results at 5th&9thgroups signals
第5、9 組反演結果對應的Ricean 分布PDF/CDF與真值分布對比如圖10 所示,其中橫縱坐標與前文定義相同。
針對上述反演統(tǒng)計結果的分析表明:
1)從表1 及圖7 結果可知,各組反演過程均成功收斂,反演結果與真值誤差按模型所定義的總體差異均在5%以內;
圖7 1~8 組信號反演結果與真值分布對比Fig.7 Comparison of inversion results with true distributions at 1thto 8thgroups signals
圖8 1~8 組反演參數變化趨勢Fig.8 Change trend of inversion parameters at 1thto 8thgroups signals
表2 1~8 組雙程信道模型參數反演結果Table 2 Parameter inversion results of go-back channel model at 1thto 8thgroups signals
圖9 第5、9 組雙程信道模型參數反演過程結果對比Fig.9 Comparison of parameter inversion results of go-back channel model at 5th&9thgroups signals
圖10 第5、9 組信號反演結果與真值分布對比Fig.10 Comparison of inversion results with true distributions at the 5th&9thgroups signals
針對水下作戰(zhàn)裝備研發(fā)對水聲環(huán)境仿真需求,基于貝葉斯全概率公式構建了水下目標回波雙程衰落信道模型,并結合針對實測時域數據的模型參數反演特點,設計了單維度串行采樣模型。模型通過M-H 核函數設計,將MCMC 方法引入全局反演算法中,以有效利用有限實測數據信息,為水下目標回波仿真模型置信度提升提供了有力而靈活的保證。文中使用仿真雙程信道數據對反演算法的有效性進行了驗證,并在最后使用多組實測數據開展反演,通過對反演結果的分析確認了雙程衰落信道模型的正確性。
值得注意的是,在使用實測數據開展模型參數反演時,除了信號SNR 條件會對反演效率和結果產生較大的影響外,準確的識別接收信號中脈沖回波所在空時區(qū)間是反演有效開展的重要保證,在算法的后續(xù)優(yōu)化中,可針對實測數據的自動選取及預處理開展更多工作,以進一步提升反演速度和結果的準確性。