基于超二次曲面模型的超聲陣列成像算法*

季 寧， 楊 濤， 李西紅

(1.西南科技大學(xué) 信息工程學(xué)院，四川 綿陽 621010；2.西南科技大學(xué) 特殊環(huán)境機(jī)器人技術(shù)四川省重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621010)

0 引 言

隨著機(jī)器人技術(shù)的發(fā)展，機(jī)器視覺的研究受到研究者的日益關(guān)注，目標(biāo)物體的三維成像對機(jī)器人獨(dú)立完成任務(wù)起關(guān)鍵作用。目前常用于成像的傳感器有相機(jī)、深度傳感器和激光掃描儀等，但是其在核輻射、強(qiáng)粉塵、黑暗、煙霧等特殊環(huán)境中發(fā)揮作用有限。雷達(dá)與聲吶雖可在這些環(huán)境中工作，但雷達(dá)成本較高，聲吶存在波束分辨率低、掃描速度慢等固有缺陷。因此基于超聲陣列的三維目標(biāo)成像具有重大研究價(jià)值。目前超聲成像技術(shù)主要用于醫(yī)學(xué)成像與無損檢測方面[1]，空氣中的超聲三維成像研究較少。

文獻(xiàn)[2]利用采集到的回波數(shù)據(jù)和圖像幾何的先驗(yàn)知識(shí)來估計(jì)和校正成像參數(shù)誤差，生成自動(dòng)聚焦圖像,但該方法得到的圖像模糊，難以分辨出目標(biāo)物。文獻(xiàn)[3]為了提高最小方差超聲成像算法的分辨率、對比度以及對噪聲的魯棒性， 提出了一種改進(jìn)的最小方差成像算法,實(shí)現(xiàn)了點(diǎn)目標(biāo)物的成像，但未對目標(biāo)物的三維成像進(jìn)行進(jìn)一步研究。文獻(xiàn)[4～6]利用被檢測目標(biāo)物的方向角來實(shí)現(xiàn)超聲陣列對目標(biāo)物的檢測與識(shí)別，本文在此基礎(chǔ)上，利用渡越時(shí)間法(time of flight,ToF)得到特征點(diǎn)與超聲陣列之間的距離信息，從而得到目標(biāo)物特征點(diǎn)的坐標(biāo)信息，對三維數(shù)據(jù)點(diǎn)進(jìn)行超二次曲面擬合以實(shí)現(xiàn)目標(biāo)物的三維成像[5,6]。

文獻(xiàn)[7]提出了基于遺傳算法的兩階段擬合方法，實(shí)現(xiàn)了超二次曲面三維模型的重構(gòu),但遺傳算法需要調(diào)整較多的參數(shù)，計(jì)算效率較低。文獻(xiàn)[8]通過激光掃描儀獲得散亂點(diǎn)云數(shù)據(jù)，根據(jù)點(diǎn)云中目標(biāo)表面各點(diǎn)到目標(biāo)的均方距離，建立了三維目標(biāo)位姿估計(jì)的非線性目標(biāo)函數(shù)，實(shí)現(xiàn)了三維目標(biāo)的定位。文獻(xiàn)[9]將超二次曲面參數(shù)擬合問題轉(zhuǎn)化成非線性最小二乘問題，用Levenberg-Marquardt 算法進(jìn)行參數(shù)擬合,但LM算法不適用于目標(biāo)函數(shù)非線性程度很高的大殘量問題。

綜上所述，本文通過超聲陣列圍繞目標(biāo)物旋轉(zhuǎn)，分別從目標(biāo)物的前、后、左、右、上、下六個(gè)方向獲得回波信號(hào)，對回波信號(hào)峰峰值時(shí)間進(jìn)行切分，采用改進(jìn)的Unitary-ESPRIT算法對每個(gè)時(shí)間片段回波信號(hào)中特征點(diǎn)的方位角進(jìn)行估計(jì)，同時(shí)根據(jù)渡越時(shí)間法得到特征點(diǎn)與發(fā)射陣列之間的距離，利用三角幾何關(guān)系計(jì)算特征點(diǎn)的坐標(biāo)，對其進(jìn)行坐標(biāo)變換以得到特征點(diǎn)在全局坐標(biāo)系下的坐標(biāo)。借助改進(jìn)的LM算法對目標(biāo)物的特征點(diǎn)數(shù)據(jù)進(jìn)行參數(shù)擬合，獲得超二次曲面參數(shù)，以實(shí)現(xiàn)目標(biāo)物的三維成像。

1 超聲陣列的分布

超聲波換能器的中心頻率會(huì)直接影響超聲波的衰減系數(shù)和吸收系數(shù)，進(jìn)而對空間中超聲波的聲場分布產(chǎn)生影響，提高超聲波換能器的中心頻率可提高橫向分辨率，但是超聲波頻率過高會(huì)造成數(shù)據(jù)量大和衰減過大，不利于信號(hào)的采集和處理。為了在保證檢測精度的同時(shí)避免超聲波能量的過快衰減，本文選擇頻率為40 kHz的波換能器。隨著超聲波換能器直徑的增大檢測的空間分辨率會(huì)提高，同時(shí)引起柵瓣的幾率也會(huì)增大。在綜合考慮了分辨率及柵瓣等因素后，選擇壓電片直徑為11 mm的超聲波換能器。

對制作成本、衰減、分辨率等因素進(jìn)行綜合考慮后[10]，本文接收陣列采用陣元間距為25 mm 、陣元大小為8×8的二維面陣，發(fā)射陣元采用單發(fā)射方式且位于接收陣列中間。超聲陣列排布如圖1所示。

圖1 超聲陣列排布示意

2 空間特征點(diǎn)坐標(biāo)信息的獲取

通過渡越時(shí)間法與改進(jìn)的Unitary-ESPRIT算法[11]得到空間特征點(diǎn)與超聲陣列之間的距離及其對應(yīng)的方位角，并將其轉(zhuǎn)換為坐標(biāo)，關(guān)系如圖2所示。

圖2 空間特征點(diǎn)投影示意

由圖2可獲得每個(gè)空間特征點(diǎn)的坐標(biāo)，xk=dtanθcosφ,yk=dtanθsinφ,zk=d。其中,θ為z軸與OP的夾角，φ為OP在xoy平面的投影與x軸的夾角，d為特征點(diǎn)與發(fā)射陣列之間的距離。

本實(shí)驗(yàn)從目標(biāo)物的前、后、左、右、上、下六個(gè)方向來獲得特征點(diǎn)的坐標(biāo)，利用坐標(biāo)變換對空間特征點(diǎn)以陣列中心為原點(diǎn)的局部坐標(biāo)系下的坐標(biāo)進(jìn)行轉(zhuǎn)換得到以物體中心為原點(diǎn)的全局坐標(biāo)系下的坐標(biāo)，從而獲得目標(biāo)物全部特征點(diǎn)在全局坐標(biāo)系下的坐標(biāo)。

2.1 特征點(diǎn)距離信息的獲取

利用TOF來計(jì)算特征點(diǎn)的空間距離的基本原理為：從發(fā)射陣列發(fā)射的超聲波在反射點(diǎn)處所形成反射回波中提取特征點(diǎn)，再根據(jù)發(fā)射與接收的時(shí)間差來計(jì)算特征點(diǎn)的空間距離。

所有特征點(diǎn)都包含在渡越時(shí)間為的回波信號(hào)中。將時(shí)間段Δt等分為多個(gè)時(shí)間段，求出這些時(shí)間片段所對應(yīng)的空間特征點(diǎn)，并用Unitary-ESPRIT算法得到相應(yīng)的仰角和方位角，當(dāng)時(shí)間片段很小時(shí)，可以認(rèn)為這段時(shí)間內(nèi)所對應(yīng)波形的所有特征點(diǎn)都在與發(fā)射陣列距離相等的等距面上。由于不同特征點(diǎn)所對應(yīng)的方位角和仰角不同則分布在等距面上的位置不同，其與發(fā)射陣列的距離可以利用這一小段時(shí)間的中心值與陣列發(fā)射信號(hào)速度的關(guān)系函數(shù)來確定。

2.2 空間特征點(diǎn)的波達(dá)方向估計(jì)

波達(dá)方向(direction of arrival,DOA)估計(jì)算法中MUSIC算法需耗時(shí)較長，Capon算法空間分辨率較低；ESPRIT算法雖然有耗時(shí)小、分辨率高的優(yōu)點(diǎn)，但其對相干源的分辨率還是較低。為解決相干源問題，將空間平滑去相干算法和Unitary-ESPRIT算法相結(jié)合。該算法基本過程如下：

1)利用去相干平滑算法獲得輸出信號(hào)矩陣X對應(yīng)的協(xié)方差矩陣Rxx。

2)計(jì)算出實(shí)矩陣對應(yīng)的協(xié)方差矩陣Ryy

Ryy=(QN?QM)Rxx(QN?QM)H

(1)

式中 ?為Kronecter積，M,N分別為x,y方向陣元數(shù)量。當(dāng)M或N為奇數(shù)時(shí),QN與QM的表達(dá)式

(2)

當(dāng)M或N為偶數(shù)時(shí)，QN與QM的表達(dá)式

(3)

3)根據(jù)[Re(Ryy)，Im(Ryy)]對相關(guān)矩陣進(jìn)行估計(jì)。

4)對相關(guān)矩陣進(jìn)行特征值分解后取較大的特征值的信號(hào)矢量，從而構(gòu)成平滑后的信號(hào)子空間Es。

5)通過特征分解獲得特征值，從而估計(jì)出信源方向角。

3 二次曲面參數(shù)擬合

二次曲面是由三元二次方程所表示的曲面，因其具有局部變化不敏感和視點(diǎn)不變性等優(yōu)點(diǎn)而被廣泛應(yīng)用于計(jì)算機(jī)圖形學(xué)、計(jì)算機(jī)視覺等領(lǐng)域。其標(biāo)準(zhǔn)方程為[12]

(4)

式中 (x，y，z)為超二次曲面上的點(diǎn)，a1,a2和a3為超二次曲面在x,y,z方向的尺寸參數(shù)，ε1和ε2為形狀參數(shù)。當(dāng)參數(shù)ε1和ε2發(fā)生改變時(shí)可以形成不同的曲面。對式(4)進(jìn)行變形可擴(kuò)大超二次曲面模型的表示范圍。超二次曲面的變形公式

(8)

式中μx,μy為由錐化參數(shù)kx,ky定義的錐化函數(shù)。

對空間三維點(diǎn)數(shù)據(jù)進(jìn)行超二次曲面參數(shù)擬合時(shí)需計(jì)算空間三維點(diǎn)與超二次曲面之間的距離

(9)

設(shè)超二次曲面中心為O，其與空間三維點(diǎn)P的連線為OP。則式(9)中d0表示OP的長度，ds表示OP與超二次曲面上的交點(diǎn)到點(diǎn)O的距離。為了實(shí)現(xiàn)全局誤差最小，可將其轉(zhuǎn)化為非線性最小二乘優(yōu)化問題，則擬合誤差f為

(10)

目前常用的超二次曲面參數(shù)擬合方法中牛頓法很難收斂即使收斂找到點(diǎn)也未必最優(yōu)，最速下降法存在鋸齒現(xiàn)象且計(jì)算效率較低。此外，考慮到本實(shí)驗(yàn)特征點(diǎn)數(shù)量較多，用遺傳算法尋找最優(yōu)解耗時(shí)較長。LM算法[13]比共扼梯度法和梯度下降法收斂速度更快，但該算法擬合效果受初值影響較大，因此，本研究利用改進(jìn)的LM算法[14]對超二次曲面參數(shù)進(jìn)行擬合。其迭代規(guī)則如下

xi+1=xi+h

(11)

式中xi為參數(shù)向量，h為迭代步長

(12)

式中λ為比例因子。目標(biāo)函數(shù)實(shí)際減小量和估計(jì)減小量之比r為

(13)

算法的實(shí)現(xiàn)過程如下：

1)設(shè)置x0，λ的初值與最小值λmin。

2)計(jì)算雅克比矩陣J，更新目標(biāo)向量。

3)計(jì)算在xi+1處目標(biāo)函數(shù)f的值。

4)計(jì)算目標(biāo)函數(shù)的值并根據(jù)式(13)對r值進(jìn)行判斷。當(dāng)r值大于0時(shí)接受該參數(shù)向量，反之則放棄。當(dāng)r值接近0時(shí)，使λ増大4倍；當(dāng)r值接近1時(shí)，則λ如式(12)；當(dāng)r值在0與1之間時(shí)，λ值不變。返回步驟(2)獲得新的xi+1，直至f的值滿足要求。

4 仿真結(jié)果分析

采用MATLAB軟件工具，模擬8×8的超聲陣列作為陣列模型，為了獲得目標(biāo)物整體的特征點(diǎn)信息，對幾何物體超聲陣列分別從物體前、后、左、右、上、下六個(gè)方向進(jìn)行超聲波發(fā)射與接收仿真，在每個(gè)方向上利用Unitary-ESPRIT算法獲取對應(yīng)的特征點(diǎn)方位角信，同時(shí)根據(jù)渡越時(shí)間法得到該方向上空間特征點(diǎn)的距離。將角度與距離轉(zhuǎn)換為以陣列中心為原點(diǎn)的特征點(diǎn)坐標(biāo)。進(jìn)行適當(dāng)?shù)淖鴺?biāo)變換得到該方向上的特征點(diǎn)在以物體中心為原點(diǎn)的坐標(biāo)系下的坐標(biāo)信息。最后將從六個(gè)方向得到的特征點(diǎn)坐標(biāo)信息進(jìn)行匯總，得到物體的全部特征點(diǎn)信息，擬合出該物體對應(yīng)的超二次曲面參數(shù)，從而得到幾何物體的三維圖像。

實(shí)驗(yàn)過程中考慮到陣列的接收范圍與可識(shí)別角度，超聲陣列與目標(biāo)物的距離設(shè)置為1 m。本問題屬于大殘量問題， 初值取0.6。

圖3(a1),(b1),(c1)為超聲陣列獲得的球體、長方體、圓柱體特征點(diǎn)坐標(biāo)信息，圖3(a2)(b2),(c2)為LM算法擬合后得到的物體圖像，圖3(a3),(b3),(c3)為LM改進(jìn)算法擬合后得到的物體圖像。

利用LM算法與LM改進(jìn)算法對球體特征點(diǎn)進(jìn)行擬合得到的參數(shù)實(shí)驗(yàn)數(shù)據(jù)分別如表1所示。

表1 球體參數(shù)LM算法改進(jìn)前后擬合實(shí)驗(yàn)數(shù)據(jù)

圖3 三種形狀體數(shù)據(jù)擬合前后對比

利用LM算法與LM改進(jìn)算法對長方體特征點(diǎn)進(jìn)行擬合得到的參數(shù)實(shí)驗(yàn)數(shù)據(jù)分別如表2所示。

表2 長方體參數(shù)LM算法改進(jìn)前后擬合實(shí)驗(yàn)數(shù)據(jù)

利用LM算法與LM改進(jìn)算法對圓柱體特征點(diǎn)進(jìn)行擬合得到的參數(shù)實(shí)驗(yàn)數(shù)據(jù)分別如表3所示。

表3 圓柱體參數(shù)LM算法前后改進(jìn)擬合實(shí)驗(yàn)數(shù)據(jù)

在實(shí)驗(yàn)過程中隨機(jī)選取了多組初值進(jìn)行擬合，LM算法擬合時(shí)初始值的選取對結(jié)果影響較大，擬合到的參數(shù)未必能得到目標(biāo)物的圖片，初始值選取不恰當(dāng)時(shí)則無法擬合出超二次曲面參數(shù)，利用改進(jìn)的LM算法擬合則較好地改善了這一問題，擬合結(jié)果優(yōu)于LM算法。表1、表2、表3是選取其中一組初值得到的擬合參數(shù)，從表中可看出兩種得到的擬出合參數(shù)值與標(biāo)準(zhǔn)參數(shù)值之間存在一定誤差，但改進(jìn)的LM算法的參數(shù)誤差普遍小于LM算法的參數(shù)誤差，且改進(jìn)后LM算法的參數(shù)誤差基本控制在10 %以內(nèi)。圖3給出了兩種算法擬合后球體、長方體以及圓柱體的圖像，從中也可看出LM改進(jìn)算法擬合后的圖像效果更好。

5 結(jié) 論

本文通過超聲陣列獲得物體的回波信號(hào)，并對其峰峰值時(shí)間進(jìn)行切分，利用Unitary-ESPRIT算法與渡越時(shí)間法得到空間特征點(diǎn)坐標(biāo)信息，通過改進(jìn)的LM算法擬合出超二次曲面參數(shù)，最終實(shí)現(xiàn)三維目標(biāo)的成像。由仿真結(jié)果可知：該算法可有效解決LM算法對初值依賴過大的問題。本方法可擴(kuò)展到多目標(biāo)對象的復(fù)雜場景三維成像，也可對方向角估計(jì)算法進(jìn)行改進(jìn)，以達(dá)到更好的效果。