CN112215830A - Method for judging impact damage characteristic types of aerospace heat-proof materials - Google Patents

Method for judging impact damage characteristic types of aerospace heat-proof materials Download PDF

Info

Publication number
CN112215830A
CN112215830A CN202011132748.9A CN202011132748A CN112215830A CN 112215830 A CN112215830 A CN 112215830A CN 202011132748 A CN202011132748 A CN 202011132748A CN 112215830 A CN112215830 A CN 112215830A
Authority
CN
China
Prior art keywords
center
transient thermal
thermal response
temperature
max
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011132748.9A
Other languages
Chinese (zh)
Other versions
CN112215830B (en
Inventor
黄雪刚
雷光钰
殷春
石安华
周浩
赵君尧
柳森
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Original Assignee
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center filed Critical Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority to CN202011132748.9A priority Critical patent/CN112215830B/en
Publication of CN112215830A publication Critical patent/CN112215830A/en
Application granted granted Critical
Publication of CN112215830B publication Critical patent/CN112215830B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N25/00—Investigating or analyzing materials by the use of thermal means
    • G01N25/72—Investigating presence of flaws
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/194—Segmentation; Edge detection involving foreground-background segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10048—Infrared image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20081—Training; Learning

Abstract

The invention discloses a method for judging impact damage characteristic types of aerospace heat-proof materials, 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 point in the three-dimensional matrix; 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; extracting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on self-adaptive weight adjustment to form a matrix Y; converting the three-dimensional matrix into a two-dimensional matrix, and performing linear transformation 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. The method makes up the defects of the traditional clustering method in setting the clustering number and selecting the initial clustering center, and ensures the accuracy of defect feature extraction.

Description

Method for judging impact damage characteristic types of aerospace heat-proof materials
Technical Field
The invention belongs to the technical field of damage detection and evaluation of spacecrafts, and particularly relates to a method for judging impact damage characteristic types of a spaceflight heat-proof material.
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 piecerIt 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 guaranteedrKeeping 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 typesmax,Lmax＞LrClustering, 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 pointsrHowever, when judging the number of defect types, L is obtained by selecting different judgment difference modesrThe result of (2) is also different, and L is reusedrAnd 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 utilizing an algorithm, can obtain an accurate clustering result for the infrared thermal image sequence data which is acquired by experiments and has high dimension and complex distribution, ensures the accuracy of further analyzing and mining the data, has better effect of selecting the thermal response of the representative temperature point by utilizing a multi-objective optimization algorithm, and has better representative temperature resultAnd (5) table property.
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 of the direction vector, when the shape of the front surface of the PF can be approximately displayed by a population through a certain iteration number, the direction vector is periodically adjusted according to the sparseness and density degree of population distribution, and then the weight vector is distributed for the subproblems according to the direction vector, so that a uniformly distributed Pareto optimal solution set is obtained to maintain the diversity of the population, and 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 method for judging impact damage characteristic types of aerospace thermal protection materials, comprising the steps of:
step one, representing a thermal image sequence of a test piece pixel point 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 Mmax,jmax,tmax) Wherein imaxIndicates the number of rows where the temperature peak point is located, jmaxIndicates the number of columns where the temperature peaks are located, tmaxRepresenting 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;
extracting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on self-adaptive weight adjustment to form a matrix Y;
step eight, changing the three-dimensional matrix M into a two-dimensional matrix, and performing linear transformation on the two-dimensional matrix by using a 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 image0And background C1And obtaining the impact damage 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
Figure BDA0002735683430000031
Calculating the temperature peak point M (i)max,jmax,tmax) To the nearest temperature point Mc(imax,jc,tmax) Is calculated by the correlation coefficient PCreWherein j iscThe number of columns where the nearest temperature point is located, if
Figure BDA0002735683430000041
Then j isc＝jc+1 calculating the temperature peak M (i) from near to farmax,jmax,tmax) Correlation coefficient PC with the remaining temperature pointsreUp to
Figure BDA0002735683430000042
When it is time, stop calculating, remember
Figure BDA0002735683430000043
Number of established pixels as line step length CLrAccording to CLrAnd 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
Figure BDA0002735683430000044
Calculating the temperature peak point M (i)max,jmax,tmax) To the nearest temperature point Mc(ic,jmax,tmax) Is calculated by the correlation coefficient PCreWherein icThe number of rows where the nearest temperature point is located, if
Figure BDA0002735683430000045
Then ic＝ic+1 calculating the temperature peak M (i) from near to farmax,jmax,tmax) Correlation coefficient PC with the remaining temperature pointsreUp to
Figure BDA0002735683430000046
At the time of stoppingCalculating and recording
Figure BDA0002735683430000047
The number of established pixels is the column step length RLqAccording to RLqAnd 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, step four divide the three-dimensional matrix M into N data blocks, and represent the nth data block as MnAnd setting a threshold CC, wherein N is R · Q, N is 1, 2.., N;
step S52, the temperature peak point M (i)max,jmax,tmax) Corresponding transient thermal response M (i)max,jmaxT), 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 calculatedn(m, n, T), T ═ 1,2, …, correlation PC between T and X (: 1)re(ii) a If PCre＜CC，Mn(m, n, T), T ═ 1,2, …, T is regarded as a new feature and stored in X (: z); otherwise, let m be m + RLqContinuously 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 + CLrAnd 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: initially 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 large difference into another class, updating the cluster number L to be L +1, after all transient thermal responses in a transient thermal response set X (: z) are subjected to iterative calculation, finally obtaining the clustered 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 11Center＝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-XmAll transient thermal responses X of | < rmWherein X ismE.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:
Figure BDA0002735683430000051
wherein, omega is a high-dimensional sphere area with the Center as the Center radius r; xmSatisfying [ Center-X ] for transient thermal response distributed in omega regionmL < r; k is transient thermal response X distributed in an omega regionmThe number of (2); g (x) is a logarithmic kernel function
Figure BDA0002735683430000052
The derivation is negative;
step S65, the Center updates along the direction of the shift mean, and the movement distance is | | shift (Center) |:
Figure BDA0002735683430000053
wherein the content of the first and second substances,tthe Center represents the updated cluster Center;
step S66, if | | | shift (Center) | | ≧ epsilon, let Center ═ epsilontCenter, t ═ t +1, return to step S63 and continue iterative computation until | | shift (Center) | < epsilon, compare beforeThe number of iterations t, find the distributiontAll transient thermal response sets in the high-dimensional sphere area with Center as the Center and radius r, namely satisfying | survivaltCenter-Xb||＜r，XbAll transient thermal responses X for e X (: z)bThese transient thermal responses are considered to belong to the same clustertNDs；
Step S67, calculating the updated cluster centertCenter and other already existing cluster centersi'Distance of Center, if anytCenter-i'If Center | < ω is true, then consider it to betCenter and the nearest one to itkThe Center belongs to the same class, wherein
Figure BDA0002735683430000054
Go to step S68; otherwise go to step S69;
step S68, collectingtTransient thermal response in NDs is ascribed tokCenter, k is 1,2, …, and L is a cluster of the cluster CenterkNDs, will be assembledtTransient thermal response markers in NDs are accessed and clusteredkThe number of occurrences in NDs is added to 1 and calculatedtCenter andkweighted average of Centert'Center，kNeutral separation of NDst'Recent transient thermal response of Centert”The Center is the cluster Center of the current cluster, i.e. the current cluster CenterkCenter＝t”Center, the number of clusters L equals L, let Center equal to Lt”The Center goes to step S610;
step S69, if nottCenter-i'Let i '═ i' +1, consider that a new cluster appears in the current iterationi'NDs of whichi'NDs＝tNDs，tThe Center is a new clusteri'Clustering centers of NDs, to be aggregatedtTransient thermal response markers in NDs are accessed and clusteredi' Adding 1 to the times of appearance in NDs, wherein the clustering number L is equal to L +1, and the current clustering centeri'Center＝tCenter, let Center ═tThe Center goes to step S610;
step S610 returns to t +1Step S63, until the loop termination condition T is satisfied, T and the number of accessed transient thermal responses marked in X (: z) reaches 90%, the algorithm is terminated, the clustering is ended, and the clustering number L and the clustering center of each cluster are outputi'Center, i ═ 1,2, …, L; the transient thermal response marked as visited in the set X (: z) goes to step S611 to classify the category of the transient thermal response; 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 occursp，XpWhen 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 clusteri'NDs, i ═ 1,2, …, L found free from XpK transient thermal responses with Euclidean distance being nearest are calculated to obtain XpToi'Crowdedness distance of k nearest transient thermal responses in NDs:
Figure BDA0002735683430000063
wherein the content of the first and second substances,
Figure BDA0002735683430000064
represents XpEuclidean distance to the j, j ═ 1,2, …, k nearest transient thermal responses in class i'; mixing XpIf 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 visitedq，XqBelongs to X (: z) to each category clustering centeri'Of CenterEuclidean distance:
Figure BDA0002735683430000061
wherein the content of the first and second substances,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 XqThe step S613 is switched to output the transient thermal response classification condition;
step S613, outputting the transient thermal response division conditioni'NDs, i ═ 1,2, …, L, and the cluster centers for each classi'Center，i'NDs represents the set of transient thermal responses of class i'.
Preferably, the seventh step of extracting a representative of each type of transient thermal response by using a decomposition multi-objective optimization algorithm based on adaptive weight adjustment, and the specific step of forming the matrix Y includes:
in step S71, when the i ', i' is 1, L classes of transient responses are represented, a multi-objective function is defined:
minimizeF(i'X)＝(f1(i'X),...,fL(i'X))T
wherein f is1(i'X) a transient response selected for the i' th class transient responsei'The intra-class Euclidean distance of X is expressed as:
Figure BDA0002735683430000062
wherein the content of the first and second substances,i'xhfor transient responsei'The pixel value of X at the h-th time, i.e. the temperature value,i'Centerhthe pixel value of the ith' type transient response cluster center at the h-th moment, namely the temperature value,j'Centerhthe pixel value of the j' th class transient response cluster center at the h moment is the temperature value;
step S72, Multi-object evolution Algorithm Based on Decomposition, MOEA/D, using the multi-objective function given in step S71 to select the representative of the i' th class transient responsei'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
Figure BDA0002735683430000071
Wherein
Figure BDA0002735683430000072
Corresponding direction vector
Figure BDA0002735683430000073
Comprises the following steps:
Figure BDA0002735683430000074
wherein the content of the first and second substances,
Figure BDA0002735683430000075
step 722, decomposing the L multi-target questions into N sub-questions by using a Tchebycheff decomposition algorithm, wherein each sub-question is as follows:
Figure BDA0002735683430000076
wherein the content of the first and second substances,
Figure BDA0002735683430000077
i'z*for the reference point corresponding to the multi-objective function of the ith' class,
Figure BDA0002735683430000078
is a function fl(i'X) a corresponding reference point;
step S723, initializing the evolution time h to 0, to obtaini'Randomly generating N values in Ω
Figure BDA0002735683430000079
Omega, wherein,i'Xn(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 seti'Initializing a reference point when EP (0) is phi;
step S725, initializing neighborhood vectors, and solving the nth weight vector of the ith class
Figure BDA00027356834300000710
The most recent TT weight vectors, and will be
Figure BDA00027356834300000711
The index set of the most recent TT weight vectors is marked as Bi'(n)＝{n1,n2,...,nTT}，n＝1,2,...,N；
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 the iteration when h is 1, terminating the iteration when h is P, and applying each weight vector
Figure BDA00027356834300000712
N is 1,2, …, and the updating operation is performed by N according to the iteration number h, and the specific steps include:
setting the population to update the weight vector and the direction vector once every iteration T times, namely if rem (h, T) is 0, turning to the step A, updating the weight vector and the direction vector, and then updating the individual operation in the step B; otherwise, directly turning to the step B to start updating the individuals;
step A, updating a weight vector and a direction vector: setting a threshold value zeta and a direction vector set W;
step A1, calculating individuals in the generated i' th class solution seti'XrR 1,2, …, N and its adjacent k solutions
Figure BDA00027356834300000713
The distance of (c):
Figure BDA00027356834300000714
wherein S is
Figure BDA0002735683430000081
Andi'Xrthe covariance matrix of (a);
step A2, if the calculated distance is larger than the threshold value, Dm(i'Xr) Zeta indicate individuali'XrThe other individuals in the dissociation set are far, the sparsity level is large, and isolated points are foundi'XrCorresponding direction vector
Figure BDA0002735683430000082
And weight vector
Figure BDA0002735683430000083
Continuing to operate in the step A3, otherwise, turning to the step B;
step A3, generating a new sub-question: according to the property of the Chebyshev decomposition model, when selecting a representative for the ith' class, a given objective function F (i'X) and a reference pointi'z*The optimal weight vector is a sub-problem of the Chebyshev decomposition
Figure BDA0002735683430000084
N is 1,2, …, N reaches a minimum value in order to allow the algorithm to converge to isolated pointsi'XrThe weight vector for constructing the subproblems is:
Figure BDA0002735683430000085
wherein the content of the first and second substances,
Figure BDA0002735683430000086
l＝1,2,…,L；
by usingi'XrGenerate new sub-questionsTitle:
Figure BDA0002735683430000087
step A4, obtaining
Figure BDA0002735683430000088
Corresponding direction vector
Figure BDA0002735683430000089
Figure BDA00027356834300000810
Step A5, calculating the direction vector by the following formula
Figure BDA00027356834300000811
Two nearest direction vectors
Figure BDA00027356834300000812
And find the corresponding individual
Figure BDA00027356834300000813
Figure BDA00027356834300000814
Wherein, Σ is
Figure BDA00027356834300000815
And
Figure BDA00027356834300000816
the covariance matrix of (a);
step A6, two individuals
Figure BDA00027356834300000817
Randomly generates a solution betweeni'XnewAs a new solution：
Figure BDA00027356834300000818
Step A7, usingi'XnewReplacement ofi'Xr；
Step A8, updating the direction vector:
Figure BDA00027356834300000819
instead of the former
Figure BDA00027356834300000820
Step a9, updating the weight vector:
Figure BDA00027356834300000821
instead of the former
Figure BDA00027356834300000822
Obtaining the updated weight vector
Figure BDA00027356834300000823
The most recent TT weight vectors are recorded into an index set
Figure BDA00027356834300000824
N is 1,2, …, N;
step B, updating the individuals, comprising the following steps:
step B1, 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 Pi'(n)：
Figure BDA0002735683430000091
Wherein, Bi'(n) is a weight vector index set;
step B2, 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 B3, calculating the adaptive value of each particle, and for each particle, matching the adaptive value with the best position p experiencedbestIf 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 gbestComparing the adaptive values, and if the adaptive values are better, taking the adaptive values as the current global best position;
step B4, 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:
vji(t+1)＝vid(t)+c1·r1·(Pji(t)-xji(t))+c2·r2·(NPji(t)-xji(t))
xji(t+1)＝xji(t)+vji(t+1)
wherein, c is1And c2Is a normal number, called a learning factor, c1Adjusting the step length of the particles flying to the direction of the best position of the particles; c. C2Adjusting the step length of the flight of the particles to the global best position; r is1And r2Is [0,1 ]]A random number in between; pji(t) are random individuals within the non-dominated solution set P _ set; NPji(t) are the particles within the dominance solution set NP _ set;
step C, updating the reference pointi'z (h): for each L ═ 1, 2., L, ifi'zl(h)＞fl(i'Y (h)) toi'zl(h)＝fl(i'Y(h))，i'zl(h) Represents the function f of the ith' class at the h evolutionhA reference point of (d);
step D, pair has n2Each particle in the NP _ subset of particles
Figure BDA0002735683430000092
k＝1,2,…,n2And has n1Each particle in a subset of particles P _ se
Figure BDA0002735683430000093
m＝1,2,…,n1Comparing one by one, calculating the fitness fit of the corresponding particles,i'Xkand all ofi'Xm，m＝1,2,…,n1After comparison, k is not present, so that fit: (i'Xm)≤fit(i'Xk) Or fit (i'Xm)≥fit(i'Xk) Is established, theni'XmMay also be a non-dominant solution, willi'XmThe P _ set subset is also added;
step E, updatei'EP (h): updating P _ set subset and its corresponding number of particles n1And the number of particles n to which the NP _ set subset corresponds2Adding particles in the P _ set subseti'EP (h); meanwhile, n is n +1, namely the next weight vector is processed;
step 728, loop determination: 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 responsei'EP(h)；
Step S729, fromi'Selecting a temperature transient thermal response optimization solution in EP (h)i'REP, representative of class i' transient responsesi'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 performing linear transformation on the two-dimensional matrix M by using the matrix Y to obtain the two-dimensional image matrix R (x, Y) includes the specific steps of: starting each frame in the three-dimensional matrix S 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 frame, 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:
Figure BDA0002735683430000101
a two-dimensional image matrix R (x, y) is obtained, wherein,
Figure BDA0002735683430000102
is a pseudo-inverse of the matrix Y, OTA transpose of the two-dimensional matrix O.
Preferably, the step nine of extracting features of the two-dimensional image matrix R (x, y) by using an improved one-dimensional Ostu segmentation algorithm to obtain defect features includes:
step S91, calculating the gray probability distribution P of the imaget：
Figure BDA0002735683430000103
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 imagetThe number of pixel points with the gray level of t in the image is counted;
step S92, calculating the target class C0And background class C1Probability of occurrence of each omega0And ω1And their corresponding mean values of gray levels mu0And mu1Thereby obtaining the spacing between the two types
Figure BDA0002735683430000104
And the degree of dispersion d of each class is determined0And d1：
Figure BDA0002735683430000105
Figure BDA0002735683430000106
Figure BDA0002735683430000107
Figure BDA0002735683430000108
Step S93, calculating a threshold value selection function h (t):
Figure BDA0002735683430000109
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, realizing segmentation of the gray level image, and obtaining the impact damage defect characteristics 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 decomposition multi-objective optimization method adopting the self-adaptive 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. After the population is set to evolve for a certain number of generations by the adaptive weight adjustment strategy, the periodic adjustment of the search direction is tried again under the condition that the solution set can show the approximate distribution of the front end of the PF, the efficiency of the whole search process is improved, a new subproblem is constructed based on the properties of the Chebyshev decomposition model, and the actual situation of the current multi-target requirement solution problem is better met.
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 an embodiment of a decomposition multi-objective optimization infrared thermal image defect feature extraction method based on adaptive 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 method for judging impact damage characteristic types of aerospace heat-proof materials, which comprises the following steps:
step one, representing a thermal image sequence of a test piece pixel point 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 Mmax,jmax,tmax) Wherein imaxIndicates the number of rows where the temperature peak point is located, jmaxIndicates the number of columns where the temperature peaks are located, tmaxRepresenting 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;
extracting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on self-adaptive weight adjustment to form a matrix Y;
step eight, changing the three-dimensional matrix M into a two-dimensional matrix, and performing linear transformation on the two-dimensional matrix by using a 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 image0And background C1And obtaining the impact damage defect characteristics of the test piece.
The invention makes up the defects of the traditional clustering method in setting the clustering number and selecting the initial clustering center, obtains the uniformly distributed optimal solution set by periodically adjusting the distribution of the direction vectors in the aspect of selecting representative temperature thermal response, can ensure the diversity of the solution set, and the selected characteristic points have the most representative characteristics so as to improve the accuracy of defect characteristic extraction. The invention relates to a decomposition multi-target optimization infrared thermal image defect feature extraction method based on self-adaptive weight adjustment, which comprises the steps of selecting transient thermal response of representative pixel points from a thermal image sequence conversion step length, clustering by adopting a mean shift clustering algorithm to obtain the category of the transient thermal response of each pixel point, obtaining a dimension reduction result of the thermal image sequence, considering the pixel value (temperature value) similarity of each category pixel point and the same category pixel points, simultaneously considering the difference between the pixel point (temperature point) and the different category pixel points (temperature points), constructing a corresponding multi-target function, obtaining a uniformly distributed non-dominated solution set which is closest to a real solution set by utilizing a decomposition multi-target evolution algorithm based on self-adaptive weight adjustment to ensure that the representative temperature points selected from each damage set are most representative, and finally utilizing an improved one-dimensional Ostu threshold segmentation algorithm to extract features, thereby extracting the defect characteristics of the infrared thermal image. The differences among different types of temperature points and the similarities among the same type of temperature points are comprehensively considered, and the uniformly distributed weight vectors cannot guarantee that a uniformly distributed Pareto optimal solution set can be obtained, so that the accurate selection of representative pixel points (temperature points) is realized, and the accuracy of defect feature extraction is guaranteed.
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: is provided withSet threshold value
Figure BDA0002735683430000121
Calculating the temperature peak point M (i)max,jmax,tmax) To the nearest temperature point Mc(imax,jc,tmax) Is calculated by the correlation coefficient PCreWherein j iscThe number of columns where the nearest temperature point is located, if
Figure BDA0002735683430000122
Then j isc＝jc+1 calculating the temperature peak M (i) from near to farmax,jmax,tmax) Correlation coefficient PC with the remaining temperature pointsreUp to
Figure BDA0002735683430000131
When it is time, stop calculating, remember
Figure BDA0002735683430000132
Number of established pixels as line step length CLrAccording to CLrAnd 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
Figure BDA0002735683430000133
Calculating the temperature peak point M (i)max,jmax,tmax) To the nearest temperature point Mc(ic,jmax,tmax) Is calculated by the correlation coefficient PCreWherein icThe number of rows where the nearest temperature point is located, if
Figure BDA0002735683430000134
Then ic＝ic+1 calculating the temperature peak M (i) from near to farmax,jmax,tmax) Correlation coefficient PC with the remaining temperature pointsreUp to
Figure BDA0002735683430000135
When it is time, stop calculating, remember
Figure BDA0002735683430000136
The number of established pixels is the column step length RLqAccording to RLqAnd 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, step four divide the three-dimensional matrix M into N data blocks, and represent the nth data block as MnAnd setting a threshold CC, wherein N is R · Q, N is 1, 2.., N;
step S52, the temperature peak point M (i)max,jmax,tmax) Corresponding transient thermal response M (i)max,jmaxT), 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 calculatedn(m, n, T), T ═ 1,2, …, correlation PC between T and X (: 1)re(ii) a If PCre＜CC，Mn(m, n, T), T ═ 1,2, …, T is regarded as a new feature and stored in X (: z); otherwise, let m be m + RLqContinuously 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 + CLrAnd 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 step six employs a mean shift algorithm, and the specific method for automatically classifying the typical transient thermal response extracted in the step five includes: initially 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 large difference into another class, updating the cluster number L to be L +1, after all transient thermal responses in a transient thermal response set X (: z) are subjected to iterative calculation, finally obtaining the clustered 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 11Center＝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-XmAll transient thermal responses X of | < rmWherein X ismE.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:
Figure BDA0002735683430000141
wherein, omega is a high-dimensional sphere area with the Center as the Center radius r; xmSatisfying [ Center-X ] for transient thermal response distributed in omega regionmL < r; k is transient thermal response X distributed in an omega regionmThe number of (2); g (x) is a logarithmic kernel function
Figure BDA0002735683430000142
The derivation is negative;
step S65, the Center updates along the direction of the shift mean, and the movement distance is | | shift (Center) |:
Figure BDA0002735683430000143
wherein the content of the first and second substances,tthe Center represents the updated cluster Center;
step S66,If | | | shift (Center) | | ≧ epsilon, let Center ═ epsilontCenter, t ═ t +1, return to step S63 and continue the iterative computation until | | shift (Center) | < epsilon, find out the previous iteration times t to distribute in order totAll transient thermal response sets in the high-dimensional sphere area with Center as the Center and radius r, namely satisfying | survivaltCenter-Xb||＜r，XbAll transient thermal responses X for e X (: z)bThese transient thermal responses are considered to belong to the same clustertNDs；
Step S67, calculating the updated cluster centertCenter and other already existing cluster centersi'Distance of Center, if anytCenter-i'If Center | < ω is true, then consider it to betCenter and the nearest one to itkThe Center belongs to the same class, wherein
Figure BDA0002735683430000144
Go to step S68; otherwise go to step S69;
step S68, collectingtTransient thermal response in NDs is ascribed tokCenter, k is 1,2, …, and L is a cluster of the cluster CenterkNDs, will be assembledtTransient thermal response markers in NDs are accessed and clusteredkThe number of occurrences in NDs is added to 1 and calculatedtCenter andkweighted average of Centert'Center，kNeutral separation of NDst'Recent transient thermal response of Centert”The Center is the cluster Center of the current cluster, i.e. the current cluster CenterkCenter＝t”Center, the number of clusters L equals L, let Center equal to Lt”The Center goes to step S610;
step S69, if nottCenter-i'Let i '═ i' +1, consider that a new cluster appears in the current iterationi'NDs of whichi'NDs＝tNDs，tThe Center is a new clusteri'Clustering centers of NDs, to be aggregatedtTransient thermal response markers in NDs are accessed and clusteredi'Adding 1 to the times of appearance in NDs, wherein the clustering number L is equal to L +1, and the current clustering centeri'Center＝tCenter, let Center ═tThe 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 percenti'Center, i ═ 1,2, …, L; the transient thermal response marked as visited in the set X (: z) goes to step S611 to classify the category of the transient thermal response; 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 occursp，XpWhen 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 clusteri'NDs, i ═ 1,2, …, L found free from XpK transient thermal responses with Euclidean distance being nearest are calculated to obtain XpToi'Crowdedness distance of k nearest transient thermal responses in NDs:
Figure BDA0002735683430000151
wherein the content of the first and second substances,
Figure BDA0002735683430000152
represents XpEuclidean distance to the j, j ═ 1,2, …, k nearest transient thermal responses in class i'; mixing XpIf the calculated congestion degree distance is the largest, the step goes to step S613 to output the transient thermal response division condition;
step S612, transient state not marked as accessedClassifying the state thermal response: calculating transient thermal response X not marked as visitedq，XqBelongs to X (: z) to each category clustering centeri'European distance from Center:
Figure BDA0002735683430000153
wherein the content of the first and second substances,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 XqThe step S613 is switched to output the transient thermal response classification condition;
step S613, outputting the transient thermal response division conditioni'NDs, i ═ 1,2, …, L, and the cluster centers for each classi'Center，i'NDs represents the set of transient thermal responses of class i'.
In the above technical solution, the seventh step of extracting a representative of each type of transient thermal response by using a decomposition multi-objective optimization algorithm based on adaptive weight adjustment, and the specific steps of forming the matrix Y include:
in step S71, when the i ', i' is 1, L classes of transient responses are represented, a multi-objective function is defined:
minimizeF(i'X)＝(f1(i'X),...,fL(i'X))T
wherein f is1(i'X) a transient response selected for the i' th class transient responsei'The intra-class Euclidean distance of X is expressed as:
Figure BDA0002735683430000154
wherein the content of the first and second substances,i'xhfor transient responsei'The pixel value of X at the h-th time, i.e. the temperature value,i'Centerhthe pixel value of the ith' type transient response cluster center at the h-th moment, namely the temperature value,j'Centerhclustering the pixel value of the j' th class transient response center at the h-th momentI.e. a 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 responsei'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
Figure BDA0002735683430000161
Wherein
Figure BDA0002735683430000162
Corresponding direction vector
Figure BDA0002735683430000163
Comprises the following steps:
Figure BDA0002735683430000164
wherein the content of the first and second substances,
Figure BDA0002735683430000165
step 722, decomposing the L multi-target questions into N sub-questions by using a Tchebycheff decomposition algorithm, wherein each sub-question is as follows:
Figure BDA0002735683430000166
wherein the content of the first and second substances,
Figure BDA0002735683430000167
i'z*for the reference point corresponding to the multi-objective function of the ith' class,
Figure BDA0002735683430000168
is a function fl(i'X) a corresponding reference point;
step S723. Initializing the evolution number h to 0 fromi'Randomly generating N values in Ω
Figure BDA0002735683430000169
Wherein the content of the first and second substances,i'Xn(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 seti'Initializing a reference point when EP (0) is phi;
step S725, initializing neighborhood vectors, and solving the nth weight vector of the ith class
Figure BDA00027356834300001610
The most recent TT weight vectors, and will be
Figure BDA00027356834300001611
The index set of the most recent TT weight vectors is marked as Bi'(n)＝{n1,n2,...,nTT}，n＝1,2,...,N；
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 the iteration when h is 1, terminating the iteration when h is P, and applying each weight vector
Figure BDA00027356834300001612
N is 1,2, …, and the updating operation is performed by N according to the iteration number h, and the specific steps include:
setting the population to update the weight vector and the direction vector once every iteration T times, namely if rem (h, T) is 0, turning to the step A, updating the weight vector and the direction vector, and then updating the individual operation in the step B; otherwise, directly turning to the step B to start updating the individuals;
step A, updating a weight vector and a direction vector: setting a threshold value zeta and a direction vector set W;
step A1, calculating individuals in the generated i' th class solution seti'XrR 1,2, …, N andadjacent k solutions
Figure BDA00027356834300001613
The distance of (c):
Figure BDA00027356834300001614
wherein S is
Figure BDA00027356834300001615
Andi'Xrthe covariance matrix of (a);
step A2, if the calculated distance is larger than the threshold value, Dm(i'Xr) Zeta indicate individuali'XrThe other individuals in the dissociation set are far, the sparsity level is large, and isolated points are foundi'XrCorresponding direction vector
Figure BDA0002735683430000171
And weight vector
Figure BDA0002735683430000172
Continuing to operate in the step A3, otherwise, turning to the step B;
step A3, generating a new sub-question: according to the property of the Chebyshev decomposition model, when selecting a representative for the ith' class, a given objective function F (i'X) and a reference pointi'z*The optimal weight vector is a sub-problem of the Chebyshev decomposition
Figure BDA0002735683430000173
N is 1,2, …, N reaches a minimum value in order to allow the algorithm to converge to isolated pointsi'XrThe weight vector for constructing the subproblems is:
Figure BDA0002735683430000174
wherein the content of the first and second substances,
Figure BDA0002735683430000175
l＝1,2,…,L；
by usingi'XrGenerating a new sub-problem:
Figure BDA0002735683430000176
step A4, obtaining
Figure BDA0002735683430000177
Corresponding direction vector
Figure BDA0002735683430000178
Figure BDA0002735683430000179
Step A5, calculating the direction vector by the following formula
Figure BDA00027356834300001710
Two nearest direction vectors
Figure BDA00027356834300001711
And find the corresponding individual
Figure BDA00027356834300001712
Figure BDA00027356834300001713
Wherein, Σ is
Figure BDA00027356834300001714
And
Figure BDA00027356834300001715
the covariance matrix of (a);
step A6, two individuals
Figure BDA00027356834300001716
Randomly generates a solution betweeni'XnewAs a new solution:
Figure BDA00027356834300001717
step A7, usingi'XnewReplacement ofi'Xr；
Step A8, updating the direction vector:
Figure BDA00027356834300001718
instead of the former
Figure BDA00027356834300001719
Step a9, updating the weight vector:
Figure BDA00027356834300001720
instead of the former
Figure BDA00027356834300001721
Obtaining the updated weight vector
Figure BDA00027356834300001722
The most recent TT weight vectors are recorded into an index set
Figure BDA00027356834300001723
N is 1,2, …, N;
step B, updating the individuals, comprising the following steps:
step B1, 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 Pi'(n)：
Figure BDA0002735683430000181
Wherein, Bi'(n) is a weight vector index set;
step B2, 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 B3, calculating the adaptive value of each particle, and for each particle, matching the adaptive value with the best position p experiencedbestIf 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 gbestComparing the adaptive values, and if the adaptive values are better, taking the adaptive values as the current global best position;
step B4, 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:
vji(t+1)＝vid(t)+c1·r1·(Pji(t)-xji(t))+c2·r2·(NPji(t)-xji(t))
xji(t+1)＝xji(t)+vji(t+1)
wherein, c is1And c2Is a normal number, called a learning factor, c1Adjusting the step length of the particles flying to the direction of the best position of the particles; c. C2Adjusting the step length of the flight of the particles to the global best position; r is1And r2Is [0,1 ]]A random number in between; pji(t) are random individuals within the non-dominated solution set P _ set; NPji(t) are the particles within the dominance solution set NP _ set;
step C, updating the reference pointi'z (h): for each L ═ 1, 2., L, ifi'zl(h)＞fl(i'Y (h)) toi'zl(h)＝fl(i'Y(h))，i'zl(h) Represents the function f of the ith' class at the h evolutionhA reference point of (d);
step D, pair has n2Each particle in the NP _ subset of particles
Figure BDA0002735683430000182
k＝1,2,…,n2And has n1Each particle in a subset of particles P _ se
Figure BDA0002735683430000183
m＝1,2,…,n1Comparing one by one, calculating the fitness fit of the corresponding particles,i'Xkand all ofi'Xm，m＝1,2,…,n1After comparison, k is not present, so that fit: (i'Xm)≤fit(i'Xk) Or fit (i'Xm)≥fit(i'Xk) Is established, theni'XmMay also be a non-dominant solution, willi'XmThe P _ set subset is also added;
step E, updatei'EP (h): updating P _ set subset and its corresponding number of particles n1And the number of particles n to which the NP _ set subset corresponds2Adding particles in the P _ set subseti'EP (h); meanwhile, n is n +1, namely the next weight vector is processed;
step 728, loop determination: 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 responsei'EP(h)；
Step S729, fromi'Selecting a temperature transient thermal response optimization solution in EP (h)i'REP, representative of class i' transient responsesi'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 transformation on the two-dimensional matrix by using the matrix Y to obtain a two-dimensional image matrix R (x, Y) specifically includes: starting each frame in the three-dimensional matrix S from a first column, connecting a next column to the tail of a previous column to form a new column to obtain T-column pixel values corresponding to the T frames, sequentially placing the T-column pixel values according to time sequence to form a two-dimensional matrix O, and using a matrix Y to form a two-dimensional matrix OThe two-dimensional matrix O is linearly transformed, i.e.:
Figure BDA0002735683430000191
a two-dimensional image matrix R (x, y) is obtained, wherein,
Figure BDA0002735683430000192
is a pseudo-inverse of the matrix Y, OTA transpose of the two-dimensional matrix O.
In the above technical solution, the ninth step of performing feature extraction on the two-dimensional image matrix R (x, y) by using an improved one-dimensional Ostu segmentation algorithm to obtain the defect feature specifically includes:
step S91, calculating the gray probability distribution P of the imaget：
Figure BDA0002735683430000193
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 imagetThe number of pixel points with the gray level of t in the image is counted;
step S92, calculating the target class C0And background class C1Probability of occurrence of each omega0And ω1And their corresponding mean values of gray levels mu0And mu1Thereby obtaining the spacing between the two types
Figure BDA0002735683430000194
And the degree of dispersion d of each class is determined0And d1：
Figure BDA0002735683430000195
Figure BDA0002735683430000196
Figure BDA0002735683430000197
Figure BDA0002735683430000198
Step S93, calculating a threshold value selection function h (t):
Figure BDA0002735683430000199
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, realizing segmentation of the gray level image, and obtaining the impact damage defect characteristics of the test piece.
Example (b):
in the present embodiment, there are two kinds of defects on the test piece, i.e., defect 1 filled with no material and defect 2 filled with a poor thermal conductive material.
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 asBacPOINT、Def1POINT andDef2POINT, as shown in fig. 8, 9, and 10.
Three transient thermal response representations are obtained by using the existing method for selecting the transient thermal response representation based on difference:ANK80、BNK138andcNK51the 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:ANK3、BNK156andcNK46the curves are shown in fig. 14, 15, and 17, 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 has the fastest rising rate and the highest temperature amplitude, and the rising rate and the temperature amplitude of the temperature point of the defect 1 are also higher than those of the material. Compared with the three characteristics, the temperature point of the defect 2 releases heat most quickly, and the temperature point of the defect 1 releases heat most slowly.
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.9894 0.9956 0.9996
The invention 0.9966 0.9932 0.9703
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 method for judging the impact damage characteristic type of a space heat-proof material is characterized by comprising the following steps:
step one, representing a thermal image sequence of a test piece pixel point 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 Mmax,jmax,tmax) Wherein imaxIndicates the number of rows where the temperature peak point is located, jmaxIndicates the number of columns where the temperature peaks are located, tmaxRepresenting 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;
extracting a representative of each type of transient thermal response by adopting a decomposition multi-objective optimization algorithm based on self-adaptive weight adjustment to form a matrix Y;
step eight, changing the three-dimensional matrix M into a two-dimensional matrix, and performing linear transformation on the two-dimensional matrix by using a 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 image0And background C1And obtaining the impact damage defect characteristics of the test piece.
2. The method for interpreting the impact damage feature types of the aerospace heat protection material as claimed in claim 1, wherein the step three is to divide the temperature of the row with the frame number of the temperature peak point, and the specific method for dividing the row with the frame number of the temperature peak point into R data blocks comprises the following steps: setting a threshold value
Figure FDA0002735683420000011
Calculating the temperature peak point M (i)max,jmax,tmax) To the nearest temperature point Mc(imax,jc,tmax) Is calculated by the correlation coefficient PCreWherein j iscThe number of columns where the nearest temperature point is located, if
Figure FDA0002735683420000021
Then j isc＝jc+1, calculating the temperature peak M (i) from near to farmax,jmax,tmax) Correlation coefficient PC with the remaining temperature pointsreUp to
Figure FDA0002735683420000022
When it is time, stop calculating, remember
Figure FDA0002735683420000023
Number of established pixels as line step length CLrAccording to CLrAnd 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 interpreting the impact damage characteristic types of the aerospace heat protection material as claimed in claim 1, wherein the fourth step is to perform temperature division on the frame number columns where the temperature peak points are located, and the specific method for dividing the frame number columns where the temperature peak points are located into Q data blocks comprises the following steps: setting a threshold value
Figure FDA0002735683420000024
Calculating the temperature peak point M (i)max,jmax,tmax) To the nearest temperature point Mc(ic,jmax,tmax) Is calculated by the correlation coefficient PCreWherein icThe number of rows where the nearest temperature point is located, if
Figure FDA0002735683420000025
Then ic＝ic+1 calculating the temperature peak M (i) from near to farmax,jmax,tmax) Correlation coefficient PC with the remaining temperature pointsreUp to
Figure FDA0002735683420000026
When it is time, stop calculating, remember
Figure FDA0002735683420000027
The number of established pixels is the column step length RLqAccording to RLqAnd 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 interpreting the impact damage characteristic types of the aerospace heat protection materials as claimed in claim 1, wherein the step of extracting the transient thermal response in five blocks and in steps comprises the following specific steps:
step S51, step three, step four divide the three-dimensional matrix M into N data blocks, and represent the nth data block as MnAnd setting a threshold CC, wherein N is R · Q, N is 1, 2.., N;
step S52, the temperature peak point M (i)max,jmax,tmax) Corresponding transient thermal response M (i)max,jmaxT), 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 calculatedn(m, n, T), T ═ 1,2, …, correlation PC between T and X (: 1)re(ii) a If PCre＜CC，Mn(m, n, T), T ═ 1,2, …, T is regarded as a new feature and stored in X (: z); otherwise, let m be m + RLqContinuously 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 + CLrAnd 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 interpreting the impact damage feature types of the aerospace heat protection materials as claimed in claim 1, wherein the step six employs a mean shift algorithm, and the method for automatically classifying the typical transient thermal responses extracted in the step five comprises the following steps: initially 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 large difference into another class, updating the cluster number L to be L +1, after all transient thermal responses in a transient thermal response set X (: z) are subjected to iterative calculation, finally obtaining the clustered 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 loop parameter T to 1, i' to 1, the loop termination condition T to T, and keeping the initial cluster when T to 1Center, note1Center＝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-XmAll transient thermal responses X of | < rmWherein X ismE.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:
Figure FDA0002735683420000031
wherein, omega is a high-dimensional sphere area with the Center as the Center radius r; xmSatisfying [ Center-X ] for transient thermal response distributed in omega regionmL < r; k is transient thermal response X distributed in an omega regionmThe number of (2); g (x) is a logarithmic kernel function
Figure FDA0002735683420000032
The derivation is negative;
step S65, the Center updates along the direction of the shift mean, and the movement distance is | | shift (Center) |:
Figure FDA0002735683420000033
wherein the content of the first and second substances,tthe Center represents the updated cluster Center;
step S66, if | | | shift (Center) | | ≧ epsilon, let Center ═ epsilontCenter, t ═ t +1, return to step S63 and continue the iterative computation until | | shift (Center) | < epsilon, find out the previous iteration times t to distribute in order totAll transient thermal response sets in the high-dimensional sphere area with Center as the Center and radius r, namely satisfying | survivaltCenter-Xb||＜r，XbAll transient thermal responses X for e X (: z)bThese transient thermal responses are considered to belong to the same clustertNDs；
Step S67, calculating the updated cluster centertCenter and other already existing cluster centersi'Distance of Center, if anytCenter-i'If Center | < ω is true, then consider it to betCenter and the nearest one to itkThe Center belongs to the same class, wherein
Figure FDA0002735683420000042
Go to step S68; otherwise go to step S69;
step S68, collectingtTransient thermal response in NDs is ascribed tokCenter, k is 1,2, …, and L is a cluster of the cluster CenterkNDs, will be assembledtTransient thermal response markers in NDs are accessed and clusteredkThe number of occurrences in NDs is added to 1 and calculatedtCenter andkweighted average of Centert'Center，kNeutral separation of NDst'Recent transient thermal response of Centert”The Center is the cluster Center of the current cluster, i.e. the current cluster CenterkCenter＝t”Center, the number of clusters L equals L, let Center equal to Lt”The Center goes to step S610;
step S69, if nottCenter-i'Let i '═ i' +1, consider that a new cluster appears in the current iterationi'NDs of whichi'NDs＝tNDs，tThe Center is a new clusteri'Clustering centers of NDs, to be aggregatedtTransient thermal response markers in NDs are accessed and clusteredi'Adding 1 to the times of appearance in NDs, wherein the clustering number L is equal to L +1, and the current clustering centeri'Center＝tCenter, let Center ═tThe 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 percenti'Center，i'＝1,2,…,L；The transient thermal response marked as visited in the set X (: z) goes to step S611 to classify the category of the transient thermal response; 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 occursp，XpWhen 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 clusteri'NDs, i ═ 1,2, …, L found free from XpK transient thermal responses with Euclidean distance being nearest are calculated to obtain XpToi'Crowdedness distance of k nearest transient thermal responses in NDs:
Figure FDA0002735683420000041
wherein the content of the first and second substances,
Figure FDA0002735683420000043
represents XpEuclidean distance to the j, j ═ 1,2, …, k nearest transient thermal responses in class i'; mixing XpIf 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 visitedq，XqBelongs to X (: z) to each category clustering centeri'European distance from Center:
Figure FDA0002735683420000051
wherein the content of the first and second substances,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 XqThe step S613 is switched to output the transient thermal response classification condition;
step S613, outputting the transient thermal response division conditioni'NDs, i ═ 1,2, …, L, and the cluster centers for each classi'Center，i'NDs represents the set of transient thermal responses of class i'.
6. The interpretation method for impact damage feature types of aerospace heat protection materials as claimed in claim 1, wherein the seventh step adopts a decomposition multi-objective optimization algorithm based on adaptive weight adjustment to extract the representation of each type of transient thermal response, and the specific step of forming the matrix Y comprises:
in step S71, when the i ', i' is 1, L classes of transient responses are represented, a multi-objective function is defined:
minimize F(i'X)＝(f1(i'X),...,fL(i'X))T
wherein f is1(i'X) a transient response selected for the i' th class transient responsei'The intra-class Euclidean distance of X is expressed as:
Figure FDA0002735683420000052
wherein the content of the first and second substances,i'xhfor transient responsei'The pixel value of X at the h-th time, i.e. the temperature value,i'Centerhthe pixel value of the ith' type transient response cluster center at the h-th moment, namely the temperature value,j'Centerhthe 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 i' th class of transient stateRepresentation of the responsei'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
Figure FDA0002735683420000053
Wherein
Figure FDA0002735683420000054
Figure FDA0002735683420000055
Corresponding direction vector
Figure FDA0002735683420000056
Comprises the following steps:
Figure FDA0002735683420000061
wherein the content of the first and second substances,
Figure FDA0002735683420000064
step 722, decomposing the L multi-target questions into N sub-questions by using a Tchebycheff decomposition algorithm, wherein each sub-question is as follows:
Figure FDA0002735683420000062
wherein the content of the first and second substances,
Figure FDA0002735683420000065
i'z*for the reference point corresponding to the multi-objective function of the ith' class,
Figure FDA0002735683420000066
is a function fl(i'X) a corresponding reference point;
step S723, initializing the evolution time h to 0, to obtaini'Randomly generating N values in Ω
Figure FDA0002735683420000067
Wherein the content of the first and second substances,i'Xn(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 seti'Initializing a reference point when EP (0) is phi;
step S725, initializing neighborhood vectors, and solving the nth weight vector of the ith class
Figure FDA0002735683420000068
The most recent TT weight vectors, and will be
Figure FDA0002735683420000069
The index set of the most recent TT weight vectors is marked as Bi'(n)＝{n1,n2,...,nTT}，n＝1,2,...,N；
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 the iteration when h is 1, terminating the iteration when h is P, and applying each weight vector
Figure FDA00027356834200000610
Updating according to the iteration number h, and specifically comprising the following steps:
setting the population to update the weight vector and the direction vector once every iteration T times, namely if rem (h, T) is 0, turning to the step A, updating the weight vector and the direction vector, and then updating the individual operation in the step B; otherwise, directly turning to the step B to start updating the individuals;
step A, updating a weight vector and a direction vector: setting a threshold value zeta and a direction vector set W;
step A1, calculating individuals in the generated i' th class solution seti'XrR 1,2, …, N and its adjacent k solutions
Figure FDA00027356834200000611
The distance of (c):
Figure FDA0002735683420000063
wherein S is
Figure FDA00027356834200000612
Andi'Xrthe covariance matrix of (a);
step A2, if the calculated distance is larger than the threshold value, Dm(i'Xr) Zeta indicate individuali'XrThe other individuals in the dissociation set are far, the sparsity level is large, and isolated points are foundi'XrCorresponding direction vector
Figure FDA0002735683420000076
And weight vector
Figure FDA0002735683420000077
Continuing to operate in the step A3, otherwise, turning to the step B;
step A3, generating a new sub-question: according to the property of the Chebyshev decomposition model, when selecting a representative for the ith' class, a given objective function F (i'X) and a reference pointi'z*The optimal weight vector is a sub-problem of the Chebyshev decomposition
Figure FDA0002735683420000078
To a minimum value in order to allow the algorithm to converge to isolated pointsi'XrThe weight vector for constructing the subproblems is:
Figure FDA0002735683420000071
wherein the content of the first and second substances,
Figure FDA0002735683420000072
by usingi'XrGenerating a new sub-problem:
Figure FDA0002735683420000073
step A4, obtaining
Figure FDA0002735683420000079
Corresponding direction vector
Figure FDA00027356834200000710
Figure FDA0002735683420000074
Step A5, calculating the direction vector by the following formula
Figure FDA00027356834200000711
Two nearest direction vectors
Figure FDA00027356834200000712
And find the corresponding individual
Figure FDA00027356834200000713
Figure FDA0002735683420000075
Wherein, Σ is
Figure FDA00027356834200000714
And
Figure FDA00027356834200000715
the covariance matrix of (a);
step A6, two individuals
Figure FDA00027356834200000716
Randomly generates a solution betweeni'XnewAs a new solution:
Figure FDA00027356834200000717
step A7, usingi'XnewReplacement ofi'Xr；
Step A8, updating the direction vector:
Figure FDA00027356834200000718
instead of the former
Figure FDA00027356834200000719
Step a9, updating the weight vector:
Figure FDA00027356834200000720
instead of the former
Figure FDA00027356834200000721
Obtaining the updated weight vector
Figure FDA00027356834200000722
The most recent TT weight vectors are recorded into an index set
Figure FDA00027356834200000723
Performing the following steps;
step B, updating the individuals, comprising the following steps:
step B1, setting the scope of evolution: setting a threshold value delta according to a random number rand epsilon (0,1) generated randomlyDefining a neighborhood region Pi'(n)：
Figure FDA0002735683420000081
Wherein, Bi'(n) is a weight vector index set;
step B2, 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 B3, calculating the adaptive value of each particle, and for each particle, matching the adaptive value with the best position p experiencedbestIf 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 gbestComparing the adaptive values, and if the adaptive values are better, taking the adaptive values as the current global best position;
step B4, 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:
vji(t+1)＝vid(t)+c1·r1·(Pji(t)-xji(t))+c2·r2·(NPji(t)-xji(t))
xji(t+1)＝xji(t)+vji(t+1)
wherein, c is1And c2Is a normal number, called a learning factor, c1Adjusting the step length of the particles flying to the direction of the best position of the particles; c. C2Adjusting the step length of the flight of the particles to the global best position; r is1And r2Is [0,1 ]]A random number in between; pji(t) are random individuals within the non-dominated solution set P _ set; NPji(t) are the particles within the dominance solution set NP _ set;
step C, updating the reference pointi'z (h): for each L ═ 1, 2., L, ifi'zl(h)＞fl(i'Y (h)) toi'zl(h)＝fl(i'Y(h))，i'zl(h) Represents the function f of the ith' class at the h evolutionhA reference point of (d);
step D, pair has n2Each particle in the NP _ set subset of particles
Figure FDA0002735683420000082
And has n1Each particle in the P _ set subset of particles
Figure FDA0002735683420000083
Comparing one by one, calculating the fitness fit of the corresponding particles,i'Xkand all ofi'Xm，m＝1,2,…,n1After comparison, k is not present, so that fit: (i'Xm)≤fit(i'Xk) Or fit (i'Xm)≥fit(i'Xk) Is established, theni'XmMay also be a non-dominant solution, willi'XmThe P _ set subset is also added;
step E, updatei'EP (h): updating P _ set subset and its corresponding number of particles n1And the number of particles n to which the NP _ set subset corresponds2Adding particles in the P _ set subseti'EP (h); meanwhile, n is n +1, namely the next weight vector is processed;
step 728, loop determination: 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 responsei'EP(h)；
Step S729, fromi'Selecting a temperature transient thermal response optimization solution in EP (h)i'REP, representative of class i' transient responsesi'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. Aerospace thermal protection material impact damage as defined in claim 1The feature type interpretation method is characterized in that the step eight of converting the three-dimensional matrix M into a two-dimensional matrix and performing linear transformation on the two-dimensional matrix by using the matrix Y to obtain a two-dimensional image matrix R (x, Y) comprises the following specific steps of: starting each frame in the three-dimensional matrix S 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 frame, 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:
Figure FDA0002735683420000094
a two-dimensional image matrix R (x, y) is obtained, wherein,
Figure FDA0002735683420000095
is a pseudo-inverse of the matrix Y, OTA transpose of the two-dimensional matrix O.
8. The method for interpreting the impact damage feature types of the aerospace heat-shielding material as claimed in claim 1, wherein the step nine of extracting the features of the two-dimensional image matrix R (x, y) by using the improved one-dimensional Ostu segmentation algorithm to obtain the defect features comprises the following specific steps:
step S91, calculating the gray probability distribution P of the imaget：
Figure FDA0002735683420000091
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 imagetThe number of pixel points with the gray level of t in the image is counted;
step S92, calculating the target class C0And background class C1Probability of occurrence of each omega0And ω1And their corresponding mean values of gray levels mu0And mu1Thereby obtaining the spacing between the two types
Figure FDA0002735683420000096
And the degree of dispersion d of each class is determined0And d1：
Figure FDA0002735683420000092
Figure FDA0002735683420000093
Figure FDA0002735683420000101
Figure FDA0002735683420000102
Step S93, calculating a threshold value selection function h (t):
Figure FDA0002735683420000103
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, realizing segmentation of the gray level image, and obtaining the impact damage defect characteristics of the test piece.
CN202011132748.9A 2020-10-21 2020-10-21 Method for judging impact damage characteristic types of aerospace heat-proof materials Active CN112215830B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011132748.9A CN112215830B (en) 2020-10-21 2020-10-21 Method for judging impact damage characteristic types of aerospace heat-proof materials

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011132748.9A CN112215830B (en) 2020-10-21 2020-10-21 Method for judging impact damage characteristic types of aerospace heat-proof materials

Publications (2)

Publication Number Publication Date
CN112215830A true CN112215830A (en) 2021-01-12
CN112215830B CN112215830B (en) 2022-03-04

Family

ID=74056325

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011132748.9A Active CN112215830B (en) 2020-10-21 2020-10-21 Method for judging impact damage characteristic types of aerospace heat-proof materials

Country Status (1)

Country Link
CN (1) CN112215830B (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112884716A (en) * 2021-01-28 2021-06-01 中国空气动力研究与发展中心超高速空气动力研究所 Method for strengthening characteristics of ultra-high-speed impact damage area
CN112906713A (en) * 2021-01-28 2021-06-04 中国空气动力研究与发展中心超高速空气动力研究所 Aerospace composite material damage visualization feature extraction method
CN113295732A (en) * 2021-04-22 2021-08-24 杭州申昊科技股份有限公司 Pipeline robot capable of detecting pipeline defects and control method and control system thereof
CN117351005A (en) * 2023-12-01 2024-01-05 四川纳拓新材料技术有限公司 Method and system for detecting coating defects of carbon-coated foil

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170294013A1 (en) * 2016-04-11 2017-10-12 North Carolina Agricultural And Technical State University Normalized Defect Characterization of Pulse Thermographic Nondestructive Evaluation
CN107256322A (en) * 2017-08-17 2017-10-17 北京航空航天大学 A kind of composite laminated plate delamination damage recognition methods that index is merged based on high sensitivity
CN108830839A (en) * 2018-05-29 2018-11-16 电子科技大学 A kind of thermal image defect inspection method of the pressure vessel based on the segmentation of ranks variable step
CN109767438A (en) * 2019-01-09 2019-05-17 电子科技大学 A kind of thermal-induced imagery defect characteristic recognition methods based on dynamic multi-objective optimization
US20190228517A1 (en) * 2018-05-29 2019-07-25 University Of Electronic Science And Technology Of China Method for separating out a defect image from a thermogram sequence based on feature extraction and multi-objective optimization
CN110211103A (en) * 2019-05-23 2019-09-06 电子科技大学 Comentropy additivity based on infrared thermal imaging obscures defect characteristic and analyzes reconstructing method
CN110294147A (en) * 2019-05-07 2019-10-01 中国空气动力研究与发展中心超高速空气动力研究所 A kind of protection of space debris configuration damping screen method for estimating damage

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170294013A1 (en) * 2016-04-11 2017-10-12 North Carolina Agricultural And Technical State University Normalized Defect Characterization of Pulse Thermographic Nondestructive Evaluation
CN107256322A (en) * 2017-08-17 2017-10-17 北京航空航天大学 A kind of composite laminated plate delamination damage recognition methods that index is merged based on high sensitivity
CN108830839A (en) * 2018-05-29 2018-11-16 电子科技大学 A kind of thermal image defect inspection method of the pressure vessel based on the segmentation of ranks variable step
US20190228517A1 (en) * 2018-05-29 2019-07-25 University Of Electronic Science And Technology Of China Method for separating out a defect image from a thermogram sequence based on feature extraction and multi-objective optimization
CN109767438A (en) * 2019-01-09 2019-05-17 电子科技大学 A kind of thermal-induced imagery defect characteristic recognition methods based on dynamic multi-objective optimization
CN110294147A (en) * 2019-05-07 2019-10-01 中国空气动力研究与发展中心超高速空气动力研究所 A kind of protection of space debris configuration damping screen method for estimating damage
CN110211103A (en) * 2019-05-23 2019-09-06 电子科技大学 Comentropy additivity based on infrared thermal imaging obscures defect characteristic and analyzes reconstructing method

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
HAONAN ZHANG等: "《Design of Hypervelocity-Impact Damage Evaluation Technique Based on Bayesian Classifier of Transient Temperature Attributes》", 《IEEE ACCESS》 *
TING XUE等: "《A New Spacecraft Impact Damage Feature Extraction Algorithm Based on Dynamic Multi-Objective Optimization Method》", 《2019 IEEE 58TH CONFERENCE ON DECISION AND CONTROL (CDC)》 *
田建辉: "《基于混合数值法的功能梯度材料板瞬态热响应分析》", 《固体力学学报》 *
黄雪刚: "《B_4C-Al基复合材料对空间碎片超高速撞击的防护应用研究》", 《稀有金属材料与工程》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112884716A (en) * 2021-01-28 2021-06-01 中国空气动力研究与发展中心超高速空气动力研究所 Method for strengthening characteristics of ultra-high-speed impact damage area
CN112906713A (en) * 2021-01-28 2021-06-04 中国空气动力研究与发展中心超高速空气动力研究所 Aerospace composite material damage visualization feature extraction method
CN112906713B (en) * 2021-01-28 2022-03-04 中国空气动力研究与发展中心超高速空气动力研究所 Aerospace composite material damage visualization feature extraction method
CN112884716B (en) * 2021-01-28 2022-03-18 中国空气动力研究与发展中心超高速空气动力研究所 Method for strengthening characteristics of ultra-high-speed impact damage area
CN113295732A (en) * 2021-04-22 2021-08-24 杭州申昊科技股份有限公司 Pipeline robot capable of detecting pipeline defects and control method and control system thereof
CN117351005A (en) * 2023-12-01 2024-01-05 四川纳拓新材料技术有限公司 Method and system for detecting coating defects of carbon-coated foil
CN117351005B (en) * 2023-12-01 2024-02-06 四川纳拓新材料技术有限公司 Method and system for detecting coating defects of carbon-coated foil

Also Published As

Publication number Publication date
CN112215830B (en) 2022-03-04

Similar Documents

Publication Publication Date Title
CN112215830B (en) Method for judging impact damage characteristic types of aerospace heat-proof materials
CN112233099B (en) Reusable spacecraft surface impact damage characteristic identification method
CN109086700B (en) Radar one-dimensional range profile target identification method based on deep convolutional neural network
CN107784320B (en) Method for identifying radar one-dimensional range profile target based on convolution support vector machine
CN110348399B (en) Hyperspectral intelligent classification method based on prototype learning mechanism and multidimensional residual error network
CN109190665B (en) Universal image classification method and device based on semi-supervised generation countermeasure network
CN109767438B (en) Infrared thermal image defect feature identification method based on dynamic multi-objective optimization
CN110298280B (en) Ocean vortex identification method based on MKL multi-feature fusion
CN109063724B (en) Enhanced generation type countermeasure network and target sample identification method
Hu et al. Deep convolutional neural networks for hyperspectral image classification
CN107103338B (en) SAR target recognition method integrating convolution features and integrated ultralimit learning machine
CN105913081B (en) SAR image classification method based on improved PCAnet
CN111914728B (en) Hyperspectral remote sensing image semi-supervised classification method and device and storage medium
CN104899549A (en) SAR target recognition method based on range profile time-frequency image identification dictionary learning
CN107944483B (en) Multispectral image classification method based on dual-channel DCGAN and feature fusion
CN108764310B (en) SAR target recognition method based on multi-scale multi-feature depth forest
CN102930300B (en) Method and system for identifying airplane target
CN104732244A (en) Wavelet transform, multi-strategy PSO (particle swarm optimization) and SVM (support vector machine) integrated based remote sensing image classification method
CN107528824B (en) Deep belief network intrusion detection method based on two-dimensional sparsification
CN112819775A (en) Segmentation and reinforcement method for damage detection image of aerospace composite material
CN107578063B (en) Image Spectral Clustering based on fast selecting landmark point
CN113095417B (en) SAR target recognition method based on fusion graph convolution and convolution neural network
CN114118362A (en) Method for detecting structural surface crack under small sample based on generating type countermeasure network
CN110070485A (en) A kind of high-spectrum image dimensionality reduction method
CN112381108A (en) Bullet trace similarity recognition method and system based on graph convolution neural network deep learning

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant