郭 凱
(南京炮兵學院 牽引火炮系,南京 211132)
槍械是我軍應用最為廣泛的武器裝備之一,其主要功能是對作戰人員等有生目標實施動能侵徹殺傷[1]。槍彈對明膠靶標的毀傷參數是衡量槍彈殺傷效能的重要指標,主要包括:彈丸侵徹最大瞬時空腔體積、窄傷道長度、動能傳遞量3個方面。目前,國內針對槍彈殺傷效能評估方面開展了很多理論和實驗方面的研究,在數值仿真計算和圖像處理等領域取得了一定的成果[1-6]。但是,在彈丸侵徹明膠靶標毀傷參數獲取方面,文獻[3-5]采用Gauss濾波、最大熵等方法和原理進行圖像增強,存在不能保護圖像輪廓線細節信息的缺陷。文獻[6]采用水平集等方法進行圖像分割,存在不能去除偽輪廓的現象。
本文在詳細分析高速攝影系統拍攝到的槍彈侵徹明膠靶標圖像基礎上,提出改進型自適應均值濾波法,實現了圖像中混合噪聲的有效去除,保護了圖像輪廓線的細節信息。采用二值化輪廓跟蹤法,對濾波后的圖像進行分割,獲取初步清晰、連續的空腔輪廓線[5]。針對偽輪廓凸性方向,采用捕捉拐點方式,提出基于曲線割線斜率的單調性去除偽輪廓線法,并綜合運用基于最小二乘的多項式擬合、基于貝濟曲線的切線斜率差法,去除了不符合要求的偽輪廓和尖角效應,實現輪廓線的精確提取。運用橢圓截面法確定了三維空腔模型參數,并最終實現了槍彈對明膠靶標的毀傷參數計算。
由于彈丸侵徹明膠靶標的高溫性、高速性以及拍攝系統性能差異性的客觀因素,圖像不可避免地存在噪聲污染、輪廓線不清晰等現象。為了提高毀傷參數獲取的準確性,必須先進行圖像預處理,主要工作包括基于改進型自適應均值濾波法的圖像增強與基于二值化輪廓跟蹤法的圖像分割。
高速攝影抓拍的瞬時圖像很容易受到“高斯-脈沖”混合噪聲污染,傳統濾波方式分為中值濾波與均值濾波[7]。中值濾波采用排序統計理論,將濾波窗口像素灰度中間值作為待處理像素的新值輸出,達到有效抵制噪聲的目的。但高斯噪聲按照正態分布污染所有像素點,無論怎樣排序得到的結果仍然是污染值,因此中值濾波對高斯噪聲抑制效果不大。均值濾波采用求取灰度平均值的方式,實現濾波窗口內所有像素的平滑處理,完成去除噪聲工作,但脈沖噪聲點的灰度值往往非常大(小),對均值影響很大,導致濾波后圖像失真以及輪廓線模糊。
本文提出了改進型自適應均值濾波法,采用自適應選擇濾波窗口以及改進濾波窗口灰度計算方式,有效地去除混合噪聲,實現了保護圖像輪廓線細節信息的目的。如圖1所示,白點為5×5濾波窗口,黑點區域表示子濾波窗口。對圖像的每一個像素用9種形狀不同的子窗口濾波,依次計算各窗口內灰度均方差,將求取的均方差按照降序法排序。依據均方差尖銳區域大于平緩區域原理,選擇均方差最小子窗口作為濾波窗口,實現窗口大小自適應選擇以及保護圖像細節信息的目的。
統計該窗口中各個灰度級像素出現的概率,求取均值濾波區域像素灰度的數學期望,并將其作為當前像素的新灰度值輸出,完成濾除噪聲工作,公式如下:
P(rk)=nk/N
(1)
(3)
式中:nk為第k級灰度的像素數,rk為第k個灰度級,N為不同形狀子窗口中對應的像素個數,P(rk)為該灰度級出現的概率;E為各濾波子窗口的數學期望,xl為像素的灰度值,Pl為該灰度級出現的概率;σ為濾波子窗口的均方差,μ為濾波子窗口灰度的平均值或數學期望。

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

圖2 3種濾波算法處理結果

濾波方法RPSNEMS中值濾波24.43234.34均值濾波27.67111.15本文算法濾波55.5330.69
彈創圖像中的空腔輪廓線是準確構建空腔三維模型的基礎。文中采用二值化輪廓跟蹤法對濾波后的圖像進行分割,從而獲取彈創空腔輪廓線[7]。
處理過程如圖3所示,對圖像進行二值化處理,提取主體特征像素區域;利用形態學腐蝕切斷區域連通性;運用區域面積測量實現空腔體的精確分割;對二值化區域進行輪廓跟蹤,實現輪廓線的初步提取。

圖3 圖像分割算法示意圖
彈創空腔正交圖像輪廓線是圖像中灰度躍變劇烈的區域,其提取質量的高低直接影響毀傷參數的計算結果。受圖像拍攝質量影響,初步提取的輪廓線存在不連續、不光滑、尖角和偽輪廓線現象,見圖3(d)。為便于后續彈創空腔模型的建立,必須開展進一步的輪廓線精確化提取工作。
為改善輪廓線質量,根據初步提取的輪廓線形狀特征,確定精確化提取方法,主要包括曲線割線單調性去偽法、最小二乘擬合法[7]、切線斜率差去尖角法[7]。
2.1.1 基于曲線割線斜率的單調性去除偽輪廓線法
針對偽輪廓凸性方向,采用捕捉拐點方式,去除不符合要求的偽輪廓,如圖4所示,具體方法如下。
圖4(a)中,設AC段弧線為偽輪廓,掃描捕捉曲線底部任意一個點設為起始點E,定步長獲取E點左側曲線點A,B,C,D,并構成割線EA,EB,EC,ED;分別計算相應的割線斜率fEA,fEB,fEC,fED。顯然拐點C的割線斜率最大,將C點捕捉。圖4(b)中,再將C點左側曲線反顯為背景白色,其它部分曲線采用相同方式處理。如此,完成圖像偽輪廓去除工作。

圖4 去除偽輪廓線算法示意圖
2.1.2 基于最小二乘法的多項式擬合輪廓線法
針對輪廓線的不連續性,采用逐行掃描方式捕獲輪廓線各像素點,形成待擬合坐標數據集,結合最小二乘法擬合出多項式函數,實現輪廓線的連續繪制工作。
具體方法:將圖像劃分為4個象限,逐行掃描獲取某一象限的待處理像素點,形成坐標數據集合(xj,yj)。設擬合函數為
p(x)=a0φ0(x)+a1φ1(x)+…+asφs(x)
(4)
2.1.3 基于貝濟曲線的切線斜率差去尖角法
輪廓線擬合時,相鄰兩曲線在交點處各自斜率的不同造成交接處產生尖角現象。如圖5(a)所示,點A為曲線AB與AC相交后形成的尖角點,f0,f1分別是AB與AC在O0與O1點處的切線斜率,ε為f0,f1的斜率差。
貝濟曲線包括4個節點:起始點、終止點和2個相互分離的中間點,當滑動中間兩點時,可以自由改變曲線形狀,繪制出光滑曲線。在圖5(b)中,令貝濟曲線的起始節點P0與O0重合、終止節點P3與O1重合、中間節點P1,P2與A點重合,并先繪制出初始貝濟曲線O1O0。令fP0與fP3分別為經過貝濟曲線節點P0與P3的切線斜率,θ為fP0,fP3的斜率差。不斷調節節點P1,P2位置,改變貝濟曲線O1O0的形狀,此時斜率差θ的值也隨之改變,當|θ|≤|ε|時,停止節點調整,得到最終去除尖角效應的貝濟曲線。

圖5 去除尖角效應的算法示意圖
由于拍攝條件所限,只能在2幅正交圖像基礎上進行三維重構,需要對提取的輪廓線進行相應的模型簡化處理。
具體處理過程:彈創空腔是彈丸在明膠里擠壓、翻滾、旋轉傳遞能量時造成明膠迅速膨脹而形成的瞬時腔體。圖6(a)、6(b)分別為正交(水平與垂直方向)空腔輪廓線投影圖。由于彈創空腔具有弧形表面,其沿x軸(子彈飛行)垂直方向的任一截面均可近似為橢圓,如圖6(c)所示。橢圓的長、短軸ap,bq可利用圖6(a)、6(b)中的2個正交輪廓線像素差值獲取,如此,可逐步確定整個空腔截面的橢圓參數。

圖6 三維彈創空腔模型的建立

圖7 彈創空腔正交圖像
在建立三維彈創空腔模型的基礎上,對圖7中某型彈丸穿過明膠時形成的2幅正交空腔圖像進行處理,給出了彈創空腔體積、彈丸侵徹靶標窄傷道長度、彈丸動能傳遞量3種毀傷參數的計算方法。
數字圖像中長度的計量單位為像素,在實際計算時,需要將其按照一定的比例關系轉換成實際的物理單位。本文采用在明膠中預先放置標定塊(邊長為15 mm的正方體型鉛塊)的方式,如圖7所示,通過一系列的圖像處理工作,提取標定塊輪廓線,處理步驟如圖8所示。以垂直方向為例,圖9為標定塊輪廓提取過程示意圖。

圖8 標定塊輪廓線獲取步驟

圖9 標定塊輪廓提取過程示意圖
掃描標定塊輪廓線,得邊長的像素值為30,由此確定標定比例尺kb為每像素0.50 mm。
2.2節中建立的三維彈創空腔可看成由若干個以橢圓截面為底的圓柱體組成,其體積可采用積分的方式進行計算,圖10為空腔體積計算模型。

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

(6)
計算結果為4 410 959.0 mm3。
窄傷道長度是彈丸侵徹明膠后明膠靶標邊緣(窄傷道起點)到腔體末端(窄傷道終點)的長度,如圖11所示。

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

圖12 窄傷道計算示意圖
彈丸侵入明膠靶標后,受明膠體擠壓、碰撞以及熱能的影響,速度逐漸遞減,能量隨深度的增加也逐漸衰減。
根據Sturdivan提出的EKE法[8],該方法將士兵一次隨機命中時喪失戰斗力的概率P(I/H)與在機體中投射物傳遞能量的期望值EKE(expected kinetic energy)相聯系。
式中:α,β,γ為與士兵承擔任務及喪失戰斗力有關的曲線擬合常數,而[8]
式中:m為殺傷元質量;vh為彈丸在侵徹深度h處的瞬時速度;v0為彈丸碰擊靶標時速度;Ph為當侵徹深度為h時彈丸存留在人體內的概率。
從式(8)可以看出,彈丸瞬時速度的求取是動能傳遞量計算的關鍵。
本文利用一組彈丸侵徹明膠靶標序列圖像,通過一系列的圖像增強與圖像分割技術,分別提取彈丸特征。采用矩理論確定彈丸的重心坐標,將每一幀圖像中彈丸重心橫坐標與起始幀彈丸重心橫坐標相減,得到彈丸位移像素值,再通過標定計算得出相應的水平位移值。
圖13展示了其中某一幀圖像中彈丸特征提取過程,圖中黑框部分為彈丸位置。

圖13 彈丸特征提取過程
彈丸水平行程與行程時間的比值為彈丸在這一行程中的水平平均速度,當行程時間非常小時,該速度可近似為彈丸瞬時速度。為更精確地計算瞬時速度,本文采用拉格朗日插值法對前面求取的一序列彈丸位移值進行了插值計算試驗。最終的彈丸瞬時速度v和動能傳遞量J計算結果匯總在表2中。

表2 彈丸瞬時速度和動能傳遞量匯總表
針對彈丸侵徹明膠靶標圖像的嚴重噪聲、不清晰輪廓線等特點,本文采用改進型自適應均值濾波法和二值化輪廓跟蹤法實現了圖像增強和分割的預處理,利用基于曲線割線斜率的單調性去除偽輪廓線法、基于最小二乘法的多項式擬合輪廓線法和基于貝濟曲線的切線斜率差去尖角法,實現了圖像的輪廓線提取,
在橢圓截面法的基礎上,確定了三維彈創空腔模型參數。運用高速攝影系統拍攝的彈丸侵徹明膠靶標圖像進行了試驗驗證,計算出了彈創空腔體積、窄傷道長度和彈丸動能傳遞量。結果表明,本文提出的彈丸毀傷參數獲取方法是可行的,為進一步研究槍械殺傷效能評估提供了有效技術途徑和手段。
[1] 劉萌秋.創傷彈道學[M].北京:解放軍出版社,1991.
LIU Meng-qiu.Wound ballistics[M].Beijing:Liberation Army Publishing House,1991.(in Chinese)
[2] 陳強華,王永娟.輕武器殺傷效能評估理論與計算[J].兵工自動化,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] 林勇,王經瑾.明膠彈創空腔閃光X射線投影圖像三維重建[J].清華大學學報,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] 賀成,王濤.創傷彈道空腔圖像邊緣檢測技術研究[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的創傷彈道三維重構[J].彈箭與制導學報,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] 汪傳忠,史俊.基于水平集的彈創瞬時空腔X射線投影圖像分割方法[J].測試技術學報,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,阮秋琦,等.數字圖像處理[M].第三版.北京:電子工業出版社,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.