張寅權(quán)， 王 寧， 李小雷， 張 爽

(中國海洋大學信息科學與工程學院，山東 青島 266100)

一種非迭代的聲學反演點散射體散射系數(shù)方法*

張寅權(quán)， 王 寧， 李小雷， 張 爽

(中國海洋大學信息科學與工程學院，山東 青島 266100)

多體散射環(huán)境下，本文提出一種非迭代的聲學反演點散射體散射系數(shù)的方法。這種方法運用非向量數(shù)學運算，從而避免測量多基地響應(yīng)矩陣，與現(xiàn)有的方法相比更容易實現(xiàn)。本文對兩種布設(shè)環(huán)境數(shù)值仿真驗證該方法的有效性，并在噪聲環(huán)境下將本文方法與現(xiàn)有的兩種非迭代反演方法進行比較。仿真結(jié)果表明：(1)反演公式中矩陣的條件數(shù)由布設(shè)環(huán)境決定，布設(shè)環(huán)境通過條件數(shù)影響反演結(jié)果；(2)噪聲對反演結(jié)果的影響程度取決于布設(shè)環(huán)境；(3)布設(shè)環(huán)境對三種方法的影響程度不同，本文方法與另外兩種方法在不同布設(shè)環(huán)境中各有優(yōu)劣。

非迭代； 反演； 散射系數(shù)； 點散射體

聲散射反問題根據(jù)測量到的聲波散射數(shù)據(jù)反演散射體位置、散射系數(shù)等信息。聲散射反問題在超聲無損探傷、醫(yī)學超聲成像、結(jié)石治療等領(lǐng)域具有廣泛的應(yīng)用價值。

時反[1](Time reversal)是聲波散射體定位和散射系數(shù)反演的基礎(chǔ)。聲學時反最早由Fink、Prada等人[2-3]提出，通過對接收的聲信號進行時間反轉(zhuǎn)并反向傳播，實現(xiàn)未知介質(zhì)環(huán)境下聲波的自適應(yīng)聚焦。如果介質(zhì)中的聲場格林函數(shù)已知，就能夠根據(jù)發(fā)射信號確定聚焦點的位置。傳統(tǒng)時反只能聚焦在聲源或散射最強的散射體位置，為了能夠聚焦在散射相對較弱的散射體上，Prada等人提出時反矩陣特征分解方法[4-6](DORT)，建立散射體與時反矩陣特征向量的一一對應(yīng)。根據(jù)某個特征向量發(fā)射聲波，聲場可以聚焦在相應(yīng)的散射體上。DORT方法要求散射體到接收換能器陣的聲場格林函數(shù)向量相互正交，這使得DORT方法的空間分辨率比較低。Devaney等人將時反矩陣信號子空間與散射體到接收陣的聲場格林函數(shù)向量聯(lián)系起來，提出時反MUSIC方法[7-9]，極大提高了散射體定位的分辨率。Lehman[10]進一步將時反 MUSIC方法推廣到收發(fā)分置、非互易介質(zhì)環(huán)境。以上方法都是建立單體散射環(huán)境下，針對多體散射環(huán)境，Gruber[11]利用Newman級數(shù)證明時反 MUSIC方法在多體散射環(huán)境下仍然適用。Devaney[12]利用Foldy-Lax方程[13-15]描述多體散射，給出時反 MUSIC方法在多體散射環(huán)境下成立的一種簡潔證明。

散射系數(shù)的反演建立在散射體定位的基礎(chǔ)上。對這個問題的研究相對較少。Devaney[12]最早提出一種多體散射環(huán)境下反演點散射體散射系數(shù)的方法。該方法利用Foldy-Lax方程建立關(guān)于點散射體散射系數(shù)的非線性迭代方程組，并通過迭代求解點散射體的散射系數(shù)。為了避免迭代運算帶來的收斂性問題，Marengo[16]給出一種非迭代的反演多體散射環(huán)境下點散射體散射系數(shù)的方法。該方法根據(jù)點散射體到接收陣或發(fā)射陣的聲場背景格林函數(shù)向量(介質(zhì)中不存在散射體時的聲場格林函數(shù))線性無關(guān)，推導(dǎo)出點散射體散射系數(shù)關(guān)于背景格林函數(shù)向量的解析公式。然而噪聲會導(dǎo)致該解析公式與背景格林函數(shù)向量的線性無關(guān)不一致，從而使得散射系數(shù)的反演結(jié)果出現(xiàn)較大偏差。為了提高噪聲環(huán)境下散射系數(shù)反演的魯棒性，Xudong Chen[17]提出一種基于最小二乘法的非迭代方法。上述散射系數(shù)反演方法均用于反演多體散射環(huán)境下點散射體散射系數(shù)，要求散射體位置和背景格林函數(shù)已知，并且需要測量散射體的多基地響應(yīng)矩陣(Multistatic responsematrix)，即多個聲源依次發(fā)射聲波，多個接收器記錄聲波與散射體作用后的散射波。

在散射體位置和背景格林函數(shù)已知的環(huán)境下，本文提出一種新的非迭代反演點散射體散射系數(shù)的方法。這種方法采用相同維度的向量對應(yīng)元素相除這種特殊運算(不屬于矩陣運算)，避免對矩陣求逆，從而只需要測量單個聲源產(chǎn)生的散射場，相對多基地響應(yīng)矩陣更容易獲得。本文首先利用Foldy-Lax方程描述點散射體環(huán)境下的多體散射聲場并闡述本文方法的基本原理；然后通過對兩種布設(shè)環(huán)境作數(shù)值仿真，驗證本文方法的有效性，并在噪聲環(huán)境下同Marengo[16]和Xudong Chen[17]兩種非迭代方法比較反演誤差。仿真結(jié)果表明：

(1)無噪聲環(huán)境下，本文方法能夠準確反演散射系數(shù)；

(2)噪聲環(huán)境下，三種方法的反演結(jié)果受布設(shè)環(huán)境的影響。不同布設(shè)環(huán)境，本文方法與另外兩種方法各有優(yōu)劣。

1 多個點散射體環(huán)境下聲場的Foldy-Lax方程

(1)

點散射體位置的頻域總聲場表示為：

(2)

Eq.1和Eq.2就是描述多個點散射體環(huán)境下聲場的Foldy-Lax方程。為了后面描述方便，將其分別寫成矩陣的形式：

(3)

(4)

其中：

(5)

(6)

(7)

(8)

(9)

G0(r,s,ω)=

(10)

G0(s,s,ω)=

(11)

2 散射系數(shù)的反演原理

同其它反演點散射體散射系數(shù)的方法[12,16-17]一樣，本文提出的反演方法要求點散射體位置和背景格林函數(shù)已知。其中，點散射體位置可以通過時反MUSIC方法或其它途徑確定。下面闡述本文方法的基本原理。

當接收器數(shù)目N大于點散射體數(shù)目M時，Eq.10矩陣G0(r,s,ω)列滿秩，由Eq.3得:

(12)

將Eq.12代入Eq.4，點散射體位置的聲場表示為，

(13)

接下來是本文方法的關(guān)鍵。根據(jù)Eq.6和Eq.9，

(14)

(15)

其中，“{ }m”表示取花括號中向量的第m個元素。Eq.12，Eq.13和Eq.15即本文反演點散射體散射系數(shù)的方法。

(16)

(17)

其中，上標“(1)，(2)，…，(M)”表示不同的發(fā)射聲源。除了極個別情形，矩陣B可逆，散射系數(shù)表示為

Λ=AB-1，

(18)

由于本文方法采用Eq.15這樣的非矩陣運算，所以避免對矩陣求逆，只需單個聲源就可以操作。參考文獻12，16和17反演散射系數(shù)時都采用嚴格的矩陣運算，因此它們都需要測量多基地響應(yīng)矩陣。

3 數(shù)值仿真

這一節(jié)利用MATLAB數(shù)值仿真兩種點散射體布設(shè)環(huán)境下的聲場，檢驗本文方法的有效性，并且在不同信噪比條件下將本文方法與Marengo[16]和XudongChen[17]兩種非迭代方法的反演結(jié)果進行比較。

3.1 仿真環(huán)境

前面提到，Marengo[16]和XudongChen[17]兩種方法需要測量多基地響應(yīng)矩陣，為了同這兩種方法進行比較，本文通過收發(fā)合置的聲學換能器陣構(gòu)造多基地響應(yīng)矩陣，本文方法的聲源選擇中間位置的換能器。

3.2 無噪聲環(huán)境仿真結(jié)果

由圖3和圖4可見，無噪聲環(huán)境下，本文方法能夠準確地反演散射系數(shù)。

3.3 噪聲環(huán)境仿真結(jié)果

噪聲通過影響多基地響應(yīng)矩陣干擾散射系數(shù)反演。被噪聲污染的多基地響應(yīng)矩陣表示為：

(19)

其中：K表示無噪聲時的多基地響應(yīng)矩陣；Noise表示噪聲矩陣。同文獻12、16、17相同，本文采用加性高斯白噪聲(Additive white Gaussian noise)，矩陣元素是相互獨立的零均值復(fù)高斯隨機變量，信噪比定義為：

(20)

其中“‖ ‖”表示矩陣的F范數(shù)。

噪聲環(huán)境下，散射系數(shù)的反演結(jié)果通過散射系數(shù)相對誤差的百分比表示：

(21)

不同布設(shè)環(huán)境，不同信噪比，散射系數(shù)反演的相對誤差如圖5、6所示。兩圖為100次快拍統(tǒng)計平均的結(jié)果。對比圖5和6可見：

(1)相同信噪比條件下，同一種方法在布設(shè)環(huán)境1中的反演誤差小于布設(shè)環(huán)境2；(2)不同布設(shè)環(huán)境，本文方法與另外兩種方法各有優(yōu)劣。

本文方法與Marengo[16]方法對應(yīng)不同的矩陣的廣義逆，與Xudong Chen[17]方法相比，矩陣的廣義逆在反演公式中的嵌入方式不同。因此改變布設(shè)環(huán)境對三種方法的產(chǎn)生不同程度的影響，本文方法與其它兩種方法在不同布設(shè)環(huán)境下各有優(yōu)劣。

4 結(jié)論

在散射體位置和背景格林函數(shù)已知的多體散射環(huán)境下，本文提出一種非迭代反演點散射體散射系數(shù)的方法。這種方法采用相同維度的向量對應(yīng)元素相除這種特殊運算，不需要測量多基地響應(yīng)矩陣。通過對兩種布設(shè)環(huán)境數(shù)值仿真發(fā)現(xiàn)：

(1)無噪聲環(huán)境下，本文方法能夠準確地反演散射系數(shù)；

(2)噪聲環(huán)境下，不同布設(shè)環(huán)境，矩陣G0(r，s，ω)和k的條件數(shù)越大，散射系數(shù)的反演誤差越大；

(3)噪聲環(huán)境下，布設(shè)環(huán)境對三種方法的影響程度不同。不同布設(shè)環(huán)境，本文方法與另外兩種非迭代方法各有優(yōu)劣。

[1] Fink M. Time-reversed acoustics [J]. Rep Prog Phys, 2000, 63:1933-1995.

[2] Fink M, Prada C. Self-focusing in inhomogeneous media with time reversal acoustic mirrors [C]. Montreal,IEEEUltrasonics Symposium Proceedings, 1989, 2: 681-686.

[3] Prada C. The iterative time reversal mirror: A solution to self-focusing in the pulse echo mode [J]. J Acoust Soc Amer, 1991, 90: 1119-1129.

[4] Prada C. The iterative time reversal process: Analysis of the convergence [J]. J Acoust Soc Amer,1995, 97: 62-71.

[5] Prada C. Decomposition of the time reversal operator: Detection and selective focusing on two scatterers [J]. J Acoust Soc Amer,1996, 99: 2067-2076.

[6] Mordant M, Prada C, Fink M.Highly resolved detection and selective focusing in a waveguide using the D.O.R.T. Method [J]. J Acoust Soc Amer, 1999, 105: 2634-2642.

[7] Devaney A J. Super-resolution Processing of Multi-static Data Using Time Reversal and MUSIC[EB/OL]. www.ece.neu.edu/faculty/devaney, 2000.

[8] Lev-AriH, Devaney A J.The time-reversal technique reinterpreted: Subspace-based signal processing for multi-static target location [C]. Cambridge, MA, Sensor Array and Multichannel Signal Processing Workshop IEEE, 2000: 509-513.

[9] Prada C,Thomas J L. Experimental subwavelength localization of scatterers by decomposition of the time reversal operator interpreted as a covariance matrix [J]. J Acoust Soc Amer, 2003, 114: 235-243.

[10] Lehman S, DevaneyA J. Transmission mode time-reversal super-resolution imaging [J]. J Acoust Soc Amer, 2003, 113: 2742-2753.

[11] Gruber F K, Marengo E A, Devaney A J. Time-reversal imaging with multiple signal classification considering multiple scattering between the targets [J]. J Acoust Soc Amer, 2004, 115: 3042-3047.

[12] Devaney A J, Marengo E A, Gruber F K. Time-reversal-based imaging and inverse scattering of multiply scattering point targets [J]. J Acoust Soc Amer, 2005, 118: 3129-3138.

[13] Foldy L L. The multiple scattering of waves: I. General theory of isotropic scattering by randomly distributed scatterers [J]. Phys Rev,1945, 67: 107-119.

[14] Lax M. Multiple scattering of waves [J]. Rev Mod Phys, 1951, 23: 287-310.

[15] Ishimaru A. Wave Propagation and Scattering in Random Media [M]. New York: IEEE Press, 1997.

[16] Marengo E A, Gruber F K.Noniterative analytical formula for inverse scattering of multiply scattering point targets [J].J Acoust Soc Amer, 2006, 120: 3782-3788.

[17] Chen X.A robust noniterative method for obtaining scattering strengths of multiply scattering point targets(L) [J]. J Acoust Soc Amer, 2007, 122: 1325-1327.

[18] 羅家洪， 方衛(wèi)東. 矩陣分析引論 [M]. 廣州: 華南理工大學出版社， 2006. Luo J, Fang W.Introduction to Matrix Analysis[M]. Guangzhou: South China University of Technology Press, 2006.

責任編輯 陳呈超

A Non-Iterativeacoustic Method to Obtain Scattering Coefficients of Multiple Point Scatterers

ZHANG Yin-Quan, WANG Ning, LI Xiao-Lei, ZHANG Shuang

(College of Information Science and Technology, Ocean University of China, Qingdao 266100, China)

In the environment of multi-body scattering, a non-iterative acoustic method is proposed in this paper to obtain the scattering coefficients of point scatters. The method takes advantage of anon-vector mathematic operation to avoid measuringmulti-static response matrix.Comparing with existing methods, the proposed method is more convenient to operate. To verify the method,two layout environmentsare simulatednumerically and the simulation results of the proposed method are compared with the results of two existing non-iterative methods. The above simulation has taken into account the effect of noise. The result of simulationdemonstrates that: (1) thecondition number of the matrixes that appear in the formulas of the aforementionedthree methods is determined by the layout environment. It is this way that the layout environment affects simulation result; (2)the influence degree of noisedepends on layout environment; (3) the influence degreeof layout environment on the condition number of the matrixesin the above three methodsis different.The performance of the three methods depends on specific layoutenvironment.

non-iterative； inversion； scattering coefficient； point scatter

國家自然科學基金項目(11374270;11374271)資助

2014-10-15；

2015-09-20

張寅權(quán)(1986-)，男，博士生。E-mail:zhyq_ouc@126.com

** 通訊作者： E-mail:zhyq_ouc@126.com

O429

A

1672-5174(2017)04-126-06

10.16441/j.cnki.hdxb.20140412

張寅權(quán)， 王寧， 李小雷， 等. 一種非迭代的聲學反演點散射體散射系數(shù)方法[J]. 中國海洋大學學報(自然科學版), 2017, 47(4): 126-131.

ZHANG Yin-Quan, WANG Ning, LI Xiao-Lei, et al. A non-iterativeacoustic method to obtain scattering coefficients of multiple point scatterers[J]. Periodical of Ocean University of China, 2017, 47(4): 126-131.

Supported by the National Natural Science Foundation of China (11374270;11374271)