一種高解析度測深側掃聲納信號處理方法
2023-05-04 05:37:56 2
專利名稱：一種高解析度測深側掃聲納信號處理方法
技術領域：
本發明涉及一種高解析度的聲納信號處理方法，特別是涉及一種利用波達方向估計(Directions of Arrival，簡稱DOA)技術，獲得海底微地形地貌的高解析度測深側掃聲納信號處理方法。
背景技術：
二十世紀五十年代後期，人們開始研製側掃聲納。它的聲納陣在垂直平面內的波束角很寬，一般寬達100°，甚至更寬。在水平面內很窄，一般為1°左右。聲納陣發射聲脈衝，海底產生反向散射聲波，依時間的先後被聲納陣接收。聲納陣在前進過程中，不斷發射，不斷接收，就獲得了海底的聲圖，它反映海底的地貌狀態。側掃聲納研製成功後獲得了廣泛的應用，產生了顯著的經濟和社會效益。實際工作中，不但需要海底的地貌，同時還需要海底地形，往往將側掃聲納與多波束測深聲納一起使用。為了簡化設備，提高效益，人們自然想起使側掃聲納在獲得海底地貌的同時，也能獲得海底地形，這就誕生了測深側掃聲納。二十世紀九十年代後期，測深側掃聲納採用多條平行線陣，測量它們之間聲回波的相位差，一次發射可以獲得數百至上千個測深點，在解析度上它優於多波束測深系統，但是它有兩個重要的缺陷，長期影響它發展。首先，在聲納陣正下方附近，它的測深精度很差；其次，不能區分不同方向同時到達的回波，因此，在存在信道多途和地形比較複雜的水域測深側掃聲納不能正常工作。
2005年3月29日授權給朱維慶、劉曉東等人的美國專利US6873570 B2「Highresolution bathymetric sonar system and measuring method for measuring thephysiognomy of the seabed」中，公開了兩種技術，解決了上述測深側掃聲納的兩個重要缺陷。第一，聲納陣由多條相互平行的等間距線陣組成，相鄰線陣間距在一個λ和λ/2之間，λ為聲納中心頻率對應的水中聲波波長，聲納正下方附近的測深精度達到了高解析度數字測深儀的精度；第二，海底自動檢測——多子陣波達方向估計方法，它採用波達方向估計(DOA)技術，由聲納陣的時空相關函數矩陣提取聲波的信息，包括聲波入射角和振幅等。
聲納陣接收的信號中除了目標信號外，還有噪聲信號。因此聲納陣時空相關函數矩陣在泛函空間裡可以分解成信號子空間和噪聲子空間，兩個子空間相互垂直。對聲納陣時空相關函數矩陣的信號處理方法一般分兩大類一類是譜基方法，它包括噪聲子空間法，又稱零空間法，在小樣本，低信噪比和高信號相干性時，此類方法的性能下降；另一類是參量法，它包括信號子空間法。參量法的性能明顯優於譜基方法。美國專利US 6873570 B2中，公開了海底自動檢測——多子陣波達方向估計方法，它屬於信號子空間法，性能優良，克服了水聲信道和複雜海底地形引起的多途信號幹擾，提取了海底直達回波信號。但是它沒有公開如何解決平行線陣間的互耦(Mutul Coupling)影響問題。
2000年10月10日授權給P.Kraeutner等人的美國專利613041「Imagingmethods and apparatus using model-based array signal processing」中採用波達方向估計(DOA)技術中的零空間(即噪聲子空間)法，對聲納陣時空相關函數矩陣進行處理，獲得了比常規波束形成技術高的解析度。但是，該專利沒有給出正下方附近測深精度，也沒有研究平行線陣間的互耦影響。
文獻1W Xu，Stewart W K.的文章「Coherent source direction estimation forthree-row bathymetric sidescan sonars」，OCEANS』99，MTS/IEEE，Seattle，Washington，299-304.文中提出了相干源方向估計(Coherent Source Direction Estimator，簡稱CSDE)方法，對於三排線陣的測深側掃聲納和2個信號源的情況，進行了模擬計算，並與ESPRIT(Estimating Signal Parameters via Rotational InvarianceTechniques，簡稱ESPRIT)方法進行了比較。模擬計算表明，對於角間隔10°的高度相關的信號源，信噪比10db以上，CSDE工作良好。同時ESPRIT對於非相干信號源有好的魯棒性。該文獻建議把兩種方法聯合應用，代替已有的差分相位法。
綜上所述，現有的技術中存在的不足主要有1.高解析度測深側掃聲納對聲納陣的相位特性有嚴格要求，它是決定解析度多高的重要參數。多條平行線性陣間的互耦影響，尤其是相鄰線陣元間的互耦影響是重要相位誤差來源。如何由信號處理的方法減小這類誤差，現有技術中，包括本專利前面提到的兩個專利和一篇文獻，均沒有提及。
2.高解析度測深測掃聲納的接收聲納陣是由多條等間距平行線陣元組成，接收陣分為幾個子陣，每個子陣包括幾條線陣元，其中有一個好的組合，它的方位估計誤差最小。現有的技術中，包括本專利前面提到的兩個專利和一篇文獻，均沒有提及獲得最好組合的方法。
發明內容
本發明的目的在於克服已有技術的上述不足，提出了多子陣子空間擬合的信號處理方法。
為了達到上述目的，本發明採取如下技術方案一種高解析度測深側掃聲納信號處理方法，包括如下步驟a)確定每個子陣選取方式中的子陣；b)計算每個子陣選取方式的子陣對於目標的入射角；c)計算每個子陣選取方式的入射角的方差；d)選取入射角標準差最小的子陣選取方式作為聲納陣的工作方式。
在上述技術方案中，進一步地，步驟a)具體包括如下子步驟1)對於互耦矩陣自由度為h的聲納陣，h≥2，將首尾陣元去掉，得聲納陣輸出信號X，由下式表示X＝AGS+N(3)其中X＝[xhxh+1… xL-(h-1)z]T(4a)S＝[s1s2… sM]T(4b) G=diagq=-(h-1)h-1b|q|ej(q1)g(1)...q=-(h-1)h-1b|q|ej(qM)g(M),b0=1---(5)]]>其中，M為目標數，即信號源數目，L為陣元數，N為零均值空間白噪聲，X和N為L-2(h-1)維向量A為空間陣列的流型矩陣；信號S是由M個獨立高斯信號源構成；2)求聲納陣相關函數的估計 即
R^=1Nn=1NXnXnH---(6)]]>3)對聲納陣輸出信號的相關函數估計 進行特徵值分解R^=AGR^SGHAH+2I=U^S^SU^SH+U^N^NU^NH---(7)]]>其中 是信號相關函數的估計值，σ2是噪聲方差的估計值， 和 分別是信號的特徵向量和特徵值的估計值， 和 分別是噪聲的特徵向量和特徵值的估計值；上標H表示共軛轉置運算；I是單位矩陣。
4)求子陣運算J1＝[Il-10](l-1)×l，J2＝
(l-1)×l(8)其中，M＜l≤L-2(h-1)；在上述技術方案中，進一步地，步驟b)具體包括如下子步驟5)求子陣對應的特徵向量估計值；U^S1=J1U^S---(9)]]>U^S2=J2U^S]]>6)由上一步驟5)可得到估計值 和 得到 存在互耦時的多子陣子空間擬合算法的解^=(U^S1HU^S1)-1U^S1HU^S2---(10)]]>7)求 運算^=C-1^C,---(11)]]>其中^=diagei^1...ei^M---(12a)]]>C=GR^SGHAHU^S^S---(12b)]]>^S=^S-2I---(12c)]]>8)求聲波入射角估計值 運算^i=sin-1(^i/kd),i=1,2,...M---(13)]]>在上述技術方案中，進一步地，步驟c)具體包括如下子步驟9)在快拍數N大於100時，均勻線陣的多子陣子空間擬合算法的估計誤差 是聯合高斯分布，其均值和方差分別為E{^i-i}=0---(14a)]]>E((^i-i)2)=22NRe{(iHi)([P-1]i,i+2[P-1(AHA)-1P-1]i,i)}---(14b)]]>其中i＝1，……，M；P=GR^sGH---(15a)]]>iH=[(A1HA1)-1A1HFi]i---(15b)]]> 式中A1＝[Il0]A，Fi是(l-1)×l矩陣，[X]i(γ)表示矩陣X的第i行， 表示信號相關函數估計值；根據方差的平方根是標準差，計算所有子陣對的入射角估計標準差，對所有子陣對的入射角估計標準差進行算術平均，此平均值作為最終估計標準差；10)由步驟8)和步驟9)計算聲入射角估計值 的最終估計標準差；11)對所有可能的子陣選取方式，重複步驟4)-步驟10)，求出所有可能的子陣選取方式的入射角最終估計標準差。
在上述技術方案中，進一步地，步驟d)中比較上一步驟11)中得出的所有可能的子陣選取方式的入射角最終估計標準差，確定入射角最終估計標準差最小的那個子陣選取方式，選取該子陣選取方式作為聲納的工作方式。
與現有技術相比，本發明的優點在於(1)與一般的只從傳感器陣，或只從DOA信號處理單方面抑制陣元間互耦不同，本發明由兩方面共同抑制陣元間的互耦，在此基礎上就可研製性能良好的聲納。本發明對於聲納設計的指導意義在於一方面儘量將聲納陣的互耦矩陣自由度h減小，這樣需要刪除的陣元少，代價小；另一方面多子陣子空間擬合技術在陣元存在互耦時仍有良好性能。
(2)本發明表明，在聲納陣元數一定的情況下，陣元間有互耦或沒有互耦時，均存在一種子陣對的構成方式，它的相位估計標準差最小，採用這種子陣對的構成方式可以使得聲納更準確地估計目標方位。
圖1入射角估計標準差(度)與子陣數的關係均勻線陣陣元間不存在互耦，信噪比20dB，陣元數8個，快拍數100，信號源數1，曲線自上而下對應的入射角為80°，60°，40°，20°和0°；圖2入射角估計標準差(度)與子陣數的關係均勻線陣陣元間存在互耦，自由度3，耦合係數b0＝1，b1＝0.2exp(jπ/6)，b2＝0.05exp(jπ/5)，信噪比20dB，陣元數8個，快拍數100，信號源數1，曲線自上而下對應的入射角為80°，60°，40°，20°和0°；圖3為應用本發明的一實施例的測試程序的流程圖。
具體實施例方式
下面結合附圖和具體實施方式
對本發明作進一步詳細描述首先，本發明是提供一種DOA信號處理技術，克服線陣元間互耦影響引起的聲納陣相位誤差。
設一窄帶平面波s(t)入射到由L個陣元組成的均勻接收線陣上，假設只考慮相鄰陣元間的相互作用，則接收陣輸入s(t)與輸出xi，i＝1…L，的關係由下式表示x1x2x2......xL-2xL-1xL=1b10000.........0b11b1000.........00b11b100.........0............................................................0...0000b11b100...00000b11b10...000000b111eiei2......ei(L-3)ei(L-2)ei(L-1)gs(t)---(1)]]>其中，g(θ)是線陣元的指向性，它們是相同的，b1是耦合因子，表示相鄰線陣元間的互耦影響。φ＝kd sinθ，k是波數，d是相鄰線陣元的間距，θ是聲波入射角。上式是討論相鄰線陣元之間的作用，所以沒有考慮噪聲的影響。式(1)右邊第一個矩陣是互耦矩B，其中第一行非零元素的個數稱為互耦矩陣的自由度h，式(1)中h＝2。將矩陣B中的第一行和最後一行去掉，得矩陣B1，與之對應的線陣元輸出為xi，i＝2，……L-1。在i＝2，……L-1中任取兩個子陣，每個子陣包含的線陣元數相同，例如取i＝2，……L-2，和i＝3，……L-1。則得x3x4......xL-2xL-1=b11b10...............0b11b1....................................................................................b11b100000......b11b11eiei2......ei(L-3)ei(L-2)gs(t)ei=x2x3......xL-1xL-2ei---(2)]]>式(2)表明，不計入首尾陣元，把接收陣分成兩個相鄰子陣，每個子陣包含相同陣元數，則子陣的輸出信號之比為eiφ，與聲納陣不存在互耦時的結果一致。也就是說，由本發明的方法對於存在互耦和不存在互耦時eiφ的均值是相同的，它是無偏估計。由此φ得到的聲波入射角θ。聲波對海底的掠射角α＝θ+θm，θm為聲納陣安裝角，它等於聲納陣平面與空間垂直平面間的夾角。
進而，本發明提供一種多子陣DOA信號處理技術，在陣元數給定的情況下，獲得子陣數和每個子陣包含的線陣元數的優化組合，提高聲納的測深精度。例如類似於(2)式，取i＝2，…L-3，或3，…L-2或4，…L-1等。
一種高解析度測深側掃聲納信號處理方法，包括如下步驟1)對於互耦矩陣自由度為h的聲納陣，h≥2，將首尾陣元去掉，得聲納陣輸出信號X，由下式表示X＝AGS+N (3)其中X＝[xhxh+1…xL-(h-1)]T(4a)S＝[s1s2…sM]T(4b) G=diagq=-(h-1)h-1b|q|ej(q1)g(1)...q=-(h-1)h-1b|q|ej(qM)g(M),b0=1---(5)]]>
其中，M為目標數，即信號源數目，L為陣元數，N為零均值空間白噪聲，X和N為L-2(h-1)維向量；A為空間陣列的流型矩陣；信號S是由M個獨立高斯信號源構成；2)求聲納陣相關函數的估計 即R^=1Nn=1NXnXnH---(6)]]>3)對聲納陣輸出信號的相關函數估計 進行特徵值分解R^=AGR^SGHAH+2I=U^S^SU^SH+U^N^NU^NH---(7)]]>其中 是信號相關函數的估計值，σ2是噪聲方差的估計值， 和 分別是信號的特徵向量和特徵值的估計值， 和 分別是噪聲的特徵向量和特徵值的估計值；上標H表示共軛轉置運算；I是單位矩陣。
4)求子陣運算J1＝[Il-10](l-1)×l，J2＝
(l-1)×l(8)其中，M＜l≤L-2(h-1)；本步驟中的求子陣運算採用文獻1中公開的技術；5)求子陣對應的特徵向量估計值；U^S1=J1U^S---(9)]]>U^S2=J2U^S]]>6)由上一步驟5)得到估計值 和 得到 存在互耦時的多子陣子空間擬合算法的解^=(U^S1HU^S1)-1U^S1HU^S2---(10)]]>7)求 運算^=C-1^C,---(11)]]>其中^=diagei^1...ei^M---(12a)]]>C=GR^SGHAHU^S^S---(12b)]]>^S=^S-2I---(12c)]]>8)求聲波入射角估計值 運算^i=sin-1(^i/kd),i=1,2,...M---(13)]]>
9)在快拍數N大於100時，均勻線陣的多子陣子空間擬合算法的估計誤差 是一個零均值的聯合高斯分布，均勻線陣的多子陣子空間擬合算法的估計誤差 是聯合高斯分布，其均值和方差分別為E{^1-i}=0---(14a)]]>E((^i-i)2)=22NRe{(iHi)([P-1]i,i+2[P-1(AHA)-1P-1]i,i)}---(14b)]]>其中i＝1，……，M；P=GR^SGH---(15a)]]>iH=[(A1HA1)-1A1HFi]i---(15b)]]> 式中A1＝[Il0]A，Fi是(l-1)×l矩陣，[X]i(γ)表示矩陣X的第i行， 表示信號相關函數估計值； 為零均值表示在考慮到陣元間存在互耦後，與不存在互耦時一樣，多子陣子空間擬合算法是無偏估計；對所有子陣對的入射角估計標準差進行算術平均，此均值作為最終估計標準差；10)由步驟8)和步驟9)計算聲入射角估計值 的最終估計標準差；11)對所有可能的子陣選取方式，重複步驟4)-步驟10)；求出所有可能的子陣選取方式的入射角最終估計標準差；12)比較上一步驟11)中得出的所有可能的子陣選取方式的入射角估計標準差，確定存在入射角估計標準差最小的那個子陣選取方式，選取該子陣選取方式作為聲納的工作方式。
在上述步驟1)中，如果目標不相干，互耦矩陣自由度h＝2，將首尾陣元去掉，得聲納陣輸出X，由下式表示X＝AGS+N其中X＝[x2x3…xL-1]TS＝[s1s2…sM]T
G=diagq=-11b|q|ej(q1)g(1)...q=-11b|q|ej(qM)g(M),b0=1]]>其中，L為陣元數，N為零均值空間白噪聲，為L-2維矢量；A為空間陣列的流型矩陣。
將本發明的技術用於湖上試驗或海上試驗，可以採取如下步驟A)發射陣向水中發射聲脈衝信號；B)接收陣接收水體和海底的反向散射回波信號；C)對接收的回波進行濾波和採樣，得到由式(3)表示的聲納陣輸出信號X；D)重複上述實施例的步驟4)-10)；求出所可能的子陣選取方式的入射角估計標準差；E)比較上一步驟D)中得出的所有可能的子陣選取方式的入射角估計標準差，選取入射角估計標準差最小的那種子陣選取方式作為聲納的工作方式。
本發明的聲納信號處理方法可以作為專用測量程序應用於聲納系統，一般將該專用測量程序裝載在聲納系統的主控計算機中。如圖3所示的程序流程圖，此處基於上述本發明的實施例，圖中的符號如 等等的意義與上述實施例中相同，程序執行過程如下步驟401是開始步驟，啟動計算機的程序，使聲納處於工作狀態。
在步驟402和403中，對聲納系統的軟體和硬體進行初始化。
在步驟404中，主控計算機生成發射信號。
在步驟405中，向海底發射聲脈衝到流體介質中，例如海水中。
在步驟406中，接收從流體介質反向散射的回波信號。
在步驟407中，對回波信號進行解調濾波。
在步驟408中，將回波信號從模擬信號轉換成數位訊號，並對每個回波信號逐一執行步驟409～步驟413。
在步驟409中，計算聲納陣相關函數估計 在步驟410中，執行 特徵值分解運算。
在步驟411中，執行子陣運算，計算子陣對應的特徵向量估計值。
在步驟412中，計算 在步驟413中，計算 在步驟414中，計算 的方差。
在步驟415中，計算 的標準差。
在步驟416中，計算聲波入射角估計值 在步驟417中，計算聲波對海底的掠射角估計值 在步驟418中，存儲 和 的標準差。同時反饋到步驟411，重複步驟411至步驟417的計算，直至所有可能的子陣方式選取完畢，並選取 和 的標準差最小的子陣選取方式作為聲納的工作方式。
需要說明的是，上述步驟401-步驟408均採用本領域技術人員公知的技術，故本發明中沒有對其進行詳細說明。其餘步驟409-步驟418採用前述實施例中的公式進行計算。
下面以應用本實施例得到的一些具體數據為例來說明本發明的效果。
執行本實施例中的步驟1)-步驟10)，求出所有可能的子陣選取方式的入射角估計標準差，典型結果圖示於圖1和圖2。圖1和圖2顯示的是入射角估計標準差與子陣數的關係。圖1中，均勻線陣陣元間不存在耦合，互耦矩陣的自由度h＝1，陣元數為8，信號源數M為1，g(θ)＝1，信源入射角分別為80°，60°，40°，20°和0°。圖1表明入射角愈大，入射角估計標準差愈大；存在某子陣數，它的入射角估計標準差最小。圖2表示的是均勻線陣陣元間存在耦合，互耦矩陣的自由度h＝3，b0＝1，b1＝0.2ejπ/6，b2＝0.05ejπ/5，其它參數與圖1相同。對G進行了歸一化，其模的最大值為1。其中陣元數8個是指實際參與計算的陣元數，不包括去掉的四個首尾陣元。由圖2可知由於存在互耦，入射角估計標準差有些增加，入射角愈大，入射角估計標準差愈大；與圖1一樣存在某子陣數，它的入射角估計標準差最小。
對於兩個目標的情況即M＝2，與一個目標的結果相似。
最後所應說明的是，以上實施例僅用以說明本發明的技術方案而非限制。儘管參照實施例對本發明進行了詳細說明，本領域的普通技術人員應當理解，對本發明的技術方案進行修改或者等同替換，都不脫離本發明技術方案的精神和範圍，其均應涵蓋在本發明的權利要求範圍當中。
權利要求
1.一種高解析度測深側掃聲納信號處理方法，包括如下步驟a)確定每個子陣選取方式中的子陣；b)計算每個子陣選取方式的子陣對於目標的入射角；c)計算每個子陣選取方式的入射角的方差；d)選取入射角標準差最小的子陣選取方式作為聲納陣的工作方式。
2.根據權利要求1所述高解析度測深側掃聲納信號處理方法，其特徵在於，步驟a)具體包括如下子步驟1)對於互耦矩陣自由度為h的聲納陣，h≥2，將首尾陣元去掉，得聲納陣輸出信號X，由下式表示X＝AGS+N其中X＝[xhxh+1…xL-(h-1)]TS＝[s1s2…sM]T G=diagq=-(h-1)h-1b|q|ej(qi)g(i)q=-(h-1)h-1b|q|ej(qM)g(M),b0=1]]>其中，M為目標數，即信號源數目，L為陣元數，N為零均值空間白噪聲，X和N為L-2(h-1)維向量；A為空間陣列的流型矩陣；信號S是由M個獨立高斯信號源構成；2)求聲納陣相關函數的估計 即R^=1Nn=1NXnXnH]]>3)對聲納陣輸出信號的相關函數估計 進行特徵值分解R^=AGR^SGHAH+2I=U^S^SU^SH+U^N^NU^NH]]>其中 是信號相關函數的估計值，σ2是噪聲方差的估計值， 和 分別是信號的特徵向量和特徵值的估計值， 和 分別是噪聲的特徵向量和特徵值的估計值；上標H表示共軛轉置運算；I是單位矩陣；4)求子陣運算J1＝[Il-10](l-1)×l，J2＝
(l-1)×l其中，M＜l≤L-2(h-1)。
3.根據權利要求2所述高解析度測深側掃聲納信號處理方法，其特徵在於，步驟b)具體包括如下子步驟5)求子陣對應的特徵向量估計值；U^S1=J1U^S]]>U^S2=J2U^S]]>6)由上一步驟5)可得到估計值 和 得到 存在互耦時的多子陣子空間擬合算法的解^=(U^S1HU^S1)-1U^S1HU^S2]]>7)求 運算^=C-1^C,]]>其中^=diag{ei^1ei^l1}]]>C=GR^SGHAHU^S^S]]>^S=^S-2I]]>8)求聲波入射角估計值 運算^i=sin-1(^i/kd)]]>i＝1，2，…M。
4.根據權利要求3所述高解析度測深側掃聲納信號處理方法，其特徵在於，所述步驟c)具體包括如下子步驟9)在快拍數N大於100時，均勻線陣的多子陣子空間擬合算法的估計誤差 是聯合高斯分布，其均值和方差分別為E{^i-i}=0]]>E((^i-i)2)=22NRe{(iHi)([P-1]i,i+2[P-1(AHA)-1P-1]l,i)}]]>其中i＝1，……，M；P=GR^sGH]]>lH=[(A1HA1)-1A1HFl]i]]> 式中A1＝[Il0]A，Fl是(l-1)×l矩陣，[X]i(γ)表示矩陣X的第i行， 表示信號相關函數估計值；對所有子陣對的入射角估計標準差進行算術平均，此平均值作為最終估計標準差；10)由步驟8)和步驟9)計算聲入射角估計值 的最終估計標準差；11)對所有可能的子陣選取方式，重複步驟4)-步驟10)，求出所有可能的子陣選取方式的入射角最終估計標準差。
5.根據權利要求4所述高解析度測深側掃聲納信號處理方法，其特徵在於，在步驟d)中比較上一步驟11)中得出的所有可能的子陣選取方式的入射角最終估計標準差，確定入射角最終估計標準差最小的那個子陣選取方式，選取該子陣選取方式作為聲納的工作方式。
全文摘要
本發明涉及一種利用波達方向估計技術，獲得海底微地形地貌的高解析度測深側掃聲納信號處理方法，稱它為多子陣子空間擬合法。本發明表明，在聲納陣元數一定的情況下，陣元間有互耦或沒有互耦時，均存在一種子陣對的構成方式，它的相位估計標準差最小。本發明克服線陣元間互耦影響引起的聲納陣相位誤差，在陣元數給定的情況下，獲得子陣數和每個子陣包含的線陣元數的優化組合，提高聲納的測深精度。
文檔編號G01S15/00GK1900738SQ200610083128
公開日2007年1月24日 申請日期2006年6月5日 優先權日2005年7月22日
發明者朱維慶, 劉曉東, 徐文 申請人:中國科學院聲學研究所