水下三維聲場仿真與可視化方法研究

孫雪海，笪良龍，李玉陽

(海軍潛艇學院 海軍水下作戰環境研究所，山東 青島266071)

0 引 言

海洋水聲環境是海軍作戰的主要環境，其涵蓋了大量抽象、復雜的信息，如何快速地從如此龐大繁雜的數據中獲取對作戰有用的信息，并以直觀、易于理解的方式呈現，對幫助作戰人員感知環境規律和做出正確的指揮決策具有至關重要的作用。圖形、圖像等視覺信息具有很強的直觀性，有利于加深人們對事物、規律的感知和理解。聲作為目前用于水下目標探測的最有效手段，水下聲場仿真及其可視化是復雜水聲環境研究的關鍵技術。

水下三維聲場數據是通過三維聲傳播模型計算所得。但三維聲傳播模型的計算量大、計算時間長，給復雜水聲戰場環境中的快速聲場分析帶來了嚴峻的挑戰。隨著高性能計算技術的發展，基于集群系統的并行計算為三維聲場的快速計算提供了可能。文獻［1］研究了三維射線－簡正波－拋物方程模型在高性能計算平臺的并行計算問題。文獻［2］實現了水平不變海洋聲道中WKBZ 簡正波方法的并行計算。

水下三維聲場可視化對幫助作戰人員感知和理解水下聲傳播規律具有重要作用。文獻［3 －4］研究了基于預積分體渲染技術的水聲環境體可視化方法;文獻［5 －6］利用基于三維紋理映射直接體繪制算法實現了水下三維聲場的體可視化。

本文在前人研究的基礎上，主要研究三維拋物方程模型N×2D 弱三維近似方法的并行計算策略及高質量水下三維聲場的實時可視化問題。

1 三維聲場計算

目前相對成熟的三維聲傳播模型主要有三維射線模型 (HARPO)、三維耦合簡正波模型(CMM3D)和三維拋物方程模型(FOR3D)等［5］。拋物方程方法較射線方法計算速度慢，但在計算精度上高于射線方法;在保證相同精度的前提下，拋物方程方法的計算速度較簡正波方法快得多。在求解低頻聲波隨距離變化的聲傳播問題上拋物方程方法具有獨特的優越性。

在柱坐標系(r，φ，z)下，簡諧點源在水平變化聲道中聲場的三維Helmholtz 方程表示為［7］:

考慮不同方位φ 之間耦合關系后的FOR3D 模型，求解過程復雜、計算量大，不利于進行并行計算。但在一般情況下可以忽略方位之間的耦合，采用N ×2D 弱三維近似計算三維聲場。此時，式(1)可簡化為標準的二維Helmholtz 方程:

式中:p(r，z)為聲壓;k0=ω/c 為參考波數;n = c0/c為介質折射率。

采用Tappert 方法，令:

其中包絡函數ψ(r，z)隨距離緩變，Hankel 函數滿足Bessel 微分方程:

當k0r ≥1 時，Hankel 函數可用其漸近形式代替:

利用方程(4)表示的Hankel 函數的性質，在遠場假設(k0r ≥1 )條件下，將式(3)和式(5)代入式(2)中，即可得到如下簡化的橢圓形波動方程:

引入“近軸近似”，小角度近似可表示為:

將近軸近似用于方程(6)，即可得到標準拋物方程:

采用Hardin 和Tappert 提出的分裂－ 步進FFT算法求解方程(8)，可得到某一方位上的二維聲場數據，通過計算N 個方位上的二維聲場即可近似地組成整個三維空間的聲場。圖1 為對三維聲場進行離散采樣的示意圖，圖中Δr 為距離分辨率，Δφ 為方位分辨率，Δz 為深度分辨率。

圖1 三維聲場離散采樣示意圖Fig.1 Sketch of discrete sampling of 3D acoustics energy field

采用N × 2D 弱三維計算三維聲場可以進行2層并行。第1 層是不同方位間二維聲場的并行計算;第2 層是求解方程內部算法的并行，如FFT算法的并行計算。利用曙光TC4000L 高性能計算平臺，結合平臺體系結構的特點，采用兩級并行策略:將方位扇面分配到不同節點，每個方位扇面內單獨執行與其相關的所有運算，計算過程中進程間沒有數據交換，直至計算結束將結果合并輸出，采用MPI 模型實現;FFT 算法的并行在節點內進行，計算結果存儲在共享內存中，采用OpenMP 模型實現。

采用N ×2D 方法計算得到N 個采樣方位切片上聲傳播損失隨r － z 的變化，將合并輸出的結果采用笛卡爾坐標的形式進行組織，不同方位上的聲場數據依次存儲到三維紋理中。

數據組織形式如圖2 所示，圖中紋理s 坐標表示水平距離r，紋理t 坐標表示深度z，紋理r 坐標表示方位φ。

圖2 原始聲場數據的三維紋理存儲示意圖Fig.2 Sketch of storage for 3D texture of original acoustics energy field data

采用三維紋理存儲聲場數據，可以通過設置三維紋理的參數，利用三次線性插值實現對任意方位上聲場數據的重采樣。由于使用了硬件加速，三次線性插值所帶來的額外開銷幾乎可以忽略不計。

2 三維聲場可視化

三維數據場體繪制的算法有很多種，其中基于GPU 的紋理映射體繪制算法利用了圖形流水線的并行處理功能，是眾多體繪制算法中速度最快的。基于紋理映射的體繪制算法主要包括以下2個階段(見圖3):第1 階段是紋理生成，將體數據載入系統內存，并將數據轉換為紋理數據格式存入顯存中;第2 階段是紋理繪制，根據一系列采樣多邊形對紋理數據進行重采樣并混合生成最終的圖像。

原始聲場數據的三維紋理中依次存儲了不同方位的聲傳播損失數據，是按照笛卡爾坐標系組織的，但原始聲場數據是在柱坐標系下計算得到的，因此，在傳統基于紋理映射體繪制算法的紋理繪制階段，直接對原始聲場數據的三維紋理進行重采樣無法還原真實的三維聲場。針對這一問題，本文在紋理繪制階段引入了渲染到紋理的步驟，基本思路是:首先定義與幀緩沖區對象(FBO，Frame Buffer Object)綁定的空的三維紋理。然后對原始聲場數據的三維紋理進行重采樣和坐標轉換，還原成真實三維聲場的結構，并將結果渲染到與FBO 綁定的三維紋理上。最后再對與FBO 綁定的三維紋理進行顏色映射和光照計算等操作。最終圖像生成階段用于進行重采樣的三維紋理是與FBO 綁定的三維紋理。改進后基于三維紋理映射體繪制的算法流程如圖4 所示。

圖3 基于紋理映射體繪制算法流程Fig.3 Flow of volume rendering algorithm based on texture mapping

圖4 基于三維紋理映射體繪制改進算法流程Fig.4 Flow of improved volume rendering algorithm based on 3D texture mapping

2.1 數據重采樣

對原始聲場數據進行重采樣及坐標轉換的具體步驟如下:

1)創建與FBO 綁定的三維紋理。為還原柱坐標系下三維聲場，與FBO 綁定的三維紋理的大小定義為:

式中:XFBO，YFBO，ZFBO分別為與FBO 綁定的三維紋理的長、寬、高;XOrig，YOrig分別為存儲原始聲場數據的三維紋理的長、寬。

2)設置視口及定義投影方式。視口大小設置如下:

glViewport (0，0，XFBO，YFBO);

定義二維正視投影矩陣，圓心位于幾何中心:

gluOrtho2D (－1，1，－1，1);

3)對存儲原始聲場數據的三維紋理進行重采樣。使用OpenGL 著色語言，在片元著色器中實現對原始聲場數據的重采樣及坐標轉換，并將結果渲染到與FBO 綁定的三維紋理上。

2.2 顏色映射

顏色映射的目的是將聲場數據映射為顏色并賦予不同的透明度值，通過傳輸函數實現。在水下目標探測過程中，當目標輻射噪聲的傳播損失大于某一臨界值時(聲吶優質因子，FOM)，將不能對目標實施可靠探測。為了在可視化過程中表達這一現象，采用非線性的顏色、透明度傳輸函數，如圖5所示，當傳播損失小于FOM 值時聲場可見，反之接近透明(顏色與光強值接近0)。

圖5 非線性傳輸函數示意圖Fig.5 Sketch of nonlinear transfer function

2.3 光照計算

為了增強三維聲場繪制的真實感，在片元著色器中引入了光照和明暗計算。采用Phong 局部光照模型，它包括環境光(Iambient)、漫反射(Idiffuse)和鏡面反射(Isepcular)3 個部分［8］。

式中，Ka，Kd，Ks分別為環境光反射系數、漫反射系數和鏡面反射系數;Ia為環境光強度;Ip為照射光源強度;n，l 和v 分別為歸一化的表面法向量、指向光源的向量和指向視的向量;h 為v 和l 的中間向量;n為光照指數。

3 仿真實驗

實驗條件:聲場半徑120 km;計算范圍內最大海深5 000 m;聲源為水面目標，深度5 m;聲源頻率100 Hz。三維聲場離散采樣的距離分辨率500 m，方位分辨率10°，深度分辨率20 m，聲場體數據的大小為241 ×36 ×251。

曙光TC4000L 高性能計算平臺分配36 個節點分別計算不同方位的聲場數據。同時，為了比較驗證并行計算的效率，在高性能計算平臺開辟了一個比較節點，用于串行計算。

圖形顯示計算機硬件性能:CPU:Inter (R)Core (TM)i7 － 4700HQ，8 核;顯卡:NVIDIA GeForce GT 750M，4GB 顯存。圖6 是分別進行5 次并行和串行計算總耗時的結果比較:

水下三維聲場數據在直角坐標空間和球面坐標空間的可視化效果分別如圖7 和圖8 所示。

圖6 并行和串行計算結果比較Fig.6 Results compare in parallel and series computation

圖7 直角坐標空間水下三維聲場體可視化效果圖Fig.7 Effect picture of Underwater 3D acoustics energy field in Cartesian space

圖8 球面坐標空間水面目標三維聲場體可視化效果圖Fig.8 Effect picture of surface target 3D acoustics energy field in spherical space

在圖形顯示計算機中，圖7 和圖8 所示場景的刷新頻率分別大于45 fps 和30 fps。從仿真實驗結果看出:在高性能計算平臺下對三維聲場進行并行計算具有明顯的加速效果，平均加速比約為20;改進后的基于三維紋理體繪制算法實現了水下三維聲場的實時體繪制，能夠生成高質量的水下三維聲場圖像。

4 結 語

本文實現了水下三維聲場的快速計算和實時可視化表達，為作戰人員感知和理解水下聲傳播規律提供了新的方法，拓展了指揮決策的方式。文中三維聲場的并行計算主要依托實驗室的高性能計算平臺，近年來，隨著高性能GPU 的飛速發展，基于CUDA 的并行計算越來越受到人們的重視和青睞，使PC 機進行復雜的科學計算成為可能。因此，下一步將重點研究基于CUDA 的三維聲場并行計算以及相應的數據組織和體繪制方法。

［1］張林，笪良龍，李玉陽．三維水聲傳播建模及并行算法研究［J］．系統仿真學報，2013，24(1):54 －57．ZHANG Lin，DA Liang-long，LI Yu-yang．3D Underwater sound propagation modeling and parallel algorithm［J］．Journal of System Simulation，2013，24(1):54 －57．

［2］笪良龍，范培勤，陳冬濱．水平不變海洋聲道中WKBZ簡正波方法的并行算法［J］．系統仿真學報，2006，18(7):1980 －1982．DA Liang-long，FAN Pei-qin，CHEN Dong-bin．Parallel programme of WKBZ mode approach to sound propagation in horizontally stratified oceans［J］．Journal of System Simulation，2006，18(7):1980 －1982．

［3］李玉陽，笪良龍，宋潔．水聲環境仿真及其可視化技術研究［J］．測繪科學技術學報，2009，26(1):8 －11．LI Yu-yang，DA Liang-long，SONG Jie．Research on underwater acoustic environment simulation and visualization［J］．Journal of Geomatics Science and Technology，2009，26(1):8 －11．

［4］王小林，笪良龍，楊廷武，等．基于預積分體渲染技術的水下三維聲場體可視化［J］．系統仿真學報，2007，19(21):4970 －4972．WANG Xiao-lin，DA Liang-long，YANG Ting-wu，et al．Underwater acoustic energy fields visualization based on technology of pre-integrated volume rendering［J］．Journal of System Simulation，2007，19(21):4970 － 4972．

［5］楊廷武，王小林，李玉陽，等．基于三維紋理的水下三維聲場直接體可視化［J］．系統仿真學報，2006，18(Suppl．1):322－323．YANG Ting-wu，WANG Xiao-lin，LI Yu-yang，et al．Research of direct volume rendering on underwater 3D acoustic energy field［J］．Journal of System Simulation，2006，18(Suppl．1):322 －323．

［6］王小林，笪良龍，楊廷武，等．水下3 維聲場體可視化兩種切平面實現方法［J］．中國圖象圖形學報，2008，13(2):312 －315．WANG Xiao-lin，DA Liang-long，YANG Ting-wu，et al．The realization of two clipping plane technologies in underwater acoustic energy fields visualization［J］．Journal of Image and Graphics，2008，13(2):312 －315．

［7］FINN B J，WILLIAM A K，MICHAEL B P，et al．Computation ocean acoustics［M］．Second Edition．New York，USA:Springer，2011:458 －459．

［8］KLAUS E，MARKUS H，JOE M K，et al．Real-time volume graphics［M］．Wellesley，Massachusetts，USA:A K Peters Ltd，2006:115 －116．