郭 凱

(南京炮兵學(xué)院 牽引火炮系,南京 211132)

槍械是我軍應(yīng)用最為廣泛的武器裝備之一,其主要功能是對(duì)作戰(zhàn)人員等有生目標(biāo)實(shí)施動(dòng)能侵徹殺傷[1]。槍彈對(duì)明膠靶標(biāo)的毀傷參數(shù)是衡量槍彈殺傷效能的重要指標(biāo),主要包括:彈丸侵徹最大瞬時(shí)空腔體積、窄傷道長(zhǎng)度、動(dòng)能傳遞量3個(gè)方面。目前,國(guó)內(nèi)針對(duì)槍彈殺傷效能評(píng)估方面開展了很多理論和實(shí)驗(yàn)方面的研究,在數(shù)值仿真計(jì)算和圖像處理等領(lǐng)域取得了一定的成果[1-6]。但是,在彈丸侵徹明膠靶標(biāo)毀傷參數(shù)獲取方面,文獻(xiàn)[3-5]采用Gauss濾波、最大熵等方法和原理進(jìn)行圖像增強(qiáng),存在不能保護(hù)圖像輪廓線細(xì)節(jié)信息的缺陷。文獻(xiàn)[6]采用水平集等方法進(jìn)行圖像分割,存在不能去除偽輪廓的現(xiàn)象。

本文在詳細(xì)分析高速攝影系統(tǒng)拍攝到的槍彈侵徹明膠靶標(biāo)圖像基礎(chǔ)上,提出改進(jìn)型自適應(yīng)均值濾波法,實(shí)現(xiàn)了圖像中混合噪聲的有效去除,保護(hù)了圖像輪廓線的細(xì)節(jié)信息。采用二值化輪廓跟蹤法,對(duì)濾波后的圖像進(jìn)行分割,獲取初步清晰、連續(xù)的空腔輪廓線[5]。針對(duì)偽輪廓凸性方向,采用捕捉拐點(diǎn)方式,提出基于曲線割線斜率的單調(diào)性去除偽輪廓線法,并綜合運(yùn)用基于最小二乘的多項(xiàng)式擬合、基于貝濟(jì)曲線的切線斜率差法,去除了不符合要求的偽輪廓和尖角效應(yīng),實(shí)現(xiàn)輪廓線的精確提取。運(yùn)用橢圓截面法確定了三維空腔模型參數(shù),并最終實(shí)現(xiàn)了槍彈對(duì)明膠靶標(biāo)的毀傷參數(shù)計(jì)算。

1 彈創(chuàng)空腔圖像的預(yù)處理

由于彈丸侵徹明膠靶標(biāo)的高溫性、高速性以及拍攝系統(tǒng)性能差異性的客觀因素,圖像不可避免地存在噪聲污染、輪廓線不清晰等現(xiàn)象。為了提高毀傷參數(shù)獲取的準(zhǔn)確性,必須先進(jìn)行圖像預(yù)處理,主要工作包括基于改進(jìn)型自適應(yīng)均值濾波法的圖像增強(qiáng)與基于二值化輪廓跟蹤法的圖像分割。

1.1 基于改進(jìn)型自適應(yīng)均值濾波法的圖像增強(qiáng)

高速攝影抓拍的瞬時(shí)圖像很容易受到“高斯-脈沖”混合噪聲污染,傳統(tǒng)濾波方式分為中值濾波與均值濾波[7]。中值濾波采用排序統(tǒng)計(jì)理論,將濾波窗口像素灰度中間值作為待處理像素的新值輸出,達(dá)到有效抵制噪聲的目的。但高斯噪聲按照正態(tài)分布污染所有像素點(diǎn),無(wú)論怎樣排序得到的結(jié)果仍然是污染值,因此中值濾波對(duì)高斯噪聲抑制效果不大。均值濾波采用求取灰度平均值的方式,實(shí)現(xiàn)濾波窗口內(nèi)所有像素的平滑處理,完成去除噪聲工作,但脈沖噪聲點(diǎn)的灰度值往往非常大(小),對(duì)均值影響很大,導(dǎo)致濾波后圖像失真以及輪廓線模糊。

本文提出了改進(jìn)型自適應(yīng)均值濾波法,采用自適應(yīng)選擇濾波窗口以及改進(jìn)濾波窗口灰度計(jì)算方式,有效地去除混合噪聲,實(shí)現(xiàn)了保護(hù)圖像輪廓線細(xì)節(jié)信息的目的。如圖1所示,白點(diǎn)為5×5濾波窗口,黑點(diǎn)區(qū)域表示子濾波窗口。對(duì)圖像的每一個(gè)像素用9種形狀不同的子窗口濾波,依次計(jì)算各窗口內(nèi)灰度均方差,將求取的均方差按照降序法排序。依據(jù)均方差尖銳區(qū)域大于平緩區(qū)域原理,選擇均方差最小子窗口作為濾波窗口,實(shí)現(xiàn)窗口大小自適應(yīng)選擇以及保護(hù)圖像細(xì)節(jié)信息的目的。

統(tǒng)計(jì)該窗口中各個(gè)灰度級(jí)像素出現(xiàn)的概率,求取均值濾波區(qū)域像素灰度的數(shù)學(xué)期望,并將其作為當(dāng)前像素的新灰度值輸出,完成濾除噪聲工作,公式如下:

P(rk)=nk/N

(1)

(3)

式中:nk為第k級(jí)灰度的像素?cái)?shù),rk為第k個(gè)灰度級(jí),N為不同形狀子窗口中對(duì)應(yīng)的像素個(gè)數(shù),P(rk)為該灰度級(jí)出現(xiàn)的概率;E為各濾波子窗口的數(shù)學(xué)期望,xl為像素的灰度值,Pl為該灰度級(jí)出現(xiàn)的概率;σ為濾波子窗口的均方差,μ為濾波子窗口灰度的平均值或數(shù)學(xué)期望。

圖1 改進(jìn)型自適應(yīng)均值濾波窗口

本文分別采用均值濾波、中值濾波和改進(jìn)型自適應(yīng)均值濾波進(jìn)行了圖像去噪試驗(yàn),結(jié)果見圖2。為定量評(píng)價(jià)3種方法的優(yōu)劣,采用RPSN值(peak signal to noise ratio,峰值信噪比)和EMS值(mean square error,均方誤差)進(jìn)行去噪性能評(píng)估,結(jié)果如表1所示。由表1可以明顯看出,改進(jìn)型自適應(yīng)均值濾波法的RPSN值最大,而EMS值最小,其去除混合噪聲能力最強(qiáng)。

圖2 3種濾波算法處理結(jié)果

濾波方法RPSNEMS中值濾波24.43234.34均值濾波27.67111.15本文算法濾波55.5330.69

1.2 基于二值化輪廓跟蹤法的圖像分割

彈創(chuàng)圖像中的空腔輪廓線是準(zhǔn)確構(gòu)建空腔三維模型的基礎(chǔ)。文中采用二值化輪廓跟蹤法對(duì)濾波后的圖像進(jìn)行分割,從而獲取彈創(chuàng)空腔輪廓線[7]。

處理過程如圖3所示,對(duì)圖像進(jìn)行二值化處理,提取主體特征像素區(qū)域;利用形態(tài)學(xué)腐蝕切斷區(qū)域連通性;運(yùn)用區(qū)域面積測(cè)量實(shí)現(xiàn)空腔體的精確分割;對(duì)二值化區(qū)域進(jìn)行輪廓跟蹤,實(shí)現(xiàn)輪廓線的初步提取。

圖3 圖像分割算法示意圖

2 彈創(chuàng)空腔圖像的三維建模

彈創(chuàng)空腔正交圖像輪廓線是圖像中灰度躍變劇烈的區(qū)域,其提取質(zhì)量的高低直接影響毀傷參數(shù)的計(jì)算結(jié)果。受圖像拍攝質(zhì)量影響,初步提取的輪廓線存在不連續(xù)、不光滑、尖角和偽輪廓線現(xiàn)象,見圖3(d)。為便于后續(xù)彈創(chuàng)空腔模型的建立,必須開展進(jìn)一步的輪廓線精確化提取工作。

2.1 輪廓線精確化提取

為改善輪廓線質(zhì)量,根據(jù)初步提取的輪廓線形狀特征,確定精確化提取方法,主要包括曲線割線單調(diào)性去偽法、最小二乘擬合法[7]、切線斜率差去尖角法[7]。

2.1.1 基于曲線割線斜率的單調(diào)性去除偽輪廓線法

針對(duì)偽輪廓凸性方向,采用捕捉拐點(diǎn)方式,去除不符合要求的偽輪廓,如圖4所示,具體方法如下。

圖4(a)中,設(shè)AC段弧線為偽輪廓,掃描捕捉曲線底部任意一個(gè)點(diǎn)設(shè)為起始點(diǎn)E,定步長(zhǎng)獲取E點(diǎn)左側(cè)曲線點(diǎn)A,B,C,D,并構(gòu)成割線EA,EB,EC,ED;分別計(jì)算相應(yīng)的割線斜率fEA,fEB,fEC,fED。顯然拐點(diǎn)C的割線斜率最大,將C點(diǎn)捕捉。圖4(b)中,再將C點(diǎn)左側(cè)曲線反顯為背景白色,其它部分曲線采用相同方式處理。如此,完成圖像偽輪廓去除工作。

圖4 去除偽輪廓線算法示意圖

2.1.2 基于最小二乘法的多項(xiàng)式擬合輪廓線法

針對(duì)輪廓線的不連續(xù)性,采用逐行掃描方式捕獲輪廓線各像素點(diǎn),形成待擬合坐標(biāo)數(shù)據(jù)集,結(jié)合最小二乘法擬合出多項(xiàng)式函數(shù),實(shí)現(xiàn)輪廓線的連續(xù)繪制工作。

具體方法:將圖像劃分為4個(gè)象限,逐行掃描獲取某一象限的待處理像素點(diǎn),形成坐標(biāo)數(shù)據(jù)集合(xj,yj)。設(shè)擬合函數(shù)為

p(x)=a0φ0(x)+a1φ1(x)+…+asφs(x)

(4)

2.1.3 基于貝濟(jì)曲線的切線斜率差去尖角法

輪廓線擬合時(shí),相鄰兩曲線在交點(diǎn)處各自斜率的不同造成交接處產(chǎn)生尖角現(xiàn)象。如圖5(a)所示,點(diǎn)A為曲線AB與AC相交后形成的尖角點(diǎn),f0,f1分別是AB與AC在O0與O1點(diǎn)處的切線斜率,ε為f0,f1的斜率差。

貝濟(jì)曲線包括4個(gè)節(jié)點(diǎn):起始點(diǎn)、終止點(diǎn)和2個(gè)相互分離的中間點(diǎn),當(dāng)滑動(dòng)中間兩點(diǎn)時(shí),可以自由改變曲線形狀,繪制出光滑曲線。在圖5(b)中,令貝濟(jì)曲線的起始節(jié)點(diǎn)P0與O0重合、終止節(jié)點(diǎn)P3與O1重合、中間節(jié)點(diǎn)P1,P2與A點(diǎn)重合,并先繪制出初始貝濟(jì)曲線O1O0。令fP0與fP3分別為經(jīng)過貝濟(jì)曲線節(jié)點(diǎn)P0與P3的切線斜率,θ為fP0,fP3的斜率差。不斷調(diào)節(jié)節(jié)點(diǎn)P1,P2位置,改變貝濟(jì)曲線O1O0的形狀,此時(shí)斜率差θ的值也隨之改變,當(dāng)|θ|≤|ε|時(shí),停止節(jié)點(diǎn)調(diào)整,得到最終去除尖角效應(yīng)的貝濟(jì)曲線。

圖5 去除尖角效應(yīng)的算法示意圖

2.2 基于橢圓截面的三維空腔模型參數(shù)確定

由于拍攝條件所限,只能在2幅正交圖像基礎(chǔ)上進(jìn)行三維重構(gòu),需要對(duì)提取的輪廓線進(jìn)行相應(yīng)的模型簡(jiǎn)化處理。

具體處理過程:彈創(chuàng)空腔是彈丸在明膠里擠壓、翻滾、旋轉(zhuǎn)傳遞能量時(shí)造成明膠迅速膨脹而形成的瞬時(shí)腔體。圖6(a)、6(b)分別為正交(水平與垂直方向)空腔輪廓線投影圖。由于彈創(chuàng)空腔具有弧形表面,其沿x軸(子彈飛行)垂直方向的任一截面均可近似為橢圓,如圖6(c)所示。橢圓的長(zhǎng)、短軸ap,bq可利用圖6(a)、6(b)中的2個(gè)正交輪廓線像素差值獲取,如此,可逐步確定整個(gè)空腔截面的橢圓參數(shù)。

圖6 三維彈創(chuàng)空腔模型的建立

圖7 彈創(chuàng)空腔正交圖像

3 彈丸毀傷參數(shù)的計(jì)算

在建立三維彈創(chuàng)空腔模型的基礎(chǔ)上,對(duì)圖7中某型彈丸穿過明膠時(shí)形成的2幅正交空腔圖像進(jìn)行處理,給出了彈創(chuàng)空腔體積、彈丸侵徹靶標(biāo)窄傷道長(zhǎng)度、彈丸動(dòng)能傳遞量3種毀傷參數(shù)的計(jì)算方法。

3.1 標(biāo)定比例尺的確定

數(shù)字圖像中長(zhǎng)度的計(jì)量單位為像素,在實(shí)際計(jì)算時(shí),需要將其按照一定的比例關(guān)系轉(zhuǎn)換成實(shí)際的物理單位。本文采用在明膠中預(yù)先放置標(biāo)定塊(邊長(zhǎng)為15 mm的正方體型鉛塊)的方式,如圖7所示,通過一系列的圖像處理工作,提取標(biāo)定塊輪廓線,處理步驟如圖8所示。以垂直方向?yàn)槔?圖9為標(biāo)定塊輪廓提取過程示意圖。

圖8 標(biāo)定塊輪廓線獲取步驟

圖9 標(biāo)定塊輪廓提取過程示意圖

掃描標(biāo)定塊輪廓線,得邊長(zhǎng)的像素值為30,由此確定標(biāo)定比例尺kb為每像素0.50 mm。

3.2 彈創(chuàng)空腔體積計(jì)算

2.2節(jié)中建立的三維彈創(chuàng)空腔可看成由若干個(gè)以橢圓截面為底的圓柱體組成,其體積可采用積分的方式進(jìn)行計(jì)算,圖10為空腔體積計(jì)算模型。

圖10 空腔體積計(jì)算模型

設(shè)x軸為定軸,空腔體位于點(diǎn)xa=a,xb=b處的2個(gè)垂直平面之間。水平與垂直視圖的輪廓在某點(diǎn)x處的縱坐標(biāo)差ap和bq分別代表空腔截面橢圓長(zhǎng)、短軸長(zhǎng)度。S(x)代表點(diǎn)x處的截面面積,S(x)=apbqπ。dV代表以S為底面積的體積元素(積分步長(zhǎng)),其數(shù)值計(jì)算公式為dV=S(x)dx。在閉區(qū)間[a,b]上作定積分,得空腔體積為

(6)

計(jì)算結(jié)果為4 410 959.0 mm3。

3.3 窄傷道長(zhǎng)度計(jì)算

窄傷道長(zhǎng)度是彈丸侵徹明膠后明膠靶標(biāo)邊緣(窄傷道起點(diǎn))到腔體末端(窄傷道終點(diǎn))的長(zhǎng)度,如圖11所示。

圖11 窄傷道示意圖

本文主要運(yùn)用圖像二值化、圖像輪廓跟蹤,提取出明膠靶標(biāo)與空腔輪廓線區(qū)域特征信息,圖12(a)、12(b)分別為捕捉到的窄傷道起點(diǎn)與終點(diǎn),計(jì)算兩者像素差值,再運(yùn)用標(biāo)定比例尺進(jìn)行換算,完成窄傷道長(zhǎng)度計(jì)算工作,結(jié)果為80.6 mm。

圖12 窄傷道計(jì)算示意圖

3.4 彈丸動(dòng)能傳遞量計(jì)算

彈丸侵入明膠靶標(biāo)后,受明膠體擠壓、碰撞以及熱能的影響,速度逐漸遞減,能量隨深度的增加也逐漸衰減。

根據(jù)Sturdivan提出的EKE法[8],該方法將士兵一次隨機(jī)命中時(shí)喪失戰(zhàn)斗力的概率P(I/H)與在機(jī)體中投射物傳遞能量的期望值EKE(expected kinetic energy)相聯(lián)系。

式中:α,β,γ為與士兵承擔(dān)任務(wù)及喪失戰(zhàn)斗力有關(guān)的曲線擬合常數(shù),而[8]

式中:m為殺傷元質(zhì)量;vh為彈丸在侵徹深度h處的瞬時(shí)速度;v0為彈丸碰擊靶標(biāo)時(shí)速度;Ph為當(dāng)侵徹深度為h時(shí)彈丸存留在人體內(nèi)的概率。

從式(8)可以看出,彈丸瞬時(shí)速度的求取是動(dòng)能傳遞量計(jì)算的關(guān)鍵。

本文利用一組彈丸侵徹明膠靶標(biāo)序列圖像,通過一系列的圖像增強(qiáng)與圖像分割技術(shù),分別提取彈丸特征。采用矩理論確定彈丸的重心坐標(biāo),將每一幀圖像中彈丸重心橫坐標(biāo)與起始幀彈丸重心橫坐標(biāo)相減,得到彈丸位移像素值,再通過標(biāo)定計(jì)算得出相應(yīng)的水平位移值。

圖13展示了其中某一幀圖像中彈丸特征提取過程,圖中黑框部分為彈丸位置。

圖13 彈丸特征提取過程

彈丸水平行程與行程時(shí)間的比值為彈丸在這一行程中的水平平均速度,當(dāng)行程時(shí)間非常小時(shí),該速度可近似為彈丸瞬時(shí)速度。為更精確地計(jì)算瞬時(shí)速度,本文采用拉格朗日插值法對(duì)前面求取的一序列彈丸位移值進(jìn)行了插值計(jì)算試驗(yàn)。最終的彈丸瞬時(shí)速度v和動(dòng)能傳遞量J計(jì)算結(jié)果匯總在表2中。

表2 彈丸瞬時(shí)速度和動(dòng)能傳遞量匯總表

4 結(jié)論

針對(duì)彈丸侵徹明膠靶標(biāo)圖像的嚴(yán)重噪聲、不清晰輪廓線等特點(diǎn),本文采用改進(jìn)型自適應(yīng)均值濾波法和二值化輪廓跟蹤法實(shí)現(xiàn)了圖像增強(qiáng)和分割的預(yù)處理,利用基于曲線割線斜率的單調(diào)性去除偽輪廓線法、基于最小二乘法的多項(xiàng)式擬合輪廓線法和基于貝濟(jì)曲線的切線斜率差去尖角法,實(shí)現(xiàn)了圖像的輪廓線提取,

在橢圓截面法的基礎(chǔ)上,確定了三維彈創(chuàng)空腔模型參數(shù)。運(yùn)用高速攝影系統(tǒng)拍攝的彈丸侵徹明膠靶標(biāo)圖像進(jìn)行了試驗(yàn)驗(yàn)證,計(jì)算出了彈創(chuàng)空腔體積、窄傷道長(zhǎng)度和彈丸動(dòng)能傳遞量。結(jié)果表明,本文提出的彈丸毀傷參數(shù)獲取方法是可行的,為進(jìn)一步研究槍械殺傷效能評(píng)估提供了有效技術(shù)途徑和手段。

[1] 劉萌秋.創(chuàng)傷彈道學(xué)[M].北京:解放軍出版社,1991.

LIU Meng-qiu.Wound ballistics[M].Beijing:Liberation Army Publishing House,1991.(in Chinese)

[2] 陳強(qiáng)華,王永娟.輕武器殺傷效能評(píng)估理論與計(jì)算[J].兵工自動(dòng)化,2011,30(7):28-30.

CHEN Qiang-hua,WANG Yong-juan.Evaluation theory and calculation of kill efficiency for small arms[J].Ordnance Industry Automation,2011,30(7):28-30.(in Chinese)

[3] 林勇,王經(jīng)瑾.明膠彈創(chuàng)空腔閃光X射線投影圖像三維重建[J].清華大學(xué)學(xué)報(bào),2002,42(12):1 576-1 578.

LIN Yong,WANG Jing-jin.Gelatin cavity 3-D reconstuction from flash X-ray images of wound ballistics[J].Journal of Tsinghua University,2002,42(12):1 576-1 578.(in Chinese)

[4] 賀成,王濤.創(chuàng)傷彈道空腔圖像邊緣檢測(cè)技術(shù)研究[J].計(jì)算機(jī)工程與設(shè)計(jì),2011,32(1):248-250.

HE Cheng,WANG Tao.Research on edge detection of wound ballistics cavity image[J].Computer Engineering and Design,2011,32(1):248-250.(in Chinese)

[5] 申彪,王濤.基于Matlab的創(chuàng)傷彈道三維重構(gòu)[J].彈箭與制導(dǎo)學(xué)報(bào),2011,31(1):140-142.

SHEN Biao,WANG Tao.Wound ballistics 3D reconstruction based on Matlab[J].Journal of Projectiles,Rockets,Missiles and Guidance,2011,31(1):140-142.(in Chinese)

[6] 汪傳忠,史俊.基于水平集的彈創(chuàng)瞬時(shí)空腔X射線投影圖像分割方法[J].測(cè)試技術(shù)學(xué)報(bào),2011,25(1):82-86.

WANG Chuan-zhong,SHI Jun.Level set-based segmentation of the transient cavity of ballistic wound in X-ray image[J].Journal of Test and Measurement Technology,2011,25(1):82-86.(in Chinese)

[7] RAFAEL C G,RICHARD E W,阮秋琦,等.數(shù)字圖像處理[M].第三版.北京:電子工業(yè)出版社,2011.

RAFAEL C G,RICHARD E W,RUAN Qiu-qi,et al.Digital image processing[M].3rd ed.Beijing:Electronic Industry Press,2011.(in Chinese)

[8] WILLIAM B,BEVERLY.A human ballistic mortality Model,ADA058947[R].1978.