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

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

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

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

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

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

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

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

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

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

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

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

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

1 多個(gè)點(diǎn)散射體環(huán)境下聲場(chǎng)的Foldy-Lax方程

(1)

點(diǎn)散射體位置的頻域總聲場(chǎng)表示為：

(2)

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

(3)

(4)

其中：

(5)

(6)

(7)

(8)

(9)

G0(r,s,ω)=

(10)

G0(s,s,ω)=

(11)

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

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

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

(12)

將Eq.12代入Eq.4，點(diǎn)散射體位置的聲場(chǎng)表示為，

(13)

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

(14)

(15)

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

(16)

(17)

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

Λ=AB-1，

(18)

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

3 數(shù)值仿真

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

3.1 仿真環(huán)境

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

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

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

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

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

(19)

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

(20)

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

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

(21)

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

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

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

4 結(jié)論

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

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

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

(3)噪聲環(huán)境下，布設(shè)環(huán)境對(duì)三種方法的影響程度不同。不同布設(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]. 廣州: 華南理工大學(xué)出版社， 2006. Luo J, Fang W.Introduction to Matrix Analysis[M]. Guangzhou: South China University of Technology Press, 2006.

責(zé)任編輯 陳呈超

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

國(guó)家自然科學(xué)基金項(xiàng)目(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)， 王寧， 李小雷， 等. 一種非迭代的聲學(xué)反演點(diǎn)散射體散射系數(shù)方法[J]. 中國(guó)海洋大學(xué)學(xué)報(bào)(自然科學(xué)版), 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)