李志明，范錦彪

(中北大學(xué) 儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051)

0 引言

彈著點(diǎn)是評(píng)價(jià)武器著靶精度的重要指標(biāo)，準(zhǔn)確獲取彈著點(diǎn)位置對(duì)于武器威力評(píng)估具有重要意義[1]。彈丸落地爆炸時(shí)，會(huì)形成破壞區(qū)、塑性帶以及彈性形變區(qū)三個(gè)區(qū)域，并最終在彈性形變區(qū)形成在地下傳播的彈性波——地震波。對(duì)動(dòng)能彈而言，當(dāng)彈丸質(zhì)量、落地時(shí)速度不同時(shí)，其對(duì)應(yīng)動(dòng)能不同，產(chǎn)生地震波的幅值和頻率也不同，動(dòng)能越大，產(chǎn)生地震波的幅值越大，傳播的距離也會(huì)越遠(yuǎn)，因此可為彈丸的落點(diǎn)定位提供有用信息。隨著微震定位技術(shù)發(fā)展及在許多領(lǐng)域的應(yīng)用[2-4]，彈著點(diǎn)定位也得到了越來(lái)越多學(xué)者的重視。張躍華等[5]改進(jìn)了Geiger方法和震源掃描算法，并進(jìn)行聯(lián)合計(jì)算，大幅提高了定位精度并縮短了計(jì)算時(shí)間；李鵬宇等[6]采用米字形陣列，預(yù)先通過(guò)來(lái)波方向(DOA)算法估計(jì)波速，然后將波速估計(jì)值代入到達(dá)時(shí)間差(TDOA)算法中計(jì)算初始位置，最后通過(guò)泰勒級(jí)數(shù)算法收斂定位；文獻(xiàn)[7-8]通過(guò)構(gòu)建方程組的方式求解彈著點(diǎn)坐標(biāo)。

目前的彈著點(diǎn)定位方法多數(shù)都需要提前測(cè)量地震波波速或?qū)ζ溥M(jìn)行反演，但由于地震波在地表介質(zhì)中傳播時(shí)在各個(gè)方向上的速度都不一致，即便是在同一方向上也會(huì)隨著距離的增大而衰減，因此會(huì)導(dǎo)致定位誤差增加。而彈著點(diǎn)定位的關(guān)鍵是得到彈丸落點(diǎn)坐標(biāo)，無(wú)需求解波速，因此本文提出一種基于波速方差的目標(biāo)函數(shù)，并利用單純形法搜索定位點(diǎn)坐標(biāo)。然而利用單純形法進(jìn)行搜索時(shí)，當(dāng)選擇的初始單純形不同時(shí)，定位結(jié)果會(huì)存在一定的偏差，個(gè)別條件下會(huì)使目標(biāo)函數(shù)收斂于局部最小值，出現(xiàn)異常的定位結(jié)果。本文針對(duì)這種情況，采用帶噪聲的密度聚類(lèi)(DBSCAN)算法對(duì)隨機(jī)生成不同的初始單純形進(jìn)行重復(fù)定位的結(jié)果進(jìn)行判別，剔除異常值，對(duì)剩余結(jié)果求平均作為最終的彈著點(diǎn)坐標(biāo)。

1 無(wú)速度參數(shù)測(cè)量模型

在定位系統(tǒng)中，可以直接獲得的數(shù)據(jù)只有各個(gè)測(cè)點(diǎn)的坐標(biāo)以及測(cè)點(diǎn)拾取到的地震波到達(dá)時(shí)間，彈丸落地的時(shí)刻、地震波走時(shí)和地震波傳播速度等參數(shù)難以在第一時(shí)間獲得，而對(duì)于定位來(lái)講，這些恰恰是獲得彈著點(diǎn)位置的重要參數(shù)，傳統(tǒng)采用基于地震波TDOA算法可以有效避免對(duì)彈丸落地時(shí)刻和地震波走時(shí)的求解，但是仍然需要提前測(cè)量波速。董隴軍等[9]提出基于TDOA的新方法雖然不需要提前測(cè)量波速，但是該方法將波速作為未知參數(shù)與震源坐標(biāo)一起反演求解，然而在實(shí)際應(yīng)用中，地震波的傳播速度是不均勻的，在地表介質(zhì)中傳播時(shí)，各個(gè)方向的速度都不一致，即便是在同一方向上也會(huì)隨著距離的增加而衰減，因此在求解時(shí)會(huì)導(dǎo)致定位誤差增加。李健等[10]提出無(wú)需測(cè)速的目標(biāo)函數(shù)，雖然不需要提前測(cè)量波速以及對(duì)波速進(jìn)行反演，但是，其目標(biāo)函數(shù)較為復(fù)雜，當(dāng)測(cè)點(diǎn)數(shù)量增加時(shí)，運(yùn)算量也會(huì)隨之增加。因此，本文提出了如下基于波速方差的目標(biāo)函數(shù)。

假設(shè)彈丸撞擊地面時(shí)，地表發(fā)生的是彈性形變，且地面無(wú)大型溝渠、無(wú)外界干擾源。在待測(cè)試區(qū)域構(gòu)建平面直角坐標(biāo)系，由于坐標(biāo)系原點(diǎn)所在位置不會(huì)對(duì)定位結(jié)果產(chǎn)生影響，所以可任意指定原點(diǎn)位置。假設(shè)彈著點(diǎn)坐標(biāo)為(x0,y0)，彈丸落地時(shí)刻，即彈丸撞擊地面的時(shí)刻為t0，第i、第j個(gè)測(cè)點(diǎn)的坐標(biāo)分別為(xi,yi)、(xj,yj)，第i、第j個(gè)測(cè)點(diǎn)拾取到的地震波初至?xí)r間，即彈丸落地時(shí)刻與地震波走時(shí)之和分別為ti、tj，地震波傳播速度為v，則有

(1)

(2)

兩式相減，可得

ti-tj=(Li-Lj)/v，

(3)

式中：

(4)

(5)

易得

(6)

(7)

因此，速度的方差為

(8)

(8)式即為無(wú)速度參數(shù)的目標(biāo)函數(shù)。

測(cè)點(diǎn)的位置坐標(biāo)是一定的，在一次彈丸落地事件中，測(cè)點(diǎn)拾取到地震波的初至?xí)r間也是一定的。如果(x0,y0)不是真實(shí)的彈著點(diǎn)坐標(biāo)，那么計(jì)算得到的各個(gè)速度值之間偏差會(huì)較大，因此計(jì)算得到的速度方差也會(huì)較大。只有真實(shí)的彈著點(diǎn)坐標(biāo)才會(huì)使速度方差最小，所以，尋找到可以使速度方差最小的點(diǎn)即為真實(shí)的彈著點(diǎn)。

2 模型仿真驗(yàn)證

2.1 模型參數(shù)

目標(biāo)函數(shù)(8)式中，未知參數(shù)僅為彈著點(diǎn)坐標(biāo)，無(wú)需求解波速，減少了誤差來(lái)源，具有更高的可靠性。用單純形法對(duì)其尋優(yōu)即可實(shí)現(xiàn)定位，為驗(yàn)證算法是否有效，構(gòu)建如圖1所示的仿真定位模型。假定監(jiān)測(cè)范圍為100 m×100 m的區(qū)域，共布置8個(gè)測(cè)點(diǎn)，編號(hào)分別為1～8. 在實(shí)際應(yīng)用中，每個(gè)測(cè)點(diǎn)由一個(gè)檢波器和一個(gè)波形記錄儀組成。A、B兩點(diǎn)為假定的彈著點(diǎn)。測(cè)點(diǎn)及假定彈著點(diǎn)的坐標(biāo)如表1所示。

表1 測(cè)點(diǎn)及假定彈著點(diǎn)坐標(biāo)

圖1 仿真定位模型

假設(shè)彈丸落地時(shí)刻t0=0 s，地震波的傳播速度v0=400 m/s，但由于地震波在地表介質(zhì)中傳播的速度是不均勻的，因此為了盡量模擬實(shí)際情況，在彈著點(diǎn)至各個(gè)測(cè)點(diǎn)的方向上引入±3%的速度誤差，引入方式如(9)式所示：

v=400±400×3%×rand，

(9)

式中：rand為0～1之間的隨機(jī)數(shù)。

v代入(1)式即可得到引入速度誤差后測(cè)點(diǎn)拾取到的地震波初至?xí)r間，8個(gè)測(cè)點(diǎn)拾取到彈著點(diǎn)A和B的地震波初至?xí)r間如表2所示。

表2 各測(cè)點(diǎn)拾取彈著點(diǎn)A和B地震波初至?xí)r間

2.2 模型驗(yàn)證

表1所示測(cè)點(diǎn)坐標(biāo)以及表2所示地震波初至?xí)r間代入目標(biāo)函數(shù)(8)式，并通過(guò)單純形法搜索其最小值。本文的彈著點(diǎn)定位是平面定位問(wèn)題，因此單純形法的初始迭代單純形為三角形，分別以三角形的各頂點(diǎn)坐標(biāo)作為彈著點(diǎn)計(jì)算3種情況下的目標(biāo)函數(shù)值，根據(jù)單純形法的迭代規(guī)則，通過(guò)拉伸、收縮、對(duì)稱(chēng)等變換，使其不斷向著目標(biāo)函數(shù)值最小的方向移動(dòng)，直到滿(mǎn)足迭代終止條件即可得到準(zhǔn)確的彈著點(diǎn)坐標(biāo)，本文中單純形法的迭代終止條件設(shè)為三角形各邊長(zhǎng)均值不大于5 cm.當(dāng)選取的初始迭代三角形不同時(shí)，單純形法的定位結(jié)果相互之間也會(huì)有一定的偏差，甚至?xí)鼓繕?biāo)函數(shù)收斂于局部最小值，導(dǎo)致出現(xiàn)定位異常的情況。鑒于上述情況，在待監(jiān)測(cè)區(qū)域內(nèi)隨機(jī)生成不同的初始迭代三角形進(jìn)行重復(fù)定位，觀察其定位結(jié)果的分布情況。模型中彈著點(diǎn)A和B的200次定位結(jié)果分布情況如圖2所示。

圖2 彈著點(diǎn)A和B的200次定位結(jié)果分布圖

由圖2可以看出，200次的定位結(jié)果在大多數(shù)情況下都比較集中，相互之間的偏差較小，只有少數(shù)定位出現(xiàn)異常的情況。

為了最終獲得彈著點(diǎn)坐標(biāo)，需首先剔除定位異常的結(jié)果。DBSCAN算法是一維或多維特征空間中的非參數(shù)、基于密度的聚類(lèi)算法。與劃分聚類(lèi)算法不同，DBSCAN算法不需要預(yù)先聲明聚類(lèi)的數(shù)量，而是通過(guò)評(píng)估樣本的緊密程度來(lái)劃分對(duì)應(yīng)的類(lèi)別，理論上可以找出任何形狀的聚類(lèi)。

DBSCAN算法將簇定義為密度相連點(diǎn)的最大集合，可以將數(shù)據(jù)集中不包含在任何簇中的對(duì)象視為噪聲，通過(guò)這樣的方法，能夠?qū)⒍ㄎ唤Y(jié)果中的異常值剔除。本文中鄰域半徑ε設(shè)置為0.05，密度閾值Minpts設(shè)置為4，DBSCAN算法運(yùn)行結(jié)果如圖3所示。

圖3 彈著點(diǎn) A和B的DBSCAN算法運(yùn)行結(jié)果

剔除圖3中的異常值后，分別對(duì)其余有效值的橫、縱坐標(biāo)分量求均值作為最終的定位結(jié)果。彈著點(diǎn)A的最終定位結(jié)果為(18.498 6 m, 61.191 6 m)，與A點(diǎn)的距離誤差為1.916 8 m；彈著點(diǎn)B的最終定位結(jié)果為(38.222 3 m,73.410 3 m)，與B點(diǎn)的距離誤差為3.845 9 m.測(cè)點(diǎn)至彈著點(diǎn)距離計(jì)算誤差絕對(duì)值平均值[11]為

(10)

式中：n為測(cè)點(diǎn)數(shù)量；li為計(jì)算彈著點(diǎn)位置到各測(cè)點(diǎn)的距離；di為實(shí)際彈著點(diǎn)位置到各測(cè)點(diǎn)的距離。

根據(jù)(10)式可得，彈著點(diǎn)A的定位誤差為1.97%，彈著點(diǎn)B的定位誤差為3.97%，定位誤差較小，所以該彈著點(diǎn)定位方法可行。

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

本文通過(guò)靜爆試驗(yàn)對(duì)該方法進(jìn)行驗(yàn)證，靜爆試驗(yàn)與彈丸侵徹產(chǎn)生的地震信號(hào)一樣，都是作用時(shí)間短、信號(hào)能量高度集中的信號(hào)，區(qū)別在于靜爆是由于超壓作用于地面所產(chǎn)生的地震波，而彈丸侵徹是通過(guò)撞擊產(chǎn)生的地震波。爆炸或者彈丸侵徹地面時(shí)都會(huì)形成具有一定能量的沖擊波，沖擊波的壓力遠(yuǎn)遠(yuǎn)大于土壤的抗壓強(qiáng)度，因此靠近震源區(qū)域的土壤結(jié)構(gòu)被破壞從而形成破壞區(qū)，之后隨著能量的消耗，沖擊波會(huì)衰減為應(yīng)力波，此時(shí)應(yīng)力波依然會(huì)對(duì)土壤造成塑性形變，最后隨著能量再次衰減，塑性波變?yōu)閺椥圆?，即為地震波?/p>

2019年11月在華北某靶場(chǎng)進(jìn)行了地面靜爆試驗(yàn)，藥柱當(dāng)量在50～60 kg梯恩梯(TNT)之間，共布置6個(gè)測(cè)點(diǎn)，選用西安石油儀器勘探廠生產(chǎn)的動(dòng)圈式速度型DFJ3A檢波器，其主要技術(shù)指標(biāo)[12]如表3所示。

表3 DFJ3A檢波器主要技術(shù)指標(biāo)

測(cè)點(diǎn)和爆點(diǎn)分布情況如圖4所示。圖4中，爆點(diǎn)所在位置坐標(biāo)為(0 m, 0 m)，為保證各個(gè)測(cè)點(diǎn)保持時(shí)鐘同步，采用了斷線觸發(fā)的方式，斷線的時(shí)刻記為0 s時(shí)刻，地面靜爆試驗(yàn)中所測(cè)得的典型地震波信號(hào)如圖5所示。

圖4 爆點(diǎn)和測(cè)點(diǎn)分布圖

圖5 實(shí)測(cè)典型地震波信號(hào)

測(cè)點(diǎn)坐標(biāo)以及各個(gè)測(cè)點(diǎn)拾取到的地震波初至?xí)r間如表4所示。

表4 測(cè)點(diǎn)坐標(biāo)及地震波初至?xí)r間

同樣，表4中數(shù)據(jù)代入(8)式并通過(guò)單純形法進(jìn)行200次定位，定位結(jié)果的分布情況如圖6所示。

圖6 實(shí)測(cè)200次定位結(jié)果分布圖

由圖6可以看出，與仿真定位模型的結(jié)果相似，也出現(xiàn)了少量定位異常的情況，DBSCAN算法的聚類(lèi)結(jié)果如圖7所示。剔除圖7中異常結(jié)果后，分別對(duì)其余結(jié)果的橫、縱坐標(biāo)分量求均值，即可得到爆點(diǎn)的定位結(jié)果為(-0.166 1 m, 0.239 0 m)，與真實(shí)爆點(diǎn)位置的距離誤差為0.291 1 m，測(cè)點(diǎn)至爆點(diǎn)距離誤差絕對(duì)值的平均值為0.38%.

圖7 實(shí)測(cè)DBSCAN算法運(yùn)行結(jié)果

4 結(jié)論

本文提出一種基于波速方差目標(biāo)函數(shù)的彈著點(diǎn)定位方法，通過(guò)仿真定位模型和地面靜爆試驗(yàn)對(duì)該方法進(jìn)行驗(yàn)證。得出以下主要結(jié)論：

1)用單純形法對(duì)基于波速方差的目標(biāo)函數(shù)尋優(yōu)并結(jié)合DBSCAN算法可以有效實(shí)現(xiàn)彈著點(diǎn)定位。

2)該彈著點(diǎn)定位方法無(wú)需提前測(cè)量地震波波速或?qū)Σㄋ龠M(jìn)行反演，可以減小地震波波速不均勻性給定位結(jié)果帶來(lái)的誤差，提高定位精度。

然而，靜爆與彈丸侵徹產(chǎn)生的地震波存在較大差異，所以不宜采用靜爆試驗(yàn)作驗(yàn)證，但在本文的目標(biāo)函數(shù)中已經(jīng)排除了波速的影響，且沒(méi)有用到幅度和頻率信息。所以仍有一定的借鑒意義，后續(xù)將進(jìn)行大當(dāng)量戰(zhàn)斗部試驗(yàn)進(jìn)一步驗(yàn)證。