朱正任，馮成凱，徐琪堯，韓 冰，付桂合，陽凡林,2

(1.山東科技大學(xué) 測(cè)繪與空間信息學(xué)院，山東 青島 266590； 2.海島(礁)測(cè)繪技術(shù)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室，山東 青島 266590)

多波束測(cè)深系統(tǒng)(multibeam echo sounder, MBES)不僅可以得到高精度的水深地形數(shù)據(jù)，還能獲取豐富的反向散射強(qiáng)度數(shù)據(jù)。反向散射強(qiáng)度在較大程度上可確定沉積物類型，是識(shí)別海底底質(zhì)類型的重要參數(shù)[1-3]。但是，在實(shí)際測(cè)量過程中，影響多波束反向散射強(qiáng)度的因素較多，除海底底質(zhì)的物理特性外，聲吶系統(tǒng)參數(shù)、海水特性(吸收和折射)、測(cè)量模式(發(fā)射距離和角度)等因素也會(huì)對(duì)反向散射強(qiáng)度產(chǎn)生影響[4-5]。因此，國(guó)內(nèi)外學(xué)者對(duì)反向散射強(qiáng)度數(shù)據(jù)處理進(jìn)行了廣泛的研究，并形成了基于主動(dòng)聲吶方程的較為完善的改正模型。反向散射強(qiáng)度數(shù)據(jù)經(jīng)過模型改正后，其中大部分影響可以消除，但角度響應(yīng)由于受到聲波散射機(jī)理的影響，表現(xiàn)為反向散射強(qiáng)度隨入射角而變化，傳統(tǒng)的角度響應(yīng)改正模型在復(fù)雜底質(zhì)環(huán)境下適應(yīng)性較差，改正效果不理想。

針對(duì)多波束聲吶圖像角度響應(yīng)改正問題，國(guó)內(nèi)外學(xué)者從兩種角度出發(fā)，針對(duì)AR改正后存在的殘差進(jìn)行改正。一種是從數(shù)理統(tǒng)計(jì)的角度出發(fā)，如唐秋華等[6]通過條帶兩側(cè)強(qiáng)度數(shù)據(jù)加權(quán)內(nèi)插得到中央波束區(qū)域強(qiáng)度數(shù)據(jù)進(jìn)行中央?yún)^(qū)域殘差改正；Yang等[7]采用等均值方差擬合的方法，通過對(duì)中央?yún)^(qū)域數(shù)據(jù)的壓縮和移動(dòng)進(jìn)行改正；王煜[8]通過構(gòu)建中央波束區(qū)與非中央波束區(qū)模型進(jìn)行殘差改正，上述方法均未考慮真實(shí)的聲波散射機(jī)理，改正效果往往與真實(shí)底質(zhì)有較大差異，造成一定的精度損失。另一種是從聲學(xué)機(jī)理角度出發(fā)，完善AR改正模型，如Hemmerstand等[9]建立了Lambertian模型進(jìn)行AR改正；Hellequin等[10]建立了綜合聲學(xué)模型進(jìn)行AR改正；嚴(yán)俊等[11]通過提取角度響應(yīng)曲線的AR參數(shù)，進(jìn)行AR改正模型的重構(gòu)。上述方法均有效減弱了角度響應(yīng)對(duì)聲吶圖像的影響，但未考慮不同底質(zhì)對(duì)AR改正模型的影響。Zhao等[12]建立了基于預(yù)分類的AR改正模型，但未給出高效的分類數(shù)的確定方法且計(jì)算量較大，存在一定局限性。AR曲線本身亦可以用于分類，但由于其以ping為單位，分類結(jié)果的精度低于聲吶圖像，故本研究有效結(jié)合兩者的優(yōu)勢(shì)，提出一種基于小波分析的多波束聲吶圖像角度響應(yīng)改正方法，通過小波分解技術(shù)，將AR曲線分解為長(zhǎng)波項(xiàng)與短波項(xiàng)，利用AR曲線長(zhǎng)波項(xiàng)三維概率密度圖先進(jìn)行同底質(zhì)區(qū)域劃分，然后在同底質(zhì)區(qū)域下完成AR改正。該方法可以在不需要確定分類數(shù)的前提下，有效消除不同底質(zhì)對(duì)AR改正的影響，獲取高質(zhì)量的聲吶圖像，為后續(xù)高精度的海底底質(zhì)分類提供基礎(chǔ)。

1 角度響應(yīng)曲線的小波分解及其三維概率密度表示

角度響應(yīng)曲線是反向散射強(qiáng)度隨入射角變化的曲線，由于海底底質(zhì)的不同，聲吶系統(tǒng)接收到的回波強(qiáng)度信號(hào)嚴(yán)格地說是非平穩(wěn)信號(hào)[13]，因此，需要用非平穩(wěn)隨機(jī)信號(hào)的方法進(jìn)行處理。小波變換在空間域和頻率域同時(shí)具有良好的分析特性[14]，可以對(duì)非平穩(wěn)信號(hào)進(jìn)行有效分析。故將角度響應(yīng)曲線進(jìn)行小波分解，分解后得到其長(zhǎng)波項(xiàng)(低頻信號(hào))和短波項(xiàng)(高頻信號(hào))。然而單ping角度響應(yīng)曲線的長(zhǎng)波項(xiàng)和短波項(xiàng)很難看出其規(guī)律性變化，為直觀地顯示長(zhǎng)波項(xiàng)和短波項(xiàng)的趨勢(shì)變化和分布情況，繪制長(zhǎng)波項(xiàng)和短波項(xiàng)的三維概率密度圖。

1.1 角度響應(yīng)曲線的小波分解

(1)

(2)

1.2 短波項(xiàng)三維概率密度圖及其特點(diǎn)

短波項(xiàng)代表角度響應(yīng)曲線的聲強(qiáng)特征變化，為直觀地顯示其反向散射強(qiáng)度隨入射角變化的整體趨勢(shì)，繪制短波項(xiàng)的三維概率密度圖。由于反射散射強(qiáng)度隨入射角的變化曲線呈對(duì)稱分布，選取[0°, 60°]進(jìn)行分析。

圖1(a)和1(b)分別是同底質(zhì)區(qū)域和不同底質(zhì)區(qū)域下相鄰50 ping角度響應(yīng)曲線短波項(xiàng)的三維概率密度圖(這里以小波參數(shù)：小波基coif5、分解層數(shù)1為例，說明短波項(xiàng)三維概率密度圖的特點(diǎn))。其中，橫坐標(biāo)代表入射角，縱坐標(biāo)表示反向散射強(qiáng)度，在指定的入射角處，反向散射強(qiáng)度出現(xiàn)的頻率以圖右側(cè)的顏色條表示?？梢钥闯?，無論在同底質(zhì)區(qū)域下還是不同底質(zhì)區(qū)域下短波項(xiàng)的反向散射強(qiáng)度均隨入射角的變化在0附近波動(dòng)，這是由于短波項(xiàng)代表角度響應(yīng)曲線變化的細(xì)節(jié)信息，不包含其整體變化趨勢(shì)。

圖1 短波項(xiàng)的三維概率密度圖

1.3 長(zhǎng)波項(xiàng)三維概率密度圖及其特點(diǎn)

長(zhǎng)波項(xiàng)代表角度響應(yīng)曲線的輪廓線，包含角度響應(yīng)曲線的主要變化趨勢(shì)。由于不同的底質(zhì)類型對(duì)應(yīng)不同的角度響應(yīng)曲線[15]，因此，長(zhǎng)波項(xiàng)三維概率密度圖可以直觀地顯示某一區(qū)域是否為同一底質(zhì)類型。

如圖2(a)和圖2(b)所示，分別是同底質(zhì)區(qū)域下和不同底質(zhì)區(qū)域下角度響應(yīng)曲線長(zhǎng)波項(xiàng)的三維概率密度圖。可以看出圖2(a)只包含1種類型的角度響應(yīng)曲線，認(rèn)為該區(qū)域是同一底質(zhì)類型；而圖2(b)中包含2種類型的角度響應(yīng)曲線，故該區(qū)域存在不同底質(zhì)類型。

2 多波束聲吶圖像角度響應(yīng)改正方法

多波束角度響應(yīng)曲線經(jīng)小波分解為長(zhǎng)波項(xiàng)和短波項(xiàng)。其中，長(zhǎng)波項(xiàng)包含角度響應(yīng)曲線的主要變化趨勢(shì)。因此，聲吶圖像的角度響應(yīng)改正即角度響應(yīng)曲線長(zhǎng)波項(xiàng)的改正，將改正后的長(zhǎng)波項(xiàng)與短波項(xiàng)進(jìn)行信號(hào)重構(gòu)，得到改正后的多波束聲吶圖像。

2.1 最優(yōu)小波分解參數(shù)的獲取

小波基和分解層數(shù)是小波分解的重要參數(shù)，不同分解參數(shù)對(duì)信號(hào)分解的影響不同[16]。因此，最優(yōu)參數(shù)的選擇是進(jìn)行角度響應(yīng)曲線小波分解的重要前提。本研究基于小波分解后的短波項(xiàng)不受角度響應(yīng)影響的特點(diǎn)，提出了利用短波項(xiàng)三維概率密度曲線的特征參數(shù)來確定最優(yōu)小波分解參數(shù)的方法(圖3)。

案例教學(xué)法還能實(shí)現(xiàn)教學(xué)相長(zhǎng)。在案例教學(xué)中，教師不僅是引導(dǎo)者而且也是學(xué)習(xí)者。一方面，教師掌握著教學(xué)進(jìn)程，引導(dǎo)學(xué)生思考、組織討論研究，并進(jìn)行總結(jié)、歸納。另一方面，收集案例的過程就是一個(gè)學(xué)習(xí)的過程，教師在教學(xué)中與學(xué)生共同討論同一話題，也可以從中獲得大量感性材料。教師在課堂上不再“獨(dú)唱”，既能更合理地利用體力和腦力，也能從中獲得一種成就感。

圖3 最優(yōu)小波參數(shù)獲取的流程圖

圖4(a)是短波項(xiàng)三維概率密度圖的剖面圖，表示在不同入射角處反向散射強(qiáng)度出現(xiàn)的頻率(橫坐標(biāo)為反向散射強(qiáng)度，縱坐標(biāo)為出現(xiàn)的頻率，這里以15°和60°為例)，由圖可以看出，在指定入射角處，角度響應(yīng)曲線的短波項(xiàng)近似服從N(0,σ2)的正態(tài)分布，其中σ2為方差。依據(jù)中心極限定理，角度響應(yīng)曲線序列數(shù)N>30時(shí)，子樣BSθ(N)、BSN(θ)序列均服從正態(tài)分布[17]，且均值在u=0處波動(dòng)，因此，其概率密度函數(shù)為

(3)

(4)

(5)

圖4 短波項(xiàng)三維概率密度的剖面圖及其概率密度曲線

圖5 統(tǒng)計(jì)參數(shù)對(duì)小波基和分解層數(shù)的評(píng)價(jià)結(jié)果

以短波項(xiàng)概率密度曲線的特征參數(shù)符合均值與方差近似為0的原則下(確保角度響應(yīng)曲線的特征趨勢(shì)不被分解到短波項(xiàng))，盡可能地保證較多的分解層數(shù)(確保角度響應(yīng)曲線分解得徹底)。由圖5(a)和圖5(b)可以看出，在不同小波基的情況下，均值絕對(duì)值與方差在分解層數(shù)為5時(shí)產(chǎn)生較大跳變，即分解層數(shù)大于5時(shí)，角度響應(yīng)曲線的特征趨勢(shì)會(huì)被分解到短波項(xiàng)中，因此，最優(yōu)分解層數(shù)為5層，結(jié)合表1，在相同分解層數(shù)為5時(shí)，小波基coif5的均值絕對(duì)值和方差分別為0.007 8和0.171 9，均小于小波基sym4(0.024 7,0.300 8)和小波基db6(0.011 5,0.317 5)的均值絕對(duì)值和方差。為此，將小波基coif5、分解層數(shù)為5層作為角度響應(yīng)曲線小波分解的最優(yōu)參數(shù)。圖6(b)和圖6(c)，為1 ping角度響應(yīng)曲線在小波基coif5和分解層數(shù)5時(shí)，小波分解得到的長(zhǎng)波項(xiàng)和短波項(xiàng)，其中，短波項(xiàng)不受角度響應(yīng)的影響，在強(qiáng)度值0附近波動(dòng)；長(zhǎng)波項(xiàng)變化趨勢(shì)明顯，包含角度響應(yīng)曲線的主要信息，驗(yàn)證了該參數(shù)選擇的有效性。

圖6 1 ping角度響應(yīng)曲線的小波分解

圖7 多波束條帶同底質(zhì)區(qū)域劃分流程圖

2.2 多波束條帶同底質(zhì)區(qū)域劃分

在獲取角度響應(yīng)曲線小波分解的最優(yōu)參數(shù)后，還需對(duì)多波束條帶進(jìn)行同底質(zhì)區(qū)域劃分，來消除不同底質(zhì)類型對(duì)AR改正的影響，進(jìn)行自適應(yīng)分區(qū)改正。由1.3節(jié)可知，長(zhǎng)波項(xiàng)三維概率密度圖可以判別某一區(qū)域是否為同一底質(zhì)類型，因此，采用二分法對(duì)多波束條帶進(jìn)行同底質(zhì)區(qū)域劃分，整體流程如圖7所示。

由于長(zhǎng)波項(xiàng)三維概率密度圖無法進(jìn)行定量分析，繪制同底質(zhì)區(qū)域下和不同底質(zhì)區(qū)域下長(zhǎng)波項(xiàng)三維概率密度圖(圖2)的剖面圖。在常見底質(zhì)下，入射角θ較小時(shí)(一般不超過15°)，換能器接收的多為鏡面反射，稱為D1區(qū)；當(dāng)θ變大時(shí)(一般在15°～60°之間)換能器接收的多為后向散射部分，稱為D2區(qū)，即漫反射區(qū)；θ進(jìn)一步增大時(shí)(一般超過60°)，稱為D3區(qū)或高入射區(qū)[12]。因此，選取繪制D1、D2、D3區(qū)的邊界角度15°和60°處的反向散射強(qiáng)度出現(xiàn)的頻率，如圖8(a)和圖8(b)所示。

圖8 長(zhǎng)波項(xiàng)三維概率密度圖的剖面圖

分析圖8可得，同底質(zhì)區(qū)域下，長(zhǎng)波項(xiàng)三維概率密度圖的剖面圖只存在1個(gè)波峰，當(dāng)某一區(qū)域存在不同底質(zhì)時(shí)，長(zhǎng)波項(xiàng)三維概率密度圖的剖面圖存在2個(gè)或多個(gè)波峰。故以長(zhǎng)波項(xiàng)三維概率密度在入射角15°和60°處的剖面圖波峰的個(gè)數(shù)為判定條件，采用二分法對(duì)多波束條帶進(jìn)行同底質(zhì)區(qū)域劃分(圖9)。

由于條帶左右兩部分可能存在不同的角度響應(yīng)曲線，將條帶劃分為左右兩部分進(jìn)行同底質(zhì)區(qū)域劃分(海底底質(zhì)復(fù)雜，以角度響應(yīng)曲線按塊劃分無法保證某一區(qū)域完全是同一底質(zhì)，這里對(duì)條帶進(jìn)行粗劃分，為后續(xù)AR改正提供同底質(zhì)區(qū)域范圍，提高改正效率)。

2.3 長(zhǎng)波項(xiàng)的AR改正及強(qiáng)度信號(hào)重構(gòu)

對(duì)于任意一條測(cè)線，采用2.1節(jié)中的方法可以獲得角度響應(yīng)曲線小波分解的最優(yōu)參數(shù)：小波基coif5和分解層數(shù)5；采用2.2節(jié)中的方法可以得到多波束條帶左右舷的同底質(zhì)區(qū)域。由上述分析可得對(duì)聲吶圖像角度響應(yīng)改正即對(duì)角度響應(yīng)曲線長(zhǎng)波項(xiàng)的改正。因此，算法整體流程如圖10所示。

圖9 二分法同底質(zhì)區(qū)域劃分流程圖

為避免單ping信號(hào)改正的偶然性，在同底質(zhì)區(qū)域下，將角度響應(yīng)曲線的長(zhǎng)波項(xiàng)疊加取平均，獲得該區(qū)域長(zhǎng)波項(xiàng)的均值曲線BSMean；然后將該區(qū)域每1 ping的長(zhǎng)波項(xiàng)減去該均值曲線，并統(tǒng)一歸化到長(zhǎng)波項(xiàng)漫反射區(qū)的平均值強(qiáng)度BSM，從而實(shí)現(xiàn)長(zhǎng)波項(xiàng)的AR改正。即：

(6)

3 實(shí)驗(yàn)與分析

為驗(yàn)證本研究算法的有效性，采用高分辨率淺水多波束系統(tǒng)R2SONIC2024在我國(guó)某附近海域獲取的相鄰4條帶數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，該次實(shí)驗(yàn)多波束的采樣模式為平均強(qiáng)度采樣，波束個(gè)數(shù)為256個(gè)，波束開角120°，聲吶頻率400 kHz，實(shí)驗(yàn)區(qū)水深20～30 m。選取4條帶中的1條測(cè)線采用本算法進(jìn)行AR改正。

首先，將原始多波束反向散射強(qiáng)度數(shù)據(jù)進(jìn)行傳播損失改正、聲照區(qū)面積改正等一系列預(yù)處理改正；其次，將條帶左舷和右舷的角度響應(yīng)曲線進(jìn)行小波分解，分解參數(shù)采用2.1節(jié)中得到的最優(yōu)分解參數(shù)：小波基coif5、分解層數(shù)5層，分解后得到角度響應(yīng)曲線的長(zhǎng)波項(xiàng)和短波項(xiàng)；然后，繪制長(zhǎng)波項(xiàng)的三維概率密度圖，并通過二分法對(duì)條帶進(jìn)行同底質(zhì)區(qū)域的劃分，劃分結(jié)果如圖11所示，在劃分的條帶小區(qū)域內(nèi)基本屬于同一底質(zhì)類型。最后，在同底質(zhì)區(qū)域下，進(jìn)行角度響應(yīng)曲線長(zhǎng)波項(xiàng)的AR改正，將改正后的長(zhǎng)波項(xiàng)與短波項(xiàng)進(jìn)行信號(hào)重構(gòu)，得到改正后的聲吶圖像(圖12(c))。

為了對(duì)比現(xiàn)有方法，結(jié)果中添加了未進(jìn)行AR改正的原始聲吶圖像(圖12(a))和傳統(tǒng)Lambert法則改正后的聲吶圖像(圖12(b))。通過對(duì)比分析圖12(a)和圖12(c)可以看出，經(jīng)過AR改正后的多波束聲吶圖像底質(zhì)分布均勻，條帶中央?yún)^(qū)域“高亮”的灰度異?，F(xiàn)象被消除，同時(shí)邊緣區(qū)域的數(shù)據(jù)也被歸化到平均回波強(qiáng)度上，表明了本算法的有效性。對(duì)比圖12(b)，由于Lambert法則改正算法只是固定的改正模型，沒有考慮不同底質(zhì)類型AR曲線的不同，因此，在復(fù)雜海底底質(zhì)情況下，無法進(jìn)行自適應(yīng)改正，會(huì)出現(xiàn)過度改正或改正不徹底的現(xiàn)象，同時(shí)也驗(yàn)證了本方法的優(yōu)越性。

圖11 二分法劃分區(qū)域后的結(jié)果

以測(cè)線為單位依次消除AR影響，并在地理編碼下進(jìn)行拼接，形成整個(gè)區(qū)域的海底聲吶圖像，拼接結(jié)果如圖13(c)所示。與原始聲吶圖像拼接后的結(jié)果(圖13(a))相比，聲吶圖像中條帶中央和邊緣區(qū)域“高亮”現(xiàn)象得到了有效改善，整個(gè)測(cè)區(qū)灰度變化均勻，底質(zhì)分界清晰。圖13(b)為L(zhǎng)ambert法則改正后地理編碼的結(jié)果，可以看出，該方法在一定程度上減弱了AR的影響，但仍存有一定殘差，影響了聲吶圖像的質(zhì)量。為進(jìn)一步量化AR改正的效果，選取同底質(zhì)區(qū)域連續(xù)50 ping數(shù)據(jù)繪制反向散射強(qiáng)度散點(diǎn)圖并擬合趨勢(shì)線，獲取其MIC(maximal information coefficient)指數(shù)[18]。

圖12 AR改正前后聲吶圖像比較

圖13 整個(gè)測(cè)區(qū)的AR改正比較

MIC指數(shù)是衡量變量之間相關(guān)性的參數(shù)，MIC指數(shù)越小，說明反向散射強(qiáng)度與入射角的相關(guān)性越小，即AR改正效果越好。如圖14所示，本算法的MIC指數(shù)為0.092，相較于AR改正前的0.433降低了0.341，小于Lambert法則改正的0.274；且改正后的強(qiáng)度點(diǎn)較為集中，趨勢(shì)線基本呈直線，證明了本算法的有效性。

圖14 AR改正效果比較

4 結(jié)論

針對(duì)多波束聲吶圖像角度響應(yīng)改正，本研究綜合考慮了不同底質(zhì)類型對(duì)AR改正的影響，提出一種基于小波分析的多波束聲吶圖像改正方法，該方法利用信號(hào)分解與重構(gòu)的思想，將多波束角度響應(yīng)曲線分解為長(zhǎng)波項(xiàng)和短波項(xiàng)，利用長(zhǎng)波項(xiàng)和短波項(xiàng)各自的特點(diǎn)，繪制三維概率密度圖，得到小波分解的最優(yōu)參數(shù)和多波束條帶同底質(zhì)區(qū)域的劃分范圍，然后在同底質(zhì)區(qū)域下進(jìn)行AR改正得到改正后的聲吶圖像。實(shí)驗(yàn)結(jié)果表明，對(duì)多波束條帶同底質(zhì)區(qū)域劃分，進(jìn)行自適應(yīng)分區(qū)改正，改善了傳統(tǒng)改正模型適應(yīng)性差、改正不精確等問題，減弱了角度響應(yīng)對(duì)多波束聲吶圖像的影響，提高了聲吶圖像的質(zhì)量，有效適用于淺水海域地形平坦和底質(zhì)較復(fù)雜的區(qū)域，為多波束聲吶圖像處理提供了參考。