有限元/邊界元耦合方法在海洋生物目標強度特性研究的應用?

于小濤 王新良 張吉昌 趙憲勇

(1 中國水產科學研究院黃海水產研究所 農業農村部極地漁業可持續利用重點實驗室 青島 266071)

(2 深藍漁業工程聯合實驗室 青島 266237)

(3 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

0 引言

目標強度(Target strength,TS)特性是海洋生物聲學識別和資源量評估的重要依據[1?4]，基于理論模型法的數值仿真是TS 特性研究的重要手段。受計算方法與能力的限制，最早將海洋生物目標近似為規則的球體、圓柱體以及橢球體等具有解析解的模型[5?6]，特別是具有紡錘體形狀的魚類，一般利用回轉橢球體模型(Prolate spheroid model,PSM)求解散射聲場。為簡化解析求解的難度，Stanton[7]進一步將回轉橢球體、彎曲圓柱體等截面為圓形的目標分割為一系列的有限長圓柱體，提出了基于模態級數的形變圓柱體模型(Modal series-based deformed cylinder model,MB-DCM)，該模型適用于高長短軸比、聲波近垂直入射的情況，被Gorska[8]應用于鯡魚的聲散射特性研究。Foote 等[9]針對鱈科(Gadidae)魚類等尺寸較大的不規則目標，基于高頻近似發展了Kirchhoff近似模型(Kirchhoffapproximation,KA)。Clay 等[10?11]提出了一種考慮魚體和魚鰾相互作用的Kirchhoff射線模型(Kirchhoffray mode,KRM)，在低頻區域結合MB-DCM研究了大西洋鱈魚(Gadus morhua)的TS 特性。對于磷蝦類(Euphausia)和橈足類(Copepoda)等自身聲阻抗與其所處媒介相近的弱散射目標，Chu等[12]和Stanton 等[13]建立了畸變波波恩近似模型(Distorted wave born approximation,DWBA)。我國學者在漁業種類TS 研究方面以實驗測量為主[14?17]，數值仿真研究主要將MB-DCM、PSM、KRM和KA等近似模型應用于我國的海水及淡水魚類[18?22]。這些模型對聲波頻率、入射角度以及目標聲阻抗、形態尺寸等各自均有不同條件的限制[23]，以致于有時同一模型難以描述不同種類或同一種類但不同尺寸的生物種類，使得這些模型各自存在其適用的局限性。

有限元方法(Finite element method,FEM)[24?25]是近些年興起的數值求解方法，相比于當前普遍使用的KRM 和MB-DCM 方法要求滿足聲波近垂直入射、目標高長短軸比及聲阻抗均勻等條件，以及DWBA 模型僅適用于與其所處媒介聲阻抗接近的弱散射目標，FEM 能夠適用于全頻段、任意形狀、任意非均勻聲阻抗目標的聲散射求解，且數值計算精度高。Forland 等[26]利用FEM 研究了玉筋魚(Ammodytes personatus)的寬頻聲散射特性，并分析了骨骼對反向散射強度的貢獻。由于TS 關注的是遠場問題，單純利用FEM 需要較大的求解域，產生額外的計算量。針對遠場輻射與散射問題，聲學領域新興的有限元/邊界元(Finite element method/boundary element method,FEM/BEM)耦合分析方法[27?33]分別利用FEM和BEM 求解目標邊界和外部聲場，既能保證模型準確性，又有效降低了單純FEM 的計算時耗。目前，FEM/BEM 耦合方法正逐步應用于艦船、潛艇等大目標的低頻聲輻射與聲散射特性研究，尚未推廣到魚類等海洋生物目標。

本文以球形生物、尾明角燈魚(Ceratoscopelus warmingii)和南極大磷蝦(Euphausia superba)等為例，利用FEM/BEM耦合方法仿真計算其TS，通過與球形目標的解析模型、MB-DCM 和DWBA 等傳統理論模型的數值計算結果對比，分析FEM/BEM耦合模型的可靠性及優勢。

1 聲散射求解方法

1.1 FEM/BEM耦合方法

FEM 用于在頻域數值求解非均勻Helmholtz方程得到聲壓的穩態解。散射體被均勻介質和完美匹配層包圍，需對整個域進行網格剖分，遠場散射聲壓由目標邊界聲場的Helmholtz-Kirchhoff積分求得。BEM 通過目標表面的聲壓和法向聲速求解邊界上的Helmholtz-Kirchhoff積分方程，僅需對目標表面進行網格剖分，一般用于輻射問題；而對于多數的散射情況，目標表面聲壓和法向聲速是未知的。FEM/BEM耦合方法將BEM看作一種對FEM的補充，對目標體的內部問題使用FEM，對輻射區域的目標體外部問題使用BEM，前者得到目標表面聲壓和法向速度，后者得到遠場輻射聲壓。該方法可以求解具有任意形狀、任意材料的非均勻目標的散射聲場。

理論方程 散射目標的內部聲場利用FEM 求解，聲壓滿足Helmholtz方程[34]：

其中，k1=ω/c1，p1為目標內部聲壓，ω為入射聲波角頻率，c1為目標體內聲速。目標外部聲場通過BEM 求解，聲壓可由Helmholtz-Kirchhoff積分計算[34]：

其中，k=ω/c，c為目標所處媒介的聲速，p2=ps+pinc，ps為散射聲壓，pinc為入射聲壓。G(r,r0)為外部自由場的格林函數，r為散射聲場場點的位置，r0為散射目標表面單元的位置，與目標的傾角θ有關。邊界上滿足p1s=p2s的聲壓連續條件。TS可由反向散射聲壓求得：

聲散射仿真過程 散射聲場的數值仿真采用COMSOL Multiphysics?5.4 軟件[35]，主要步驟如下：

TS基于穩態聲場計算，魚鰾和魚體簡化為流體介質(只有縱波)[8]，采用壓力聲學(頻域)和壓力聲學(邊界元)的物理場接口進行FEM/BEM 耦合求解。前者計算目標內部及表面聲場，求解域為目標實體；后者根據表面聲場計算外部輻射聲場，求解域為外部無限域。在邊界元物理場中添加入射聲壓，包括入射聲壓的頻率與方向。

目標幾何形態方面，COMSOL 既支持簡單的幾何實體建模，也支持其他建模工具的幾何導入。一般通過X 射線或CT 成像得到其三維形態結構，通過直接繪制或外部導入以重構目標。

定義目標及周圍水體的密度與聲速、目標的傾角等變量。

目標網格剖分采用自由四面體網格自動剖分方法基本滿足計算需求，網格單元尺寸不超過波長的1/6。此外，需通過對計算結果進行收斂性測試，優化網格劃分。

頻域求解中頻率是必要的參數，如果求解頻率響應，需添加參數化掃描。海洋生物在水下的體長、傾角均具有不確定性，一般需要對其TS 進行統計平均，通過添加參數化掃描，可以求解不同體長、傾角時的散射聲場。

1.2 球形目標的解析模型

解析模型主要基于分離變量法求解散射聲場，球體目標反向散射函數的表達式為[5,23]

其中，An=?1/(1+iCn)。漁業相關生物種類的組織器官主要呈流體態(氣態和液態)，根據聲壓和聲速連續的邊界條件求解得到[5,23]：

其中，jn和yn分別為第一類和第二類n階球貝塞爾函數，分別是jn和yn的一階偏導數，g和h分別為目標與其所處媒介的密度比和聲速比，k是聲波在海水中的波數，k1是聲波在目標體內的波數，a為球體半徑。目標的反向散射截面為

TS由反向散射截面求得：

1.3 MB-DCM模型

MB-DCM 模型主要用于求解圓形截面、高長寬比目標在近法向入射時的散射聲場。以有鰾魚類為例，將魚體和魚鰾分別簡化為充滿氣體和液體的回轉橢球體，其反向散射函數可表示為[8]

其中，l是目標(魚體或魚鰾)的長軸長度，u=x/(l/2)，x是沿長軸的圓柱散射體微元相對于長軸中心的距離，D(θ)是反向散射的傾角指向性函數，函數bn的表達式為[8]

其中，εn是諾伊曼系數，當n=0時，εn=1；當n>0時，εn=2。Jn和Nn分別是第一類和第二類n階柱貝塞爾函數，分別是Jn和Nn的一階偏導數，是圓柱微元的橫截面半徑，w為短軸長度。魚鰾和魚體的反向散射指向性函數表達式分別為[8]

其中，下標sb和b分別代表魚鰾和魚體，θ是入射聲波相對于目標長軸的夾角，?θ是魚鰾相對于魚體的夾角。因樣品均系出水后采集，即相當于樣品取自表層，因此文中數值計算中均將z設置為0 m。

1.4 DWBA模型

DWBA 模型主要用于目標與其所處媒介聲阻抗相近的聲散射(如浮游動物)求解，對于截面為圓形的細長形目標，可簡化為沿其中心軸線的一系列充滿流體介質的離散圓柱體微元，Stanton 等[13]求得DWBA模型的一維線性解，反向散射函數表達為

式(14)中，k1i為入射波在目標體內的矢量波數，r0表示圓柱體微元中心的矢量位置，βtilt表示圓柱體微元截面與入射波之間的夾角，aj表示每個圓柱體微元的半徑，J1表示第一類1 階柱貝塞爾函數。壓縮系數γk和γρ分別為

其中，ρ和c分別表示目標所處媒介的密度和聲速，ρ1和c1表示目標體內的密度和聲速。對于圓弧形彎圓柱體近似目標，公式(14)可簡化為[13]

其中，ρc表示彎圓柱的曲率半徑。

2 TS仿真結果

海洋生物目標與其所處媒介(周圍水體)的密度比和聲速比采用前人的測量結果[5,36]，如表1所示。

表1 目標與其所處媒介(周圍水體)的密度和聲速比[5,36]Table 1 Density and sound speed ratio between the target and surrounding water[5,36]

2.1 球形生物

橈足類浮游動物一般可以被近似為液態球形目標，而受數值計算方法與能力的限制，早期將浮游動物和魚類也近似為球形目標。由于球形目標的散射聲場具有精確的解析解，本文將其用于驗證數值模擬方法。以充氣(氣泡)和液態球體為例，分別利用球形目標的解析模型和FEM/BEM 模型計算其TS。球體半徑假設為10 mm，充氣球體和液態球體的密度比和聲速比見表1。TS隨頻率變化的計算結果如圖1所示。圖1顯示，FEM/BEM模型計算結果與解析解高度吻合，充氣球體的TS 遠大于液態球體。

圖1 充氣和液態球體TS 的頻率響應Fig.1 Frequency response of TS for a gas-filled and liquid sphere

2.2 紡錘形魚類

對于紡錘形魚類，可將其形態近似為回轉橢球體，以南海尾明角燈魚(圖2(a))為例，通過X 射線掃描(圖2(b))并測量得到魚鰾和魚體的形態學參數見表2，近似的回轉橢球體如圖2(c)所示。有鰾魚類的聲散射主要來自魚鰾，因此以魚鰾TS 代表尾明角燈魚；以魚體TS 代表同類體型的無鰾魚類。魚鰾和魚體的物理參數分別采用表1 中充氣和液態橢球體的參數。魚鰾的TS 計算采用MBDCM、FEM/BEM 以及Ye[37]的低頻修正模型，魚體TS 的數值計算采用MB-DCM、FEM/BEM 以及DWBA 模型，數值求解得到魚鰾和魚體TS 隨頻率和傾角的變化分別如圖3～5所示。

圖2 尾明角燈魚的形態學建模Fig.2 Morphological modeling for a Ceratoscopelus warmingii

圖3 尾明角燈魚魚鰾和魚體TS 的頻率響應Fig.3 Frequency response of TS for a Ceratoscopelus warmingii with a swimbladder

表2 7 條尾明角燈魚樣品的平均形態參數Table 2 Mean morphological parameters of 7 Ceratoscopelus warmingii samples

圖3 為垂直魚體長軸方向(傾角為0?)入射時，魚體和魚鰾TS隨頻率的變化。對于魚體，MB-DCM和FEM/BEM 模型得到的結果除共振頻率附近基本一致；對于魚鰾，FEM/BEM模型和MB-DCM 模型在共振附近的結果差別非常大，這是由于魚鰾長短軸比值約3.5 : 1，不足以滿足其高長短軸比的適用條件。Ye[37]建立的充氣橢球體低頻聲散射模型對MB-DCM 模型共振頻率附近的結果進行了修正，并通過解析模型進行了驗證；同時，他發現MB-DCM 模型的準確性隨著回轉橢球體長軸與短軸比值的升高而提高。本文中，FEM/BEM 模型在共振頻率附近的結果與Ye的結果基本吻合，高頻計算結果與MB-DCM模型計算結果相近。

圖4為38 kHz和120 kHz頻率入射下，魚鰾TS隨魚類游泳傾角的變化。結果顯示，FEM/BEM和MB-DCM 模型的計算結果差別非常大，FEM/BEM 模型計算的TS 從垂直入射(傾角為0?)向兩端入射隨傾角的增加TS 減小的速度更慢，且曲線相對平滑，沒有峰值和谷值的震蕩起伏。圖5 為38 kHz和120 kHz入射頻率下，魚體TS隨傾角的變化。由圖5可知，在垂直入射方向附近，FEM/BEM、DWBA 和MB-DCM 模型的計算結果一致，但隨著傾角的變化，3 種模型計算結果的差異逐漸增加；FEM/BEM 和MB-DCM 模型吻合的傾角范圍極窄，FEM/BEM 模型和DWBA 模型吻合的傾角范圍較寬。

圖4 魚鰾TS 的傾角分布Fig.4 TS as a function of tilt angle for the swimbladder

圖5 魚體TS 的傾角分布Fig.5 TS as a function of tilt angle for the fish body

2.3 細長形浮游動物

對于細長形浮游動物，以南極大磷蝦(圖6(a))為例，可將其近似為呈圓弧形彎曲的形變圓柱體，如圖6(b)所示。假設圓柱體長度L=38.35 mm[36]，截面半徑為a=L/10.5，曲率半徑ρc=10L，聲學參數采用表1 中彎曲圓柱體的參數。采用FEM/BEM和DWBA 模型數值求解得到TS 隨頻率和傾角的變化分別如圖7 和圖8 所示。圖7 表明，FEM/BEM模型得到的垂直入射時TS頻率響應曲線與DWBA模型計算結果比較一致，在曲線峰值和谷值處兩者的結果有所差異。圖8 表明，FEM/BEM 模型得到的TS 隨傾角的變化曲線與DWBA 模型計算結果整體比較吻合，但隨著傾角增加，兩者計算結果的差異逐漸增加，尤其曲線峰值和谷值處。

圖6 南極大磷蝦的形態學建模Fig.6 Morphological modeling of Euphausia superba

圖7 南極大磷蝦TS 的頻率響應Fig.7 Frequency response of TS for an Euphausia superba

圖8 南極大磷蝦TS 的傾角分布Fig.8 TS as a function of angle of incidence for an Euphausia superba

3 結論與討論

文中利用國際上已廣泛應用的MB-DCM 和DWBA 散射模型驗證FEM/BEM 耦合方法在海洋生物聲散射研究的適用性，尾明角燈魚的形態學參數為本研究團隊前期基于X 射線掃描的實測結果[38]，南極大磷蝦的形態學參數和聲學參數采用南極海洋生物資源養護委員會(CCAMLR)推薦的標準[36]?；贔EM/BEM 耦合方法的球形生物聲散射模型與解析模型TS 的計算結果完全吻合，驗證了FEM/BEM 耦合模型的可靠性。對于紡錘形魚類的魚鰾，FEM/BEM 模型既解決了MB-DCM模型在中低頻段的準確性問題，又彌補了Ye 低頻散射模型在高頻段時的不足；對于魚體(無鰾魚類)，FEM/BEM 與DWBA 模型的仿真結果相近；在傾角變化時，FEM/BEM 和MB-DCM、DWBA 模型得到的TS 峰值和谷值出現在不同角度，且傾角越大，兩者偏差越大(如圖4 和圖5 所示)。對于細長形浮游動物，FEM/BEM 與DWBA 模型的仿真結果高度一致。但DWBA 模型的基礎是玻恩近似理論，僅僅適用于密度和聲速與其所處媒介(周圍水體)接近的弱散射目標(如浮游動物)，即目標與周圍水體的密度、聲速比接近1[23]，不適用于有鰾魚類等含氣囊海洋生物的TS 計算；而FEM/BEM 耦合方法可仿真具有不規則形態和不均勻參數的任意目標(有鰾魚類、無鰾魚類、浮游動物等)，能夠更精確地模擬生物各部分的聲散射，且適用于所有入射頻率和入射方向。因此，FEM/BEM 耦合方法在具有不同復雜形態的海洋生物TS 仿真研究方面的應用前景廣闊。

關于FEM/BEM 與MB-DCM 模型仿真結果隨傾角增加產生的差異，主要由于MB-DCM 模型的TS 表達式是在聲波近似垂直入射、目標高長短軸比的假設下推導得到的[7]，當傾角增大時模型自身的適用性降低；FEM/BEM 和DWBA 模型理論上均適用于全入射角度，但兩者隨傾角增加也產生一定差異。因此，FEM/BEM 模型在聲波端向入射時的準確性尚需進一步驗證。由于FEM/BEM耦合方法具有可仿真非均勻目標散射的優點，后續應用于海洋生物聲散射建模時將對目標的形態學參數和聲學參數作精細測量，充分發揮其可構建更精準物理模型的優點。對于無其他理論模型可仿真對比的非均勻目標的散射，需結合實驗測量進一步驗證FEM/BEM耦合方法的可靠性。

海洋生物普遍以大量尺寸不一、游泳傾角不同個體構成的集群形式混合棲息于海洋環境中，因此，對不同生物集群的分類識別、資源量評估有實際應用價值的是魚群的散射特性(如頻率響應特征、平均單體TS)。后續在單體TS 特性研究的基礎上，將考慮生物集群中大量單體的體長、游泳傾角的統計特征，分析生物集群的平均散射頻率響應特征和平均單體TS，推動FEM/BEM 耦合方法實際應用于海洋生物資源聲學資源評估。

由于FEM/BEM 耦合方法基于對目標幾何的網格化離散求解波動方程，計算結果的準確性很大程度上取決于網格的疏密，一般網格最大尺寸不超過聲波波長的1/6。因此，頻率越高(目標尺寸與入射聲波波長的比值越大)，網格剖分越密，FEM/BEM 耦合方法對計算能力的需求呈指數增長?？紤]到未來研究生物集群的散射特性需在體長、游泳傾角的二維空間進行統計平均，FEM/BEM耦合方法計算過程中對頻率、體長、游泳傾角3 個維度進行參數化掃描，計算量較大。因參數化掃描產生的計算量，常用的解決方案是增加計算節點、利用并行計算等方式。同時，可結合幾何模型的優化，充分利用目標的對稱性以降低計算規模；對于具有不規則幾何形狀的生物，可將其細分為更小的子域，合理設置網格大小，通過收斂性試驗，保證一定計算精度的前提下，提高單次計算效率。