基于直線檢測算法的矢量波導(dǎo)不變量估計

莫世奇，王博，李韜哲

（1.哈爾濱工程大學(xué) 水聲技術(shù)重點實驗室，黑龍江 哈爾濱 150001；2.哈爾濱工程大學(xué) 水聲工程學(xué)院，黑龍江 哈爾濱 150001；3.杭州應(yīng)用聲學(xué)研究所，浙江 杭州 310000）

海洋波導(dǎo)的重要特點是在低頻聲場中存在穩(wěn)定可觀測的干涉結(jié)構(gòu)［1-3］，Kuperman［4］認為近代水聲學(xué)以及水聲信號處理的主要研究方向就是海洋波導(dǎo)中低頻聲場的干涉結(jié)構(gòu)，隨著矢量水聽器技術(shù)的蓬勃發(fā)展，聲場中的振速信號受到越來越多人的關(guān)注，越來越多的水聲工作者投入到對聲壓信號和振速信號聯(lián)合處理方法的研究上，共同促使聲場干涉結(jié)構(gòu)逐漸成為國內(nèi)外水聲研究的重點。

對聲場干涉結(jié)構(gòu)研究的另一個里程碑式的節(jié)點是波導(dǎo)不變量理論的提出，通過波導(dǎo)不變量描述海洋波導(dǎo)中的干涉條紋［5-6］，波導(dǎo)不變量的分布情況與波導(dǎo)環(huán)境、模態(tài)階數(shù)有關(guān)，反映到干涉圖像上就是其干涉條紋的斜率變化［7］，因此需要借助圖像處理技術(shù)提取干涉圖中干涉條紋的信息。

用于干涉條紋提取的較為常見的圖像處理算法為Hough變換和Radon變換［8-10］，二者的原理都是將像素從圖像空間映射到參數(shù)空間，這2 種變換的提取效果普遍受遍歷角度步長設(shè)置的影響較大，特別是想?yún)^(qū)分多條斜率相近的直線時，對遍歷角度的步長要求更加嚴格，會導(dǎo)致需要一個非常龐大的參數(shù)空間去承載變換后的結(jié)果，而在圖像空間中的信息較少時，在參數(shù)空間很難形成比較尖銳的峰值。而Hough變換的下一步操作就是通過設(shè)置閾值將符合門限要求的點定義為直線，如果無法累積出足量像素點，很容易發(fā)生漏檢。同時，由于空間變換類算法都需要對整個圖像空間進行處理再累積，其處理速度會受圖像空間大小的影響，而且結(jié)果必然是表征圖像空間的整體特性，無法分別表征圖像空間中各條直線的特征。此外，Hough 變換要求輸入為二值圖像，所以在提取干涉條紋斜率前需要對圖像進行二值化操作，而且二值化精度對于變換最終結(jié)果影響較大；Radon 變換是依靠線積分實現(xiàn)的，雖然不要求輸入為二值化矩陣，但是輸出結(jié)果會受輸入圖像像素點密度分布影響，而且如果直線結(jié)構(gòu)的強度分布不均勻，其提取結(jié)果也將不準確。

因此需要尋找其他算法進行干涉條紋提取。直線檢測（line segmen detection,LSD）算法由von 等［11］提出。從圖像梯度角度入手，經(jīng)過精細化處理就能夠得到一條線段。它具有不需要額外的參數(shù)輸入、能夠控制自身的誤檢數(shù)、在較短時間內(nèi)獲得較高精度的直線段檢測結(jié)果等優(yōu)點。本文提出利用LSD處理算法進行矢量聲場干涉條紋的提取，并對典型淺海環(huán)境聲場的波導(dǎo)不變量進行估計，仿真分析表明LSD直線檢測算法的檢測效果令人滿意。

1 矢量聲場和波導(dǎo)不變量理論

1.1 矢量聲場干涉結(jié)構(gòu)

理想流體介質(zhì)中聲矢量場涵蓋質(zhì)點位移re-iωt、質(zhì)點振速ve-iωt、質(zhì)點加速度ae-iωt等，對于諧和波場，它們之間可以相互導(dǎo)出。聲波傳播過程必然伴隨能量的傳遞，分析聲場能量的干涉特性，從聲壓譜入手，根據(jù)簡正波理論［12-13］，得到柱坐標系下聲壓譜密度函數(shù)為：

式中：Δknm(ω)=krn(ω)-krm(ω)為第n階與第m階簡正波水平波數(shù)差；p(r,z,ω)=為有限階簡正波和的形式下，聲場中某點聲壓的表達式，聲壓譜即為此點聲壓與其共軛復(fù)數(shù)的乘積。

在平面波情況下，分層介質(zhì)中的動能譜密度和勢能譜密度相等，聲強與聲能量密度在水平方向上可看作正比關(guān)系，所以聲壓譜的干涉特性可以推廣到振速譜中。

振速自譜為：

式中：第1部分為各階模態(tài)自身作用的非相干項；第2 部分為模態(tài)間的相干項，引起譜強度振蕩變化，直接導(dǎo)致明暗相間條紋的形成。

1.2 波導(dǎo)不變量表征干涉條紋方程

簡正波的相互干涉，淺海波導(dǎo)中的聲場會形成穩(wěn)定的條紋結(jié)構(gòu)，在其距離-頻率平面上，如果Δr?r0，Δω?ω0，r0和ω0分別為中心距離和中心頻率，可認為在此干涉條紋上聲強不變，則其全微分為0，可以表示為：

式(2)轉(zhuǎn)化為：

式中vp和vg分別為水平相速度和群速度。

干涉條紋的斜率與聲源距離、頻率之間的關(guān)系用一個參數(shù)——波導(dǎo)不變量β來描述，其表達式為：

由式(3)可以求過寬帶聲源距離-頻率域中某點(r0,ω0)的干涉條紋方程為：

即：

可見，干涉條紋上頻點與距離的β次方有關(guān)，當波導(dǎo)不變量β ≈1 時，此式描述為一族經(jīng)過頻點0 的直線方程。

2 矢量聲場干涉條紋提取

2.1 矢量聲場干涉結(jié)構(gòu)融合

聲壓譜和水平振速自譜的干涉條紋雖然在譜強度上存在差異，但是都存在穩(wěn)定的干涉結(jié)構(gòu)，且條紋具有一致性。聲壓和振速在時空上具有相關(guān)性，不同頻率及位置處的干涉條紋強度不同，信息上具有互補性，采用圖像融合能在聲壓、振速干涉結(jié)構(gòu)中相關(guān)信息最大合并的基礎(chǔ)上減少輸出的不確定度和冗余度，擴大干涉圖所含有的時間-空間信息，減少不確定性，增加可靠性，使融合后的干涉圖對聲場干涉結(jié)構(gòu)會有更全面、清晰的描述，從而更有利于干涉條紋的提取。

本文從數(shù)據(jù)級圖像融合入手，并以圖像的標準差、平均梯度和平均空間頻率為標準，對圖像融合結(jié)果進行評價。標準差反映圖像內(nèi)像素點的離散程度，數(shù)值越大表明圖像質(zhì)量越優(yōu)良，平均梯度和空間頻率反映圖像的清晰度，數(shù)值越大表示圖像越清晰。

選擇基于小波變換的圖像融合方法，將圖像分解為4 層，對各分解層從高到低分別進行融合。常規(guī)的融合方式為在高頻部分采用加權(quán)平均法或基于區(qū)域能量的融合準則等方法進行融合，在低頻部分采用平均的線性融合方法。但是在許多情況下，單用取平均的方法會影響最終的融合效果，所以本文采用平均與選擇相結(jié)合的融合規(guī)則，融合結(jié)果如表1所示。

表1 融合結(jié)果Table 1 Fusion results

將2 幅圖像的相關(guān)性作為選擇融合方式的標準時，設(shè)定2 幅圖像的相似程度系數(shù)，當2 幅圖像相關(guān)性程度較強時，即認為2幅圖像的能量相差不大，此時低頻部分采用平均的線性融合方法，高頻部分采用加權(quán)平均法；當2幅圖像相關(guān)性程度較弱時，則認為2 幅圖像局部能量相差較大，此時低頻部分仍然采用平均的線性融合方法，高頻部分基于區(qū)域能量的融合準則，選取局部區(qū)域能量較大的小波系數(shù)作為融合圖像的小波系數(shù)。

從融合結(jié)果的評價指標中可以看出，將2 幅不同質(zhì)量的圖像進行融合后，平均與選擇相結(jié)合的融合方法可以保持圖像的質(zhì)量與清晰度不會下降，并能夠得到一定的增強效果。

2.2 LSD算法提取干涉條紋

2.2.1 LSD算法原理

LSD 算法是一種直線檢測分割算法，它能在線性的時間內(nèi)得出亞像素級精度的檢測結(jié)果，并且可以自己控制誤檢的數(shù)量。其算法流程如圖1所示。

圖1 LSD算法提取直線段流程Fig.1 Flow chart of line segment extraction by LSD algorithm

1)對輸入圖像各像素點梯度值進行計算。盡可能少使用其他像素，以減少對彼此的依賴性。按照圖2 所示選取每個像素點右下的4 個像素點計算梯度，分為x和y方向，并計算相應(yīng)的梯度幅值，梯度方向的垂直方向構(gòu)成水平線場（level-line）。

圖2 像素點計算范圍Fig.2 Calculation range of pixels

2) 對梯度幅值進行排序。在圖像中具有較大梯度幅值的像素點往往意味著該像素點表征較強的邊緣，而較小梯度幅值的像素點意味著梯度變換緩慢的區(qū)域，即圖像中平坦的區(qū)域。由于像素值的量化問題，利用較小幅值的像素點進行運算會引入誤差，采取設(shè)置閾值來拒絕此類像素點參與計算，閾值為：

式中τ為角度允許誤差，一般選擇為π 8。

將過檢后的梯度幅值進行降序排列并分配到列表中，從具有最高梯度幅值的像素點入手，依次以列表中的像素作為種子點，進行區(qū)域生長。區(qū)域生長原則是將種子點鄰域內(nèi)誤差小于τ且未被使用的像素點納入生長區(qū)域，再以生長區(qū)域內(nèi)滿足條件的像素點為中心向外生長。整個區(qū)域的初始角度為種子點的角度，之后每添加一個新的像素到該區(qū)域，區(qū)域的角度更新為：

當區(qū)域中所有新增點的八鄰域都不滿足偏差條件時，停止生長。

3)對生長區(qū)域進行矩形近似。將整個區(qū)域當作一個實體，每個像素點的梯度幅值大小作為點的質(zhì)量，以此計算矩形的重心。則質(zhì)量為G(j)的點重心為(cx,cy)，該矩形的主慣性軸方向為矩陣M的最小特征值對應(yīng)的特征向量角度。

4)對矩形近似結(jié)果進行評價。評價過程分為2個階段：第1 階段的依據(jù)為近似矩形內(nèi)與矩形方向相同的像素點密度是否滿足閾值。同向點不滿足占全部像素點70%的條件，則需要以矩形的短邊為標準裁剪并重新進行區(qū)域生長，直到滿足閾值條件；第2 階段為采用文獻［17-18］提出的Helmholtz 原理，計算矩形中的誤檢數(shù)量(number of false alarms,NFA)值來確定是否接受或拒絕提取的直線，不滿足NFA值則改變矩形R區(qū)域使其滿足閾值為止。

2.2.2 干涉條紋提取

聲壓和振速在時空上具有相關(guān)性，信息上具有互補性，采用圖像融合方法將二者進行融合處理，使干涉圖對聲場干涉結(jié)構(gòu)具有更全面、清晰、多樣性的描述。采用小波變換的方法，從數(shù)據(jù)級圖像融合入手，將聲壓譜、水平振速譜的干涉圖像進行融合處理，并將處理后得到的圖像作為LSD算法的輸入圖像。

仿真環(huán)境：聲速梯度選取等聲速梯度，聲源設(shè)置位于水下4 m 處，接收器位于水下30 m 處；相應(yīng)的環(huán)境參數(shù)設(shè)置如圖3所示。

圖3 等聲速梯度下的環(huán)境參數(shù)Fig.3 Environmental parameters under constant sound velocity gradient

在等聲速梯度下，使用LSD 算法對干涉條紋進行提取，結(jié)果如圖4所示，可以得到較為清晰的干涉條紋線段；采用Hough 變換及Radon 變換對干涉條紋進行提取，參數(shù)空間中峰值尖銳，提取效果較好。

圖4 等聲速梯度LSD算法提取結(jié)果Fig.4 Extraction results of LSD algorithm under constant sound velocity gradient

圖5 等聲速梯度空間變換類算法提取結(jié)果Fig.5 Extraction results of spatial transformation algorithm under constant sound velocity gradient

仿真環(huán)境：聲速梯度選取SWellEx-96 試驗實測的聲速梯度，其表現(xiàn)為淺海負躍層波導(dǎo)環(huán)境，聲源設(shè)置位于水下4 m 處，接收器位于水下100 m 處，相應(yīng)的環(huán)境參數(shù)設(shè)置如圖6 所示。圖7 為LSD 算法提取結(jié)果，圖8為空間變換類算法提取結(jié)果。

圖6 SWellEx-96實驗聲速梯度下的環(huán)境參數(shù)Fig.6 Environmental parameters in SWellEx-96 experi‐ment under sound velocity gradient

圖7 LSD算法提取結(jié)果Fig.7 Extraction results of LSD algorithm

圖8 空間變換類算法提取結(jié)果Fig.8 Extraction results of spatial transformation algorithm

對比圖5和圖8的提取結(jié)果，當波導(dǎo)環(huán)境較為復(fù)雜時，空間變換類算法的參數(shù)空間中出現(xiàn)多個亮點，峰值出現(xiàn)偏移，因此從參數(shù)空間中提取到合適的最值點變得較為困難。

通過比較不同聲速梯度下兩類算法的提取結(jié)果，可見LSD算法檢測效果令人滿意。并且相較于空間變換類算法有其優(yōu)勢所在：能夠分別表征圖像空間中各條直線的特征，這是空間變換類算法所沒有的。

2.3 融合算法對直線提取的增強

為了更貼合實際的海洋波導(dǎo)環(huán)境，設(shè)置SWel‐lEx-96實驗聲速梯度和沉積層等海底環(huán)境共同作用下，利用LSD 算法進行干涉條紋提取時，在近程運動部分會出現(xiàn)許多較短的誤檢線段如圖9 所示，為了減少這些較短線段的干擾，先采用圖像融合的方式，將聲壓自譜與水平振速自譜進行融合，再使用LSD算法進行提取。

圖9 融合前后干涉條紋提取結(jié)果Fig.9 Extraction result of interference fringes before and after fusion

對比圖9 與圖10，從提取結(jié)果中線段長度的統(tǒng)計數(shù)量可以看出，采用圖像融合的方式，能夠在一定程度上減少未融合前的較短檢測線段數(shù)量，增加檢測結(jié)果的有效性。

圖10 融合前后干涉條紋統(tǒng)計結(jié)果Fig.10 Extraction result of interference fringes before and after fusion

3 波導(dǎo)不變量的估計

由式(4)、(5)可以看到，給出了2 種計算波導(dǎo)不變量β的方法，分別為通過建模得到簡正波的相速度和群速度，利用相慢度和群慢度計算和利用圖像處理方法提取距離-頻率域干涉圖中的干涉條紋斜率，從而計算β［16］，仿真環(huán)境設(shè)置與第2節(jié)相同。

取觀測頻帶100～500 Hz，設(shè)置勻速運動目標速度2.5 m/s，總觀測時間為8 000 s，定義目標朝向接收器運動時時間為負，遠離接收器運動時時間為正。取其中-4 000～-1 000 s 的中遠程運動部分，此時運動方式可以定義為靠近運動，其運動方式示意圖及聲場干涉結(jié)構(gòu)仿真結(jié)果如圖11所示。

利用相慢度和群慢度計算各相鄰簡正波決定的波導(dǎo)不變量結(jié)果如圖12 所示，其數(shù)值維持在1～1.1。在相同環(huán)境下，將聲壓自譜與水平振速自譜進行融合，并利用LSD直線檢測算法提取聲場干涉結(jié)構(gòu)如圖13所示，并將其中所有值計算波導(dǎo)不變量如圖14所示。

圖12 波導(dǎo)不變量計算結(jié)果Fig.12 Computes the result of waveguide invariant

圖13 融合后LSD算提取結(jié)果Fig.13 LSD algorithm extraction result after fusion

圖14 線段特征統(tǒng)計直方圖Fig.14 Statistical histogram of line segments features

按照長度排序，選取前100 的直線段計算波導(dǎo)不變量。從圖14中可以看出，提取結(jié)果中長度小于100 的占所有線段數(shù)目中58%以上。這些較短的直線中，有很大一部分是誤檢結(jié)果，如果這部分直線參與到波導(dǎo)不變量的計算中，會拉低最終計算的波導(dǎo)不變量的可信程度。

選擇線段提取結(jié)果中長度較長的部分進行計算如圖15所示，在計算中以長度作為權(quán)重有利于降低誤差，提高計算結(jié)果的可信度。

圖15 線段重新排序Fig.15 Line segments are re-ordered

從圖16 中可以看出，仿真的波導(dǎo)不變量數(shù)值大約在0.92～1.04。以直線長度為權(quán)重計算波導(dǎo)不變量的平均值，得到波導(dǎo)不變量值為0.981 4。

圖16 波導(dǎo)不變量計算結(jié)果Fig.16 Computes the result of waveguide invariant

4 SWellEx-96試驗數(shù)據(jù)處理

4.1 試驗環(huán)境

試驗數(shù)據(jù)選自SWellEx-96 試驗期間的事件S5。在此事件中，大部分聲源拖曳海域深度為180～220 m，另一半試驗的聲源是沿180 m等深線進行拖曳的。試驗區(qū)域不存在較大的干擾，聲源艦船在所有陣列的南側(cè)開始向北航行，航行速度保持在5 kn（2.5 m/s），航行軌跡及接收陣布放情況見圖17。

圖17 SWellEx-96試驗中S5事件情況Fig.17 Experimental situation of S5 event in SWellEx-96

艦船航行區(qū)域海深約225 m，共布放有3組聲學(xué)傳感器，分別為垂直陣列（VLA）、傾斜陣列（TLA）和水平陣列（HLA），處理數(shù)據(jù)選自VLA 的測量記錄，21個陣元布放在深度94.124～212.5 m。

艦船牽引了2 個聲源，且聲源拖曳在一條等深線上，“淺源（J-13）”大約在9 m 深，發(fā)射109～385 Hz的9 個頻率，編號為C-109-9S；“深源（J-15）”大約在54 m 深，發(fā)射49～400 Hz 的多條線譜。

4.2 試驗數(shù)據(jù)處理

選取位于水下139.12 m 處水聽器的接收數(shù)據(jù)進行處理，得到時頻圖，如圖18所示。

圖18 接收陣元LOFARFig.18 LOFAR of receiving element

從處理結(jié)果可以看到，此次試驗線譜和干涉條紋都很明顯。線譜分布在49、64、94、112、130、148、166、201、235、283、338、388 Hz 附近，在400～-700 Hz可以看到清晰的干涉條紋。

選取400～700 Hz、300～900 s 范圍進行LSD 算法提取干涉條紋。

4.3 試驗數(shù)據(jù)處理結(jié)論

圖19 為事件S5 LSD 算法提取結(jié)果。干涉條紋進行Hough 變換，如圖20 所示會在變換域形成較為寬大的亮斑，無法在參數(shù)空間體現(xiàn)圖像的細節(jié)。同樣進行Radon 變換，無法在參數(shù)空間中形成可靠的亮斑，不能提取出參考線段的有效角度。使用LSD算法可以檢測出較為清晰的線譜線段，剔除線譜可得干涉條紋，但是受到環(huán)境干擾等影響，加之水聽器接收到的信號強度較弱，導(dǎo)致條紋強度不夠，因此檢測條紋需要進一步處理。

圖19 事件S5 LSD算法提取結(jié)果Fig.19 Extraction results of LSD algorithm of S5

圖20 事件S5空間變換類算法提取結(jié)果Fig.20 Extraction results of spatial transformation algorithm of S5

5 結(jié)論

1）對于大多數(shù)淺海條件，聲場中存在較為穩(wěn)定的干涉結(jié)構(gòu)，受到聲速梯度、聲源和接收器深度等因素的影響，空間變換類算法隨著干涉條紋偏移，提取結(jié)果也發(fā)生偏移。使用LSD 算法仍能較為詳細地提取干涉條紋。

2）利用干涉條紋斜率計算波導(dǎo)不變量相較于簡正波群慢度、相慢度的方法更為直觀，在對波導(dǎo)不變量進行估計時可結(jié)合二者進行使用。

在本文仿真部分中對于矢量場應(yīng)用的設(shè)計是利用其共有的干涉聲場信息進行相同結(jié)構(gòu)增強，但是在實驗數(shù)據(jù)處理中發(fā)現(xiàn)聲壓信號振速信號往往不夠強，導(dǎo)致增強效果有限且會引入新的誤差，在后續(xù)的工作中，考慮研究結(jié)合弱信號提取的圖像處理技術(shù)。