Kuo ZHANG ,Jianliang HUO ,Shengzhe WANG ,Xiao ZHANG ,Yiting FENG

1School of Automation Engineering,University of Electronic Science and Technology of China,Chengdu 611731,China

2Norla Institute of Technical Physics,Chengdu 610041,China

Abstract:To ensure the safety and reliability of spacecraft during multiple space missions,it is necessary to conduct in-situ nondestructive detection of the spacecraft to judge the damage caused by the hypervelocity impact of micrometeoroids and orbital debris (MMOD).In this paper,we propose an innovative quantitative assessment method based on damage reconstructed image mosaic technology.First,a Gaussian mixture model clustering algorithm is applied to extract images that highlight damage characteristics.Then,a mosaicking scheme based on the ORB feature extraction algorithm and an improved M-estimator SAmple Consensus (MSAC) algorithm with an adaptive threshold selection method is proposed which can create large-scale mosaicked images for damage detection.Eventually,to create the mosaicked images,the damage characteristic regions are segmented and extracted.The location of the damage area is determined and the degree of damage is judged by calculating the centroid position and the perimeter quantitative parameters.The efficiency and applicability of the proposed method are verified by the experimental results.

Key words:Hypervelocity impact;Damage information extraction;Image mosaicking;Damage localization;Quantitative assessment

1 Introduction

With increasing human activities in space,more and more satellites,rockets,and probes are launched into the Earth orbit.Consequently,the risk of hypervelocity impact (HVI) caused by meteoroid/orbital debris (M/OD) continues to rise and it has become a major threat to space activities of spacecraft.This will also have a great impact on spacecraft,including the damage caused by surface craters and embedded impurities (Lamb,2018;Adushkin et al.,2020;Huang et al.,2020).Considering the uncertainty of the occurrence of a spacecraft HVI event,the impact damage has high complexity.Therefore,it is necessary to adopt effective nondestructive testing techniques and to evaluate the potential damage to the spacecraft.

The optical pulsed thermography (OPT)method has been widely used due to its high effi-ciency,high safety,and low loss (Cheng et al.,2018;Kaur and Mulaveesala,2020;Dua et al.,2021).In the OPT detection system,Smurov et al.(2016)and Karnati and Liou (2020) used infrared camera (IR camera) to record the distributions of Joule heat which present different transient thermal responses(TTRs) in the detection area as infrared thermal data.Recently,algorithms for extracting characteristic information from thermal data have been studied,such as independent component analysis(ICA)and principal component analysis(PCA)(Rajic,2002;Khan et al.,2008;Gao et al.,2014;Liang et al.,2016).Huang et al.(2018)introduced a rapid adaptive method of removing redundant data.Yin et al.(2021) proposed an improved method of automatically identifying the defect.Zhang HN et al.(2020) discussed the use of the variable interval to search for thermal feature components.

Many scholars have proposed quantitative methods for assessing damage information in different platforms and solution contexts.Li et al.(2019) characterized the impact damage types using passive thermography and facilitated the identification of impact damage mode and the evaluation of damage degree based on hot spot characteristics.Yin et al.(2019) proposed a novel HVI damage evaluation method based on an active infrared thermal wave image detection technology with multiobjective feature extraction optimization to achieve quantitative evaluation of M/OD HVI damages.

However,none of the above methods can meet the requirements of large-scale detection,and they have strict requirements for detection conditions.Because spacecraft are so large,it is impossible to completely understand all the damage information(such as distribution and extent) due to the limited range of a single shot.To extract damage information and reconstruct it in two-dimensional images which are easier to process,we refer to large-scale remote sensing image processing methods (Surya Kumari and Bhavani,2017;Laraqui et al.,2018;Zhang WP et al.,2018),and apply a mosaicking method to detect a wide area of spacecraft surface.

The method based on feature points is widely used and has better performance because of its high detection speed and low computational complexity.Among the diverse feature-based methods,scale invariant feature transform(SIFT)and speeded-up robust feature (SURF) are two of the most popular feature extraction algorithms.These algorithms can maintain the scaling and rotation invariance of images and are robust to noise,displacements,geometric deformation,and illumination (Lowe,2004;Bay et al.,2008).In conventional mosaicking schemes,for feature matching the distance ratio method and the RANdom SAmple Consensus (RANSAC) or Mestimator SAmple Consensus(MSAC)algorithm are often adopted to find the matched point pairs after screening and to estimate the geometric transformation matrix model parameters to realize the image mosaicking process (Richter et al.,2014;Magri and Fusiello,2017).Traditional schemes rely on the experience of researchers to select the iterative threshold for screening the correct matching point combination,which is inaccurate and not universal.It is of great significance to use the optimization algorithm to calculate the appropriate threshold value to ensure matching point accuracy and save manpower.

Color features directly reflect different damage characteristics,so we can preliminarily extract damage regions by the color image segmentation method.The limitation of the classic RGB color space is that the three colors have significant linear correlation.In contrast,theL?a?b?(Lab for simplification)color space is the most uniform color space independent of equipment because of its uniformity on the perception level.Therefore,the Euclidean distance can be used as the metric criterion in the Lab color space,to further realize color feature segmentation using a classification algorithm (Chandrakanth and Sandhya,2015;Fida et al.,2017).

Based on the clustering idea,in this study a Gaussian mixture model (GMM) is introduced to classify and divide data information according to the TTRs,and to separate images that highlight damage information(Chen et al.,2015;Qiu et al.,2019).The damage characteristics are extracted from the infrared thermal video to reconstruct the damage image.Then the ORB algorithm,which has high feature point detection efficiency,is adopted,and an optimized MSAC feature-matching algorithm based on an adaptive threshold selection method is proposed to obtain an accurate set of matching points to ensure the correctness of the mosaicking results(Rublee et al.,2011).Finally,the damage area is further segmented and extracted,and the location and quantitative calculations are performed on it.The efficiency and applicability of the quantitative assessment scheme are verified by the results.

2 Problem statement and presentation

In this study,we aim at achieving a rapid and effective quantitative assessment of spacecraft damage while meeting the needs of large-scale detection,and propose three main processing modules to deal with the three key issues shown in Fig.1.They are based on damage information extraction,mosaic technology for damage image reconstruction,and quantitative assessment.

Fig.1 Schematic of the key issues

The main innovations and contributions of this paper are as follows:(1) Based on the idea of clustering,we classify and extract damage information from the infrared thermal data,and reconstruct it in a two-dimensional image highlighting the damaged area;(2) We introduce an innovative use of mosaic technology based on damage reconstructed images,which can satisfy the large-scale detection requirement;(3) A new quantitative assessment method is effectively applied to spacecraft damage detection.

3 Damage reconstruction module

3.1 Damage information extraction

The thermal video information is sampled as the image sequence Is(f) with a certain number of frames(f=1,2,...,F1).The frame-by-frame change of each pixel is denoted as a TTR.Then,we vectorize each image frame and combine them into a reconstruction matrixOttrwith the dimension ofF1×F2,whereF2=1,2,...,p×qandp×qis determined by the infrared camera resolution.Figs.2a—2d show the data conversion.Each column of the reconstruction matrixOttris the corresponding TTR of a pixel,soOttrcan be regarded as the TTR data composed of multiple characteristic regions.Based on the TTR characteristics,the typical TTRs with the same patterns are similar and can be classified.

The thermal data can be approximated using a GMM composed of TTRs in various characteristic regions,so we can use the GMM algorithm which can describe the probability distribution of several Gaussian components to classifyOttrinto the category with the highest probability.Then,the mean vectors of each Gaussian model that can most typically reflect the data distribution characteristics of each Gaussian model are selected to constitute the damage characteristic matrixGd,which is shown in Fig.2.The flowchart of theOttrclassification implemented by the GMM algorithm (Algorithm 1) is as follows:

Fig.2 Schematic of the data conversion and clustering:(a) image sequence Is(f);(b) fth frame of Is(f);(c) vec[Is(fth)]T;(d) reconstruction matrix Ottr;(e) damage characteristic matrix construction

As shown in Fig.2e,the GMM algorithm is used to divide the reconstruction matrixOttrinto three clusters.The clustering results areG1,G2,andG3,which are combined into the damage characteristic matrixGd(F1×3).

3.2 Damage image reconstruction

When introducing the pseudo-inverse matrixwe obtain the observation result matrixOr(n×M)ofGdunder the reconstruction matrixOttrconditions bywhich can highlight all kinds of characteristic information.Similarly,matrixOris converted back to the original image size to reconstruct various images of typical characteristics.As shown in Fig.3,matrixOris reconstructed into three types of characteristic images(Or1,Or2,Or3).Each type represents a different detection area location.

The essence of the damage reconstructed image(DRI) is a matrix that contains temperature information.Based on the reconstructed temperature

whereαimis the mixing coefficient and=1(i=1,2,...,n),oooi ∈{Ottr(:,1),Ottr(:,2),···,Ottr(:,n)},andθ={μ1,μ2,...,μM,Σ1,Σ2,...,ΣM}denotes the parametric vector set of mean vectorμand covariance matrixΣ.3:Calculate the posterior probability ofoooicoming from themthGaussian mixture distribution.Use Eq.(A2) to find the posterior probability ofoooicoming from themthGaussian mixture distribution:

4:The logarithmic likelihood functionL(θ,θv) is used to carry out iterative calculation ofθ.It is calculated by Eq.(A3):

Update parameterθv+1using the following formula:

5:The cluster mark of each sampleoooiis determined byGm=and is divided into the corresponding clusters.

6:The mean vectors of eachGmcluster are combined to form the damage characteristic matrix.

information represented by each element in the matrix,a mapping is formed to the color space.In DRIs,the damage characteristic area is reconstructed in red and the background area is blue.

4 Mosaicking module

We propose a mosaicking scheme for DRIs,despite a wide range of detecting requirements,which can realize complete integration of globally detected information on a large scale.First,considering its detection speed and robustness,the ORB algorithm is applied.Then,after the initial rough matching method with Euclidean distance calculation,numerous mismatches,which will seriously affect the precision and effect of mosaicking,need to be removed.To accomplish this task,an optimization algorithm,such as the MSAC algorithm,is used to search for inliers.

However,in the conventional MSAC algorithm,the threshold is set empirically by researchers,which is extremely uncertain and time-consuming.Therefore,we propose an adaptive threshold-setting method to obtain the optimized feature matching results to ensure the accuracy of mosaicking.In this section,a detailed algorithm description is provided and the schematic of DRI mosaicking is shown in Fig.4.

Fig.3 Diagram sketch of image reconstruction (References to color refer to the online version of this figure)

Fig.4 Schematic of the image mosaicking scheme

4.1 ORB feature extraction algorithm for DRIs

The ORB method is improved on the basis of the FAST detector and BRIEF descriptor.Conventional FAST-16 with a circular radius of 16 is adopted here;it is favorable due to the performance boost that it can offer.For a local area of damage in a DRI,the absolute pixel value difference between pointPand the pixel point in the dashed line in Fig.5 is greater than a given threshold;that is,Eq.(1)is satisfied:

whereI(x) represents the grayscale value of a pixel point on the circumference,I(p) represents the grayscale value of target pixelp,andεis a predetermined threshold.

Fig.5 Feature detection vs.feature point detection

An intensity centroid is employed to make FAST robust against orientations.An assumption is made that the intensity of the corner is offset from its center and that it subsequently uses this vector to add an orientation.The momentsMpqof a patch used to compute the centroid are represented as follows:The centroid can be obtained asThen,a vector can be constructed fromO(center of the corner)to the centroid,and the orientation of the patch becomesθ=arctanIn the aforementioned equations,I(x,y)is the grayscale value of point(x,y)and(Cx,Cy)is the centroid.

Consider the illustration of conventional BRIEF operations.Suppose that there is a smoothed image patchp.A binary testτon this patch is represented as follows:

wherep(x) denotes the intensity at a given pointx.Hence,the feature can be written as a vector ofnbinary tests as follows:

Consider any given feature set ofnbinary tests at a particular location,a 2×nmatrix can be represented as follows:

Using the patch orientationθand the corresponding rotation matrixRθ,a steered versionSθofScan be written asSθ=RθS.Hence,the steered BRIEF operator can be represented as follows:

The angle is discretized such that every angle is a multiple of 2π/30 (12?).A lookup table of precomputed brief is constructed.The accurate set of pointsSθwill be used to compute the key point descriptor as long as key-point orientationθis constant across all the directions.

4.2 Adaptive threshold selection method for MSAC

In this subsection,an adaptive threshold selection method is developed based on the MSAC algorithm.The MSAC algorithm model can be expressed as where

Detailed steps of the adaptive threshold selection method are as follows:

Step 1:For all the initial matching point pairs,calculate the cost function.

Step 2:Calculate and record the minimum,maximum,and mean values of the cost at step 1 asEmin,Emax,andEmean,respectively.

Step 3:K,the number of thresholdsT,is considered in the range fromEmintoEmax,where the value ofKdetermines the range of the threshold.The step size of the adjacent threshold is calculated asand in accordance withT(k)=(?×k)+Emin(k=1 :?:K),the value of each threshold is determined.

Step 4:ForKthreshold values,two classes are determined,and for each matching point pair,cost is calculated.According to the selected threshold value and the calculated cost,the matching point pairs are divided into two categories (the correct matching is the first category and the incorrect matching is the second category).

Step 5:To decrease the variance in the first category(correct matches)and to increase the variance in the second category,distance values must be close to each other in the first category.The mean valueEmeanis also taken into account,where the threshold value should be close to the average value.The objective function is implemented as

Finally,the threshold with the lowest value minf(k) is selected as the final threshold according to Eq.(9).The optimal inlier point setNinlierand candidate homographyHinlierare selected,and the infrared reconstruction image mosaicking process is further implemented based on the homographyHinlier.

5 Damage quantitative assessment module

In this section,a two-stage damage information segmentation method is proposed to achieve a complete extraction of the damage areas.The folwchart is shown in Fig.6.

5.1 Segmentation and extraction of damage characteristic regions

The color feature difference of the damage reconstructed image directly reflects the different characteristic areas,so we can apply a clustering segmentation algorithm to segment the regions with damage characteristic color features(orange to red).

Fig.6 Flowchart of the quantitative assessment method (References to color refer to the online version of this figure)

First,the input DRI is converted to the Lab color space,and it is recorded as a color feature object setD(ai,bi) (i=1,2,...,n),wherenis the number of DRI pixels.The goal ofK-means clustering is to divideD(ai,bi) intoKcategories and form the damage feature data setCd={ck|k=1,2,...,K},where the cluster center ofckisok.Define the Euclidean distance measurement color difference criterion as

Based on Eq.(10),there arenkcolor feature object sample points that are divided into categoryck.Take the average value of all sample points in the class subsetckas the class centerok.The iterative update calculation method is as follows:

Calculation is repeated until the color difference takes the minimum value,which means that the sum of the distances from each sample data to its cluster center point takes the minimum value.Ifokin Eq.(11) goes to the minimum,the judgment algorithm is terminated and the damage characteristic data setCdis taken as the final clustering segmentation result.

The cluster segmentation image reflects the different DRI area information,in which the red damage characteristic region is separately extracted,removing the interference caused by the background and the lateral heat diffusion.

Moreover,we propose a binary segmentation extraction algorithm that is based on dual-threshold segmentation processing to ensure the segmentation accuracy and obtain complete segmentation results of damage characteristic regions.Suppose that the gray level set of its pixels isg ∈(0,1),that the number of all pixels in gray leveliis denoted asgi,and that the total number of pixels isThen,the probability that gray leveliappears in the image can be expressed asAfter dual-threshold (T1,T2) binarization segmentation,the pixels in the image are divided into three categories,The probability distribution of various occurrence types is

The average values of various gray levels are as follows:

The expression for the variance between classes is as follows:

From the best threshold combinationthe damage characteristic region is segmented and marked as bright white,which achieves segmentation and extraction of damage information.

5.2 Localization and quantification of damage characteristic regions

By counting the number of pixel points in the connected domain of each damage characteristic region as “1,” we compare the “1” pixels statistically and record the minimum and maximum horizontal and vertical coordinates as the upper-left corner coordinates ul(x,y) and the lower-right corner coordinates lr(x,y),respectively.We set the upper-left corner of the shooting range and the damage segmentation image as the starting point(1,1)of the spatial position coordinates and the pixel coordinates.

Assume that the detection area isM ×Nand that the DRI with resolutionm×nis segmented by the clustering segmentation algorithm and binarization segmentation algorithm.The equivalent proportional relations are expressed asFig.7 shows the damage location.The actual spatial location information of the damage characteristic regions of the test piece is estimated based on the scale conversion of the position information of image pixel points.

Fig.7 Schematic of damage characteristic area localization

The perimeter conversion coefficientλbetween the tested region and the DRI is the ratio of the actual perimeter Pmrof the tested region to the perimeter value Pmp.The calculation formula is as follows:

where Pmp_dand Pmr_drepresent the calculated and actual perimeter values of a damage characteristic region,respectively.

6 Experimental results and discussion

6.1 Experiment of damage reconstruction

In the experiment,for sampleAwith HVI and sampleBwith artificial carbon fiber damage,we collected the infrared thermal data of local areas A1—A4 and B1—B4.We set the mixing coefficient of GMM as 3,divided the data,and recorded the mean vectors of the Gaussian distributions asμ1,μ2andμ3.The mean vector curves of different detection areas reflecting the same damage region were approximately the same and can be clearly distinguished.As shown in Tables 1 and 2,we further applied the reconstruction algorithm to obtain the DRIs of samplesAandB,respectively.

6.2 Mosaicking results of reconstructed damage images

We carried out an experiment of the three methods (SIFT,SURF,and ORB) based on the Matlab2020a platform.By comparing the number of detection feature points and the time spent,the applicability and advantages of the algorithms can be determined.

Fig.8 shows the experimental results of the feature extraction process of A1-A2 mosaicking.The three methods have extracted a considerable number of feature points.The experimental data statistics of the feature extraction algorithm is shown in Table 3.Among them,the feature points detected by the SIFT algorithm were distributed mainly at the edge feature contour of the impact damage and the number of feature points was considerable,but the detection time was too long.Meanwhile,compared with the ORB algorithm,the SURF algorithm had fewer feature points detected per unit time.In summary,for the mosaicking process of damage reconstructed images,the detection speed and efficiency of the ORB algorithm were better than their counterparts of the traditional SIFT and SURF methods.

Fig.8 Experimental result comparison of feature extraction algorithms:(a) SIFT for A1;(b) SIFT for A2;(c) SURF for A1;(d) SURF for A2;(e) ORB for A1;(f) ORB for A2

Table 1 Damage characteristic mean vector curves and the corresponding DRIs of sample A

Table 2 Damage reconstructed images of sample B

The root mean square error (RMSE) was selected to quantify the accuracy of homographyHinlierestimation,and it is defined as follows:

The initial matching set obtained with the number 142 from the mosaicking process of A1-A2 was used as the input,and we found that the value of the adaptive threshold was 13.43.Next,as shown in Table 4,we set the threshold values to 1,5,10,50,100,and 13.43.It can be seen from Table 4 that whenthe threshold was set to 1,the inlier ratio was very small.When the threshold was set from 5 to 50,the ratios and RMSE values were almost the same and when the threshold value was set to 100,mismatches occurred and the RMSE value became large,which means that the mosaicking precision was low.

Table 3 Experimental results of the feature extraction algorithms

Table 4 Comparison of mosaicking results for different thresholds

Finally,the specific mosaicking process of sampleAis shown in Table 5.Similarly,the mosaicking results of sampleBare shown in Fig.9.From the final global mosaicking results of samplesAandB,it can be seen that the mosaicking results were explicit,without deformation or distortion,reflecting the global damage characteristic distribution of the test samples.The experimental results proved the efficiency and applicability of the proposed method.

Fig.9 Global mosaicking results of sample B

Table 5 Specific display of image mosaicking of sample A

6.3 Damage quantitative assessment results

Fig.10 shows the damage characteristic region segmentation extraction results of the global mosaicking based on the clustering segmentation algorithm in the Lab color space.Through the segmentation algorithm,the background noise of the test sample was removed and the lateral thermal diffusion interference caused by thermal radiation was reduced.The damage characteristic images were segmented to carry out further parameter calculation to achieve damage quantitative assessment.

Then,for the results obtained by the clustering segmentation algorithm,the binarized segmentation results of the damage characteristic area were further extracted.Because compared with single threshold processing,multi-threshold processing increases only the computational complexity and will involve more parameter adjustments,to facilitate verification and comparison of advanced methods,we introduced only the single-threshold principle.We listed the following state-of-the-art threshold-based segmentation algorithms as a comparison,namely,Otsu,the minimum error threshold method,onedimensional maximum entropy threshold segmentation method,and two-dimensional maximum entropy threshold segmentation method.

The segmentation results of the four methods are shown in Fig.11,and it can be seen that all the four methods can achieve segmentation and extraction of partial damage characteristic regions.The threshold results are shown after calculation and optimization of each threshold segmentation algorithm and the time spent on image segmentation is also provided.Compared with other threshold segmentation algorithms,the Otsu method had a great advantage in time consumption and the calculated threshold parameters can achieve partial effective segmentation of the damage characteristic areas.Therefore,based on the Otsu threshold segmentation algorithm,we used dual-threshold processing to achieve complete segmentation and extraction of the damage characteristic regions.

Fig.10 Segmentation results of the clustering segmentation algorithm:(a) damage characteristic regions of sample A;(b) damage characteristic regions of sample B

As shown in Fig.12,segmentation results of each damage characteristic region were surrounded by the smallest circumscribed rectangle delimited by the upper-left corner coordinates ul(x,y) and the lower-right corner coordinates lr(x,y).Both samplesAandBwere square plates of 1 m×1 m,and the perimeter was denoted as 4 m.The reconstructed image mosaicking results of samplesAandBshowed that the characteristic parameters of perimeter were calculated as 3296.077 and 3200.316 respectively,and the resolutions were 855×845 and 883×752 respectively.Based on this,the conversion coefficients of perimeter were calculated asλA=1.21 andλB=1.25.Finally,the image damage characteristic regional coordinates and perimeter parameter Pmp_dwere converted into the actual regional spatial position coordinates lr(x′,y′) and ul(x′,y′) and the actual damage perimeter Pmr_d,respectively.A1—A4 and B1—B9 damage characteristic regions were extracted from samplesAandBand the statistics is shown in Tables 6 and 7.

Fig.11 Experimental result comparison of threshold segmentation algorithms:(a) Otsu;(b) minimum error threshold;(c)one-dimensional maximum entropy;(d) two-dimensional maximum entropy

Fig.12 Segmentation marking results of damage characteristic regions of samples AAA (a) and BBB (b)

The actual regional spatial position coordinates correspond to the position information of each damage characteristic region in the test sample and reflect the damage distribution information.The actual damage perimeter value Pmr_dof each damage characteristic region obtained by conversion directly reflects the damage degree,thus realizing quantitative assessment of the damage information.

7 Conclusions

In this study,we proposed a quantitative damage assessment method for spacecraft based on the damage reconstructed image mosaic technology.Based on the GMM clustering algorithm,the damage characteristics were extracted from the infrared thermal video information to reconstruct the image.On this basis,an improved mosaicking scheme was proposed.Based on the global mosaicking results,the damage characteristic region was further segmented and extracted.Damage localization and quantitative calculation of the damage information were realized.The experimental results verified the efficiency and applicability of the quantitative assessment scheme.

Table 6 Quantitative parameter statistics for damage characteristic regions of sample A

Table 7 Quantitative parameter statistics for damage characteristics regions of sample B

Contributors

Kuo ZHANG and Jianliang HUO designed the research.Kuo ZHANG,Shengzhe WANG,and Xiao ZHANG processed the data.Kuo ZHANG and Jianliang HUO drafted the paper.Shengzhe WANG helped organize the paper.Xiao ZHANG and Yiting FENG revised and finalized the paper.

Compliance with ethics guidelines

Kuo ZHANG,Jianliang HUO,Shengzhe WANG,Xiao ZHANG,and Yiting FENG declare that they have no conflict of interest.