非均勻流體介質內部散射聲場重建的逐層計算方法*

唐少杰 向 宇 石梓玉

(1 廣西科技大學 廣西汽車零部件與整車技術重點實驗室 柳州 545006)

(2 廣西科技大學機械與汽車工程學院 柳州 545006)

0 引言

超聲層析成像廣泛應用于生物醫學、無損檢測、地球物理等領域[1-5]，其基本原理是利用超聲波照射物體，由外部超聲散射場求解物體內部結構，亦可稱為逆散射問題。對于流體介質內部的散射聲場的求解，即已知介質內部結構，由對比度函數求解介質內部散射聲場分布，以往采用矩量法(Method of moment,MoM)求解上述問題。MoM 的相關研究最早可追溯至20 世紀60年代，Richmond[6-7]采用MoM 精確描述了微波通過二維介質時的散射場。隨后，Johnson 等[8]將MoM 應用于超聲前向散射問題中。自此MoM 在超聲技術中獲得廣泛應用[9-12]。但在超聲問題中，為獲得較高分辨率，入射聲波頻率常取200 kHz 以上，且為得到良好的計算精度，需在一個波長內劃分3～10個單元，由此所得網格密度極大，若此時采用MoM 離散整個求解域將形成大規模矩陣，普通計算機算力難以求解，限制了其在實際中的應用。

國內外學者為使MoM 適用于大規模問題提出了一系列改進方法，通常分為兩類。第一類以快速多極子算法及多層快速多級子算法為代表[13-16]，該類算法將MoM中“點對點”的直接運算通過中間變換和傳遞分解劃分為“塊對塊”的間接計算，并采用迭代求解器求解線性方程組以避免顯式生成系數矩陣，以此降低內存使用量與計算復雜度。第二類以亞全域基函數法[17]、特征基函數法[18]為代表，該類算法使用高階基本函數離散向量積分方程，產生一個高度壓縮的線性標量矩陣方程，大幅度減少了矩陣方程中的變量個數。

針對MoM 計算非均勻流體介質內部散射聲場形成矩陣規模大進而對算力要求高的問題，本文結合聲散射基本理論與近場聲全息技術推導出逐層離散、逐層計算非均勻介質內部散射聲場的理論公式，并給出對應的幾何離散模型。在相同模型、條件下應用MoM 與逐層算法進行數值仿真，結果表明逐層算法可以有效重建流體介質內部散射聲場并減小矩陣求解規模。

1 基本理論

當超聲波入射非均勻流體介質時，由于散射作用，該介質成為次級聲源，向四周輻射散射波。此時，空間任意點r的全場聲壓可由式(1)描述：

其中，p(r)、pin(r)和ps(r)分別為r處的全場聲壓、入射聲壓和散射聲壓。

式(1)滿足非齊次亥姆霍茲方程[5]，可用Lippmann-Schwinger積分方程表示：

其中，S為非均勻流體介質存在區域；G(r,r′)是自由空間Green 函數，在二維問題中為(kd)[19-20]，其中的(kd)表示0 階第2類Hankel 函數，d=‖r-r′‖，k為波數；o(r′)為像函數，是描述物體內部介質聲學特性參數的函數，o(r′)=k2(n2(r′)-1)，n2(r′)-1 為對比度函數，n(r′)=c0(r′)/c(r′)，c0(r′)為背景介質聲速，c(r′)為介質聲速；散射場ps(r)表示為

以往使用MoM 離散計算區域簡化求解過程。MoM 采用脈沖基和點匹配，將式(3)轉化為代數方程組的形式，若介質內部結構已知，即各離散單元的像函數o(r′)已知，可求出空間任意點的全場與散射場。MoM 計算思路簡單，易于編程實現，但更適用于低激勵頻率、單元大且少的問題，隨著激勵頻率增高，網格劃分逐漸密集，離散整個求解域S將形成極大規模的復數滿秩矩陣，普通計算機已無法滿足算力要求，不利于實際應用。

2 逐層算法

為克服MoM 上述缺陷，提出一種結合聲散射基本理論與近場聲全息技術逐層求解流體介質內部聲場的計算方法。為便于計算，將介質存在區域拓展為規則形狀(Region of interest,ROI)，若對該區域逐層離散、逐層求解，可逐步得到整個ROI 區域的聲場分布。通過這種逐層計算方法，對比于MoM的矩陣求解規模明顯減小。

假設ROI 區域可劃分N層單元，由第1 節理論可知，計算第n層(n=1,2,···,N)任意單元處的全場，應由入射聲源、前n-1 層已離散計算得到聲場分布的單元聲源、第n層所有單元聲源及第n+1～N層暫未離散區域聲源共同作用得到，求解計算層聲場的關鍵在于表示未離散區域聲源對計算層單元處的聲場。因近場聲全息技術可以匹配檢測面上的信息重建聲源外部聲場，且無論是邊界元法(Boundary element method,BEM)還是波疊加法(Wave superposition method,WSM)均可將全部已知量和待求量集中在邊界或虛擬邊界上，達到降維、簡化計算的效果[19,21]，所以近場聲全息技術可以應用于求解未離散區域聲源對第n層單元處的聲場計算問題。又因WSM 具有算法簡單、計算效率高的優點，在近場聲全息中得到了廣泛應用，所以本文采用WSM求解上述問題。

2.1 WSM及其離散形式

WSM 基本思想是：在計算聲源外部散射場時，次級聲源向外散射的聲場可由連續分布于其內部的虛擬源的聲場疊加替代，虛擬源強可以通過匹配檢測面的相關信息獲取。如圖1 所示，O為坐標原點；S0為散射聲源真實邊界；S′為置于散射聲源內部的虛擬邊界；Ωout為散射聲源外部散射域；r′為虛擬邊界上虛擬源的位置向量；r為Ωout區域內任意場點的位置向量。

圖1 WSM 示意圖Fig.1 Diagram of WSM

聲源外任意一點的全場為

其中，Q(r′)為虛擬源強；G(r,r′)為關于虛擬源和場點位置矢量的傳遞矩陣；散射場為

2.2 離散模型描述

逐層算法采用聲全息的聲場重構原理，為保證聲場重建的穩定性，應將虛擬邊界S′、聲源邊界S0與檢測邊界設置H為共形，又因檢測邊界常為圓形，所以本文將ROI區域、虛擬邊界與檢測邊界均設置為圓形，如2 圖所示。ROI 區域等寬劃分N層，前N-1 層每一層劃分M個矩形單元，第N層為1 個圓形單元。

如圖2所示，將ROI區域劃分為3個部分，S2與S1邊界之間的區域為前n-1 層聲源(已通過逐層縮進計算得到聲場分布的區域)，記為A區域；S1與S0邊界之間為第n層聲源，記為B區域；S0邊界以內為未離散區域，記為C區域，該區域聲源對S0邊界以外場點的散射聲場可用式(5)表示；在ROI 區域外的圓形邊界為聲壓檢測邊界。

圖2 逐層算法幾何離散模型示意圖Fig.2 Diagram of layer-by-layer geometric discrete model

2.3 逐層計算算法

場點r處的全場應由A、B、C區域所有散射聲源以及入射聲源共同作用得到，式(6)可重寫為

其中，pA(r)、pB(r)、pC(r)分別為A、B、C三個區域散射聲源作用在C區域外任意點r處的聲壓。

A區域包含M×(n-1)個單元，在場點r處的聲壓為

B區域包含M個單元，在場點r處的聲壓為

為簡化計算，將C區域中虛擬邊界S′同樣離散為M個單元，在場點r處的聲壓為

oA=diag(o())，表示由A區域單元像函數形成的(M×(n-1))× (M×(n-1))維對角陣；oB=diag(o())，表示由B區域單元像函數形成的M×M維對角陣；

整理式(11)可得C區域虛擬源的源強列向量QC表達式為

QC為僅關于pB的函數，其他參數均為已知。為得到pB，需將場點配于B區域每個單元中心：

聯立式(12)、式(13)，整理得到pB表達式：

式(14)中，E為M×M維單位矩陣。

由此便得到了B區域所有單元的全場聲壓，計算出B區域聲壓分布后納入區域A作為已知量，隨著n增大。逐層內縮，重復上述過程即可得到聲源內部即ROI區域所有單元的聲場分布情況。需注意因WSM計算時聲源邊界S0內部存在虛擬面S′，因此計算層最大值nmax不能取到ROI 區域最內層，nmax大小為當虛擬面取到ROI 區域最內層時的計算層n。此時還未求解的總規模相比整個ROI區域規模非常小，可直接使用MoM 計算，具體計算過程不加以贅述。

特殊的，當n=1 時，ROI 區域僅由B、C區域組成，此時在檢測面H上配點可得pH為

對B區域(第n層)單元中心配點得到全場聲壓pB為

聯立式(15)、式(16)解得：

逐層算法優勢不僅在于可以減小矩陣規模，而且能減小數據檢測規模。在MoM 中，聲源外部散射方程組通常表示發射器從某一個方向上發射超聲波所對應的方程組，一般情況下檢測陣列點數H小于ROI 區域總單元數M×(N-1)+1，此時式為關于像函數o的欠定方程組，為使其可解，需從Φ個方向上發射超聲波，得到Φ×H個散射方程，使[Φ×H]≥M×(N-1)+1，此時方程組雖然可解，但也極大地擴充了矩陣規模與檢測規模；逐層算法則不需通過上述方案補全不充分的檢測陣列數據，因一次僅計算一層單元，輕易滿足H≥M。

3 仿真

為驗證“逐層算法”是否能夠正確重建流體介質內部聲場分布，本文對8 種模型進行仿真，將“逐層算法”與MoM得到的計算結果進行對比。在仿真時，仿真軟件采用MATLAB R2021a；計算機處理器為Intel(R) Core(TM) i7-10700 CPU@2.90 GHz，運行內存為16 GB。為得到更直觀的結果，直接對比ROI區域的散射場分布。因逐層算法需獲取檢測面信息作為已知量用于求解，檢測面上測點的散射場由MoM得到。

入射聲壓、ROI區域相關參數如表1所示。

表1 入射聲壓、ROI 區域相關參數Table 1 Relevant parameters of incident sound pressure and ROI area

表1 中所給ROI 區域單元劃分基于極坐標系，第n層任意單元節點周向坐標為r=a-nl，n=1,2,3,···；第n層中第m個單元節點徑向坐標為θ=2mπa/l,m=1,2,3,···；結合上述ROI區域劃分方法及表1 給出的入射聲壓相關參數，任意場點處入射聲壓可表示為

式(18)中，r為場點的位置向量，θ為r與x軸正方向的夾角。

對上述8 種模型設置不同的介質分布進行仿真，繪制ROI區域散射聲壓云圖，如圖3所示。

圖3 ROI 區域散射聲壓云圖Fig.3 Scattered Sound Pressure Cloud Chart in ROI Area

相較于MoM 一次性求解全部單元的未知量，逐層算法的優勢主要在于將求解域分層后逐次計算，進而可減小計算機內存消耗，并降低對計算機的算力要求。如表2 所示，以模型7、模型8 為例，給出了逐層算法計算一層占用內存與MoM 計算占用內存。

表2 MoM 與逐層算法求解一層時所形成矩陣規模與內存占用情況Table 2 Matrix size and memory occupation of moment method and layer by layer algorithm when solving one layer

在模型7 中，ROI 區域內包含2395 個單元，每層單元數有126個單元，共劃分20層。當使用MoM求解時，一次需要計算的矩陣規模為2395×2395，即需存儲及處理5,736,025個元素，生成并處理矩陣占用了577 MB 內存；而采用逐層算法求解時，求解域劃分為20 層，除最內層包含一個圓形單元外，其余每層均僅包含126 個矩形單元，一次運算時形成的矩陣規模為126×126，僅需存儲及處理15,876個元素，生成并處理矩陣僅占用20 MB 內存，相較于MoM 大幅度減小。同理，在模型8 中，使用MoM 求解，生成并處理矩陣占用了424 MB 內存；采用逐層算法求解，生成并處理矩陣僅占用17 MB內存。

需要說明的是，文中為了采用MoM 的計算結果作為對比，因而將模型的規模設置得較小。在實際應用中，求解域的尺寸和入射波的頻率可能會比算例中大得多。例如，若求解域的半徑為0.5 m，入射波頻率為1 MHz，同樣在一個波長λ內劃分5 個單元，那么總共就需要劃分1668 層，每層10,472 個單元。此時，若采用逐層算法，每一層的矩陣規模均為10,472×10,472，經計算發現，在不考慮系統及軟件需消耗的基礎內存情況下，內存占用已經達到8.15 GB；而若采用MoM，其生成矩陣的規模為17,456,824×17,456,824，遠遠超過逐層算法，按照逐層算法一層計算所占內存與層數相乘估計，需消耗運行內存8.15 GB×1668=13,586 GB，遠超普通計算機算力限額。由此可見，所需求解問題的規模越大，逐層算法相對于MoM的優勢就越明顯。

4 結論

針對超聲散射求解非均勻流體介質內部散射聲場的問題，本文提出了一種計算方案——逐層算法，結合近場聲全息理論與聲散射理論，實現逐層對非均勻流體介質內部散射聲場求解。在數值仿真上以MoM 計算結果為參考，仿真結果表明，本文提出的逐層算法可以有效地重建非均勻流體介質內部散射聲場，且由于一次僅計算一層單元，大幅度減小了矩陣的運算規模。

雖然本文中“逐層算法”只給出了在非均勻流體介質內部散射聲場重建問題上的應用，但當介質對比度未知時，通過在式(11)、式(13)之間反復迭代，即首先給出介質對比度的初始估計，將其帶入式(13)中求得較近似的全場，再將求得的全場帶入式(11)中求解像函數，如此反復迭代便可逼近每一層單元像函數的真值。因此使用“逐層算法”重建非均勻流體介質內部散射聲場問題是超聲層析成像的第一步，本文提出并驗證了“逐層算法”的可行性，為“逐層算法”在超聲層析成像中的應用奠定了基礎。