低秩及幾何約束的嵌入式陣列膛口波測(cè)向方法

王旭， 趙兆， 許志勇

(南京理工大學(xué) 電子工程與光電技術(shù)學(xué)院, 江蘇 南京 210094)

0 引言

單兵便攜式聲探測(cè)系統(tǒng)、機(jī)器人、無(wú)人平臺(tái)等現(xiàn)代偵察裝備對(duì)小型麥克風(fēng)陣列大氣被動(dòng)聲探測(cè)系統(tǒng)的需求日益增長(zhǎng)，相關(guān)功能涵蓋對(duì)戰(zhàn)場(chǎng)聲源目標(biāo)的檢測(cè)、定位、識(shí)別等[1-3]。其中，便攜式槍聲探測(cè)系統(tǒng)因能快速捕獲突發(fā)槍擊事件并確定其射擊點(diǎn)位置信息[4]，在近年國(guó)際反恐戰(zhàn)爭(zhēng)中戰(zhàn)果顯著，已被各國(guó)重視，具有重要的軍事國(guó)防、公共安全價(jià)值。

彈頭出膛時(shí)形成的膛口沖擊波[5]是主要的槍聲定位信息。在眾多的麥克風(fēng)陣列聲源測(cè)向方法中[6-8]，基于陣元間波達(dá)時(shí)間差(TDOA)測(cè)量的最小二乘(LS)估計(jì)方法因?qū)崿F(xiàn)復(fù)雜度低、環(huán)境適應(yīng)性好等優(yōu)點(diǎn)而廣泛應(yīng)用于包括槍聲定位在內(nèi)的各種實(shí)時(shí)聲探測(cè)任務(wù)中[8-9]。然而，受外形結(jié)構(gòu)、空間尺寸、使用運(yùn)輸?shù)葪l件限制，便攜式槍聲探測(cè)系統(tǒng)所用的麥克風(fēng)陣列通常采用嵌入式安裝方式，導(dǎo)致部分陣元由于外殼遮擋而與目標(biāo)處于非通視(NLOS)狀態(tài)?？紤]到遠(yuǎn)距離傳播的槍聲膛口波信號(hào)能量主要集中在數(shù)百赫茲的低頻范圍，NLOS陣元相關(guān)的實(shí)際TDOA主要受殼體衍射傳播效應(yīng)影響，不再與陣元間距和目標(biāo)來(lái)向保持理想全通視條件下的理論對(duì)應(yīng)關(guān)系，進(jìn)而造成基于LS的膛口波方向估計(jì)結(jié)果出現(xiàn)明顯偏差。

有關(guān)NLOS聲源目標(biāo)定位問(wèn)題，許多文獻(xiàn)主要關(guān)注聲探測(cè)系統(tǒng)整體因障礙阻隔而無(wú)法直視目標(biāo)時(shí)的解決方案[10-12]，一般不涉及殼體遮擋對(duì)嵌入式陣列目標(biāo)方向估計(jì)的影響及其處理方法。而現(xiàn)有針對(duì)僅有部分NLOS陣元的研究通常采用球麥克風(fēng)陣列聲場(chǎng)模型來(lái)獲得精確的聲源目標(biāo)方位估計(jì)[13]，但這種基于模型的衍射傳播效應(yīng)抑制方法一般較為復(fù)雜，難以實(shí)時(shí)實(shí)現(xiàn)；同時(shí)受成本、結(jié)構(gòu)、使用方式等制約，便攜式槍聲探測(cè)系統(tǒng)通常陣元數(shù)很少，很難嚴(yán)格滿足均勻球形陣列的模型要求，實(shí)踐中也難以通過(guò)建模查表的方式來(lái)解決，而現(xiàn)有基于數(shù)據(jù)的方案雖然可以通過(guò)求解基于低秩和異常值個(gè)數(shù)約束的非凸優(yōu)化問(wèn)題有效校正陣元間TDOA測(cè)量噪聲及個(gè)別異常值[14]，但卻無(wú)法校正由于殼體遮擋衍射傳播所造成的實(shí)際TDOA偏差，因?yàn)檫@種偏差既非測(cè)量誤差也非異常值。

針對(duì)上述問(wèn)題，本文以一種小型嵌入式四元平面麥克風(fēng)陣列系統(tǒng)作為研究對(duì)象，通過(guò)在施加低秩約束的同時(shí)再增加幾何約束，實(shí)現(xiàn)對(duì)實(shí)測(cè)TDOA矩陣中隨機(jī)測(cè)量誤差和殼體遮擋衍射傳播偏差的聯(lián)合抑制，同時(shí)選取距離目標(biāo)最近的通視麥克風(fēng)作為參考陣元，形成了一種新的槍聲膛口波測(cè)向算法，明顯減小了LS測(cè)向偏差的統(tǒng)計(jì)均值與方差。

1 問(wèn)題描述

便攜式槍聲探測(cè)系統(tǒng)采用360°全方位遠(yuǎn)距離探測(cè)方式，目標(biāo)通常遠(yuǎn)至數(shù)百米甚至上千米，相應(yīng)的嵌入式四元平面麥克風(fēng)陣列實(shí)物及其俯視二維測(cè)向示意圖如圖1所示。圖1中S為遠(yuǎn)場(chǎng)目標(biāo)，4個(gè)陣元M1、M2、M3、M4以嵌入式方式分別安裝于扁平長(zhǎng)方形殼體的4個(gè)側(cè)面，頂端均與殼體側(cè)面平齊。令陣列平面中心為本地坐標(biāo)系原點(diǎn)O，各陣元的坐標(biāo)向量分別為m1=[47 mm, 39.5 mm]T、m2=[34 mm, -52.5 mm]T、m3=[-47 mm, -39.5 mm]T、m4=[-34 mm, 52.5 mm]T。實(shí)際槍聲探測(cè)任務(wù)中的目標(biāo)俯仰角通常較小，以200 m外射擊點(diǎn)目標(biāo)為例，即使目標(biāo)高度達(dá)到30 m(約10層樓高)，目標(biāo)相對(duì)于陣列平面的俯仰角也僅為8.5°. 本文僅討論槍聲膛口波的方位角估計(jì)問(wèn)題，較小的目標(biāo)俯仰角對(duì)其方位角測(cè)量性能的影響可被忽略[9]。

由圖1(b)可以看出，總是存在至少1個(gè)陣元與目標(biāo)之間處于NLOS狀態(tài)?？紤]到遠(yuǎn)距離傳播的膛口波能量主要集中在低頻范圍，NLOS陣元接收到的槍聲信號(hào)主要是沿殼體衍射傳播而來(lái)，傳播延時(shí)增大，使得相關(guān)的實(shí)際TDOA偏離理想全通視直線傳播模型下的理論值。以來(lái)自第二象限的遠(yuǎn)場(chǎng)目標(biāo)為例，為盡量減小衍射傳播效應(yīng)的不利影響，可選擇波達(dá)時(shí)間(TOA)最小的通視陣元M4作為參考麥克風(fēng)，則在理想全通視條件下，目標(biāo)方向?qū)?yīng)的理想陣元間TDOA向量b應(yīng)滿足[8-9]：

(1)

式中：s=[sx,sy]T=[cosφcosθ, cosφsinθ]T為遠(yuǎn)場(chǎng)目標(biāo)方向矢量，θ和φ分別為目標(biāo)方位角(逆時(shí)針為正)和俯仰角(向上為正)；τij為陣列Mi和Mj之間的TDOA，i、j=1,2,3,4,i≠j；c為聲速；A為陣列系數(shù)矩陣。

(1) 式所示的超定方程存在唯一的線性LS解為

s=A?b=(ATA)-1ATb=[x,y]T,

(2)

式中：A?為矩陣A的Moore-Penrose逆；x和y分別為cosφcosθ和cosφsinθ的估計(jì)，從而可得遠(yuǎn)場(chǎng)目標(biāo)方位角的LS估計(jì)為

θ=arctan (y/x),

(3)

arctan(·)表示四象限反正切函數(shù)。

假設(shè)200 m處某個(gè)膛口波目標(biāo)的方位角和俯仰角分別為120°和0°，結(jié)合4個(gè)陣元位置坐標(biāo)可計(jì)算出理想全通視直線傳播模型下各陣元的TOA理論值，進(jìn)而可得以最小TOA陣元M4作為參考麥克風(fēng)時(shí)的TDOA理想值向量b=[-0.151 0 ms, -0.364 3 ms, -0.213 4 ms]T，代入(2)式和(3)式可得目標(biāo)方位角的LS估計(jì)為120.000 5°；即使目標(biāo)距離近至本文驗(yàn)證試驗(yàn)所用的15 m，目標(biāo)方位角的LS估計(jì)也為120.006 6°，估計(jì)誤差仍可忽略。但實(shí)際陣元間TDOA因存在殼體遮擋引起的衍射傳播延時(shí)而將大于相應(yīng)的理想值，再加上不可避免的隨機(jī)測(cè)量噪聲和異常值等實(shí)際因素，基于理想全通視模型的目標(biāo)方位角LS估計(jì)結(jié)果將出現(xiàn)明顯偏差。例如前期試驗(yàn)中在120°方位角處的5次TDOA實(shí)測(cè)值對(duì)應(yīng)的目標(biāo)方位角估計(jì)分別為125.59°、127.07°、125.45°、130.47°、126.73°，誤差高達(dá)5.45°～10.47°. 更為棘手的是，這種測(cè)向誤差是隨目標(biāo)來(lái)向和環(huán)境因素的不同而變化，造成在目標(biāo)位置未知的實(shí)際應(yīng)用條件下，殼體遮擋衍射傳播效應(yīng)引起的測(cè)向誤差方向和誤差程度都難以實(shí)時(shí)預(yù)測(cè)，因而槍聲探測(cè)系統(tǒng)的膛口波測(cè)向性能必須根據(jù)實(shí)際測(cè)量數(shù)據(jù)進(jìn)行校正。

2 TDOA矩陣校正

考慮到槍聲探測(cè)系統(tǒng)采用的麥克風(fēng)陣列形狀和陣元位置信息均為已知的先驗(yàn)知識(shí)且通常固定不變，本文以圖1所示的四元嵌入式平面麥克風(fēng)陣列為例，提出一種TDOA矩陣校正方法，通過(guò)同時(shí)施加低秩及幾何約束分別抑制TDOA測(cè)量值中的隨機(jī)測(cè)量噪聲和衍射傳播影響，從而明顯減小嵌入式平面陣列的膛口波測(cè)向偏差。

令遠(yuǎn)場(chǎng)目標(biāo)信號(hào)傳播至陣元Mi的波達(dá)時(shí)間為ti，i=1,2,3,4. 定義陣列的TOA向量為x=[t1,t2,t3,t4]T，以及陣元Mi和Mj之間的τij=ti-tj，i、j=1,2,3,4,i≠j，結(jié)合全1向量1=[1, 1, 1, 1]T可得相應(yīng)的反對(duì)稱TDOA矩陣

M=x1T-1xT.

(4)

顯然，(1)式中TDOA向量b的所有信息均已包含在矩陣M中。當(dāng)同時(shí)存在TDOA測(cè)量噪聲以及少量異常值(含丟失值賦值)時(shí)，實(shí)測(cè)TDOA矩陣可表示為

=M+N+S,

(5)

式中：N為反對(duì)稱的測(cè)量噪聲矩陣；S為T(mén)DOA矩陣中的加性異常值。對(duì)于俯仰角較小的遠(yuǎn)場(chǎng)目標(biāo)，(4)式所示理想全通視TDOA矩陣M的秩為2，不會(huì)出現(xiàn)所有陣元與目標(biāo)距離相同的平凡解情況(此時(shí)M為全零矩陣)?；谶@一事實(shí)，文獻(xiàn)[14]采用低秩和異常值個(gè)數(shù)約束將(5)式中對(duì)矩陣M的估計(jì)任務(wù)轉(zhuǎn)化為一種非凸優(yōu)化問(wèn)題并求解。值得注意的是，文獻(xiàn)[14]中方法有效的前提是存在一個(gè)隱含假設(shè)，就是(5)式右式中的M+S不能分解為(4)式所示形式，即不存在一種實(shí)際陣列構(gòu)型，它的無(wú)噪TOA向量x能通過(guò)(4)式得到秩2的反對(duì)稱矩陣M+S.

然而對(duì)于嵌入式平面陣列，當(dāng)部分NLOS陣元受殼體遮擋衍射傳播效應(yīng)影響時(shí)，情況將與之不同。若不考慮其他任何TDOA異常值和隨機(jī)測(cè)量誤差，實(shí)際的無(wú)噪TDOA矩陣仍然可以表示為M+S形式，但此時(shí)的M和S將分別由理想全通視TDOA值和衍射傳播TDOA偏差構(gòu)成。由于嵌入式陣列中NLOS陣元因其對(duì)應(yīng)的實(shí)際TOA大于理想全通視情況而可等效于陣列形狀發(fā)生了改變，此時(shí)(5)式中的M+S將與受衍射傳播影響的實(shí)際無(wú)噪TOA向量對(duì)應(yīng)，即無(wú)論是M和S本身還是M+S作為一個(gè)整體，都可以分解為(4)式所示的反對(duì)稱秩2矩陣構(gòu)成方式。由于衍射傳播TDOA偏差矩陣S還隨目標(biāo)來(lái)向變化且難以實(shí)時(shí)預(yù)測(cè)，因此僅采用秩2約束無(wú)法分辨M和S，使得文獻(xiàn)[14]所提方法因其隱含假設(shè)不再成立而無(wú)法準(zhǔn)確重構(gòu)理想全通視TDOA矩陣M，即不能校正衍射傳播引起的TDOA偏差。

為了同時(shí)抑制TDOA實(shí)際測(cè)量值中的隨機(jī)測(cè)量噪聲與衍射傳播偏差，本文在現(xiàn)有低秩約束基礎(chǔ)上又施加了對(duì)應(yīng)陣列形狀特征的幾何約束，即將(5)式中對(duì)理想全通視TDOA矩陣M的估計(jì)任務(wù)轉(zhuǎn)化為如下有約束凸優(yōu)化問(wèn)題：

(6)

式中：‖·‖F(xiàn)為Frobenius范數(shù)；約束矩陣D=[1,c]，c為幾何約束向量，包括1Tx=0和cTx=0兩部分，前者與理想全通視TDOA矩陣M應(yīng)為秩2的低秩要求密切相關(guān)，用于抑制高斯隨機(jī)測(cè)量噪聲[14]，后者則體現(xiàn)了圖1所示平行四邊形麥克風(fēng)陣列所對(duì)應(yīng)的形狀特征幾何約束。由于τ12=τ43，因而(6)式中

c=[1, -1, 1, -1]T.

(7)

針對(duì)(6)式應(yīng)用Lagrange乘數(shù)法構(gòu)造函數(shù)

(8)

式中：λ1、λ2為拉格朗日乘子。

考慮到

則令(8)式分別對(duì)x、λ1、λ2的1階導(dǎo)數(shù)為0，可得

(9)

(10)

最后可得校正后的理想全通視TDOA矩陣估計(jì)為

(11)

3 嵌入式平面陣膛口波測(cè)向方法

本文提出的嵌入式平面陣膛口波測(cè)向算法主要包括3個(gè)主要步驟：參考麥克風(fēng)選取、TDOA矩陣校正、LS方位角估計(jì)。

首先，考慮到與目標(biāo)距離最小的陣元處于通視狀態(tài)，選擇該陣元作為參考麥克風(fēng)有利于減小(1)式中TDOA向量b受殼體遮擋衍射傳播效應(yīng)的影響程度，從而為降低LS方位角估計(jì)偏差提供良好基礎(chǔ)。由(4)式可知，TDOA矩陣M中每一列都是TOA向量x加上不同的時(shí)間參考點(diǎn)，因而理論上每一列中與最近通視陣元序號(hào)對(duì)應(yīng)的元素均應(yīng)在該列中數(shù)值最小，即TOA最短。但需要強(qiáng)調(diào)的是，由于TDOA估計(jì)過(guò)程中不可避免地會(huì)存在測(cè)量噪聲甚至出現(xiàn)個(gè)別異常值，實(shí)測(cè)TDOA矩陣中每列各元素之間的大小關(guān)系有可能發(fā)生偶然隨機(jī)變化。因此采用多數(shù)表決方式確定參考麥克風(fēng)序號(hào)，即首先記錄實(shí)測(cè)TDOA矩陣中每列最小值所在的行序號(hào)，然后以最小值出現(xiàn)記錄最多的行序號(hào)作為距離目標(biāo)最近的通視陣元序號(hào)，并選取該陣元作為參考麥克風(fēng)。

需要說(shuō)明的是，通過(guò)改變各矩陣或向量的維數(shù)以及幾何約束的表現(xiàn)形式，上述膛口波測(cè)向算法就可直接推廣到陣元數(shù)量或陣列形狀各不相同的其他嵌入式麥克風(fēng)陣列系統(tǒng)。

4 試驗(yàn)驗(yàn)證

驗(yàn)證試驗(yàn)在較空曠的室外環(huán)境進(jìn)行，試驗(yàn)布設(shè)示意圖及試驗(yàn)現(xiàn)場(chǎng)照片如圖2所示。圖2中的嵌入式四元平面麥克風(fēng)陣列系統(tǒng)被水平固定在三腳架上，以陣列中心為坐標(biāo)原點(diǎn)，通過(guò)水平順時(shí)針旋轉(zhuǎn)陣列系統(tǒng)來(lái)等效從小到大地改變目標(biāo)方位角。作為模擬射擊點(diǎn)目標(biāo)的揚(yáng)聲器與陣列高度基本相同，放置在距離陣列中心15 m處，該距離遠(yuǎn)大于陣列孔徑，因而可模擬遠(yuǎn)場(chǎng)情況。試驗(yàn)中，目標(biāo)方位角從0°開(kāi)始以30°步進(jìn)遞增變化一周，在每個(gè)目標(biāo)方位角位置進(jìn)行26次獨(dú)立試驗(yàn)，12個(gè)方位角位置共計(jì)進(jìn)行312次獨(dú)立試驗(yàn)。每次試驗(yàn)中揚(yáng)聲器播放一段含有外場(chǎng)實(shí)測(cè)的槍聲膛口波信號(hào)錄音。

圖2 試驗(yàn)布設(shè)示意圖及試驗(yàn)現(xiàn)場(chǎng)Fig.2 Experimental layout and experimental site

采用理想全通視直線傳播模型下基于TDOA的常規(guī)LS測(cè)向算法(下文簡(jiǎn)稱常規(guī)方法)作為性能對(duì)比，并定義測(cè)向誤差為目標(biāo)方位角估計(jì)值與真實(shí)值之差。

圖3給出了本文方法和常規(guī)方法在12個(gè)測(cè)試方位角上各26次獨(dú)立試驗(yàn)的測(cè)向誤差分布。圖3中，常規(guī)方法的測(cè)角誤差在一些方位角處散布較大，而經(jīng)過(guò)本文方法校正后，各方位角的測(cè)向誤差分布則已趨向集中。

圖3 本文方法和常規(guī)方法的測(cè)向誤差分布Fig.3 Distribution of azimuth estimation errors of the proposed and conventional methods

圖4比較了兩種方法的測(cè)向誤差總體分布直方圖。對(duì)比圖4(a)、圖4(b)可以看出，常規(guī)方法的測(cè)向誤差主要集中分布在-2.03°～5.14°范圍內(nèi)(上下四分位數(shù)誤差范圍)，分布均勻程度較高，且誤差絕對(duì)值的最大值高達(dá)23°. 而本文方法測(cè)向誤差則主要集中分布在-1.38°～2.84°范圍內(nèi)(上下四分位數(shù)誤差范圍)，相應(yīng)頻度隨著誤差絕對(duì)值增大而很快下降，誤差絕對(duì)值的最大值也降為10°.

圖4 本文方法和常規(guī)方法的測(cè)向誤差分布直方圖Fig.4 Histograms of azimuth estimation errors of the proposed and conventional methods

圖4中，常規(guī)方法和本文方法的測(cè)向誤差統(tǒng)計(jì)均值正負(fù)標(biāo)準(zhǔn)差分別為1.73°±5.87°和0.90°±3.55°，統(tǒng)計(jì)中值分別為1.88°和0.87°. 實(shí)際目標(biāo)方位角通常未知且隨機(jī)，通常假設(shè)實(shí)際目標(biāo)方位角在0°～360°之間均勻分布。為使全方位的總體測(cè)角性能相對(duì)均衡，同時(shí)考慮到兩種方法測(cè)向誤差的統(tǒng)計(jì)均值與統(tǒng)計(jì)中值基本相等，實(shí)際應(yīng)用時(shí)可預(yù)先補(bǔ)償上述統(tǒng)計(jì)均值。兩種方法在均值補(bǔ)償后的測(cè)向誤差絕對(duì)值所對(duì)應(yīng)的累積概率分布對(duì)比如圖5所示。需要強(qiáng)調(diào)的是，就槍聲探測(cè)系統(tǒng)實(shí)戰(zhàn)環(huán)境下的技術(shù)指標(biāo)而言，當(dāng)前以英國(guó)Ears SWATS II代、美國(guó)Boomerang Warrior-X系統(tǒng)為代表，這些國(guó)外同類便攜式產(chǎn)品的槍聲膛口波測(cè)向性能通常以±7.5°為測(cè)向精度考核標(biāo)準(zhǔn)[4]。通過(guò)考察圖5可知：針對(duì)測(cè)向誤差絕對(duì)值在7.5°以內(nèi)的累積概率，常規(guī)方法僅為0.86；而本文方法達(dá)到0.94，測(cè)向性能明顯提升。

圖5 均值補(bǔ)償后兩種方法的絕對(duì)測(cè)向誤差累積概率分布Fig.5 Cumulative probability distributions of absolute azimuth estimation errors after mean value compensation

5 結(jié)論

對(duì)于小型嵌入式平面陣槍聲探測(cè)系統(tǒng)，部分陣元將因殼體遮擋而不能直視目標(biāo)。低頻能量集中的槍聲膛口波信號(hào)傳播至NLOS陣元會(huì)因衍射傳播而出現(xiàn)額外延時(shí)，導(dǎo)致實(shí)際陣元間TDOA不再滿足與陣元間距和目標(biāo)方向之間的理論函數(shù)關(guān)系，使得基于理想全通視直線傳播模型的常規(guī)LS測(cè)向結(jié)果發(fā)生較大誤差，該誤差在實(shí)際中很難預(yù)測(cè)和建模。通過(guò)在低秩約束的基礎(chǔ)上增加陣列形狀幾何約束來(lái)構(gòu)建凸優(yōu)化問(wèn)題，同時(shí)選取距離目標(biāo)最近的通視麥克風(fēng)作為參考陣元，能夠有效抑制TDOA隨機(jī)測(cè)量噪聲與殼體遮擋衍射傳播效應(yīng)，明顯提高基于TDOA的LS測(cè)向精度。試驗(yàn)結(jié)果表明，所提方法的膛口波測(cè)向誤差絕對(duì)值在7.5°以內(nèi)的累積概率由常規(guī)方法的0.86提升至0.94，誤差分布也更為聚集，能夠滿足實(shí)際便攜式槍聲探測(cè)系統(tǒng)的測(cè)向指標(biāo)要求。