基于GPU的海底混響信號(hào)快速仿真方法研究

李　利，劉興華

（大連測控技術(shù)研究所，遼寧　大連 116013）

基于GPU的海底混響信號(hào)快速仿真方法研究

李利，劉興華

（大連測控技術(shù)研究所，遼寧大連 116013）

在計(jì)算海底混響信號(hào)時(shí)，根據(jù)混響產(chǎn)生的物理機(jī)理，以射線聲學(xué)為基礎(chǔ)，用 Lambert 散射定律計(jì)算海底反向散射強(qiáng)度，采用單元散射模型建立海底混響信號(hào)模型。利用 GPU 相比于 CPU 具有更高的浮點(diǎn)運(yùn)算能力和內(nèi)存帶寬的特點(diǎn)，采用 GPU 進(jìn)行計(jì)算海底混響信號(hào)。通過對(duì)仿真混響信號(hào)的處理分析，在散射點(diǎn)較少時(shí)，混響信號(hào)包絡(luò)更接近于 K 分布，而隨著散射點(diǎn)的增多，混響信號(hào)的包絡(luò)接近于瑞利分布。符合混響信號(hào)的一般統(tǒng)計(jì)特性。該方法能快速仿真出混響信號(hào)，達(dá)到高效的目的，為以后混響信號(hào)的實(shí)時(shí)演示驗(yàn)證提供一條可供選擇的途徑。

射線聲學(xué)；混響信號(hào)；GPU運(yùn)算

0　引　言

混響是影響主動(dòng)聲吶工作性能的重要干擾因素之一，混響的預(yù)報(bào)對(duì)聲吶的使用和設(shè)計(jì)有著重要意義。對(duì)混響信號(hào)仿真的研究主要包括混響強(qiáng)度仿真［1-3］和混響時(shí)間序列仿真。當(dāng)聲吶系統(tǒng)所采用的信號(hào)處理方法主要依賴能量時(shí)，常用混響強(qiáng)度對(duì)其進(jìn)行估計(jì)和評(píng)價(jià)。隨著聲吶系統(tǒng)的發(fā)展，如多波束和復(fù)雜相干信號(hào)處理技術(shù)的應(yīng)用，就需要對(duì)混響時(shí)間序列進(jìn)行仿真?；祉憰r(shí)間序列仿真的方法主要有 2 種，一種是根據(jù)混響信號(hào)的概率分布，仿真所需分布的混響信號(hào)［4-6］；另一種是以點(diǎn)散射模型為基礎(chǔ)的仿真方法，即通過計(jì)算散射體或通過網(wǎng)格劃分海底為小散射元的散射信號(hào)在接收點(diǎn)的疊加得到混響信號(hào)，很多文獻(xiàn)中混響仿真方法都是以此模型為基礎(chǔ)的改進(jìn)和擴(kuò)展［7-13］，該方法能夠比較準(zhǔn)確地仿真混響信號(hào)，但由于海洋中存在大量散射體，在混響仿真時(shí)存在計(jì)算量大的問題。一般來說，同一時(shí)期的 GPU 其浮點(diǎn)運(yùn)算能力和內(nèi)存帶寬要比 CPU 高一個(gè)數(shù)量級(jí)左右，并且顯卡廠商 NVIDIA 推出的運(yùn)算平臺(tái) CUDA 簡單易用，使得科研工作者能夠把更多的精力放在算法上，將已有算法進(jìn)行簡單的修改就可以得到十分可觀的加速比，因此使用 GPU 并行計(jì)算非常適合科學(xué)計(jì)算。基于這個(gè)問題，采用 GPU 運(yùn)算的方法則能很好地解決由計(jì)算量大引起的運(yùn)算時(shí)間長的缺點(diǎn)，這樣就能快速準(zhǔn)確地仿真混響信號(hào)。

1　海底混響信號(hào)模型

在淺海近程，由于海底的一次散射對(duì)海底混響的貢獻(xiàn)最大，因此對(duì)海底混響信號(hào)模型的建立，也只考慮海底 1 次散射，共 4 類聲線，聲線如圖1 所示。

圖1　聲線示意圖Fig.1　Schematic plan of sound ray

混響信號(hào)模型不考慮海面的散射；海水均勻、無吸收、恒定聲速；入射聲波及散射聲波按球面波擴(kuò)張。海底水平，海水深度為 H，聲源與接收水聽器合置且聲源到海面的距離為 h，恒定聲速為 c，海面反射系數(shù)為 m。海底散射單元的大小依據(jù)海底散射系數(shù)空間相關(guān)半徑選取，海底的聲壓散射系數(shù)幅值和相位在此空間相關(guān)半徑之內(nèi)為同一數(shù)值，而在此空間相關(guān)半徑之外，則為服從相同分布的其他數(shù)值。

粗糙海底散射強(qiáng)度由 Lambert 定律確定，聲強(qiáng) Ii為入射波以掠射角 θi入射到粗糙單位面元上，各方向上的散射聲強(qiáng)度為：

式中：μ 為比例常數(shù)；θs為散射掠射角。

由于聲強(qiáng)與聲壓幅值的平方成比例，則散射聲壓幅值為：

式中：μ' 為比例常數(shù)；Pi和 Ps分別為入射聲壓幅值和散射聲壓幅值。

根據(jù)混響的統(tǒng)計(jì)特性：混響的瞬時(shí)幅值為正態(tài)分布；瞬時(shí)相位為（0～2π）均勻分布，令海底聲壓散射系數(shù)為：

式中：μ'p服從高斯分布 N（a，σ2）（a 為均值，σ 為方差），φ 服從（0～2π）均勻分布。

1）第 1 類聲線 tk時(shí)刻混響信號(hào)

若發(fā)射信號(hào)為脈寬為 τ 的窄帶脈沖信號(hào) s（t），tk時(shí)刻散射聲的散射區(qū)域?yàn)橐煌膱A環(huán)，將散射區(qū)域依據(jù)海底散射系數(shù)相關(guān)半徑劃分為 N1個(gè)散射元，各散射元散射系數(shù)的幅值相互獨(dú)立，且皆服從相同高斯分布，相位也相互獨(dú)立并服從（0～2π）均勻分布，因此散射信號(hào)在接收點(diǎn)迭加后為：

由于聲源與接收水聽器合置，因此sin（θi1）=sin（θs1）=（H-h）/r1，式（4）可化為：

2）第 2 類聲線 tk時(shí)刻混響信號(hào)

同理，tk時(shí)刻此類聲線的海底散射元為 N2，則第2 類聲線的混響信號(hào)為：

式中：tk=（r1+r2+r3）/c；sin（θi2）=（H-h）/r1；其他參數(shù)同前。

3）第 3 類聲線 tk時(shí)刻混響信號(hào)

同理，tk時(shí)刻此類聲線的海底散射元為N3，則第 2類聲線的混響信號(hào)為：

式中：tk=（r1+r2+r3）/c；sin（θi3）= （H+h）/（r2+r3）；sin（θs3）=（H-h）/r1，其他參數(shù)同前。

4）第 4 類聲線 tk時(shí)刻混響信號(hào)

同理，tk時(shí)刻此類聲線的海底散射元為 N4，則第2 類聲線的混響信號(hào)為：

式中：tk=2（r2+r3）/c；sin（θi4）= sin（θs4）=（H+h）/（r2+r3）；其他參數(shù)同前。

將這四類聲線的混響信號(hào)按時(shí)間累加即得總混響信號(hào)。

2　基于 GPU 的快速計(jì)算方法

CPU 為提高分支指令的處理速度，其很多部件都用于做分支預(yù)測，以及在分支預(yù)測錯(cuò)誤的時(shí)候修正和恢復(fù)算術(shù)邏輯單元的結(jié)果，這將大大增加器件的復(fù)雜度，使其更側(cè)重于靈活高效的處理速度而非計(jì)算能力。相比之下，GPU 在設(shè)計(jì)之初就定位于大量數(shù)據(jù)的并行計(jì)算，發(fā)展到現(xiàn)在已經(jīng)具有非常強(qiáng)大的計(jì)算能力，其核心如圖2 所示。以市面上最強(qiáng)桌面 CPU—Intel Core i7 5960 X 來說，它具有 8 個(gè)核心 16 個(gè)線程，默認(rèn)頻率 3.3 GHz，配合 AVX2 指令集其運(yùn)算能力可以達(dá)到 422 Gflop/s，其內(nèi)存帶寬大約為 68 GB/s；而同一時(shí)期的最強(qiáng)桌面 GPU—GTX TITAN X 具有 3 072 個(gè)CUDA 核心，默認(rèn)頻率 1 006 MHz，配合 FMA 指令集其運(yùn)算能力可以達(dá)到 6 181 Gflop/s，其內(nèi)存帶寬大約為336 GB/s，并且 1 臺(tái)計(jì)算機(jī)中可以搭配多個(gè) GPU，由此可見 GPU 并行計(jì)算的巨大優(yōu)勢。

CUDA（Compute Unified Device Architecture）是顯卡廠商 NVIDIA 推出的運(yùn)算平臺(tái)，在 CUDA 的架構(gòu)下，程序分為 host 端和 device 端 2 個(gè)部分。Host 端是指在 CPU 上執(zhí)行的部分，而 device 端則是在 GPU 上執(zhí)行的部分。Device 端的程序又稱為“kernel”。通常host 端程序會(huì)將數(shù)據(jù)準(zhǔn)備好后，復(fù)制到顯卡的內(nèi)存中，再由顯示芯片執(zhí)行 device 端程序，完成后再由host 端程序?qū)⒔Y(jié)果從顯卡的內(nèi)存中取回。

圖2　CPU 與 GPU 運(yùn)算單元示意圖Fig.2　Schematic plan of CPU and GPU arithmetic element

基于此原理，將仿真混響信號(hào)的程序進(jìn)行 GPU 編譯，流程如圖3 所示。

同樣仿真 1.1 s 的海底混響信號(hào)，用 GPU 進(jìn)行運(yùn)算與用 CPU 進(jìn)行運(yùn)算的時(shí)間對(duì)比如表1 所示。由表1 可看出，運(yùn)用 GPU 進(jìn)行運(yùn)算能極大幅度提高運(yùn)算速度，減少運(yùn)算時(shí)間，這為以后混響信號(hào)的實(shí)時(shí)演示驗(yàn)證提供一條途徑。

3　仿真混響信號(hào)分析

當(dāng)發(fā)射信號(hào)為 20 kHz 的 CW 信號(hào)，幅值為 1 V，脈寬 τ=5 ms，水中聲速 c=1 500 ms，海深 50 m，海面反射系數(shù) m=0.9，收發(fā)換能器距水面 5 m 時(shí)，仿真得到如圖4 所示的混響信號(hào)波形。

圖3　海底混響信號(hào) GPU 運(yùn)算流程圖Fig.3　Flow chart of reverberation signal calculation by using GPU

表1　CPU 與 GPU 運(yùn)算時(shí)間對(duì)比Tab.1　The comparison of operation time between Cpu and Gpu

圖4　仿真混響信號(hào)波形圖Fig.4　The simulation oscillogram of reverberation signal

取 0.15 s 處的一段混響信號(hào)進(jìn)行頻數(shù)統(tǒng)計(jì)，并用直方圖表示，這種直方圖可以估計(jì)總體的概率密度，而為便于觀察，可提取直方圖高度是頻率組距的情況的包絡(luò)，用曲線形式近似反映樣本總體的概率密度。用此方法對(duì)此段仿真信號(hào)得到的樣本分布形式進(jìn)行估計(jì)，并分別對(duì)用 K 分布和瑞利分布曲線擬合與仿真樣本數(shù)據(jù)統(tǒng)計(jì)得到的結(jié)果進(jìn)行比較，如圖5 所示。

同樣的方法對(duì) 0.5 s 處的一段混響信號(hào)進(jìn)行處理，得到的結(jié)果如圖6 所示。

由圖5 可知，當(dāng)在近程時(shí)，各個(gè)散射單元散射聲信號(hào)疊加后形成的混響信號(hào)的包絡(luò)趨于 K 分布，與傳統(tǒng)理論中混響信號(hào)的包絡(luò)呈瑞利分布不相符，這是由于傳統(tǒng)理論中，只有當(dāng)散射波的數(shù)量足夠多時(shí)，其混響信號(hào)的包絡(luò)才會(huì)服從瑞利分布，而在近程，散射單元的個(gè)數(shù)較少，不滿足大數(shù)定理，因此會(huì)有此現(xiàn)象產(chǎn)生。而圖6 則說明在相對(duì)遠(yuǎn)程處，散射單元的個(gè)數(shù)增多，滿足大數(shù)定理，則此時(shí)混響信號(hào)的包絡(luò)就會(huì)趨于瑞利分布而非呈 K 分布。進(jìn)而驗(yàn)證混響信號(hào)的合理性。

圖5　混響信號(hào)包絡(luò)統(tǒng)計(jì)概率Fig.5　The statistical probability of reverberation signal envelope

圖6　混響信號(hào)包絡(luò)統(tǒng)計(jì)概率Fig.6　The statistical probability of reverberation signal envelope

4　結(jié)　語

根據(jù)混響產(chǎn)生的物理機(jī)理，以射線聲學(xué)為基礎(chǔ)，用 Lambert 散射定律計(jì)算海底反向散射強(qiáng)度，采用單元散射模型建立海底混響信號(hào)模型。通過對(duì)仿真混響信號(hào)的處理分析，在散射點(diǎn)較少時(shí)，混響信號(hào)包絡(luò)更接近于 K 分布，而隨著散射點(diǎn)的增多，混響信號(hào)的包絡(luò)接近于瑞利分布，符合混響信號(hào)的一般統(tǒng)計(jì)特性，驗(yàn)證了混響信號(hào)模型的合理性。由于 GPU 相比于 CPU具有更高的浮點(diǎn)運(yùn)算能力和內(nèi)存帶寬，因此采用 GPU進(jìn)行計(jì)算，能快速仿真出混響信號(hào)，達(dá)到高效的目的，為今后混響信號(hào)的實(shí)時(shí)演示驗(yàn)證提供一條可供選擇的途徑。

［1］姚萬軍，蔡志明.淺海混響的建模與預(yù)報(bào)［J］.海軍工程大學(xué)學(xué)報(bào)，2009，21（2）: 88-91，112.

［2］袁駿，肖卉，張明敏，等.雙基地聲吶海底混響面元?jiǎng)澐峙c建模［J］.探測與控制學(xué)報(bào)，2009，31（4）: 1-4.

［3］張明輝.三維環(huán)境海洋混響強(qiáng)度衰減規(guī)律研究［D］.哈爾濱:哈爾濱工程大學(xué)，2005: 5-7.

［4］方世良.海洋混響信號(hào)的序貫仿真［J］.聲學(xué)技術(shù)，1996，15（3）: 101-104.

［5］ABRAHAM D A，LYONS A P.Reverberation envelope statistics and their dependence on sonar bandwidth and scattering patch size［J］.IEEE journal of Oceanic Engineering，2004，29（1）:126-137.

［6］ABRAHAM D A，LYONS A P.Simulation of non-rayleigh reverberation and clutter［J］.IEEE Journal of Oceanic Engineering，2004，29（2）: 347-362.

［7］蔡平，梁國龍，葛鳳翔，等.界面混響信號(hào)的仿真研究［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2000，21（4）: 31-35.

［8］陳文劍，孫輝，朱建軍，等.基于分?jǐn)?shù)階傅里葉變換混響抑制的目標(biāo)回波檢測方法［J］.聲學(xué)學(xué)報(bào)，2009，34（5）: 408-415.

［9］郭熙業(yè)，蘇紹璟，王躍科.海底混響統(tǒng)計(jì)建模與仿真方法研究［J］.兵工學(xué)報(bào)，2009，30（7）: 940-944.

［10］郭熙業(yè)，蘇紹璟，王躍科.基于射線理論的海底混響建模研究［J］.聲學(xué)技術(shù)，2009，28（3）: 203-207.

［11］姚萬軍，蔡志明.基于簡正波的淺?；祉懶蛄蟹抡妫跩］.聲學(xué)技術(shù)，2009，28（1）: 25-28.

［12］姚萬軍，蔡志明，衛(wèi)紅凱.淺?；祉懡５穆暿櫪碚摚跩］.聲學(xué)學(xué)報(bào)，2009，34（3）: 223-228.

［13］ZHANG M H，SUN H，CHEN W J.Simulation model of bottom reverberation signals for horizontal bistatic receiving array［C］//Proceedings of 2008 IEEE International ultrasonics symposium.Beijing: IEEE，2008: 1437-1440.

［14］ABRAHAM D A，LYONS A P.Exponential scattering and K-distributed reverberation［C］//Proceedings of 2001 MTS/IEEE Conference and Exhibition.Honolulu，HI: IEEE，2001，3: 1622-1628.

Research on the reverberation signal fast simulation method based on GPU

LI Li，LIU Xing-hua
（DaLian Scientific Test and Control Technology Institute，Dalian 116013，China）

According to the physical process of reverberation and the theory of ray acoustics，calculation of sea bottom back scattering intensity with Lambert scattering law is carried out and a model of reverberation signal model is set up by using the cell scattering model.Because GPU has a higher floating-point computing performance and memory bandwidth compared to CPU，fast calculation of the reverberation signal simulate is achieved by using GPU.The results show that the envelope of the reverberation signal is more close to the K distribution with few scatter points and it is close to the Rayleigh distribution with the increase of scattered points which is accord with general statistical characteristics of reverberation signal.

ray acoustics；reverberation signals；GPU

TB566

A

1672-7619（2016）06-0128-04

10.3404/j.issn.1672-7619.2016.06.026

2015-12-14；

2016-02-29

李利（1985-），女，碩士，工程師，主要從事艦船噪聲測量與分析及水聲信號(hào)處理方面的工作。