CN112233099A - Reusable spacecraft surface impact damage characteristic identification method - Google Patents

Abstract

The invention discloses a reusable spacecraft surface impact damage characteristic identification method, which comprises the following steps: representing a thermal image sequence of the pixel points acquired by the thermal infrared imager by using a three-dimensional matrix; finding a temperature peak value in the three-dimensional matrix M; dividing the temperature of the line of the frame number of the temperature peak point; dividing the temperature of the frame number row of the temperature peak point; extracting transient thermal response in a long time by blocks step by step; automatically classifying the extracted typical transient thermal response by adopting a mean shift algorithm; selecting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on gradient search weight adjustment to form a matrix Y; changing the three-dimensional matrix into a two-dimensional matrix, and performing linear change on the two-dimensional matrix by using a matrix Y to obtain a two-dimensional image matrix R (x, Y); and (5) performing feature extraction on the two-dimensional image matrix R (x, y) by using a one-dimensional Ostu segmentation algorithm. According to the invention, through the uniformly distributed solution set, the difference and the similarity are comprehensively considered, and the defect feature extraction accuracy is improved.

Description

Reusable spacecraft surface impact damage characteristic identification method

Technical Field

The invention belongs to the technical field of spacecraft damage detection and evaluation, and particularly relates to a reusable spacecraft surface impact damage feature identification method.

Background

For the reusable shuttle spacecraft, in order to ensure that the spacecraft cannot generate serious aerodynamic heat damage when entering/reentering the atmosphere, various heat-proof materials need to be paved on the surface of the spacecraft to ensure the operation safety of the spacecraft. However, the physical damage of the surface heat-proof material can be caused by the high-speed impact of the coating-stripping fragments of the spacecraft itself or the ultra-high-speed impact of the space fragments on the earth orbit. Particularly during the performance of a mission in a space debris environment by a spacecraft, the probability of being struck by a tiny space debris is very high. It is worth noting that the ultra-high speed impact can directly cause the surface of the aerospace material to generate various damages such as impact pits, cracks, internal layer cracks and the like, and accumulated damage defects caused by multiple impacts. Meanwhile, secondary projectiles generated by space debris collision are larger in quantity and wider in coverage area, and the secondary projectiles can cause local dense multi-defect damage on the surface of the material as secondary pollution. Due to the very frequent occurrence of space debris impacting the spacecraft and the complexity of ultra-high speed impact damage, various types of spacecraft exposed to the space debris environment must be rapidly and effectively detected and evaluated for various types of impact damage.

The infrared thermal imaging detection technology is that a thermal excitation source is used for actively or passively heating a tested piece, a high-speed high-resolution thermal infrared imager is used for obtaining the distribution difference of the surface temperature field of the tested piece caused by damage defects, and the material defects are detected through the analysis and the processing of an infrared thermal image sequence. The technology comprehensively utilizes the thermal imaging technology and the image sequence processing technology, has the advantages of large single detection area, visual detection result, high detection speed, non-contact, simple use and the like, and is very suitable for the automatic on-orbit spacecraft surface impact damage detection.

Because the thermal image sequence obtained by the thermal infrared imager has huge data volume and strong noise interference, the thermal image sequence needs to be subjected to feature extraction in order to obtain a better detection effect. When processing a thermal image sequence, there are methods based on single-frame image processing and also methods based on image sequence processing. The method based on single-frame image processing only considers the temperature distribution information of the test piece at a certain moment, and cannot reflect the temperature conditions of the test piece at different moments, so that the obtained processing result is incomplete and one-sided. And the infrared image sequence data is analyzed and processed, so that the characteristic information of complex impact damage can be quickly and effectively extracted, and the impact damage degree and damage type of the space debris can be effectively judged.

In the Chinese invention patent application with publication number CN108717069A and name of 'a high-pressure container thermal imaging defect detection method based on line variable step length segmentation' published in 2018, 10, month and 30, for the infrared thermal image sequence with huge data volume, the transient thermal response of similar pixel points is removed through variable step interval search, after effective characteristic information is extracted, clustering dimension reduction processing is needed, the method comprises the steps of dividing pixel points into different damage sets according to transient thermal response temperature characteristics of the pixel points, selecting representative thermal response temperature points capable of representing characteristics of the damage sets from the damage sets according to differences among different types of temperature points and similarities among the same type of temperature points, forming a two-dimensional matrix by temperature data of the representative thermal response temperature points of the types, extracting defect characteristics through linear change, and finally displaying defect information in a visual mode.

The clustering algorithm used in the prior infrared nondestructive testing needs to manually set the clustering number, and the algorithm cannot be used for automatically dividing the defect types by depending on the temperature characteristics of pixel points in the space.

The number of clusters represents the division condition of the defect types, so the selection of the number of clusters can influence the clustering effect to a greater extent. In the actual detection of the defects, the laser lamp is used for irradiating the surface of the material, the structural information of the surface and the subsurface of the material is obtained according to the temperature field change information recorded by the thermal infrared imager, the defect condition in the detection test piece cannot be observed, and only the temperature information of the infrared thermal image sequence with huge data quantity is acquired, so that the actual defect type number L of the detection test piece_rIt is not predictable in advance and therefore it is not predictable in advance that the classification of defects into classes is most appropriate. Initially, the number of defect types L 'that may occur is set in advance empirically from the practical point of view of the material being inspected, but L' and L cannot be guaranteed_rKeeping consistent, so that the classification result is not good; in order to solve the problem that the number of classified defect types is inconsistent with the actual number of defect types, the cluster number L is set in consideration of all possible defect types_max,L_max＞L_rClustering, and determining the actual defect type number L according to the transient thermal response of representative temperature points selected from each damage set and the difference of the temperature points_rHowever, when judging the number of defect types, L is obtained by selecting different judgment difference modes_rThe result of (2) is also different, and L is reused_rAnd re-clustering, wherein the obtained results including the clustering centers are different, and the real clustering result cannot be obtained. In order to solve the problem, the invention automatically sets the clustering number according to the difference of the data by using an algorithm, and can obtain an accurate clustering result for the high-dimensional infrared thermal image sequence data with complex distribution acquired by an experiment, thereby ensuring that the next step is carried outThe accuracy of analyzing and mining the data is improved, the effect of selecting the representative temperature point thermal response by utilizing the multi-objective optimization algorithm is better, and the selected result is more representative.

The clustering center only considers the local relation in the iteration process, so the selection of the clustering center also influences the optimization of the objective function and influences the convergence speed of the algorithm and the clustering result. Under high-speed striking, the material surface and inside can produce the defect of multiple different grade type, from the experiment of doing before, because the defect type is more, in local space, the temperature variation of pixel is complicated various and be the irregular distribution, when iterating current clustering center, all the other categorised defects probably have a certain or more pixel very near with clustering center for the objective function reduces rapidly, forms local minimum, causes the objective function to become the multimodal function that has a plurality of extreme points. If the selected initial clustering center is located near a certain minimum point in the local space, when iterative updating is carried out, a minimum value point may be selected, the objective function derivative corresponding to the minimum value point is 0, so that the membership matrix with the derivative of 0 and the clustering center are both current values, at this time, the clustering result is not updated, the objective function representing the connection among all damage features cannot be continuously iterated to find the optimal classification mode, and the algorithm converges to the local minimum value. And the difference of the clustering centers is not large, the difference of defect damage sets obtained by corresponding classification is not large, so that the defect classification is not clear, the defect region can not be well distinguished from other regions, the classification effect is not ideal, and the accuracy of selecting each damage set to represent a thermal response temperature point by establishing a multi-target problem by using the clustering centers is influenced. In order to solve the problem, the size of a search area is determined by setting a bandwidth, the density of pixel points in the space is estimated by using a kernel function, the area with the densest pixel points in the space is searched in an iteration process, a clustering center drifts towards a local density maximum value point along the direction of increasing the density of the pixel points, the clustering center gradually approaches to an optimal position by continuously moving, and the target function is prevented from being trapped in a local extremum.

In the Chinese patent application of the multi-objective optimization infrared thermal image defect feature extraction method based on uniform evolution, when the representative temperature points of each damage set are selected by using a multi-objective optimization algorithm to represent the features of the damage set, the optimal solution is searched by uniformly distributed direction vectors corresponding to uniformly generated weight vectors. However, the uniformly distributed search direction vectors can only improve the diversity of solution set distribution as much as possible, the front surface of the PF is unknown in the actual solving process, the PF is possibly nonuniform for a complex high-dimensional multi-target optimization problem, and in a nonuniform area, the subproblems corresponding to the part of the direction vectors do not have Pareto optimal solutions, so that the diversity of the solution set solved by the algorithm is reduced finally. In order to solve the problem, the invention designs a redistribution strategy for adjusting the direction vector based on gradient search, the direction vector is adjusted based on the gradient search according to the distribution sparse and dense degree of the population after each iteration, then the weight vector is distributed for the subproblems according to the direction vector, and a uniformly distributed Pareto optimal solution set is obtained to maintain the population diversity, so that the selected representative thermal response temperature point is most representative. However, when the weight vectors are redistributed, the distribution situation of the real solution needs to be known in advance, different weight vectors are distributed according to different solution set density degrees, and the adjustment of the weight vectors can influence the distribution situation of the solution sets in the space and the dominance relation between objective function solutions, so that when the weight vectors are adjusted, if the adjustment is not enough, the algorithm cannot achieve the expected effect, and the ideal Pareto optimal solution set with uniform distribution cannot be obtained; too much adjustment may result in the generated solution set having no solution to make the resolved subproblem approach a stable minimum value, thereby affecting the convergence of the algorithm.

Disclosure of Invention

An object of the present invention is to solve at least the above problems and/or disadvantages and to provide at least the advantages described hereinafter.

To achieve these objects and other advantages in accordance with the purpose of the invention, there is provided a reusable spacecraft surface impact damage feature identification method, comprising the steps of:

step one, representing a thermal image sequence of a pixel point on the surface of a test piece acquired by an infrared thermal imager by using a three-dimensional matrix M, wherein M (i, j, T) represents a pixel in the ith row and the jth column, the third dimension T represents acquisition time, and T is 1,2, … and T;

step two, finding a temperature peak point M (i) in the three-dimensional matrix M_max,j_max,t_max) Wherein i_maxIndicates the number of rows where the temperature peak point is located, j_maxIndicates the number of columns where the temperature peaks are located, t_maxRepresenting the frame number of the temperature peak point;

step three, dividing the temperature of the line of the frame number of the temperature peak point, dividing the line of the frame number of the temperature peak point into R data blocks, and recording the line step length;

step four, carrying out temperature division on the row of the frame number of the temperature peak point, dividing the row of the frame number of the temperature peak point into Q data blocks, and recording the row step length;

step five, extracting the transient thermal response in the data block obtained in the step three and the step four in a blocking and step-length manner, and finishing the extraction of the typical transient thermal response;

step six, adopting a mean shift algorithm to automatically classify the typical transient thermal responses extracted in the step five to obtain the number of clustered damage sets and the category of each transient thermal response;

selecting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on gradient search weight adjustment to form a matrix Y;

step eight, changing the three-dimensional matrix M into a two-dimensional matrix, and performing linear change on the two-dimensional matrix by using the matrix Y to obtain a two-dimensional image matrix R (x, Y);

ninthly, extracting features of the two-dimensional image matrix R (x, y) by using an improved one-dimensional Ostu segmentation algorithm, and dividing the image into targets C by setting a threshold k according to the gray features of the image₀And background C₁And obtaining the defect characteristics of the test piece.

Preferably, in the third step, the temperature of the line of the frame number where the temperature peak point is located is divided, and the specific method for dividing the line of the frame number where the temperature peak point is located into R data blocks includes: setting a threshold value

Calculating the temperature peak point M (i)_max,j_max,t_max) To the nearest temperature point M^c(i_max,j_c,t_max) Is calculated by the correlation coefficient PC_reWherein j is_cThe number of columns where the nearest temperature point is located, if

Then j is_c＝j_c+1, calculating the temperature peak M (i) from near to far_max,j_max,t_max) Correlation coefficient PC with the remaining temperature points_reUp to

When it is time, stop calculating, remember

Number of established pixels as line step length CL^rAccording to CL^rAnd dividing the temperature into lines, and dividing the lines of the frame number of the temperature peak point into R data blocks.

Preferably, the step four of dividing the temperature of the column of the frame number where the temperature peak point is located into the temperature, and the specific method of dividing the column of the frame number where the temperature peak point is located into the Q data blocks includes: setting a threshold value

Calculating the temperature peak point M (i)_max,j_max,t_max) To the nearest temperature point M^c(i_c,j_max,t_max) Is calculated by the correlation coefficient PC_reWherein i_cThe number of rows where the nearest temperature point is located, if

Then i_c＝i_c+1 calculating the temperature peak M (i) from near to far_max,j_max,t_max) Correlation coefficient PC with the remaining temperature points_reUp to

When it is time, stop calculating, remember

The number of established pixels is the column step length RL^qAccording to RL^qAnd dividing the temperature into columns, and dividing the columns of the frame number of the temperature peak point into Q data blocks.

Preferably, the step of extracting the transient thermal response in five blocks and steps comprises the following specific steps:

step S51, step three, and step four divide the three-dimensional matrix M into N data blocks, where N is R · Q, and denote the nth data block as MⁿAnd setting a threshold CC, where N is 1,2, …, N;

step S52, the temperature peak point M (i)_max,j_max,t_max) Corresponding transient thermal response M (i)_max,j_maxT), T is 1,2, …, T is stored in the set X (: 1), and the transient thermal response M of the mth row and nth column of each data block is calculatedⁿ(m, n, T), T ═ 1,2, …, correlation PC between T and X (: 1)_re(ii) a If PC_re＜CC，Mⁿ(m, n, T), T ═ 1,2, …, T is regarded as a new feature and stored in X (: z); otherwise, let m be m + RL^qContinuously calculating the correlation degree of the next temperature point and X (: 1); if m > k, the column calculation is performed: let m be m-k, n be n + CL^rAnd k is the row number of the nth data block, and if n is more than l, wherein l is the column number of the nth data block, the typical transient thermal response extraction is finished.

Preferably, wherein the step six employs a mean shift algorithm, the method for automatically classifying the typical transient thermal response extracted in the step five includes: defaulting that all data belong to one class, namely the cluster number L is 1, calculating by an iterative algorithm, dividing the transient thermal response of the pixel points with larger difference into the other class, updating the cluster number L to be L +1, after all the transient thermal responses in a transient thermal response set X (: z) are subjected to iterative calculation, finally obtaining the cluster damage set number L and the class to which each transient thermal response belongs, and the specific steps comprise:

step S61, setting a bandwidth range r, an iteration termination condition epsilon and a threshold omega;

step S62, randomly selecting a transient thermal response in the set X (: z) as the initial cluster Center point Center, setting the cluster number L to 1, the cycle parameter T to 1, i' to 1, the cycle end condition T to T, keeping the initial cluster Center when T is 1, and recording the initial cluster Center when T is 1₁Center＝Center；

Step S63, find out all transient thermal response sets distributed in the high-dimensional sphere area with Center as the Center and radius r, namely, satisfy | Center-X^mAll transient thermal responses X of | < r^mWherein X is^mE.g., X (: z), and the transient thermal responses are considered to belong to the same cluster NDs;

step S64, taking the Center as the Center, calculates the vector of each transient thermal response in the Center to cluster set NDs, and adds these vectors to obtain the shifted mean:

wherein, omega is a high-dimensional sphere area with the Center as the Center radius r; x^mSatisfying [ Center-X ] for transient thermal response distributed in omega region^mL < r; k is transient thermal response X distributed in an omega region^mThe number of (2); g (x) is a binuclear function

The derivation is negative;

step S65, the Center updates along the direction of the shift mean, and the movement distance is | | shift (Center) |:

^tthe Center represents the updated cluster Center;

step S66, if | | | shift (Center) | | | > epsilon, letCenter＝^tThe Center, t is t +1, the step S63 is returned to continue the iterative computation until | | shift (Center) | < epsilon, and the previous iteration times t are found out to be distributed so as to^tAll transient thermal response sets in the high-dimensional sphere area with Center as the Center and radius r, namely satisfying | survival^tCenter-X^b||＜r，X^bAll transient thermal responses X for e X (: z)^bThese transient thermal responses are considered to belong to the same cluster^tNDs；

Step S67, calculating the updated cluster center^tCenter and other already existing cluster centers_i'Distance of Center, if any^tCenter-_i'If Center | < ω is true, then consider it to be^tCenter and the nearest one to it_kThe Center belongs to the same class, wherein

Go to step S68; otherwise go to step S69;

step S68, collecting^tTransient thermal response in NDs is ascribed to_kCenter, k is 1,2, …, and L is a cluster of the cluster Center_kNDs, will be assembled^tTransient thermal response markers in NDs are accessed and clustered_kThe number of occurrences in NDs is added to 1 and calculated^tCenter and_kweighted average of Center^t'Center，_kNeutral separation of NDs^t'Recent transient thermal response of Center^t″The Center is the cluster Center of the current cluster, i.e. the current cluster Center_kCenter＝^t″Center, the number of clusters L equals L, let Center equal to L^t″The Center goes to step S610;

step S69, collecting^tTransient thermal response in NDs is ascribed to_kCenter, k is 1,2, …, and L is a cluster of the cluster Center_kNDs, will be assembled^tTransient thermal response markers in NDs are accessed and clustered_kThe number of occurrences in NDs is added to 1 and calculated^tCenter and_kweighted average of Center^t'Center，_kNeutral separation of NDs^t'Recent transient thermal response of Center^t″Center is the cluster of the current clusterCentres, i.e. current cluster centres_kCenter＝^t″Center, the number of clusters L equals L, let Center equal to L^t″The Center goes to step S610;

and step S610, making T equal to T +1, returning to step S63, ending the algorithm and clustering, and outputting the clustering number L and the clustering center of each cluster when the loop ending condition T equal to T is met and the number of accessed transient thermal responses marked in X (: z) reaches 90 percent_i'Center, i ═ 1,2, …, L; the transient thermal response marked as visited in the set X (: z) goes to step 611 to classify the category to which it belongs; the transient thermal response which is not marked as visited is transferred to step S612 to divide the category to which the transient thermal response belongs; otherwise, returning to step S63 to continue the iterative computation until the number of transient thermal responses marked as visited in X (: z) reaches 90%;

step S611, performing category classification on the accessed transient thermal response: according to the occurrence frequency of each transient thermal response in each cluster, taking the class with the largest occurrence frequency as the class to which the current transient thermal response belongs, and dividing the current transient thermal response into the classes; when a certain transient thermal response X occurs^p，X^pWhen the times of appearance of epsilon X (: z) in several clustering clusters are the same, the transient thermal response set divided into the class is determined from each clustering cluster_i'NDs, i ═ 1,2, …, L found free from X^pK transient thermal responses with Euclidean distance being nearest are calculated to obtain X^pTo_i'Crowdedness distance of k nearest transient thermal responses in NDs:

wherein,

represents X^pEuclidean distance to the j, j ═ 1,2, …, k nearest transient thermal responses in class i'; mixing X^pIf the calculated congestion degree distance is the largest, the step goes to step S613 to output the transient thermal response division condition;

step S612, classifying the transient thermal response which is not marked as accessed: calculating transient thermal response X not marked as visited^q，X^qBelongs to X (: z) to each category clustering center_i'European distance from Center: categorizing transient thermal responses that are not marked as visited: calculating transient thermal response X not marked as visited^q，X^qBelongs to X (: z) to each category clustering center_i'European distance from Center:

wherein,_i'center, i ═ 1,2, …, L is the cluster Center for each cluster. Calculating the class with the nearest Euclidean distance as the current transient thermal response X^qThe step S613 is switched to output the transient thermal response classification condition;

step S613, outputting the transient thermal response division condition_i'NDs, i ═ 1,2, …, L, and the cluster centers for each class_i'Center，_i'NDs represents the set of transient thermal responses of class i'.

Preferably, the seventh step adopts a decomposition multi-objective optimization algorithm based on gradient search weight adjustment to select a representative of each type of transient thermal response, and the formation of the matrix Y includes the following steps:

in step S71, when the i '(i' ═ 1., L) th class transient response is represented, a multi-objective function is defined:

minimizeF(_i'X)＝(f₁(_i'X),...,f_L(_i'X))^T

wherein f is₁(_i'X) a transient response selected for the i' th class transient response_i'The intra-class Euclidean distance of X is expressed as:

f_l(_i'X),l＝2,3,...,La transient response selected for the i' th class transient response_i'The Euclidean distance between L-1 classes of X is calculated according to the calculated Euclidean distance between L-1 classes

The components are renumbered and the components are,

expressed as:

_i'x_hfor transient response_i'The pixel value of X at the h-th time, i.e. the temperature value,_i'Center_hthe pixel value of the ith' type transient response cluster center at the h-th moment, namely the temperature value,_j'Center_hthe pixel value of the j' th class transient response cluster center at the h moment is the temperature value;

step S72, Multi-objective evolution Algorithm Based on Decomposition, MOEA/D, step S71 to give Multi-objective function, and select the representative of the i' th class transient response_i'REP, i' belongs to (1, 2.. gtorel), and the specific steps comprise:

step S721, initializing the population size to be N, and randomly generating N weight vectors uniformly distributed by a Tchebycheff decomposition algorithm

Wherein

Corresponding direction vector

Comprises the following steps:

wherein,

step S722, decomposing the L multi-target questions into N sub-questions by utilizing a Tchebycheff decomposition algorithm, wherein each sub-question is as follows:

wherein,

_i'z^*for the reference point corresponding to the multi-objective function of the ith' class,

is a function f_l(_i'X) a corresponding reference point;

step S723, initializing the evolution time h to 0, to obtain_i'Randomly generating N values in Ω

_i'Xⁿ(0) Indicating the value of the nth individual in the ith' class at initialization,_i'omega is a value range determined by transient thermal response of the ith class T time dimension;

step S724, initializing Pareto leading edge solution set_i'Initializing reference point for EP (0) < phi >

Wherein,

step S725, initializing neighborhood vectors, and solving the nth weight vector of the ith class

The most recent TT weight vectors;

step S726, initializing position information and speed information of the particle swarm according to the independent variable space dimension;

step S727, setting the maximum number of iterations P, starting iteration when h is 1, terminating iteration when h is P, setting a threshold value δ and a parameter r, and applying each weight vector to each weight vector

And performing the following updating operation according to the iteration number h:

step A, updating individuals, and specifically comprising the following steps:

step A1, setting the scope of evolution: setting a threshold value delta according to a random number rand epsilon (0,1) generated randomly, and determining a neighborhood range P_i'(n)：

Wherein, B_i'(n) is a weight vector index set;

step A2, calculating the particle fitness according to a fitness algorithm, and dividing the initial population into two subgroups which are a non-dominant subset P _ set and a dominant subset NP _ set respectively;

step A3, calculating the adaptive value of each particle, and for each particle, matching the adaptive value with the best position p experienced^bestIf the adaptive value is better, the adaptive value is taken as the current best position; for each particle, its adaptation value is compared with the global experienced best position g^bestComparing the adaptive values, and if the adaptive values are better, taking the adaptive values as the current global best position;

step a4, updating its position information and velocity information by means of the individual extremum and the global extremum for each particle of the dominant subset NP _ set:

v_ji(t+1)＝v_id(t)+c₁·r₁·(P_ji(t)-x_ji(t))+c₂·r₂·(NP_ji(t)-x_ji(t))

x_ji(t+1)＝x_ji(t)+v_ji(t+1)

wherein, c₁And c₂Is a normal number, called acceleration factor, c₁Adjusting the step length of the particles flying to the direction of the best position of the particles; c. C₂Adjusting the step length of the flight of the particles to the global best position; r is₁And r₂Is [0,1 ]]A random number in between; p_ji(t) are random individuals within the non-dominated solution set P _ set; NP_ji(t) are the particles within the dominance solution set NP _ set;

step B, updating the weight vector and the direction vector, and specifically comprising the following steps:

step B1, on the h-th iteration, for all individuals in class i

According to the preference, considering the similarity of the objects, the i' th class population of the h iteration is obtained according to the following formula

Reference points:

wherein,

is an individual

The Euclidean distance in class from the current i' th class center is used for measuring similarity, sum represents the sum of Euclidean distances in class, and the sum which enables sum to be minimum is selected

The individual is taken as a reference point, and the point is recorded

The individual is composed of a set of

Step B2, for each individual reference point

Find out in

High dimensional sphere region omega with radius r as center_hAll individuals in the table, calculate Ω_hStandard deviation of all individuals within

Wherein,

is distributed at omega_hAverage of all individuals within;

represents a distribution in Ω_hK individuals within;

step B3, calculating the standard deviation satisfying

The number of reference point individuals satisfying the condition is recorded as s, and the set of the s reference point individuals is recorded as

To the collection

Each individual reference point in

The following update operations are performed:

step B31, calculating basis weight vector:

and converted into the corresponding direction vector:

wherein,

for the reference point corresponding to the multi-objective function of the ith' class,

is a function f corresponding to a subproblem decomposed by the Chebyshev algorithm_l(_i'X) of the reference points of the two,

step B32, finding the reference point in the generated population

The most distant European individuals

Wherein,

step B33, mixing the individuals

Corresponding weight vector

Conversion to the corresponding directional vector:

step B34, calculating the generated direction vector:

wherein,

step B35, utilizing reference point

And

generate a new solution

Step B36, using the new solution

Replacement of

Direction vector

Replacement of

Corresponding weight vector

Replacement of

Wherein,

step 728, loop judgment: if n is>N, h is h +1, i.e. the next evolution is carried out; if the termination condition is not met, repeating the step S727, otherwise, obtaining a final Pareto leading edge solution set of the ith temperature transient thermal response_i'EP(h)；

Step S729, from_i'Selecting a temperature transient thermal response optimization solution in EP (h)_i'REP, representative of class i' transient responses_i'REP，i'∈(1,2,...,L)；

Step S73, the L-type transient response representatives are arranged in columns, and one column is the temperature values, which are the pixel values at T moments, to form a T × L matrix Y.

Preferably, the step eight of changing the three-dimensional matrix M into a two-dimensional matrix, and the specific method of obtaining the two-dimensional image matrix R (x, Y) by linearly changing the three-dimensional matrix M with the matrix Y includes: starting each frame in the three-dimensional matrix M from a first column, connecting a next column to the tail of a previous column to form a new column, obtaining T-column pixel values corresponding to the T frames, then sequentially placing the T-column pixel values according to time sequence to form a two-dimensional matrix O, and performing linear transformation on the two-dimensional matrix O by using a matrix Y, namely:

a two-dimensional image matrix R is obtained, wherein,

is a pseudo-inverse of the matrix Y, O^TA transpose of the two-dimensional matrix O.

Preferably, the step nine utilizes a modified one-dimensional Ostu segmentation algorithm to perform feature extraction on the two-dimensional image matrix R (x, y),dividing the image into targets C by setting a threshold k according to the gray scale characteristics of the image₀And background C₁The method comprises the following two steps:

step S91, calculating the gray probability distribution P of the image_t：

Where, t is 1,2, …, L is the number of gray levels of the image, N is the total number of pixels in the image, and N is the total number of pixels in the image_tThe number of pixel points with the gray level of t in the image is counted;

step S92, calculating the target class C₀And background class C₁Probability of occurrence of each omega₀And ω₁And their corresponding mean values of gray levels mu₀And mu₁Thereby obtaining the spacing between the two types

And the degree of dispersion d of each class is determined₀And d₁：

Step S93, calculating a threshold value selection function H (t);

step S94, extracting the maximum value in H (t), wherein the optimal threshold k is the corresponding gray level t;

and step S95, binarizing the image according to k to realize segmentation of the gray level image, thereby completing identification of the impact damage characteristics of the surface of the test piece.

The invention at least comprises the following beneficial effects:

1. the clustering number of the clustering algorithm adopted by the invention is not required to be known in advance, and the accuracy rate of the clustering result is higher; the clustering center is not based on the initial assumption, the clustering division result is relatively stable, and the defect classification unobvious defect classification caused by the fact that the clustering center is easy to fall into a local extreme value in the iteration process in the traditional clustering method is overcome.

2. The multi-objective optimization method based on the gradient search weight adjustment not only realizes the consideration of difference and comprehensiveness, but also considers the condition that the spatial solution set is not uniformly distributed in the actual solving process by using the multi-objective optimization algorithm, can effectively extract the most representative pixel points of each damage set, and realizes the effective extraction of defect characteristics and the accurate depiction of the defect outline.

3. The method adopts a weight adjustment strategy based on gradient search, can simply and effectively obtain an optimized weight vector set according to the preference of a decision maker, and is also suitable for solving the problem of high-dimensional complex multi-target optimization.

Additional advantages, objects, and features of the invention will be set forth in part in the description which follows and in part will become apparent to those having ordinary skill in the art upon examination of the following or may be learned from practice of the invention.

Description of the drawings:

FIG. 1 is a flowchart of one embodiment of a method for multi-objective optimized infrared thermal image defect feature extraction based on gradient search weight adjustment according to the present invention;

FIG. 2 is a block-wise selection exemplary transient thermal response flow diagram employed by the present invention;

FIG. 3 is a flow chart of a mean-shift clustering algorithm employed by the present invention to classify selected transient thermal responses;

FIG. 4 is a graph of results of a mean-shift clustering algorithm after classification of selected transient thermal responses;

FIG. 5 is a multi-target Pareto front surface of the temperature points of the material itself;

FIG. 6 is a multi-target Pareto front at defect 1 temperature points;

FIG. 7 is a multi-target Pareto front at defect 2 temperature points;

FIG. 8 is a graph of transient thermal response of the temperature point of the material itself;

FIG. 9 is a graph of transient thermal response at defect 1 temperature point;

FIG. 10 is a graph of transient thermal response at defect 2 temperature points;

FIG. 11 is a graph of transient thermal response for corresponding temperature points of the material itself, selected based on differences;

FIG. 12 is a graph of transient thermal response for corresponding defect 1 temperature points selected based on variability;

FIG. 13 is a graph of transient thermal response for corresponding defect 2 temperature points selected based on variability;

FIG. 14 is a graph of transient thermal response for a corresponding temperature point of the material itself, selected in accordance with the present invention;

FIG. 15 is a graph of transient thermal response for the corresponding defect 1 temperature point selected based on the present invention;

FIG. 16 is a graph of transient thermal response for corresponding defect 2 temperature points selected based on the present invention;

fig. 17 is a defect feature map extracted based on the present invention.

The specific implementation mode is as follows:

the present invention is further described in detail below with reference to the attached drawings so that those skilled in the art can implement the invention by referring to the description text.

It will be understood that terms such as "having," "including," and "comprising," as used herein, do not preclude the presence or addition of one or more other elements or groups thereof.

As shown in fig. 1-3: the invention discloses a reusable spacecraft surface impact damage characteristic identification method, which comprises the following steps:

step one, representing a thermal image sequence of a pixel point on the surface of a test piece acquired by an infrared thermal imager by using a three-dimensional matrix M, wherein M (i, j, T) represents a pixel in the ith row and the jth column, the third dimension T represents acquisition time, and T is 1,2, … and T;

step two, finding a temperature peak point M (i) in the three-dimensional matrix M_max,j_max,t_max) Wherein i_maxIndicates the number of rows where the temperature peak point is located, j_maxIndicates the number of columns where the temperature peaks are located, t_maxRepresenting the frame number of the temperature peak point;

step three, dividing the temperature of the line of the frame number of the temperature peak point, dividing the line of the frame number of the temperature peak point into R data blocks, and recording the line step length;

step four, carrying out temperature division on the row of the frame number of the temperature peak point, dividing the row of the frame number of the temperature peak point into Q data blocks, and recording the row step length;

step five, extracting the transient thermal response in the data block obtained in the step three and the step four in a blocking and step-length manner, and finishing the extraction of the typical transient thermal response;

step six, adopting a mean shift algorithm to automatically classify the typical transient thermal responses extracted in the step five to obtain the number of clustered damage sets and the category of each transient thermal response;

selecting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on gradient search weight adjustment to form a matrix Y;

step eight, changing the three-dimensional matrix M into a two-dimensional matrix, and performing linear change on the two-dimensional matrix by using the matrix Y to obtain a two-dimensional image matrix R (x, Y);

ninthly, extracting features of the two-dimensional image matrix R (x, y) by using an improved one-dimensional Ostu segmentation algorithm, and dividing the image into targets C by setting a threshold k according to the gray features of the image₀And background C₁And obtaining the defect characteristics of the test piece.

The invention relates to a method for extracting defect characteristics of multi-target optimization infrared thermal images based on gradient search weight adjustment, which comprises the steps of selecting transient thermal responses of pixel points from thermal image sequence conversion step length, classifying by adopting a mean shift clustering algorithm to obtain the category of the transient thermal response of each pixel point, considering the pixel value (temperature value) similarity of each category pixel point and the same category pixel point, considering the difference between the pixel point (temperature point) and different category pixel points (temperature points), constructing a corresponding multi-target function, obtaining a dimension reduction result of the thermal image sequence by using a method for extracting the defect characteristics of the multi-target optimization infrared thermal images based on the gradient search weight adjustment, and finally extracting the defect characteristics of the infrared thermal images by using an improved one-dimensional Ostu threshold segmentation algorithm. By the uniform evolution direction of the solution, the difference and the similarity are comprehensively considered, the accurate selection of the representative pixel point (temperature point) is realized, and the accuracy of defect feature extraction is ensured.

In the above technical solution, the third step of dividing the temperature of the line of the frame number where the temperature peak point is located, and the specific method of dividing the line of the frame number where the temperature peak point is located into R data blocks includes: setting a threshold value

Calculating the temperature peak point M (i)_max,j_max,t_max) To the nearest temperature point M^c(i_max,j_c,t_max) Is calculated by the correlation coefficient PC_reWherein j is_cThe number of columns where the nearest temperature point is located, if

Then j is_c＝j_c+1, calculating the temperature peak M (i) from near to far_max,j_max,t_max) Correlation coefficient PC with the remaining temperature points_reUp to

When it is time, stop calculating, remember

Number of established pixels as line step length CL^rAccording to CL^rAnd dividing the temperature into lines, and dividing the lines of the frame number of the temperature peak point into R data blocks.

In the above technical solution, the fourth step of performing temperature division on the column of the frame number where the temperature peak point is located, and the specific method of dividing the column of the frame number where the temperature peak point is located into Q data blocks includes: setting a threshold value

Calculating the temperature peak point M (i)_max,j_max,t_max) To the nearest temperature point M^c(i_c,j_max,t_max) Is calculated by the correlation coefficient PC_reWherein i_cThe number of rows where the nearest temperature point is located, if

Then i_c＝i_c+1 calculating the temperature peak M (i) from near to far_max,j_max,t_max) Correlation coefficient PC with the remaining temperature points_reUp to

When it is time, stop calculating, remember

The number of established pixels is the column step length RL^qAccording to RL^qAnd dividing the temperature into columns, and dividing the columns of the frame number of the temperature peak point into Q data blocks.

In the above technical solution, the step of extracting the transient thermal response in five blocks and in steps includes:

step S51, step three, and step four divide the three-dimensional matrix M into N data blocks, where N is R · Q, and denote the nth data block as MⁿAnd setting a threshold CC, where N is 1,2, …, N;

step S52, the temperature peak point M (i)_max,j_max,t_max) Corresponding transient thermal response M (i)_max,j_maxT), T1, 2, …, T being stored in the set X (: 1), countCalculating the transient thermal response M of the mth row and the nth column of each data blockⁿ(m, n, T), T ═ 1,2, …, correlation PC between T and X (: 1)_re(ii) a If PC_re＜CC，Mⁿ(m, n, T), T ═ 1,2, …, T is regarded as a new feature and stored in X (: z); otherwise, let m be m + RL^qContinuously calculating the correlation degree of the next temperature point and X (: 1); if m > k, the column calculation is performed: let m be m-k, n be n + CL^rAnd k is the row number of the nth data block, and if n is more than l, wherein l is the column number of the nth data block, the typical transient thermal response extraction is finished.

In the above technical solution, the method for automatically classifying the typical transient thermal response extracted in the step five by using a mean shift algorithm in the step six includes: defaulting that all data belong to one class, namely the cluster number L is 1, calculating by an iterative algorithm, dividing the transient thermal response of the pixel points with larger difference into the other class, updating the cluster number L to be L +1, after all the transient thermal responses in a transient thermal response set X (: z) are subjected to iterative calculation, finally obtaining the cluster damage set number L and the class to which each transient thermal response belongs, and the specific steps comprise:

step S61, setting a bandwidth range r, an iteration termination condition epsilon and a threshold omega;

step S62, randomly selecting a transient thermal response in the set X (: z) as the initial cluster Center point Center, setting the cluster number L to 1, the cycle parameter T to 1, i' to 1, the cycle end condition T to T, keeping the initial cluster Center when T is 1, and recording the initial cluster Center when T is 1₁Center＝Center；

Step S63, find out all transient thermal response sets distributed in the high-dimensional sphere area with Center as the Center and radius r, namely, satisfy | Center-X^mAll transient thermal responses X of | < r^mWherein X is^mE.g., X (: z), and the transient thermal responses are considered to belong to the same cluster NDs;

step S64, taking the Center as the Center, calculates the vector of each transient thermal response in the Center to cluster set NDs, and adds these vectors to obtain the shifted mean:

wherein, omega is a high-dimensional sphere area with the Center as the Center radius r; x^mSatisfying [ Center-X ] for transient thermal response distributed in omega region^mL < r; k is transient thermal response X distributed in an omega region^mThe number of (2); g (x) is a binuclear function

The derivation is negative;

step S65, the Center updates along the direction of the shift mean, and the movement distance is | | shift (Center) |:

^tthe Center represents the updated cluster Center;

step S66, if | | | shift (Center) | | | ≧ epsilon, let Center | |^tThe Center, t is t +1, the step S63 is returned to continue the iterative computation until | | shift (Center) | < epsilon, and the previous iteration times t are found out to be distributed so as to^tAll transient thermal response sets in the high-dimensional sphere area with Center as the Center and radius r, namely satisfying | survival^tCenter-X^b||＜r，X^bAll transient thermal responses X for e X (: z)^bThese transient thermal responses are considered to belong to the same cluster^tNDs；

Step S67, calculating the updated cluster center^tCenter and other already existing cluster centers_i'Distance of Center, if any^tCenter-_i'If Center | < ω is true, then consider it to be^tCenter and the nearest one to it_kThe Center belongs to the same class, wherein

Go to step S68; otherwise go to step S69;

step S68, collecting^tIn NDsThe transient thermal response is ascribed to_kCenter, k is 1,2, …, and L is a cluster of the cluster Center_kNDs, will be assembled^tTransient thermal response markers in NDs are accessed and clustered_kThe number of occurrences in NDs is added to 1 and calculated^tCenter and_kweighted average of Center^t'Center，_kNeutral separation of NDs^t'Recent transient thermal response of Center^t″The Center is the cluster Center of the current cluster, i.e. the current cluster Center_kCenter＝^t″Center, the number of clusters L equals L, let Center equal to L^t″The Center goes to step S610;

step S69, collecting^tTransient thermal response in NDs is ascribed to_kCenter, k is 1,2, …, and L is a cluster of the cluster Center_kNDs, will be assembled^tTransient thermal response markers in NDs are accessed and clustered_kThe number of occurrences in NDs is added to 1 and calculated^tCenter and_kweighted average of Center^t'Center，_kNeutral separation of NDs^t'Recent transient thermal response of Center^t″The Center is the cluster Center of the current cluster, i.e. the current cluster Center_kCenter＝^t″Center, the number of clusters L equals L, let Center equal to L^t″The Center goes to step S610;

and step S610, making T equal to T +1, returning to step S63, ending the algorithm and clustering, and outputting the clustering number L and the clustering center of each cluster when the loop ending condition T equal to T is met and the number of accessed transient thermal responses marked in X (: z) reaches 90 percent_i'Center, i ═ 1,2, …, L; the transient thermal response marked as visited in the set X (: z) goes to step 611 to classify the category to which it belongs; the transient thermal response which is not marked as visited is transferred to step S612 to divide the category to which the transient thermal response belongs; otherwise, returning to step S63 to continue the iterative computation until the number of transient thermal responses marked as visited in X (: z) reaches 90%;

step S611, performing category classification on the accessed transient thermal response: according to the occurrence frequency of each transient thermal response in each cluster, the class with the most occurrence frequency is taken out to be used as the current classThe class of the transient thermal response divides the current transient thermal response into the classes; when a certain transient thermal response X occurs^p，X^pWhen the times of appearance of epsilon X (: z) in several clustering clusters are the same, the transient thermal response set divided into the class is determined from each clustering cluster_i'NDs, i ═ 1,2, …, L found free from X^pK transient thermal responses with Euclidean distance being nearest are calculated to obtain X^pTo_i'Crowdedness distance of k nearest transient thermal responses in NDs:

wherein,

represents X^pEuclidean distance to the j, j ═ 1,2, …, k nearest transient thermal responses in class i'; mixing X^pIf the calculated congestion degree distance is the largest, the step goes to step S613 to output the transient thermal response division condition;

step S612, classifying the transient thermal response which is not marked as accessed: calculating transient thermal response X not marked as visited^q，X^qBelongs to X (: z) to each category clustering center_i'European distance from Center: categorizing transient thermal responses that are not marked as visited: calculating transient thermal response X not marked as visited^q，X^qBelongs to X (: z) to each category clustering center_i'European distance from Center:

wherein,_i'center, i ═ 1,2, …, L is the cluster Center for each cluster. Calculating the class with the nearest Euclidean distance as the current transient thermal response X^qThe step S613 is switched to output the transient thermal response classification condition;

step S613, outputting the transient thermal response division condition_i'NDs, i ═ 1,2, …, L, and the cluster centers corresponding to each class_i'Center，_i'NDs represents the set of transient thermal responses of class i'.

In the above technical solution, the seventh step of selecting a representative of each type of transient thermal response by using a decomposition multi-objective optimization algorithm based on gradient search weight adjustment, and forming the matrix Y includes the following steps:

in step S71, when the i '(i' ═ 1., L) th class transient response is represented, a multi-objective function is defined:

minimize F(_i'X)＝(f₁(_i'X),...,f_L(_i'X))^T

wherein f is₁(_i'X) a transient response selected for the i' th class transient response_i'The intra-class Euclidean distance of X is expressed as:

f_l(_i'x), L2, 3, L is a transient response selected from the i' th class of transient responses_i'The Euclidean distance between L-1 classes of X is calculated according to the calculated Euclidean distance between L-1 classes

The components are renumbered and the components are,

expressed as:

_i'x_hfor transient response_i'The pixel value of X at the h-th time, i.e. the temperature value,_i'Center_hthe pixel value of the ith' type transient response cluster center at the h-th moment, namely the temperature value,_j'Center_hthe pixel value of the j' th class transient response cluster center at the h moment is the temperature value;

step S72, Multi-objective evolution Algorithm Based on Decomposition, MOEA/D, step S71 to give Multi-objective function, and select the representative of the i' th class transient response_i'REP, i' belongs to (1, 2.. gtorel), and the specific steps comprise:

step S721, initializing the population size to be N, and randomly generating N weight vectors uniformly distributed by a Tchebycheff decomposition algorithm

Wherein

Corresponding direction vector

Comprises the following steps:

wherein,

step S722, decomposing the L multi-target questions into N sub-questions by utilizing a Tchebycheff decomposition algorithm, wherein each sub-question is as follows:

wherein,

_i'z^*for the reference point corresponding to the multi-objective function of the ith' class,

is a function f_l(_i'X) a corresponding reference point;

step S723, initializing the evolution time h to 0, to obtain_i'Randomly generating N values in Ω

_i'Xⁿ(0) Indicating the value of the nth individual in the ith' class at initialization,_i'omega is a value range determined by transient thermal response of the ith class T time dimension;

step S724, initializing Pareto leading edge solution set_i'Initializing reference point for EP (0) < phi >

Wherein,

step S725, initializing neighborhood vectors, and solving the nth weight vector of the ith class

The most recent TT weight vectors;

step S726, initializing position information and speed information of the particle swarm according to the independent variable space dimension;

step S727, setting the maximum number of iterations P, starting iteration when h is 1, terminating iteration when h is P, setting a threshold value δ and a parameter r, and applying each weight vector to each weight vector

And performing the following updating operation according to the iteration number h:

step A, updating individuals, and specifically comprising the following steps:

step A1, setting the scope of evolution: setting a threshold value delta according to a random number rand epsilon (0,1) generated randomly, and determining a neighborhood range P_i'(n)：

Wherein, B_i'(n) is a weight vector index set;

step A2, calculating the particle fitness according to a fitness algorithm, and dividing the initial population into two subgroups which are a non-dominant subset P _ set and a dominant subset NP _ set respectively;

step A3, calculating the adaptive value of each particle, and for each particle, matching the adaptive value with the best position p experienced^bestIf the adaptive value is better, the adaptive value is taken as the current best position; for each particle, its adaptation value is compared with the global experienced best position g^bestComparing the adaptive values, and if the adaptive values are better, taking the adaptive values as the current global best position;

step a4, updating its position information and velocity information by means of the individual extremum and the global extremum for each particle of the dominant subset NP _ set:

v_ji(t+1)＝v_id(t)+c₁·r₁·(P_ji(t)-x_ji(t))+c₂·r₂·(NP_ji(t)-x_ji(t))

x_ji(t+1)＝x_ji(t)+v_ji(t+1)

wherein, c₁And c₂Is a normal number, called acceleration factor, c₁Adjusting the step length of the particles flying to the direction of the best position of the particles; c. C₂Adjusting the step length of the flight of the particles to the global best position; r is₁And r₂Is [0,1 ]]A random number in between; p_ji(t) are random individuals within the non-dominated solution set P _ set; NP_ji(t) are the particles within the dominance solution set NP _ set;

step B, updating the weight vector and the direction vector, and specifically comprising the following steps:

step B1, on the h-th iteration, for all individuals in class i

According to the preference, considering the similarity, the first iteration is obtained according to the following formulaOf group i

Reference points:

wherein,

is an individual

The Euclidean distance in class from the current i' th class center is used for measuring similarity, sum represents the sum of Euclidean distances in class, and the sum which enables sum to be minimum is selected

The individual is taken as a reference point, and the point is recorded

The individual is composed of a set of

Step B2, for each individual reference point

Find out in

High dimensional sphere region omega with radius r as center_hAll individuals in the table, calculate Ω_hStandard deviation of all individuals within

Wherein,

is distributed at omega_hAverage of all individuals within;

represents a distribution in Ω_hK individuals within;

step B3, calculating the standard deviation satisfying

The number of reference point individuals satisfying the condition is recorded as s, and the set of the s reference point individuals is recorded as

To the collection

Each individual reference point in

The following update operations are performed:

step B31, calculating basis weight vector:

and converted into the corresponding direction vector:

wherein,

for the reference point corresponding to the multi-objective function of the ith' class,

is a function f corresponding to a subproblem decomposed by the Chebyshev algorithm_l(_i'X) of the reference points of the two,

step B32, finding the reference point in the generated population

The most distant European individuals

Wherein,

step B33, mixing the individuals

Corresponding weight vector

Conversion to the corresponding directional vector:

step B34, calculating the generated direction vector:

wherein,

step B35, utilizing reference point

And

generate a new solution

Step B36, using the new solution

Replacement of

Direction vector

Replacement of

Corresponding weight vector

Replacement of

Wherein,

step 728, loop judgment: if n is>N, h is h +1, i.e. the next evolution is carried out; if the termination condition is not met, repeating the step S727, otherwise, obtaining a final Pareto leading edge solution set of the ith temperature transient thermal response_i'EP(h)；

Step S729, from_i'Selecting a temperature transient thermal response optimization solution in EP (h)_i'REP, representative of class i' transient responses_i'REP，i'∈(1,2,...,L)；

Step S73, the L-type transient response representatives are arranged in columns, and one column is the temperature values, which are the pixel values at T moments, to form a T × L matrix Y.

In the above technical solution, the step eight of converting the three-dimensional matrix M into a two-dimensional matrix, and performing linear change on the two-dimensional matrix by using the matrix Y to obtain the two-dimensional image matrix R (x, Y) specifically includes: starting each frame in the three-dimensional matrix M from a first column, connecting a next column to the tail of a previous column to form a new column, obtaining T-column pixel values corresponding to the T frames, then sequentially placing the T-column pixel values according to time sequence to form a two-dimensional matrix O, and performing linear transformation on the two-dimensional matrix O by using a matrix Y, namely:

a two-dimensional image matrix R is obtained, wherein,

is a pseudo-inverse of the matrix Y, O^TA transpose of the two-dimensional matrix O.

In the above technical solution, the ninth step utilizes an improved one-dimensional Ostu segmentation algorithm to perform feature extraction on a two-dimensional image matrix R (x, y), and divides the image into the target C by setting a threshold k according to the gray features of the image₀And background C₁The method comprises the following two steps:

step S91, calculating the gray probability distribution P of the image_t：

Where, t is 1,2, …, L is the number of gray levels of the image, N is the total number of pixels in the image, and N is the total number of pixels in the image_tThe number of pixel points with the gray level of t in the image is counted;

step S92, calculating the target class C₀And background class C₁Probability of occurrence of each omega₀And ω₁And their corresponding mean values of gray levels mu₀And mu₁Thereby obtaining the spacing between the two types

And the degree of dispersion d of each class is determined₀And d₁：

Step S93, calculating a threshold value selection function H (t);

step S94, extracting the maximum value in H (t), wherein the optimal threshold k is the corresponding gray level t;

and step S95, binarizing the image according to k to realize segmentation of the gray level image, thereby completing identification of the impact damage characteristics of the surface of the test piece.

Example (b):

in this example, there are two defects on the test piece, namely a flat-bottom hole defect 1 whose surface is visible and an internal defect 2 whose surface is not visible.

In this embodiment, a result graph obtained by classifying the selected transient thermal response by using a mean shift clustering algorithm is shown in fig. 4.

Three known temperature points, namely transient thermal response curves of a material temperature point, a defect 1 temperature point and a defect 2 temperature point are directly extracted from a thermal image sequence of the test piece and are respectively recorded as_BacPOINT、_Def1POINT and_Def2POINT, as shown in FIGS. 8, 9,As shown in fig. 10.

Three transient thermal response representations are obtained by using the existing method for selecting the transient thermal response representation based on difference:_ANK₄、_BNK₅₀and_cNK₅the curves are shown in fig. 11, 12, and 13, and correspond to the material temperature point, the defect 1 temperature point, and the defect 2 temperature point, respectively.

By using the method for selecting the transient thermal response representatives through multi-objective optimization, three transient thermal response representatives are obtained:_ANK₂₀、_BNK₂and_cNK₂₁the curves are shown in fig. 14, 15, and 16, and correspond to the material temperature point, the defect 1 temperature point, and the defect 2 temperature point, respectively.

From the thermal response curves, it can be seen that: the temperature point of the defect 2 rises at the fastest rate and has the highest peak value. After the heat source is removed, the temperature rises reversely due to the waste heat of the heat source and the thermal diffusion phenomenon of the material, and compared with the three characteristics, the reverse rising rate of the temperature point of the defect 1 is the fastest;

the correlation between the transient thermal response curves under the two methods and the corresponding transient thermal response curves extracted directly from the thermographic sequence is shown in table 1.

TABLE 1

	Temperature point of itself	Temperature point of defect 1	Temperature point of defect 2
				Based on a difference method	0.9634	0.9701	0.9941
The invention	0.8470	0.9051	0.8777

From table 1, it can be seen that the transient thermal response curves selected by the method of the present invention are more correlated.

In the present embodiment, the defect features extracted are as shown in fig. 17.

The number of apparatuses and the scale of the process described herein are intended to simplify the description of the present invention. Applications, modifications and variations of the present invention will be apparent to those skilled in the art.

While embodiments of the invention have been described above, it is not limited to the applications set forth in the description and the embodiments, which are fully applicable in various fields of endeavor to which the invention pertains, and further modifications may readily be made by those skilled in the art, it being understood that the invention is not limited to the details shown and described herein without departing from the general concept defined by the appended claims and their equivalents.

Claims (8)

1. A reusable spacecraft surface impact damage feature identification method is characterized by comprising the following steps:

step one, representing a thermal image sequence of a pixel point on the surface of a test piece acquired by an infrared thermal imager by using a three-dimensional matrix M, wherein M (i, j, T) represents a pixel in the ith row and the jth column, the third dimension T represents acquisition time, and T is 1,2, … and T;

step two, finding a temperature peak point M (i) in the three-dimensional matrix M_max,j_max,t_max) Wherein i_maxIndicates the number of rows where the temperature peak point is located, j_maxIndicates the number of columns where the temperature peaks are located, t_maxRepresenting the frame number of the temperature peak point;

step three, dividing the temperature of the line of the frame number of the temperature peak point, dividing the line of the frame number of the temperature peak point into R data blocks, and recording the line step length;

step four, carrying out temperature division on the row of the frame number of the temperature peak point, dividing the row of the frame number of the temperature peak point into Q data blocks, and recording the row step length;

step five, extracting the transient thermal response in the data block obtained in the step three and the step four in a blocking and step-length manner, and finishing the extraction of the typical transient thermal response;

step six, adopting a mean shift algorithm to automatically classify the typical transient thermal responses extracted in the step five to obtain the number of clustered damage sets and the category of each transient thermal response;

selecting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on gradient search weight adjustment to form a matrix Y;

step eight, changing the three-dimensional matrix M into a two-dimensional matrix, and performing linear change on the two-dimensional matrix by using the matrix Y to obtain a two-dimensional image matrix R (x, Y);

ninthly, extracting features of the two-dimensional image matrix R (x, y) by using an improved one-dimensional Ostu segmentation algorithm, and dividing the image into targets C by setting a threshold k according to the gray features of the image₀And background C₁And obtaining the defect characteristics of the test piece.

2. The method for identifying the surface impact damage characteristics of the reusable spacecraft of claim 1, wherein the step three is to perform temperature division on the row of the frame number of the temperature peak point, and the specific method for dividing the row of the frame number of the temperature peak point into R data blocks comprises the following steps: setting a threshold value

Calculating the temperature peak point M (i)_max,j_max,t_max) To the nearest temperature point M^c(i_max,j_c,t_max) Is calculated by the correlation coefficient PC_reWherein j is_cThe number of columns where the nearest temperature point is located, if

Then j is_c＝j_c+1, calculating the temperature peak M (i) from near to far_max,j_max,t_max) Correlation coefficient PC with the remaining temperature points_reUp to

When it is time, stop calculating, remember

Number of established pixels as line step length CL^rAccording to CL^rAnd dividing the temperature into lines, and dividing the lines of the frame number of the temperature peak point into R data blocks.

3. The method for identifying the surface impact damage characteristics of the reusable spacecraft as claimed in claim 1, wherein the step four includes the step of performing temperature division on the frame number columns where the temperature peak points are located, and the specific method of dividing the frame number columns where the temperature peak points are located into Q data blocks includes: setting a threshold value

Calculating the temperature peak point M (i)_max,j_max,t_max) To the nearest temperature point M^c(i_c,j_max,t_max) Is calculated by the correlation coefficient PC_reWherein i_cThe number of rows where the nearest temperature point is located, if

Then i_c＝i_c+1 calculating the temperature peak M (i) from near to far_max,j_max,t_max) Correlation coefficient PC with the remaining temperature points_reUp to

When it is time, stop calculating, remember

The number of established pixels is the column step length RL^qAccording to RL^qAnd dividing the temperature into columns, and dividing the columns of the frame number of the temperature peak point into Q data blocks.

4. The method for identifying the surface impact damage characteristics of the reusable spacecraft of claim 1, wherein the step of extracting the transient thermal response in five blocks and steps comprises the following specific steps:

step S51, step three, and step four divide the three-dimensional matrix M into N data blocks, where N is R · Q, and denote the nth data block as MⁿAnd setting a threshold CC, where N is 1,2, …, N;

step S52, the temperature peak point M (i)_max,j_max,t_max) Corresponding transient thermal response M (i)_max,j_maxT), T is 1,2, …, T is stored in the set X (: 1), and the transient thermal response M of the mth row and nth column of each data block is calculatedⁿ(m, n, T), T ═ 1,2, …, correlation PC between T and X (: 1)_re(ii) a If PC_re＜CC，Mⁿ(m, n, T), T ═ 1,2, …, T is regarded as a new feature and stored in X (: z); otherwise, let m be m + RL^qContinuously calculating the correlation degree of the next temperature point and X (: 1); if m > k, the column calculation is performed: let m be m-k, n be n + CL^rAnd k is the row number of the nth data block, and if n is more than l, wherein l is the column number of the nth data block, the typical transient thermal response extraction is finished.

5. The method for identifying surface impact damage features of a reusable spacecraft of claim 1, wherein the step six employs a mean shift algorithm, and the method for automatically classifying the typical transient thermal response extracted in the step five comprises: defaulting that all data belong to one class, namely the cluster number L is 1, calculating by an iterative algorithm, dividing the transient thermal response of the pixel points with larger difference into the other class, updating the cluster number L to be L +1, after all the transient thermal responses in a transient thermal response set X (: z) are subjected to iterative calculation, finally obtaining the cluster damage set number L and the class to which each transient thermal response belongs, and the specific steps comprise:

step S61, setting a bandwidth range r, an iteration termination condition epsilon and a threshold omega;

step S62, randomly selecting a transient thermal response in the set X (: z) as the initial cluster Center point Center, setting the cluster number L to 1, the cycle parameter T to 1, i' to 1, the cycle end condition T to T, keeping the initial cluster Center when T is 1, and recording the initial cluster Center when T is 1₁Center＝Center；

Step S63, find out all transient thermal response sets distributed in the high-dimensional sphere area with Center as the Center and radius r, namely, satisfy | Center-X^mAll transient thermal responses Xm of | < r, where X^mE.g., X (: z), and the transient thermal responses are considered to belong to the same cluster NDs;

step S64, taking the Center as the Center, calculates the vector of each transient thermal response in the Center to cluster set NDs, and adds these vectors to obtain the shifted mean:

wherein, omega is a high-dimensional sphere area with the Center as the Center radius r; x^mSatisfying [ Center-X ] for transient thermal response distributed in omega region^mL < r; k is transient thermal response X distributed in an omega region^mThe number of (2); g (x) is a binuclear function

The derivation is negative;

step S65, the Center updates along the direction of the shift mean, and the movement distance is | | shift (Center) |:

^tthe Center represents the updated cluster Center;

step S66, if | | | shift (Center) | | | ≧ epsilon, let Center | |^tThe Center, t is t +1, the step S63 is returned to continue the iterative computation until | | shift (Center) | < epsilon, and the previous iteration times t are found out to be distributed so as to^tAll transient thermal response sets in the high-dimensional sphere area with Center as the Center and radius r, namely satisfying | survival^tCenter-X^b||＜r，X^bAll transient thermal responses X for e X (: z)^bThese transient thermal responses are considered to belong to the same cluster^tNDs；

Step S67, calculating the updated cluster center^tCenter and other already existing cluster centers_i'Distance of Center, if any^tCenter-_i'If Center | < ω is true, then consider it to be^tCenter and the nearest one to it_kThe Center belongs to the same class, wherein

Go to step S68; otherwise go to step S69;

step S68, collecting^tTransient thermal response in NDs is ascribed to_kCenter, k is 1,2, …, and L is a cluster of the cluster Center_kNDs, will be assembled^tTransient thermal response markers in NDs are accessed and clustered_kThe number of occurrences in NDs is added to 1 and calculated^tCenter and_kweighted average of Center^t'Center，_kNeutral separation of NDs^t'Recent transient thermal response of Center^t”The Center is the cluster Center of the current cluster, i.e. the current cluster Center_kCenter＝^t”Center, the number of clusters L equals L, let Center equal to L^t”The Center goes to step S610;

step S69, collecting^tTransient thermal response in NDs is ascribed to_kCenter, k is 1,2, …, and L is a cluster of the cluster Center_kNDs, will be assembled^tTransient thermal response marking in NDs as visitedAnd cluster it in clusters_kThe number of occurrences in NDs is added to 1 and calculated^tCenter and_kweighted average of Center^t'Center，_kNeutral separation of NDs^t'Recent transient thermal response of Center^t”The Center is the cluster Center of the current cluster, i.e. the current cluster Center_kCenter＝^t”Center, the number of clusters L equals L, let Center equal to L^t”The Center goes to step S610;

and step S610, making T equal to T +1, returning to step S63, ending the algorithm and clustering, and outputting the clustering number L and the clustering center of each cluster when the loop ending condition T equal to T is met and the number of accessed transient thermal responses marked in X (: z) reaches 90 percent_i'Center, i ═ 1,2, …, L; the transient thermal response marked as visited in the set X (: z) goes to step 611 to classify the category to which it belongs; the transient thermal response which is not marked as visited is transferred to step S612 to divide the category to which the transient thermal response belongs; otherwise, returning to step S63 to continue the iterative computation until the number of transient thermal responses marked as visited in X (: z) reaches 90%;

step S611, performing category classification on the accessed transient thermal response: according to the occurrence frequency of each transient thermal response in each cluster, taking the class with the largest occurrence frequency as the class to which the current transient thermal response belongs, and dividing the current transient thermal response into the classes; when a certain transient thermal response X occurs^p，X^pWhen the times of appearance of epsilon X (: z) in several clustering clusters are the same, the transient thermal response set divided into the class is determined from each clustering cluster_i'NDs, i ═ 1,2, …, L found free from X^pK transient thermal responses with Euclidean distance being nearest are calculated to obtain X^pTo_i'Crowdedness distance of k nearest transient thermal responses in NDs:

wherein,

represents X^pEuclidean distance to the j, j ═ 1,2, …, k nearest transient thermal responses in class i'; mixing X^pIf the calculated congestion degree distance is the largest, the step goes to step S613 to output the transient thermal response division condition;

step S612, classifying the transient thermal response which is not marked as accessed: calculating transient thermal response X not marked as visited^q，X^qBelongs to X (: z) to each category clustering center_i'European distance from Center: categorizing transient thermal responses that are not marked as visited: calculating transient thermal response X not marked as visited^q，X^qBelongs to X (: z) to each category clustering center_i'European distance from Center:

wherein,_i'center, i ═ 1,2, …, L is the cluster Center for each cluster; calculating the class with the nearest Euclidean distance as the current transient thermal response X^qThe step S613 is switched to output the transient thermal response classification condition;

step S613, outputting the transient thermal response division condition_i'NDs, i ═ 1,2, …, L, and the cluster centers for each class_i'Center，_i'NDs represents the set of transient thermal responses of class i'.

6. The method for identifying the surface impact damage characteristics of the reusable spacecraft of claim 1, wherein the seventh step adopts a decomposition multi-objective optimization algorithm based on gradient search weight adjustment to select a representative of each type of transient thermal response, and the forming of the matrix Y comprises the following steps:

in step S71, when the i '(i' ═ 1., L) th class transient response is represented, a multi-objective function is defined:

minimize F(_i'X)＝(f₁(_i'X),...,f_L(_i'X))^T

wherein f is₁(_i'X) a transient response selected for the i' th class transient response_i'The intra-class Euclidean distance of X is expressed as:

f_l(_i'x), L2, 3, L is a transient response selected from the i' th class of transient responses_i'The Euclidean distance between L-1 classes of X is calculated according to the calculated Euclidean distance between L-1 classes

The components are renumbered and the components are,

expressed as:

_i'x_hfor transient response_i'The pixel value of X at the h-th time, i.e. the temperature value,_i'Center_hthe pixel value of the ith' type transient response cluster center at the h-th moment, namely the temperature value,_j'Center_hthe pixel value of the j' th class transient response cluster center at the h moment is the temperature value;

step S72, Multi-objective evolution Algorithm Based on Decomposition, MOEA/D, step S71 to give Multi-objective function, and select the representative of the i' th class transient response_i'REP, i' belongs to (1, 2.. gtorel), and the specific steps comprise:

step S721, initializing the population size to be N, and randomly generating N weight vectors uniformly distributed by a Tchebycheff decomposition algorithm

Wherein

Corresponding direction vector

Comprises the following steps:

wherein,

step S722, decomposing the L multi-target questions into N sub-questions by utilizing a Tchebycheff decomposition algorithm, wherein each sub-question is as follows:

wherein,

_i'z^*for the reference point corresponding to the multi-objective function of the ith' class,

is a function f_l(_i'X) a corresponding reference point;

step S723, initializing the evolution time h to 0, to obtain_i'Randomly generating N values in Ω

_i'Xⁿ(0) Indicating the value of the nth individual in the ith' class at initialization,_i'when Ω is T of the i' th classThe value range determined by transient thermal response of the inter-dimension;

step S724, initializing Pareto leading edge solution set_i'Initializing reference point for EP (0) < phi >

Wherein,

step S725, initializing neighborhood vectors, and solving the nth weight vector of the ith class

The most recent TT weight vectors;

step S726, initializing position information and speed information of the particle swarm according to the independent variable space dimension;

step S727, setting the maximum number of iterations P, starting iteration when h is 1, terminating iteration when h is P, setting a threshold value δ and a parameter r, and applying each weight vector to each weight vector

And performing the following updating operation according to the iteration number h:

step A, updating individuals, and specifically comprising the following steps:

step A1, setting the scope of evolution: setting a threshold value delta according to a random number rand epsilon (0,1) generated randomly, and determining a neighborhood range P_i'(n)：

Wherein, B_i'(n) is a weight vector index set;

step A2, calculating the particle fitness according to a fitness algorithm, and dividing the initial population into two subgroups which are a non-dominant subset P _ set and a dominant subset NP _ set respectively;

step A3, calculating an adaptive value for each particle, for eachThe particles adapt their values to the best positions p experienced^bestIf the adaptive value is better, the adaptive value is taken as the current best position; for each particle, its adaptation value is compared with the global experienced best position g^bestComparing the adaptive values, and if the adaptive values are better, taking the adaptive values as the current global best position;

step a4, updating its position information and velocity information by means of the individual extremum and the global extremum for each particle of the dominant subset NP _ set:

v_ji(t+1)＝v_id(t)+c₁·r₁·(P_ji(t)-x_ji(t))+c₂·r₂·(NP_ji(t)-x_ji(t))

x_ji(t+1)＝x_ji(t)+v_ji(t+1)

wherein, c₁And c₂Is a normal number, called acceleration factor, c₁Adjusting the step length of the particles flying to the direction of the best position of the particles; c. C₂Adjusting the step length of the flight of the particles to the global best position; r is₁And r₂Is [0,1 ]]A random number in between; p_ji(t) are random individuals within the non-dominated solution set P _ set; NP_ji(t) are the particles within the dominance solution set NP _ set;

step B, updating the weight vector and the direction vector, and specifically comprising the following steps:

step B1, on the h-th iteration, for all individuals in class i

N is 1,2, …, N, and the target in the aspect of considering similarity according to preference is obtained according to the following formula

Reference points:

wherein,

is an individual

The Euclidean distance in class from the current i' th class center is used for measuring similarity, sum represents the sum of Euclidean distances in class, and the sum which enables sum to be minimum is selected

The individual is taken as a reference point, and the point is recorded

The individual is composed of a set of

Step B2, for each individual reference point

Find out in

High dimensional sphere region omega with radius r as center_hAll individuals in the table, calculate Ω_hStandard deviation of all individuals within

Wherein,

is distributed at omega_hAverage of all individuals within;

represents a scoreDistributed at omega_hK individuals within;

step B3, calculating the standard deviation satisfying

The number of reference point individuals satisfying the condition is recorded as s, and the set of the s reference point individuals is recorded as

To the collection

Each individual reference point in

The following update operations are performed:

step B31, calculating basis weight vector:

and converted into the corresponding direction vector:

wherein,

for the reference point corresponding to the multi-objective function of the ith' class,

is a function f corresponding to a subproblem decomposed by the Chebyshev algorithm_l(_i'X) of the reference points of the two,

step B32, finding the reference point in the generated population

The most distant European individuals

Wherein,

step B33, mixing the individuals

Corresponding weight vector

Conversion to the corresponding directional vector:

step B34, calculating the generated direction vector:

wherein,

step B35, utilizing reference point

And

generate a new solution

Step B36, using the new solution

Replacement of

Direction vector

Replacement of

Corresponding weight vector

Replacement of

Wherein,

step 728, loop judgment: if n is>N, h is h +1, i.e. the next evolution is carried out; if the termination condition is not met, repeating the step S727, otherwise, obtaining a final Pareto leading edge solution set of the ith temperature transient thermal response_i'EP(h)；

Step S729, from_i'Selecting a temperature transient thermal response optimization solution in EP (h)_i'REP, representative of class i' transient responses_i'REP，i'∈(1,2,...,L)；

Step S73, the L-type transient response representatives are arranged in columns, and one column is the temperature values, which are the pixel values at T moments, to form a T × L matrix Y.

7. As claimed in claim 1The reusable spacecraft surface impact damage characteristic identification method is characterized in that the step eight is to change a three-dimensional matrix M into a two-dimensional matrix, and the specific method for obtaining a two-dimensional image matrix R (x, Y) by utilizing a matrix Y to carry out linear change on the matrix comprises the following steps: starting each frame in the three-dimensional matrix M from a first column, connecting a next column to the tail of a previous column to form a new column, obtaining T-column pixel values corresponding to the T frames, then sequentially placing the T-column pixel values according to time sequence to form a two-dimensional matrix O, and performing linear transformation on the two-dimensional matrix O by using a matrix Y, namely:

a two-dimensional image matrix R is obtained, wherein,

is a pseudo-inverse of the matrix Y, O^TA transpose of the two-dimensional matrix O.

8. The method for recognizing the surface impact damage characteristics of the reusable spacecraft as claimed in claim 1, wherein the nine steps utilize a modified one-dimensional Ostu segmentation algorithm to perform feature extraction on a two-dimensional image matrix R (x, y), and divide the image into the targets C by setting a threshold k according to the gray features of the image₀And background C₁The method comprises the following two steps:

step S91, calculating the gray probability distribution P of the image_t：

Where, t is 1,2, …, L is the number of gray levels of the image, N is the total number of pixels in the image, and N is the total number of pixels in the image_tThe number of pixel points with the gray level of t in the image is counted;

step S92, calculating the target class C₀And background class C₁Probability of occurrence of each omega₀And ω₁And their corresponding mean values of gray levels mu₀And mu₁Thereby obtaining the spacing between the two types

And the degree of dispersion d of each class is determined₀And d₁：

Step S93, calculating a threshold value selection function H (t);

step S94, extracting the maximum value in H (t), wherein the optimal threshold k is the corresponding gray level t;

and step S95, binarizing the image according to k to realize segmentation of the gray level image, thereby completing identification of the impact damage characteristics of the surface of the test piece.

Images