基于智能算法的超高速撞擊碎片云圖像分析

覃俞璋， 宋 燕， 王冶天， 趙 丹

（上海理工大學 光電信息與計算機工程學院， 上海200093）

0 引 言

隨著人類空間活動的日漸增加，殘留在地球軌道的廢棄物也越來越多，產生了大量的空間撞擊碎片在地球軌道上高速運動著，對在軌航天器的安全構成極大威脅。 超高速撞擊過程中的彈靶材料會發生大形變、碎裂以及應力波作用導致的層裂破壞現象。 當撞擊薄板結構時通常會形成碎片云；當碰撞速度很高時，彈靶材料的壓力和溫度極高，碎片云會發生烙化、汽化甚至變成離子體等物理現象。 因此，根據碎片云的結構演變、質量／動量分布、碎片云熱力學狀態和相分布等碎片云特性展開研究，成為超高速碰撞研究中的重要方向。 因此，采用相應的手段對超高速撞擊碎片云圖像特征（如數量、質量和速度分布等）進行分析研究，對于指導航天器防護設計，保障航天器可靠的安全在軌運轉，提升航天器在空間碎片環境中的適應能力具有重要戰略意義。

本文采用2 mm 彈丸，以3.06 km／s 的速度撞擊厚度為2 mm 的靶板（板間距為80 mm），模擬太空中航天器的損傷情況。 對采集的圖像（如圖1）進行預處理、自適應閾值分割及canny 邊緣檢測，并對碎片云進行分析，同時根據后板的損傷情況估計碎片云中產生的有效損傷碎片數量比例。

圖1 彈丸撞擊靶板產生的碎片云圖像Fig.1 The debris cloud image produced by the projectile hitting the target

1 碎片云圖像處理

本文運用opencv 圖像分析技術對撞擊產生的碎片云圖像進行預處理（包含圖像對比度拉伸、圖像銳化、自適應平滑濾波、小波去噪等）、分割（先kmeans 聚類，然后基于區域進行圖像分割）和特征提取（canny 邊緣檢測等），分析碎片云特征分布的動態變化過程。

1.1 高斯濾波

首先，需要對碎片云圖像進行高斯濾波，一方面為了去除噪點，另一方面為圖像分割和邊緣檢測做準備。

高斯濾波就是一種線性平滑濾波，實用于消除高斯噪聲，并廣泛應用于圖像處理的減噪過程。 通俗地講，高斯濾波是對整幅圖像展開加權平均的過程，每一個像素點的數值，均應由其本身及鄰域之內的其它像素值經加權平均之后獲得。 高斯濾波的具體實施是：用一個模板（或稱卷積、掩模）掃描圖像中的每一個像素，用模板設定的鄰域內像素的加權平均灰度值去替代模板中心像素點的值。

本質上，高斯模糊就是將（灰度）圖像I 和一個高斯核進行卷積操作：

其中： ?表示卷積操作，Gσ是標準差為σ 的二維高斯核，定義為：

計算平均值時，只需把“中心點”作為原點，其余點根據其在正態曲線上的位置，分配權重，便能夠得到一個加權平均值，而這便是上述與二維高斯核進行卷積的過程。 獲得權重矩陣之后，中心點以及附近n 個點，每一點乘以自身的權重值并將這些數值相加，就是中心點的高斯模糊的數值。 對所有點重復這個過程，就得到了高斯濾波后的圖像。 高斯濾波后的超高速撞擊碎片云圖像如圖2 所示。

圖2 高斯濾波處理后的碎片云圖像Fig.2 Gaussian filtered image of debris cloud

1.2 對比度拉伸

高斯濾波后，對碎片云圖像進行對比度拉伸。對比度拉伸是圖像增強的一種方法，本質上屬于灰度變換操作。 當圖像所有像素的灰度值大部分集中在0－255 中某個很小的區間內時，整個圖像就會顯得很暗淡，對比度不高。 如果通過灰度變換，將灰度值拉伸到整個0－255 的區間，那么其對比度顯然會大幅增強，圖片也會顯得更加明亮清晰。 可以用如下公式將某個像素的灰度值映射到更大的灰度空間

其中： Imin，Imax是原始圖像的最小灰度值和最大灰度值，MIN 和MAX 是要拉伸到的灰度空間的灰度最小值和最大值。

除上述方法，對比度拉伸還可以使用直方圖位移法，即

在每個像素位置的灰度值增加一個偏移量offset。 注意，這個offset 可以是正數，也可以是負數。 若為正數，整體亮度變亮；若為負數，整體亮度變暗。 需要注意的是控制offset 值的大小，不要越界。

碎片云圖像處理結果如圖3 所示。 與處理前的圖像相比，經過對比度拉伸后的圖像更加清晰，且圖像中的無關信息、雜點等都被消除。

圖3 進行對比度拉伸處理后的碎片云圖像Fig.3 Debris cloud image after contrast stretching

1.3 自適應閾值圖像分割

為了減少數據量，同時提取感興趣區域（ROI），即碎片云的主體部分，需要對圖像進行自適應閾值圖像分割。

根據圖像上灰度值的分布，把圖像分為背景和前景兩部分，前景是按照閾值分割出來的部分，背景和前景的分界值即為待求的閾值。 遍歷不同的閾值，計算各個閾值下對應的背景和前景之間的類內方差，當類內方差取到極大值時，這時對應的閾值即為所求的閾值。

對于圖像I（x，y）， 前景（即目標）和背景的分割閾值記作T，屬于前景的像素點數占整幅圖像的比例記為ω0，其平均灰度μ0；背景像素點數占整幅圖像的比例為ω1， 其平均灰度為μ1。 圖像的總平均灰度記為μ，類間方差記為g。 假設圖像的背景較暗，并且圖像的大小為M ×N，圖像中像素的灰度值小于閾值T 的像素個數記作N0，像素灰度大于閾值T 的像素個數記作N1，則有：

將式（9）代入式（10），得到等價公式：

上述即為類間方差的公式表述。 采用遍歷的方法得到使類間方差g 最大的閾值T，即為所求的閾值。

將分割后的圖像二值化后，便得到圖4 的效果。

圖4 進行自適應閾值分割處理后的碎片云圖像Fig.4 Fragmentation cloud image after adaptive threshold segmentation

1.4 Canny 邊緣檢測

為了便于構建碎片云的數學模型，進一步除去無關信息，減少數據量，需要對圖像進行canny 邊緣檢測。

進行圖像邊緣檢測必須滿足2 個條件：一是能有效地抑制噪聲；二必須盡量精確確定邊緣的位置。根據信噪比及定位乘積進行測度，獲得最優化逼近算子。 這便是Canny 邊緣算子。 類似于Marr（LoG）邊緣檢測方法，都屬于先平滑后求導數的方法。

該算法首先使用高斯濾波器平滑圖像，原理如下：

令g（x，y） 為平滑后的圖像，用h（x，y，σ） 對圖像f（x，y） 的平滑可表示為：

其中，?代表卷積，用一階偏導的有限差分來計算梯度的幅值和方向。 已平滑的梯度可以使用2×2一階有限差分近似式來計算x 與y 偏導數的2 個陣列和

幅值和方位角可用直角坐標系到極坐標的坐標轉化公式來計算：

M［x，y］ 反映了圖像的邊緣強度；θ［x，y］ 反應了邊緣的方向。 使M［x，y］ 取得局部最大值的方向角θ［x，y］，就反映了邊緣的方向。

隨后對梯度幅值進行非極大值抑制。 因僅僅得到全局的梯度并不足以確定邊緣，因此為確定邊緣，必須保留局部梯度最大的點，而抑制非極大值。

圖5 邊緣梯度示意圖Fig.5 Edge gradient

圖6 圓周圖和3×3 窗口Fig.6 Circumference and 3?3 windows

利用梯度的方向將梯度角離散為圓周的4 個扇區之一，以便用3×3 的窗口作抑制運算。 4 個扇區的標號為0 到3，對應3×3 鄰域的4 種可能組合。在每一點上，鄰域的中心像素M［x，y］ 與沿著梯度線的2 個像素相比。 如果M［x，y］ 的梯度值小于沿梯度線的2 個相鄰像素梯度值，則令M［x，y］ ＝0。邊緣梯度示意圖如圖5 所示。 圓周圖和3×3 窗口如圖6 所示。

最后用雙閾值算法檢測和連接邊緣，如圖7 所示。 圖像既保留了碎片云的結構屬性，又最大程度地減少了數據量，便于后續的分析計算。

圖7 進行Canny 算法處理后的碎片云圖像Fig.7 Debris cloud image processed by Canny algorithm

2 后板圖像處理

本文實驗采用的是一塊19.8 cm×19.7 cm 后板。 經過超高速撞擊靶板后產生的碎片云撞擊后板而形成的彈坑，用相機拍照得到如圖8 所示的圖像。

為保證圖像中后板彈坑部分和未受撞擊部分顏色灰度值具有明顯的區別，需要對后板圖片進行預處理（包括小波去噪、分割ROI 圖像、增加圖像對比度等）。 將后板彈坑部分和未撞擊部分進行圖像二值化，對比可知彈坑部分與未受撞擊部分顏色具有明顯的區別，如圖9 所示。 為了突出后板損傷部分，同時減小數據量，對后板圖像進行自適應閾值分割。本文實驗自適應閾值為154，處理后圖像如圖10 所示。

圖8 碎片云撞擊后板形成的彈坑圖像Fig.8 The crater image formed by the debris cloud impacting the back plate

圖9 進行自適應閾值分割后的后板圖像Fig.9 Back－plate image after adaptive threshold segmentation

圖10 選擇后板圖像的感興趣區域Fig.10 Select the ROI of theback panel image

3 實驗分析

3.1 碎片云長徑比分析

對于碎片云來說，長徑比與長寬比的概念相同，即經過碎片云內部的最長徑，和與其相垂直的最長徑之比。 此參數常用來表述碎片云形貌，在碎片云的判斷中具有實用價值。 由于碎片云的顆粒較細數量又多，為了更好地觀察和分析碎片云的輪廓。 根據圖7 以及計算該碎片云的長徑比數值，可以準確構建如圖11 所示的數學模型。

圖11 碎片云運動模型Fig.11 Model of debris cloud movement

圖11 中，D1為碎片云的橫向長，D2為碎片云的徑長。 長徑比公式為：

經過測量得到D1＝1 310，D2＝648，所以碎片云徑長比D ＝2.02。

3.2 后板損傷面積及損傷率

后板內部有彈坑，經圖像處理后，彈坑與未撞擊后板部分的顏色不同， 因此對于圖像中的后板部分像素， 只需將該圖像所有的像素個數計算出來即可。 計算公式如下：

其中，M 表示后板像素個數， n 和m 分別表示后板圖像像素的行數和列數。

后板彈坑面積與后板面積的比應等于各自像素個數之比， 因此后板彈坑面積為N?S／M。 其中，S表示后板的面積，h1和h2分別表示圖像中彈坑部分和未撞擊部分的像素值。

對圖10 中的后板彈坑面積以及后板面積的像素值總數進行測量，得到一個1 854×1 858 的灰度矩陣， 且h1＝0，h2＝255 則后板像素個數： M ＝1 854 × 1 858 ＝3 444 732。 同理計算彈坑像素個數：N ＝22 731 而S ＝19.8 cm×19.7 cm＝390.06 cm2， 因此彈坑的面積為：

后板損傷率：

3.3 后板有效損傷面積

結合實際情況，本文認為后板集中損傷部分為實際的有效損傷。 如圖12 所示，則有效損傷面積計算方式如下：后板尺寸為S ＝19.8 cm×19.7 cm＝390.06cm2，經測量圖12 中有效損傷面積的縱深長度（即紅色矩形長度） R1＝2.07 cm，縱深寬度（即紅色矩形寬度） R2＝1.39 cm，則有效損傷面積S′ ＝2.07 cm×1.39 cm＝2.877 3 cm。

圖12 紅色矩形為后板實際有效損傷面積Fig.12 The red rectangle is the actual effective damage area of the rear plate

3.4 有效損傷碎片比例

如果忽略碎片云的飛行角度、速度、運動軌跡等物理特性，則碎片云中的每個碎片將水平飛行撞擊后板。 因此想要計算對后板產生有效損傷的碎片的區域面積，則需要將圖12 中有效損傷面積的縱深長度R1 映射到圖7 中的碎片云圖像，將兩者結合得到圖13，圖13 中紅線所圍矩陣部分即為碎片云中產生有效損傷的部分。 經測量圖13 中碎片云的面積S1＝620 819（像素值），產生有效損傷的碎片云面積S2＝261 811 （像素值），則有效損傷碎片比例rate ≈42.17%。

圖13 有效面積的縱向長映射在碎片云圖像中Fig.13 The longitudinal length of the effective area is mapped in the fragment cloud image

4 結束語

本文利用智能算法分析超高速撞擊碎片云圖像分析的初步探索，同時結合碎片云的物理特性、飛行角度、飛行軌跡等特性，建立了一種全新的航天器損傷估計模型。 這種模型對于指導航天器防護設計，保障航天器可靠的安全在軌運行，提高航天器在空間碎片環境中的生存能力具有重要意義。 首先將圖像進行預處理消除了大量冗余信息與噪聲，大大減少了數據處理量；其次對圖像進行自適應閾值分割和邊緣檢測操作，很好地提取碎片云的輪廓，有利于對碎片云的建模。 根據實驗分析結果表明了本文所提算法的有效性。 這種損傷估計模型的建立不僅考慮了彈丸的物理特征，還考慮了不同時刻碎片云的圖像特征，可以預見，其將更符合實際工程需要。 由于實驗結果還存在一些誤差，尚需對算法等進行進一步調整和優化，還需要更進一步的學習研究。