孫雪海，笪良龍，李玉陽

(海軍潛艇學(xué)院 海軍水下作戰(zhàn)環(huán)境研究所，山東 青島266071)

0 引 言

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

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

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

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

1 三維聲場計(jì)算

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

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

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

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

采用Tappert 方法，令:

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

當(dāng)k0r ≥1 時(shí)，Hankel 函數(shù)可用其漸近形式代替:

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

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

將近軸近似用于方程(6)，即可得到標(biāo)準(zhǔn)拋物方程:

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

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

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

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

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

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

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

2 三維聲場可視化

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

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

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

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

2.1 數(shù)據(jù)重采樣

對原始聲場數(shù)據(jù)進(jìn)行重采樣及坐標(biāo)轉(zhuǎn)換的具體步驟如下:

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

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

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

glViewport (0，0，XFBO，YFBO);

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

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

3)對存儲原始聲場數(shù)據(jù)的三維紋理進(jìn)行重采樣。使用OpenGL 著色語言，在片元著色器中實(shí)現(xiàn)對原始聲場數(shù)據(jù)的重采樣及坐標(biāo)轉(zhuǎn)換，并將結(jié)果渲染到與FBO 綁定的三維紋理上。

2.2 顏色映射

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

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

2.3 光照計(jì)算

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

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

3 仿真實(shí)驗(yàn)

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

曙光TC4000L 高性能計(jì)算平臺分配36 個(gè)節(jié)點(diǎn)分別計(jì)算不同方位的聲場數(shù)據(jù)。同時(shí)，為了比較驗(yàn)證并行計(jì)算的效率，在高性能計(jì)算平臺開辟了一個(gè)比較節(jié)點(diǎn)，用于串行計(jì)算。

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

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

圖6 并行和串行計(jì)算結(jié)果比較Fig.6 Results compare in parallel and series computation

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

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

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

4 結(jié) 語

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

［1］張林，笪良龍，李玉陽．三維水聲傳播建模及并行算法研究［J］．系統(tǒng)仿真學(xué)報(bào)，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］．系統(tǒng)仿真學(xué)報(bào)，2006，18(7):1980 －1982．DA Liang-long，F(xiàn)AN 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］李玉陽，笪良龍，宋潔．水聲環(huán)境仿真及其可視化技術(shù)研究［J］．測繪科學(xué)技術(shù)學(xué)報(bào)，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］王小林，笪良龍，楊廷武，等．基于預(yù)積分體渲染技術(shù)的水下三維聲場體可視化［J］．系統(tǒng)仿真學(xué)報(bào)，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］．系統(tǒng)仿真學(xué)報(bào)，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 維聲場體可視化兩種切平面實(shí)現(xiàn)方法［J］．中國圖象圖形學(xué)報(bào)，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．