一種基於互譜優化的互譜廣義逆波束形成方法與流程
2023-09-27 06:41:05
本發明涉及一種聲源識別算法，主要涉及運用一種基於傳聲器陣列的互譜優化廣義逆波束形成算法，在消除背景噪聲的基礎上來精確識別定位聲源。
背景技術：
現在基於傳聲器陣列的聲源識別方法主要有兩種，波束形成與聲全息。兩者各有優劣勢，其中波束形成的中遠距離聲源識別性能相對於聲全息較好，而低頻性能相對於聲全息較差。
波束形成是近年來快速發展的一種聲源識別方法，它基於規則或者非規則傳聲器陣列(如十字軸和輪型陣列)進行信號採集，然後利用如簡單延時求和算法進行聲場重構。對比傳統方法如主觀評價法、選擇隔離法等，基於多通道信號處理技術的波束形成具有測量方便、重構聲場直觀可視化等優點。相比聲強測量法，波束形成具有測量效率高和應用環境寬等優點。相對於近年來的發展的聲全息方法，波束形成在中高頻具有較高的空間解析度，而且能夠實現中遠距離的聲源重構，並且其能夠適應不同的運行工況。由於基于波束形成的噪聲源識別方法具有靈活的動態顯示範圍、幹擾能力強、中高頻高解析度高等優點，因此近年受到各國研究人員和學者的高度關注。使得其應用範圍也不斷地拓展，從航空、聲吶等軍事領域，逐漸發展到車輛噪聲測試等民用領域。與此同時衍生出了各種不同的波束形成算法，如傳統的互譜波束形成、反卷積算法以及廣義逆波束形成算法，各種不斷地推出的各色的波束形成算法各具特色又相互聯繫，不斷推進該方法的應用便利性與實用性。
現有技術中，存在一種基於傳聲器陣列的聲源識別方法，在傳統的等效源算法的基礎上公開了一種寬帶聲全息算法。其在靠近聲源面較近的等效源平面構造一系列等效源，讓這些等效源來替代實際聲源分布。通過聲源識別理論，建立一個等效源到傳聲器陣列的聲學傳遞方程。並通過傳統的梯度下降算法求解該傳遞方程，從而迭代求解等效源強度。然後通過該等效源強度，基於聲傳播過程重構在等效源面與陣列面之間的聲源分布。而且該方法在等效源求解過程中，引入一個迭代濾波過程，從而保證迭代過程的收斂效率與精準度。該算法由於引入梯度下降算法和濾波過程，能夠提高傳統等效源方法在中低頻的解析度。同時該算法也能夠實現與其他聲全息算法之間的結合或者切換，從而實現在中低頻聲源的重構。但該算法基於傳統的等效源算法重構，其在中低頻聲源以及較近的測量距離中具有較高的空間解析度，但是對於高頻聲源或者較遠距離測量，其聲源識別的空間解析度將較低。
技術實現要素：
本發明目的在於提供一種聲源重構精度高、魯棒性較好的聲源識別算法，其建構在一般廣義逆和互譜波束形成的基礎上，具有較高的重建精度，並充分結合聲壓互譜矩陣的優化，抗噪性能較強。
為實現本發明目的而採用的技術方案是這樣的，一種基於互譜優化的互譜廣義逆波束形成方法，其特徵在於，
搭建一個聲源重構系統：
該系統包括位於聲源附近的攝像頭、傳聲器陣列、多通道的信號採集器以及後處理計算機。所述攝像頭用於採集整個聲場的空間影像；所述傳聲器陣列用於測量聲場信號；所述多通道的信號採集器用於測量時域模擬信號，並將其轉化成數位訊號傳輸給計算機進行後處理；所述計算機的主要功能主要是存儲以及後處理採集到的聲場信號，並能高效顯示重構聲場分布情況；
波束形成過程包括以下步驟：
1)基於廣義逆波束形成來計算轉向矩陣：
聲源面布置n個單位聲源，則在傳聲器陣列中的傳聲器ri點處的測量得到的聲壓信號如下：
上式中，為聲場中ri處的單位激勵在r0處的複數響應，q(ri)為假設的單極子聲源強度。式(1)表達的是由聲源到傳聲器陣列之間的傳遞關係，為後續分析必要，可以將其轉化成向量—矩陣形式：
p＝aq(2)
其中a為n×m聲場傳遞矩陣，其元素成分由式(1)決定。p為m維測量聲壓向量，其元素為對應傳感器在單一頻率f對應的頻域響應。q為n維聲源強度列向量，其元素成分代表聲源平面中對應點處的強度分布。
為了求解重構聲場分布，可以將等式(2)轉化為該最小化問題：
其中β為吉鴻諾夫正則化參數，l為n×n正則化對角矩陣。通過迭代算法求解該最小化問題(3)，最終可以得到第(n)步下的轉向矩陣的表達式：
w(n)＝l(n)b(n)(4)
式(4)中b(n)為n×m矩陣，其行向量可以表達成如下形式：
b(i,:)(n+1)＝g(:,i)half/(g(:,i)halfg(:,i))(5)
式(5)其主要功能是對傳遞矩陣進行歸一化，從而保證輸出的準確性。上式中上標h表示共軛轉置，g(:,i)為歸一化後的傳遞矩陣g的第i列向量，其中g的表達式如下：
g＝al(n)(6)
式(5)中alf為可逆的m×m維矩陣，其表達式有如下形式：
alf＝(ggh+βi)-1(7)
上式(7)中正則化參數β，可以通過l曲線或者gvc法得到。
l(n)＝diag(|sp|/||sp||∞)(8)
上式(8)中diag為將向量轉化成對角矩陣形式，|*|為取向量的絕對值，||*||∞為求向量的無窮範數。sp為歸一化後的廣義逆波束形成輸出：
sp＝w(n)p(9)
在求解波束輸出sp的求解過程中也加入如下濾波過程：
式(10)中，(|sp|)max取向量|sp|中的最大元素，其中dk隨迭代次數(n)的增加而增加，從而使得濾波門限值隨迭代次數增加而逐漸減小：
dk＝0.2+0.5n(11)
求解得到基於廣義逆波束形成的一個轉向矩陣w(n)；
2)對聲壓互譜矩陣進行優化：
基於傳感器陣列測量得到的m維聲壓向量p，求得一個m×m維聲壓互譜矩陣，其表達式如下：
c＝pph(12)
優化後的聲壓互譜矩陣形式表達如下：
c優＝c-diag(x)(13)
其中，c優為消除掉對角元素的矩陣，diag(x)為一對角陣，x為噪聲分量；
x為的元素x(i,1)被限定在了[0,c(i,i)]之間，這就保持了去對角化元素中的一些優勢的同時又能夠降低通道自噪聲。為得到噪聲分量x，可以通過迭代優化算法來優化聲壓互譜矩陣。假設一系列歸一化單位向量uk＝ek(k＝1,2,...,m),則有如下不等式成立：
ukh(c-diag(x))uk≥0(14)
由上式可以推導出如下線性不等式：
上式中d(k,i)＝|u(i,k)|2(其中u(i,k)為向量uk的第i號元素)，d(k,1)＝ukhcuk。通過上式可得到一個關於噪聲分量x的一個約束條件dx≤d，由於聲壓互譜矩陣的優化的核心思想就是最小化噪聲分量-x，所以可得到一個目標函數minf＝-x。將以上最小化噪聲分量過程總結如下，即可以一個線性規劃問題：
求解最小化噪聲分量-x，並輸出優化後的聲壓互譜矩陣c優。
3)依照所述轉向矩陣和優化後的聲壓互譜矩陣重新構造波束輸出：波束輸出形式如下：
s優＝w(n)c優(w(n))h(18)
上式中優化後的聲壓互譜矩陣c優可以替代成如下形式：
s優＝w(n)(c-diag(x))(w(n))h(19)
將上式(19)取對角元素並進行開方，最終得到基於優化互譜矩陣的互譜廣義逆波束形成的輸出如下：
s輸出＝(diag(s優))1/2(20)。
本發明的技術效果是毋庸置疑的，首先基於傳統廣義逆波束形成構造轉向矩陣。其次是結合線性優化方法對聲壓互譜矩陣進行優化。最後通過求解得到的轉向矩陣和優化的聲壓互譜矩陣，結合傳統的互譜矩陣波束形成思想，構造新的基於優化互譜矩陣的互譜廣義逆波束形成算法。該算法在運用廣義逆波束形成高解析度優勢的同時，能夠充分利用基於互譜優化的互譜波束形成的抗噪性能。
附圖說明
圖1聲源重構示意圖；
圖2聲壓互譜矩陣優化流程示意圖。
具體實施方式
下面結合實施例對本發明作進一步說明，但不應該理解為本發明上述主題範圍僅限於下述實施例。在不脫離本發明上述技術思想的情況下，根據本領域普通技術知識和慣用手段，做出各種替換和變更，均應包括在本發明的保護範圍內。
由於本實施例涉及的是一個基於聲源識別的波束形成算法，其建立在聲源識別系統的基礎上，所以首先將簡要介紹該聲源識別系統。其中圖1描述了一個聲源重構系統示意圖，該系統包括位於聲源附近的攝像頭、傳聲器陣列、多通道的信號採集器以及後處理計算機。
所述攝像頭用於採集整個聲場的空間影像；，通過波束形成算法計算得到的波束輸出與空間影響相重合即可以得到整個聲場的聲源分布情況。攝像頭採用一般的數字攝像頭即可，可通過一般的usb接口傳輸到計算機。
所述傳聲器陣列用於測量聲場信號；，其結構包括了聲壓或者聲速傳感器、傳聲器支架以及相關線束。在傳聲器陣列中，傳感器的分布可以是規則化的或者是非規則化一、二、三維空間分布，這個根據實際測量需求而選擇。傳感器數目的選擇和傳聲器間距的選擇則根據需要測量的頻率範圍而決定。
所述多通道的信號採集器用於測量時域模擬信號，並將其轉化成數位訊號傳輸給計算機進行後處理；其應該包括多通道的前置放大器以及多通道的具有高精度的模擬/數字轉換器，甚至應該包括相關的濾波器，從而保證輸出到計算機的數位訊號足夠精確。
所述計算機的主要功能主要是存儲以及後處理採集到的聲場信號，並能高效顯示重構聲場分布情況；。其中計算機要求的內存512mb以上，建議是2g，cpu主頻2.0g以上比較好，這樣能夠充分保證計算的高效性。
波束形成過程包括以下步驟：
假設以上聲源重構系統中包括m個聲壓傳聲器，通過多通道信號採集器採集到的時域信號傳輸並儲存到計算機，並將採集得到的時域信號通過快速傅立葉變換轉化到頻域中。實施例中，波束形成算法主要涉及中高頻聲源的重構問題，所以我們將得到的陣列頻域信號中的某個單一頻率f成分所對應的複數信號提取出來，並將其轉化為後面算法所要處理的測量聲壓信號p。
由於該算法是基於廣義逆波束形成，互譜矩陣優化以及互譜波束形成所開發。
1)基於廣義逆波束形成來計算轉向矩陣：
，其中的轉向矩陣為歸一化之後的轉向矩陣，其基本作用於傳統互譜波束形成中的轉向向量類似。為求解該歸一化的轉向矩陣，聲源面布置n個單位聲源，則在傳聲器陣列中的傳聲器ri點處的測量得到的聲壓信號(頻域單頻成分)如下：
上式中，為聲場中ri處的單位激勵在r0處的複數響應，q(ri)為假設的單極子聲源強度。式(1)表達的是由聲源到傳聲器陣列之間的傳遞關係，為後續分析必要，可以將其轉化成向量—矩陣形式：
p＝aq(2)
其中a為n×m聲場傳遞矩陣，其元素成分由式(1)決定。p為m維測量聲壓向量，其元素為對應傳感器在單一頻率f對應的頻域響應。q為n維聲源強度列向量，其元素成分代表聲源平面中對應點處的強度分布。
為了求解重構聲場分布，可以將等式(2)轉化為該最小化問題：
其中β為吉鴻諾夫正則化參數，l為n×n正則化對角矩陣。通過迭代算法求解該最小化問題(3)，最終可以得到第(n)步下的轉向矩陣的表達式：
w(n)＝l(n)b(n)(4)
式(4)中b(n)為n×m矩陣，其行向量可以表達成如下形式：
b(i,:)(n+1)＝g(:,i)half/(g(:,i)halfg(:,i))(5)
式(5)其主要功能是對傳遞矩陣進行歸一化，從而保證輸出的準確性。上式中上標h表示共軛轉置，g(:,i)為歸一化後的傳遞矩陣g的第i列向量，其中g的表達式如下：
g＝al(n)(6)
式(5)中alf為可逆的m×m維矩陣，其表達式有如下形式：
alf＝(ggh+βi)-1(7)
上式(7)中正則化參數β，可以通過l曲線或者gvc法得到。本專利中推薦使用的正則化參數β為ggh最大特徵值的倍(其中f為聚焦聲源頻率)。
式(6)中的正則化矩陣l(n)為一重要參數，在迭代過程中隨廣義逆波束形成的輸出的變化而變化，從而保證重構過程的精確性。一般可以將其定義為如下形式：
l(n)＝diag(|sp|/||sp||∞)(8)
上式(8)中diag為將向量轉化成對角矩陣形式，|*|為取向量的絕對值，||*||∞為求向量的無窮範數。sp為歸一化後的廣義逆波束形成輸出：
sp＝w(n)p(9)
在傳統的廣義逆波束形成輸出中，由於噪聲的影響，將在非聲源方向出現幹擾旁瓣。為進一步抑制這些幹擾旁瓣的影響，從而實現更為精確聲源定位，因此在求解波束輸出sp的求解過程中也加入如下濾波過程：
該濾波過程的核心思想主要是通過假定一個隨迭代變化門限值，在迭代過程中不斷地將波束輸出向量中小於該閾值的元素強制置為零，進而抑制噪聲的影響。式(10)中(|sp|)max取向量|sp|中的最大元素，其中dk隨迭代次數(n)的增加而增加，從而使得濾波門限值隨迭代次數增加而逐漸減小：
dk＝0.2+0.5n(11)
值得說明的是，步驟1)最終目的不是為了求解波束輸出，而是為了求解得到基於廣義逆波束形成的一個轉向矩陣w(n)，即為求解到如式(2)的轉向矩陣。
2)對聲壓互譜矩陣進行優化：
基於傳感器陣列測量得到的m維聲壓向量p，求得一個m×m維聲壓互譜矩陣，其表達式如下：
c＝pph(12)
傳統基於聲壓互譜矩陣的波束形成為了去除通道自噪聲的影響，通常通過消除掉互譜矩陣的對角元素，即令c＝c-(c)對角。由於聲壓互譜矩陣中的對角元素是同通道的功率譜，同通道中的信號的功率譜可能引入通道自噪聲，所以通過去去除對角元素，確實能夠提高波束輸出的抗噪能力。不過該方法也具有較大的副作用，如導致波束輸出譜中會出現負功率，甚至可能導致聲源定位出現較大的偏差。當然實際測量中我們當然想要在保持去對角化方法抗噪優勢的同時，又能夠減小定位偏差等負面情況的影響。針對這種情況，我們基於互譜矩陣的廣義逆波束形成輸出，對聲壓互譜矩陣進行優化改進。
假設聲壓互譜矩陣中的對角元素中包括了有功分量和噪聲分量，如果在聲壓互譜矩陣的計算過程中能夠在保存有功分量的基礎上剔除噪聲分量，則此時就能夠充分運用優化後的聲壓互譜矩陣進行聲場重構。根據以上假設，我們可以將優化後的聲壓互譜矩陣形式表達如下：
c優＝c-diag(x)(13)
所以可通過優化方法剔除噪聲分量而保持有功分量。上式類似於聲壓互譜矩陣去掉對角元素，其中，c為消除掉對角元素的矩陣，diag(x)為一對角陣，其代表噪聲分量，是需要在對角元素中去除的部分，該對角陣為其對角元素對應向量x中的元素，例如x為的元素x(i,1)被限定在了[0,c(i,i)]之間，c(i,i)指的是矩陣c第i個對角元素，這就保持了去對角化元素中的一些優勢的同時又能夠降低通道自噪聲。為得到噪聲分量x，我們可以通過迭代優化算法來優化聲壓互譜矩陣。假設一系列歸一化單位向量
uk＝ek(k＝1,2,...,m),則有如下不等式成立：
ukh(c-diag(x))uk≥0(14)
由上式可以推導出如下線性不等式：
上式中d(k,i)＝|u(i,k)|2(其中u(i,k)為向量uk的第i號元素)，d(k,1)＝ukhcuk。通過上式可得到一個關於噪聲分量x的一個約束條件dx≤d，(d為一矩陣，其元素值由d(k,i)＝|u(i,k)|2決定；d為一向量，其元素值由d(k,1)＝ukhcuk決定；)由於聲壓互譜矩陣的優化的核心思想就是最小化噪聲分量-x，所以可得到一個目標函數minf＝-x。將以上最小化噪聲分量過程總結如下，即可以一個線性規劃問題：
當然求解線性規劃問題(17)的方法很多，我們可以根據不同的需要選擇不同的優化方法。最後通過選取一般的線性迭代算法求解最小化噪聲分量-x，並輸出優化後的聲壓互譜矩陣c優。
其中將聲壓互譜矩陣的優化過程總結如圖2所示，其主要包括四個過程，初始化、第一循環、第二循環以及第三循環。首先開始進行的初始化過程主要是為了給一些變量賦值。對於優化的聲壓互譜矩陣c優的初始值我們假定為沒有進行任何處理的聲壓互譜矩陣，並令u為由一系列的歸一化單位列向量組成的矩陣。其中b為向量x的幅值因子，x下和x上分別代表向量x在迭代求解過程中的下上限值，主要是為了約束其收斂的範圍，保證求解的準確性，其中t1和t2為約束變量，主要是為了提高該優化過程的收斂速度，其中0≤t1<t2≤1，最常規的做法可以將t1令為0且將t2設為1。通過設置t1和t2數值，能夠有效地在優化聲壓互譜矩陣和去對角化方法之間切換。其中迭代步數n1是一個非常重要的參數，其取值隨不同的測量環境的不同而不同，理論上其取值越大越好，不過其取值過大也會造成計算效率的降低，所以實際測量中必須要選擇可以接受的合理數值。
其中優化過程的第二和第三循環主要是為了計算約束條件中的相關參數，其中第三循環通過u去計算矩陣d。第二循環主要是通過矩陣u和聲壓互譜矩陣c去計算向量d。其中最重要的屬於第一循環，第一循環的開始部分主要是計算第二循環的循環次數，並根據矩陣u的列的維數重新賦予矩陣d以及向量d應有的維度以及初始值。當第二和第三循環完成時，已經計算出新的矩陣d以及向量d。此時第一循環的後半部分將通過常規的線性規劃方法，在約束住噪聲分量x的上下限的基礎上進行優化求解，從而計算出新的噪聲分量x。再將噪聲分量向量x轉化成對角矩陣，然後通過聲壓互譜矩陣減去該對角矩陣就可以得到優化後的聲壓互譜矩陣，這相當於在聲壓互譜矩陣的對角元素中減去噪聲分量，從而提高波束輸出結果的抗噪性能。當然在第一循環的最後還有一個重要的步奏，就是將上一步得到的矩陣u的列向量與優化後的聲壓互譜矩陣的特徵向量做併集，並將得到併集重新賦值給新矩陣u，這樣隨著迭代次數的增加，矩陣u的列數將逐漸增加。
完成以上循環過程以後，則可以輸出優化後聲壓互譜矩陣c優。該優化後的聲壓互譜矩陣可以運用到基於互譜的波束形成，從而保證輸出的精準性，以下步驟3)將通過基於廣義逆的互譜波束形成進行輸出，從而進行聲源定位識別。
3)依照所述轉向矩陣和優化後的聲壓互譜矩陣重新構造波束輸出：依照傳統互譜波束形成理論，其波束輸出可以表達成如下形式：
sout＝wcwh(17)
其中的w為由轉向向量組成的轉向矩陣，其列向量為轉向向量。其中轉向向量主要作用是通過將波束聚焦到聲源方向，從而自聲源方向的波束輸出得到聚焦而增加，在非聲源方向的波束輸出從而受到抑制，這樣通過觀察聚焦的主瓣方向即可以識別定位聲源方向。所以基於以上思想，在得到的基於廣義逆波束形成轉向矩陣與優化的聲壓互譜矩陣基礎上，可以重新構造一個基於廣義逆波束形成的互譜波束形成，其波束輸出形式如下：
s優＝w(n)c優(w(n))h(18)
上式中優化後的聲壓互譜矩陣c優可以替代成如下形式：
s優＝w(n)(c-diag(x))(w(n))h(19)
上式w(n)就是式(4)中求得的基於廣義逆波束形成求得的轉向矩陣，通過該轉向矩陣將能夠聚焦波束輸出方向，從而實現聲源識別定位。上式中輸出的s優形式上類似於傳統互譜波束形成的輸出形式，但是由於其轉向矩陣經過了廣義逆波束形成進行了優化，進一步提高了該算法的重建精度。同時由於引入了優化後的聲壓互譜矩陣，進而提高了該聲源重構方法的抗噪性能。為簡化上式的輸出，將上式(19)取對角元素並進行開方，最終得到基於優化互譜矩陣的互譜廣義逆波束形成的輸出如下：
s輸出＝(diag(s優))1/2(20)。