一種聲學(xué)測試彈著點(diǎn)坐標(biāo)的穩(wěn)健估計(jì)方法

王玉龍 盧小汐 張亞輝 李宏凱 蘇磊

摘? 要： 針對火炮彈著點(diǎn)聲檢靶系統(tǒng)中彈丸著靶坐標(biāo)精度不穩(wěn)定的問題，建立抗差最小二乘法。該方法通過等價(jià)權(quán)函數(shù)將抗差估計(jì)與最小二乘法結(jié)合起來，利用殘差對異常值進(jìn)行降權(quán)處理，從而消除異常值對精度的不良影響。最后，對火炮彈著點(diǎn)聲檢靶系統(tǒng)進(jìn)行仿真分析。分析和試驗(yàn)表明，抗差最小二乘法對異常值是穩(wěn)健的，且精度優(yōu)于一般最小二乘法，當(dāng)靶幅為10 m×10 m且風(fēng)速測試誤差為3 m/s時(shí)，火炮彈著點(diǎn)聲檢靶系統(tǒng)針對火炮彈丸彈著點(diǎn)坐標(biāo)測試精度優(yōu)于5 cm。

關(guān)鍵詞： 彈著點(diǎn)坐標(biāo); 穩(wěn)健估計(jì); 抗差最小二乘法; 抗差估計(jì); 殘差變化; 聲定位

中圖分類號： TN911.1?34; TJ306? ? ? ? ? ? ? ? ? ?文獻(xiàn)標(biāo)識碼： A? ? ? ? ? ? ? ? ? ? 文章編號： 1004?373X（2020）09?0025?04

A robust estimation method of acoustic test for impact point coordinate

WANG Yulong， LU Xiaoxi， ZHANG Yahui， LI Hongkai， SU Lei

（Huayin Ordnance Test Center， Huayin 714200， China）

Abstract： A robustified least square method is established for the instability coordinate precision of the sound detection target system for the projectile impact point. In the method， the robust estimation is combined with the least square method by the equivalent weight function， and the residual is used to reduce the weight of the outliers to eliminate the adverse effects of the outliers on the accuracy. The simulation analysis of sound detection target system for the projectile impact point was carried out. The analysis and test results show that the robustified least square method is robust to outliers and its accuracy is higher than that obtained with the general least square method. When the target area is 10 m×10 m and the measurement error of wind speed is 3 m/s， the tested coordinate accuracy of the impact point obtained with the sound detection target system for projectile impact point is within 5 cm.

Keywords： impact point coordinate; robust estimation; robustified least square method; robust estimation; residual variation; acoustic positioning

0? 引? 言

聲學(xué)立靶彈著點(diǎn)測試法屬于基于到達(dá)時(shí)間差的無源定位，一般可得到一組非線性雙曲方程組。一般采用最小二乘法（Least Square，LS）進(jìn)行求解[1?3]。在彈著點(diǎn)聲學(xué)靶系統(tǒng)試驗(yàn)測試過程中發(fā)現(xiàn)，有時(shí)傳感器測量的時(shí)差存在粗差，即異常值或野值。由于測試系統(tǒng)中每個(gè)傳感器所處位置不同，根據(jù)其測試波形處理得到的波達(dá)時(shí)刻也不同，很難從表面上區(qū)分哪些測試數(shù)據(jù)屬于異常值或野值。在采用最小二乘法確定著靶坐標(biāo)時(shí)，因?yàn)樽钚《朔▽Ξ惓Ｖ档拿舾行?，使得系統(tǒng)測試精度變低。

本文將抗差估計(jì)[4?5]與最小二乘法結(jié)合，建立了穩(wěn)健估計(jì)方法（Robust Least Square，RLS），消除了異常值對定位結(jié)果精度的影響，且計(jì)算流程與最小二乘法一致。最后，給出了仿真與試驗(yàn)結(jié)果。

1? 定位原理

聲學(xué)靶系統(tǒng)主要采用傳聲器陣列，采集超音速彈丸產(chǎn)生的激波信號，再根據(jù)傳聲器坐標(biāo)與激波時(shí)間差，確定彈丸通過空中虛擬平面的坐標(biāo)。聲學(xué)靶及立靶坐標(biāo)系示意圖如圖1所示。其中，傳聲器采用雙圓環(huán)布陣，每個(gè)圓環(huán)上均勻分布12個(gè)傳聲器，[xz]平面為水平面，[z]軸指向炮口方向，[xy]平面為豎直平面，[xy]平面即為虛擬靶平面，雙圓環(huán)陣列位于靶平面內(nèi)。

計(jì)算分析和實(shí)驗(yàn)表明，風(fēng)對定位精度的影響不可忽略。假定風(fēng)平行地面運(yùn)動(dòng)，且在彈丸著靶時(shí)刻靶面周圍的局部風(fēng)場是均勻的。因?yàn)轱L(fēng)在[z]方向的分量[vz]影響較小[6?8]，所以主要考慮風(fēng)在[x]方向的分量[vx]。利用矢量疊加原理可得到定位方程組如下：

[（x1-x-vx（t0+t1））2+（y1-y）2=v0（t0+t1）（x2-x-vx（t0+t2））2+（y2-y）2=v0（t0+t2）? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??（xn-x-vx（t0+tn））2+（yn-y）2=v0（t0+tn）] （1）

式中：[ti，i=1，2，…，n]是激波到達(dá)各個(gè)傳感器的相對時(shí)間（即時(shí)延或時(shí)差）;[t0]是激波視速度在靶面上傳播至首先觸發(fā)傳感器的時(shí)間;[v0]是視速度。

方程組（1）中的未知參量為[x，y，v0，t0]。當(dāng)[n≥4]時(shí)方程組（1）有解;當(dāng)[n=4]時(shí)，方程組（1）是非線性方程組;[n>4]時(shí)，方程組（1）是超定非線性方程組（矛盾方程組）。求解方程組（1）即得[x，y，v0，t0]。

2? 彈著點(diǎn)穩(wěn)健估計(jì)方法

2.1? 迭代公式

最小二乘法理論完善，算法成熟，得到了廣泛應(yīng)用。但是，當(dāng)觀測數(shù)據(jù)含有粗差（或異常值）時(shí)，經(jīng)典最小二乘法的結(jié)果可能因異常值的存在而不可靠。在聲靶試驗(yàn)過程中，由于傳感器偶然抖動(dòng)、風(fēng)吹等不確定因素的影響，傳感器測量的時(shí)差有時(shí)會(huì)出現(xiàn)異常值。此時(shí)，應(yīng)用最小二乘法計(jì)算時(shí)，導(dǎo)致著靶測試精度低。穩(wěn)健估計(jì)（也叫抗差估計(jì)）是指在存在異常值的情況下，通過選擇適當(dāng)?shù)墓烙?jì)方法，使估計(jì)結(jié)果盡可能少地受到異常值干擾。穩(wěn)健估計(jì)的原則是充分利用有效信息，限制利用可用信息，排除有害信息。

抗差最小二乘法通過等價(jià)權(quán)將抗差估計(jì)與最小二乘法結(jié)合起來，其實(shí)質(zhì)是等價(jià)權(quán)函數(shù)的設(shè)計(jì)。當(dāng)實(shí)測時(shí)差不含異常值時(shí)，抗差最小二乘法的結(jié)果與最小二乘法結(jié)果一致;當(dāng)實(shí)測時(shí)差含有異常值時(shí)，抗差最小二乘法能夠消除異常值的影響。

下面詳細(xì)推導(dǎo)抗差最小二乘法的迭代公式，記：

[fi=（x-xi+vx（t0+ti））2+（y-yi）2-v20（ti+t0）2，? ? ? ? ? ? ? ? ? ? ? ? ? ? ?i=1，2，…，n]? ?（2）

設(shè)[f： D?R4→Rn]，[f=（f1，f2，…，fn）T]。求解非線性方程組（1）等價(jià)于求解如下非線性最小二乘問題。

[minx∈D12fT（x）f（x）]? ?（3）

利用Taylor級數(shù)展開將非線性問題轉(zhuǎn)化為線性問題，在[x0=（x0，y0，v00，t00）]展開得：

[f（x0）+DfT（x0）（x-x0）=0]? ?（4）

式中：[Df（x0）=Df1DxDf1DyDf1Dv0Df1Dt0Df2DxDf2DyDf2Dv0Df2Dt0????DfnDxDfnDyDfnDv0DfnDt0x0]是Jacobi矩陣;[DfiDx=2（x-xi+vx（t0+ti））];[DfiDy=2（y-yi）];[DfiDv0=][-2v0（ti+t0）2];[DfiDt0=-2v20（ti+t0）];[i=1，2，…，n]。

記[b=-f（x0）]，[A=DfT（x0）]，[X=（x-x0）]，則：

[AX=b] （5）

進(jìn)而，有：[x1=x0+（ATP0A）-1ATP0b]。

一般地，抗差最小二乘迭代公式如下：

[xk+1=xk+（ATkPkAk）-1ATkPkbk，k=0，1，2，…] （6）

式中[Pk]稱為等價(jià)權(quán)矩陣或等價(jià)權(quán)函數(shù)。特別地，當(dāng)取[Pk=I]時(shí)，抗差最小二乘法退化為最小二乘法。

就數(shù)學(xué)表達(dá)式而言，抗差最小二乘法與加權(quán)最小二乘法基本相同，其區(qū)別是權(quán)函數(shù)的內(nèi)容不同。加權(quán)最小二乘法的權(quán)是先驗(yàn)的，而抗差最小二乘法的權(quán)是殘差的函數(shù)。

2.2? 等價(jià)權(quán)矩陣的確定

權(quán)函數(shù)的作用在于：將所有的觀測數(shù)據(jù)劃分為正常觀測值、可利用觀測值和粗差觀測值（異常值）三個(gè)部分。對于正常觀測值，使其保持原有的權(quán)不變，對余差較大的可利用觀測值進(jìn)行降權(quán)處理，而那些粗差觀測值，則使其權(quán)為零。經(jīng)過這樣的處理，就可以根據(jù)各個(gè)觀測值及“后驗(yàn)權(quán)”，求得最接近于正常情況下的結(jié)果。本文中權(quán)函數(shù)選為IGGⅢ權(quán)函數(shù)（IGG是中國科學(xué)院測量與地球物理研究所的英文縮寫），其表達(dá)式如下：

[Pi（vi）=1，? ? vi≤k0k0vik1-vik1-k02，? ? k0k1] （7）

式中：[vi=δiσ]，[δi=fi]是殘差，[σ]是標(biāo)準(zhǔn)化殘差，[σ=i=1nδi（n-1）]，[k0∈（1，1.5）]，[k1∈（2.5，8.0）]，本文取[k0=1.5]，[k1=3]。

IGGⅢ權(quán)函數(shù)充分考慮了正常階段、可疑階段和淘汰階段三個(gè)部分，公式簡單，易于實(shí)現(xiàn)，如圖2所示。

2.3? 算法流程圖

抗差最小二乘法的計(jì)算流程如圖3所示。

2.4? 等價(jià)權(quán)和殘差變化的算例

以某型榴彈試驗(yàn)為例進(jìn)行分析，圖4給出了迭代過程中不同傳感器對應(yīng)的殘差和等價(jià)權(quán)，其中橫坐標(biāo)代表傳感器編號，縱坐標(biāo)分別代表殘差、標(biāo)準(zhǔn)化殘差和等價(jià)權(quán)。

由圖4可知，在迭代的過程中，有多個(gè)傳聲器的等價(jià)權(quán)發(fā)生變化，隨著迭代的不斷進(jìn)行，最后6號傳聲器的等價(jià)權(quán)減小，并降為0，即最終判斷6號傳感器存在異常。圖5給出了對應(yīng)的彈著點(diǎn)誤差隨迭代次數(shù)的變化曲線。

由圖5可知，隨著迭代次數(shù)的增加，水平誤差基本不變，而高低誤差由42.3 cm降低到1.0 cm，即抗差最小二乘法消除了異常值的影響。

3? 仿真分析

影響聲學(xué)靶測量誤差的因素主要有：傳感器坐標(biāo)測試誤差、傳感器時(shí)延誤差、溫度測試誤差和風(fēng)速測試誤差。傳感器坐標(biāo)測試誤差和時(shí)延誤差由設(shè)備決定，溫度變化比較穩(wěn)定。聲靶計(jì)算修正時(shí)用到的環(huán)境風(fēng)速是激波在靶面附近傳播幾毫秒間的瞬態(tài)風(fēng)速，測試誤差較大。下面以靶幅10 m×10 m，傳感器坐標(biāo)測試誤差為1 mm，傳感器時(shí)延誤差為5 μs，溫度測試誤差為5 ℃，風(fēng)速測試誤差分別為1 m/s，2 m/s，3 m/s和4 m/s時(shí)，給出誤差分布圖，如圖6，圖7所示。

由圖6和圖7知，隨著風(fēng)速測試誤差的增加，水平誤差和高低誤差也隨之增加。表1給出了在不同風(fēng)速測試誤差時(shí)的均方根誤差。

4? 試驗(yàn)數(shù)據(jù)結(jié)果與分析

4.1? 鞭炮模擬試驗(yàn)

采用鞭炮爆炸激發(fā)傳聲器，事先測量出鞭炮炸點(diǎn)和傳聲器的相對關(guān)系，作為真值，檢測聲學(xué)靶的精度。分別采用LS法、RLS法和交會(huì)法進(jìn)行計(jì)算，對比不同算法的精度。其中，交會(huì)法采用文獻(xiàn)[8]中的交會(huì)計(jì)算方法。試驗(yàn)結(jié)果如圖8所示。

由圖8可知，交會(huì)法的水平均方根誤差為5.0 cm，高低均方根誤差為6.0 cm;LS法的水平均方根誤差為2.8 cm，高低均方根誤差為8.1 cm;RLS法的水平均方根誤差為2.8 cm，高低均方根誤差為2.9 cm。從而可知，LS法容易受到異常值的影響，而RLS法能夠保證測量的穩(wěn)健性，并且RLS法的精度優(yōu)于交會(huì)法。

4.2? 某型榴彈試驗(yàn)

某型榴彈試驗(yàn)同樣分別采用LS法、RLS法和交會(huì)法進(jìn)行計(jì)算，對比不同算法的精度。結(jié)果如圖9所示。

由圖9知，交會(huì)法的水平均方根誤差為5.5 cm，高低均方根誤差為7.6 cm;LS法的水平均方根誤差為3.2 cm，高低均方根誤差為15.4 cm;RLS法的水平均方根誤差為3.3 cm，高低均方根誤差為4.5 cm。與鞭炮模擬試驗(yàn)結(jié)果類似，RLS法精度最高。

5? 結(jié)? 論

本文針對火炮彈著點(diǎn)聲學(xué)靶系統(tǒng)定位計(jì)算方法進(jìn)行研究，建立了基于抗差最小二乘法的穩(wěn)健估計(jì)方法，提高了彈丸著靶坐標(biāo)的測試精度。對火炮彈著點(diǎn)聲學(xué)靶系統(tǒng)進(jìn)行仿真分析，風(fēng)速測試誤差越大，著靶坐標(biāo)測試精度越低。當(dāng)風(fēng)速測試誤差為3 m/s時(shí)，著靶坐標(biāo)測試精度為5 cm。

參考文獻(xiàn)

[1] 王磊，陳昭男.基于被動(dòng)聲學(xué)的超聲速武器脫靶量解算方法[J].應(yīng)用聲學(xué)，2018，37（3）：385?390.

[2] 苗晟，周維，唐浩，等.一種聲源定位系統(tǒng)設(shè)計(jì)[J].計(jì)算機(jī)科學(xué)，2013，40（11A）：398?400.

[3] 劉哲，陳日林，滕鵬曉，等.基于平面?zhèn)髀暺麝嚵械穆曉炊ㄎ幌到y(tǒng)[J].聲學(xué)技術(shù)，2011，30（2）：123?128.

[4] 張洋，張志剛，錢棟，等.基于IGGⅢ模型的高精度抗差坐標(biāo)轉(zhuǎn)換算法研究[J].全球定位系統(tǒng)，2017，42（5）：29?32.

[5] 李蟬，張士峰，張力軍.應(yīng)用抗差無跡卡爾曼濾波的再入彈道處理技術(shù)[J].國防科技大學(xué)學(xué)報(bào)，2017，39（1）：1?5.

[6] MSI. Strafe scoring target type 590 air?to?ground training [EB/OL]. [2014?07?20]. http：//www.msinstruments.com.

[7] 蔣東東，狄長安.八點(diǎn)線式聲學(xué)立靶彈著點(diǎn)檢測系統(tǒng)研究[J].電子測量技術(shù)，2010，33（4）：19?21.

[8] 馮斌，石秀華.雙三角陣聲靶測試系統(tǒng)研究[J].應(yīng)用聲學(xué)，2012，31（2）：140?144.

[9] 董杰，高樓.彈丸任意角度入射彈著點(diǎn)聲學(xué)檢測模型[J].計(jì)算機(jī)仿真，2015，32（1）：10?14.