CN112016627B - 一种在轨航天器微小撞击损伤可视化检测评估方法 - Google Patents

一种在轨航天器微小撞击损伤可视化检测评估方法 Download PDF

Info

Publication number
CN112016627B
CN112016627B CN202010920381.0A CN202010920381A CN112016627B CN 112016627 B CN112016627 B CN 112016627B CN 202010920381 A CN202010920381 A CN 202010920381A CN 112016627 B CN112016627 B CN 112016627B
Authority
CN
China
Prior art keywords
transient thermal
thermal response
matrix
column
dimensional
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.)
Active
Application number
CN202010920381.0A
Other languages
English (en)
Other versions
CN112016627A (zh
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 CN202010920381.0A priority Critical patent/CN112016627B/zh
Publication of CN112016627A publication Critical patent/CN112016627A/zh
Application granted granted Critical
Publication of CN112016627B publication Critical patent/CN112016627B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • 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
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Mathematical Physics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Biophysics (AREA)
  • Molecular Biology (AREA)
  • Computational Linguistics (AREA)
  • Biomedical Technology (AREA)
  • Computing Systems (AREA)
  • Software Systems (AREA)
  • Multimedia (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

本发明公开了一种在轨航天器微小撞击损伤可视化检测评估方法，包括以下步骤：将红外热像仪获取的航天器撞击损伤试件的热图像序列用三维矩阵表示；从三维矩阵中选出上升率最小的瞬态热响应对应的像素点；基于最小上升率瞬态热响应确定分块大小尺寸；进行分块并确定块内搜索步长；分块分步长选取瞬态热响应；对选取的瞬态热响应进行分类；基于SC‑SD谱聚类算法实现动态预测并进行多目标优化选取每类瞬态热响应的代表；基于斯皮尔曼相关系数进行实际缺陷类别数判别，并根据实际缺陷类别数矩阵变化得到二维图像；利用区域卷积神经网络R‑CNN对二维图像进行特征提取，得到航天器撞击损伤试件的缺陷特征图像，从而实现对微小撞击损伤的可视化检测评估。

Description

一种在轨航天器微小撞击损伤可视化检测评估方法
技术领域
本发明属于航天器空间碎片撞击损伤检测评估技术领域，更具体地说，本发明涉及一种在轨航天器微小撞击损伤可视化检测评估方法。
背景技术
随着人类空间活动的不断开展，日益增多的空间碎片给人类的航天活动和在轨航天器的安全造成了极大的威胁，空间碎片问题已成为人类空间开发和实践中逐渐发展和严重的一个现实问题。特别是地球轨道上存在的数量巨大的厘米级以下微小空间碎片，由于无法进行有效地监测预警和主动规避，已经成为威胁在轨运行航天器与航天员的潜在杀手，由微小空间碎片撞击导致的航天事故屡见不鲜。当各类航天器遭受空间碎片撞击时，如何实现对撞击损伤信息的采集获取和分析评估是非常重要的。考虑到空间碎片撞击事件是在一个动态环境中随机发生的，长期在轨航天器在服役期间通常会遭受不定期多次撞击，其撞击损伤数量、损伤程度、损伤位置和损伤类型等都具有不可预知性。因此，对微小空间碎片撞击损伤进行快速而精准的在轨检测、识别和分析，可以方便在太空的航天员和地面工作人员及时做出正确的操作和采取必要的措施，为航天器空间碎片撞击事件的在轨评估与决策提供重要技术支持，从而保障航天员和航天器的安全以及航天任务的顺利完成。
红外热成像装置及其相关技术广泛应用于航空航天领域，并在航天器损伤检测与评估方面发挥了重要作用。基于红外热成像原理的损伤检测评估技术是以红外辐射特性为基础，利用不同结构或材料的热辐射物理特性不同，对材料表面及内部的不均匀性或异常进行检测，其具有快速、非接触、无污染、单次检测面积大、结果可视化、适用材料种类广等优点，非常适合于对空间碎片撞击造成的复杂损伤进行原位在役检测。基于红外热成像装置采集被测对象在外部热激励环境下的表面温度场变化数据(既红外热图像序列)，可以获得被测对象不同损伤区域在空间和时间维度上的瞬态热响应信息，进而利用相应的特征提取处理算法可实现碎空间碎片撞击造成的复杂损伤缺陷的可视化检测和评估。同时，考虑到长期在轨航天器会多次遭受微小空间碎片撞击，因此通过对定期采集的红外热图像序列数据进行分析处理，可以实现对微小空间碎片撞击损伤进行在轨监测。可见，无论是对于微小空间碎片撞击损伤的检测评估，还是对于撞击损伤变化过程的在轨监测而言，如何从大量的红外热图像序列数据中自动、快速、精确地提取和分离出损伤特征信息都是至关重要的。实践表明，在处理红外热图像序列数据过程中，利用合适的多目标优化方法可以提高数据分析处理的效率和精度，从而实现航天器在空间碎片环境下产生微小撞击损伤的可视化评估方法。
在2018年11月30日申请的名称为“一种基于动态多目标优化的红外热图像缺陷特征提取方法”(公布号为201811451744.X)和“基于均匀进化的多目标优化红外热图像缺陷特征提取方法”(公布号为201811451866.9)的中国发明专利申请中，当多目标优化问题的环境发生了变化时，采用了基于预测的多方向预测策略对每一次环境改变后的理想PS位置进行了预测，在预测过程中，算法在m-1和m次外部环境时分别存储能够充分描述m-1和m次时刻的PS，
Figure BDA0002666530010000021
和
Figure BDA0002666530010000022
的形状和多样性的W个代表瞬态热响应构成m-1和m次外部环境的多向预测集
Figure BDA0002666530010000023
和
Figure BDA0002666530010000024
在m+1次多目标环境，利用m-1和m次的多向预测集
Figure BDA0002666530010000025
和
Figure BDA0002666530010000026
预估出新环境的
Figure BDA0002666530010000027
点进化轨迹，并在预测的PS附近生成新的初始种群解以加速新环境下的多目标优化算法收敛，从而提高动态多目标优化算法的运行效率。但在
Figure BDA0002666530010000028
中寻找构成m-1次外部环境的多向预测集
Figure BDA0002666530010000029
和在
Figure BDA00026665300100000210
中寻找构成m次外部环境的多向预测集
Figure BDA00026665300100000211
的W个多向预测集元素时，它利用
Figure BDA00026665300100000212
中的瞬态热响应与多向预测集
Figure BDA00026665300100000213
中的初始元素的距离，
Figure DA00026665300132302986
中各瞬态热响应与多向预测集
Figure DA00026665300132308507
中的初始元素的距离作为选择新增的多向预测集元素的度量尺度，但实际应用中单一的用距离作为度量尺度，新增的代表瞬态热响应不能充分描述
Figure BDA0002666530010000031
和
Figure BDA0002666530010000032
的形状和多样性信息，从而影响预测的精准度，故需要增强多向预测集和
Figure BDA0002666530010000033
代表
Figure BDA0002666530010000034
和
Figure BDA0002666530010000035
的能力。本发明在该基础上，采取完全自适应的SC-SD谱聚类算法，基于谱图理论，把所有的PS数据看做空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的，将聚类问题转化为图划分问题，避免了利用单个瞬态热响应点与其他点的距离大小去寻找每次环境下PS形状和多样性的过程，计算量大大减少。利用谱聚类操作，发现每次环境中PS解集中瞬态热响应数据的真实分布情况，基于PS数据的真实分布情况找到能够充分描述当前PS的形状和多样性的瞬态热响应数据作为代表瞬态热响应。并且基于样本邻域标准差的完全自适应的SC-SD谱聚类算法能够自适应地确定样本的尺度参数，消除人为给定参数的主观性，增加聚类结果稳定性。
同时，在多目标优化问题的分解方面，它采取了切比雪夫分解方法，使得各个瞬态热响应类别的前沿近似解集
Figure BDA0002666530010000036
中的解受权重向量的指导沿着一定方向朝实际前沿解集PF进化。但切比雪夫分解方法在面对高于二维的多目标优化问题时难以得到均匀的近似前沿面解集
Figure BDA0002666530010000037
以权重向量λ＝(λ1,λ2,λ3)T为例，由于切比雪夫聚合函数形式本身的原因，指导属于权重向量λ＝(λ1,λ2，λ3)T的解的进化方向为方向向量为λ'＝(1/λ1,1/λ2,1/λ3)T的一条直线。由于解的进化方向不是沿着权重向量本身所在的直线进化，因此就算权重向量是均匀的，得到的解也是不均匀的，造成获得各个瞬态热响应类别的前沿近似解集
Figure BDA0002666530010000038
中解的疏密程度不一，对于航天航空领域多目标环境下的实际前沿解集PF上本身解集稀疏的区域可能收敛不到最优解，也就无法准确找到表征缺陷信息的瞬态热响应，造成缺陷检测失败和漏检的情况出现，而对于本身解过于密集的区域，其实际解集的多样性有限，故只需要用少量的近似前沿解就能够描述此区域的解集多样性情况，若仍然采用相同数量的均匀的权重向量去指导解的进化，容易重复找到表征同一类缺陷信息的瞬态热响应，不仅会造成资源的浪费，还因为解集密集和解集稀疏的区域实际上完成收敛所需的迭代次数不一致，造成算法整体时间增加，使得对动态多目标环境的响应变慢。本发明采取基于惩罚项的边界交叉法进行多目标优化问题的分解，从聚合函数本身入手，将聚合函数的形式改进为解与权重向量λ＝(λ1,λ2,λ3)T的垂直距离和沿着权重向量的平行方向与参考点的距离两者的加权和，从分解形式上将解的进化方向限制在权重向量本身身上，使得多目标优化算法的解沿着权重向量的方向进化，使得在面对航空航天领域中存在的大量的高于二维的多目标优化问题时也能够得到分布均匀的前沿近似解集
Figure BDA0002666530010000041
保证了在进化过程中位于各类瞬态热响应实际前沿PF上本身最优解集稀疏区域的解也能够进化到最优解所在的位置，提高检测准确度，减少得到进化到实际前沿最优解集密集区域所需的权重向量，避免资源浪费，提升算法速度。并引入惩罚因子，平衡种群解之间的收敛性与多样性，减少算法整体运行时间，从而应对航天器的复杂多目标优化问题环境。
此外，原专利中针对红外热图像数据块，基于整体温度最大值点进行了小尺寸的行、列块分割，造成所需处理的数据块个数较多，面对航空航天过程中需要频繁和实时对目标进行缺陷检测的需求，小尺寸分块的计算量大，处理效率慢，难以对缺陷产生后进行及时、迅速的发现和检测，同时，过度的数据块分割会造成表征缺陷部分的瞬态热响应数据的搜素遗漏和误去除，影响缺陷检测精度。本发明改进了针对红外热图像数据块的分块操作部分，首先基于瞬态热响应上升率找到三维数据块中代表背景区域的部分，根据背景区域的尺寸将数据块分成较大尺寸的数据块，在不同数据块中采取不同的搜索步长进行数据搜索，有效避免了基于温度最大值进行小尺寸分割块所造成的过度去冗余现象，提高数据处理速度和检测精确度。原专利中根据前沿最优解集
Figure BDA0002666530010000042
选择得到各类的瞬态热响应代表时，采取随机选择的方式，使得算法检测精度和准确性不高。本发明提出从前沿近似解集
Figure BDA0002666530010000043
中选出第i'类瞬态热响应的代表i'REP的基于加权隶属度的具体方案，避免随机选择带来的不确定性。原专利中在设定缺陷类别总数时采取固定缺陷类别数然后从红外数据中寻找各类别缺陷对应的热数据进行检测，造成缺陷类别的错误识别和漏检测。本发明提出基于斯皮尔曼相关系数的缺陷类别数判别方法，避免固定缺陷类别数检测带来的错检和漏检现象，提高检测精确度。原专利中表征动态多目标环境变化程度的指标采取了简单的算术平均值，造成对极端的函数值变化过于敏感。本发明改进动态多目标环境下每次环境变化剧烈程度的判别指标公式，避免环境变化后极端函数值对判别值产生过大的影响，提升应对动态环境变化的响应速度。
发明内容
本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。
为了实现根据本发明的这些目的和其它优点，提供了一种在轨航天器微小撞击损伤可视化检测评估方法，包括以下步骤：
步骤一、将红外热像仪获取的航天器撞击损伤试件的热图像序列用三维矩阵S表示，其中三维矩阵S的元素S(i,j,t)表示热图像序列的t帧热图像的第i行、第j列的像素值；
步骤二、从三维矩阵S的所有瞬态热响应中，选取上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)，其中，izz、jzz分别表示上升率最小的瞬态热响应对应的像素点所在行的行数、所在列的列数；
步骤三、基于最小上升率瞬态热响应确定分块大小尺寸；
步骤四、根据分块尺寸进行分块并确定块内搜索步长；
步骤五、分块分步长选取瞬态热响应；
步骤六、采用无监督聚类算法对选取的瞬态热响应进行分类；
步骤七、基于完全自适应SC-SD谱聚类算法实现动态预测并利用基于惩罚项的边界交叉法进行多目标优化选取每类瞬态热响应的代表，构成矩阵Y；
步骤八、基于斯皮尔曼相关系数进行实际缺陷类别数判别，并根据实际缺陷类别数矩阵变化得到二维图像；
步骤九、利用区域卷积神经网络R-CNN对二维图像进行特征提取，得到航天器撞击损伤试件的缺陷特征。
优选的是，其中，所述步骤三基于最小上升率瞬态热响应确定分块大小尺寸的具体步骤包括：
设置分块行阈值K_THVr，分块列阈值K_THVc；依次计算上升率最小的瞬态热响应对应的像素点所在列的其他温度点S(iN,jzz,t)与上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)的相关度
Figure BDA0002666530010000061
其中iN表示距离上升率最小的瞬态热响应对应的像素点最近的第N个点；一直寻找到第一个
Figure BDA0002666530010000062
小于分块行阈值K_THVr的像素点，统计像素点个数N；
依次计算上升率最小的瞬态热响应对应的像素点所在行的其他温度点S(izz,jM,t)与上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)的相关度
Figure BDA0002666530010000063
其中jM表示上升率最小的瞬态热响应对应的像素点最近的第M个点；一直寻找到第一个
Figure BDA0002666530010000064
小于分块列阈值K_THVc的像素点，统计像素点个数M。
优选的是，其中，所述步骤四根据分块尺寸进行分块并确定块内搜索步长的具体步骤包括：
根据基于分块行、列阈值得到的像素点个数M、N，将三维矩阵依次分解为个数为K个，大小为N×M的子三维矩阵块kS(in,jm,t)，其中k表示第k个子三维矩阵块，in、jm、t分别表示第k个子三维矩阵块的第in行、第jm列、第t帧的像素值，n＝1,2,…,N，m＝1,2,…,M，t＝1,2,…T，T为三维矩阵S帧的总数量；
在第k个子三维数据块kS(in,jm,t)中，以子三维数据块中心点为圆心，向四周搜索，找到子三维数据块内中心点，记为kS(kiN/2,kjM/2,kt)，其中，kiN/2,kjM/2,kt分别表示第k个子三维数据块中块内极大值像素点所在行数、列数以及帧数；
设置第k个子三维数据块的块内行阈值R_THVk，依次计算子三维数据块内中心点所在帧、所在行的温度点kS(kin',kjM/2,kt)与子三维数据块内中心点kS(kiN/2,kjM/2,kt)的相关度
Figure BDA0002666530010000065
其中kin'表示在第k个子三维数据块中，子三维数据块内中心点最近的第n'个像素点；计算并统计
Figure BDA0002666530010000066
的温度点个数，记为kRSS，作为第k个子三维数据块的块内行步长；
设置第k个子三维数据块的块内列阈值C_THVk，依次计算子三维数据块内中心点所在帧、所在列的温度点kS(kiN/2,kjm',kt)与子三维数据块内中心点kS(kiN/2,kjM/2,kt)的相关度
Figure BDA0002666530010000067
其中kjm'表示在第k个子三维数据块中，距子三维数据块内中心点kS(kiN/2,kjM/2,kt)最近的第m'个像素点；计算并统计
Figure BDA0002666530010000071
的温度点个数，记为kCSS，作为第k个子三维数据块的块内列步长。
优选的是，其中，所述步骤五分块分步长选取瞬态热响应的具体步骤包括：
步骤S51、根据步骤三统计的M、N个像素值对三维矩阵进行分块，得到K个个大小为N×M×T的子三维数据块，其中kS(in,jm,t)表示第k个子三维数据块内的第in行、第jm列像素的瞬态热响应，其中t＝1,2,…,T，T为原三维矩阵S帧的总数量；
步骤S52、对于每个子三维数据块内的像素点kS(in,jm,t)，设置阈值DD，初始化集合编号g＝1，初始化像素点位置in＝1,jm＝1，并将块内极大值kS(kizz,kjzz,ktzz)对应的瞬态热响应kS(kizz,kjzz,t)，t＝1,2,…,T，存储在集合X(g)中；然后计算子三维数据块内的像素点kS(in,jm,t)中的位于in行，jm列的瞬态热响应与集合X(g)间的相关度Rei,j，并判断：
如果Rei,j<DD，则g＝g+1，并将瞬态热响应kS(in,jm,t)作为一个新特征存储在集合X(g)中；否则，令in＝in+kRSS，继续计算下一个瞬态热响应kS(in,jm,t)，t＝1,2,…,T与集合X(g)的相关度；如果in＞N，则令in＝in-N，jm＝jm+kCSS，即变化到第jm+kCSS列进行计算，如果jm＞M,则第k个子三维矩阵的瞬态热响应选取完毕，k＝k+1；其中，N、M分别为第k个子三维矩阵数据块kS(:,:,:)的行数、列数。
优选的是，其中，所述步骤六采用无监督聚类算法对选取的瞬态热响应进行分类的具体方法为：将步骤五选取的所有K个数据块的所有集合X(g)即瞬态热响应采用模糊C均值聚类FCM算法分为L类，得到每个瞬态热响应所属的类别，具体包括以下步骤：
步骤S61、设置聚类数目L，初始化迭代次数c＝0，设定终止迭代条件阈值ε；
步骤S62、利用公式
Figure BDA0002666530010000072
计算隶属度矩阵U；
其中，i'＝1,2,…,L，c∈L，n'dk'＝||xk'-i′V||，n'＝i',j'，n'dk'表示第k'个像素点与第i'聚类中心i'V的欧氏距离，xk'表示第k'个像素点的坐标；τ为常数；i'uk'表示第k'个像素点隶属于第i'类的程度；
步骤S63、更新聚类中心
Figure BDA0002666530010000081
其中，
Figure BDA0002666530010000082
表示第k'个像素点的热响应值；
步骤S64、如果迭代次数到达最大值L或者前后两次聚类中心之差绝对值小于ε，则算法结束，并输出隶属度矩阵U和聚类中心V，再进入步骤S65；否则，令c＝c+1，返回步骤S62；
步骤S65、利用隶属度最大化准则对所有像素点去模糊化，得到每个像素点所属类别，即Mk'＝argi'max(i'uk')。
优选的是，其中，所述步骤七基于动态多目标选取每类瞬态热响应的代表构成矩阵Y的具体步骤包括：
步骤S71、在第m+1次外部环境下，对第i'(i'＝1,…,L)个类瞬态热响应选代表时，定义多目标函数：
Figure BDA0002666530010000083
其中，
Figure BDA0002666530010000084
为第m+1次外部环境下第i'类瞬态热响应选出的一个瞬态热响应
Figure BDA0002666530010000085
的类内欧氏距离，表示为：
Figure BDA0002666530010000086
为第i'类瞬态热响应选出的一个瞬态热响应
Figure BDA0002666530010000087
的L-1个类间欧氏距离，由计算出的L-1个类间欧氏距离
Figure BDA0002666530010000088
重新编号组成，
Figure BDA0002666530010000089
表示为：
Figure BDA00026665300100000810
Figure BDA00026665300100000811
为瞬态热响应
Figure BDA00026665300100000812
在第t时刻的像素值即温度值，
Figure BDA00026665300100000813
为第i'类瞬态热响应聚类中心在第t时刻的像素值即温度值，
Figure BDA00026665300100000814
为第j'类瞬态热响应聚类中心在第t时刻的像素值即温度值；
步骤S72、第m-1次和第m次环境下得到的多目标函数近似前沿解集分别为
Figure BDA0002666530010000091
和
Figure BDA0002666530010000092
对应的种群瞬态热响应，解集分别为
Figure BDA0002666530010000093
和
Figure BDA0002666530010000094
其数目分别为
Figure BDA0002666530010000095
和
Figure BDA0002666530010000096
在环境变化后，根据第m-1次和第m次环境的历史信息，预测计算第m+1次环境下的近似前沿解集的初始化种群瞬态热响应，步骤如下：
步骤S721、
Figure BDA0002666530010000097
是从
Figure BDA0002666530010000098
解集中随机选择NE个瞬态热响应
Figure BDA0002666530010000099
构成的瞬态热响应集，n'＝1,2,..,NE，计算
Figure BDA00026665300100000910
集中代表瞬态热响应的个数W，用于获得第m+1次环境下多方向预测集：
Figure BDA00026665300100000911
Figure BDA00026665300100000912
Figure BDA00026665300100000913
其中，
Figure BDA00026665300100000914
是从
Figure BDA00026665300100000915
解集中随机选择NE个瞬态热响应
Figure BDA00026665300100000916
构成的瞬态热响应集，n'＝1,2,..,NE；
步骤S722、选择W+1个代表瞬态热响应，其具体步骤包括：
步骤A、代表瞬态热响应构成的PS多方向预测集由两部分组成：
一是
Figure BDA00026665300100000917
解集瞬态热响应的中心，记为
Figure BDA00026665300100000918
Figure BDA00026665300100000919
其中，
Figure BDA00026665300100000920
为解集
Figure BDA00026665300100000921
中第n个瞬态热响应；
二是基于完全自适应谱聚类算法得到的能够充分描述当前PS形状和多样性的W个代表瞬态热响应集
Figure BDA00026665300100000922
此时，集合
Figure BDA00026665300100000923
中代表瞬态热响应的个数
Figure BDA00026665300100000924
为W个；
步骤B、用完全自适应谱聚类算法SC_SD聚类算法将解集
Figure BDA00026665300100000925
中的瞬态热响应聚类成聚类集合
Figure BDA00026665300100000926
其具体步骤包括：
步骤B1、为消除瞬态热响应不同属性的不同量纲对聚类结果的影响，采用最大最小标准化方法对瞬态热响应点进行标准化：
Figure BDA0002666530010000101
其中
Figure BDA0002666530010000102
表示解集
Figure BDA0002666530010000103
中第n个瞬态热响应的第d个时刻的温度值取值d∈{1,2,…,T}，
Figure BDA0002666530010000104
表示解集
Figure BDA0002666530010000105
中所有瞬态热响应样本在第d个个时刻的温度值取值；
步骤B2、计算标准化后每个瞬态热响应样本的全局标准差，以瞬态热响应n为例，n＝1,2,…,NE，其全局标准差公式如下：
Figure BDA0002666530010000106
其中
Figure BDA0002666530010000107
为标准化后瞬态热响应样本n和标准化后瞬态热响应样本j之间的欧氏距离
Figure BDA0002666530010000108
NE为解集
Figure BDA0002666530010000109
中瞬态热响应样本总数；
步骤B3、以stdn为瞬态热响应样本n的领域半径，n＝1,2,…,NE，统计瞬态热响应样本n邻域内的瞬态热响应样本数Cn；利用瞬态热响应样本n在相应邻域的局部标准差获得瞬态热响应样本n对应的自适应局部尺度参数σn：
Figure BDA00026665300100001010
其中Cn为瞬态热响应样本n对应领域半径内瞬态热响应样本的数量；
步骤B4、计算基于自适应的局部尺度参数的亲和矩阵M，亲和矩阵M为由其第n行，j列的元素mnj,n,j∈{1,2,…,NE}构成的NE×NE的方阵，其元素mnj由以下公式确定：
Figure BDA00026665300100001011
步骤B5、计算亲和矩阵M的度矩阵D：
Figure BDA00026665300100001012
度矩阵D为
Figure BDA00026665300100001013
是由dn,n∈{1,2,…,NE}组成的大小为NE×NE对角矩阵，mnj为亲和矩阵的第n行、第j列元素；
步骤B6、计算拉普拉斯矩阵
Figure BDA0002666530010000111
步骤B7、计算拉普拉斯矩阵L的特征值，令λ1≥λ2≥…≥λW≥0为拉普拉斯矩阵L的前W个最大特征值，找到其对应的特征向量ξ1,ξ2,…,ξW，构造矩阵
Figure BDA0002666530010000112
其中ξn为列向量，n∈{1,2,…,W}，令vij表示矩阵V中的第i行、第j列个元素；
步骤B8、按行对矩阵V进行归一化，得到矩阵U，U的第i行，第j列元素值uij由以下公式确定：
Figure BDA0002666530010000113
其中，uij、vik分别表示矩阵U中的第i行、第j列个元素和矩阵V中的第i行、第k列个元素；
步骤B9、将矩阵U每行uN，N∈{1,2,…,NE}，当作一个样本，采用SD_K-medoids算法聚类，其具体步骤包括：
步骤B91、根据下面公式计算矩阵U中每个样本的密度：
Figure BDA0002666530010000115
其中uN,uj,ut分别表示矩阵U的第N，第j，第t行，选取样本密度值最小的样本作为初始聚类中心，将剩余样本作为集合LU，需聚类成的总簇数为W，设置聚类标记数w，置w＝1；
步骤B92、根据公式：
Figure BDA0002666530010000116
计算各个样本的方差VaN,N∈{1,2,…,NE}；
步骤B93、在剩余样本集LU中选择方差最小的样本，记为uv_s(w)，加入到初始聚类中心集合C中，即C＝C∪{uv_s(w)}；
步骤B94、根据公式：
Figure BDA0002666530010000117
计算样本uv_s(w)的标准差Sv_s(w)，以Sv_s(w)为样本uv_s(w)的领域半径rv_s(w)，即
Figure BDA0002666530010000121
计算样本uv_s(w)的领域NBv_s(w)，使其满足：
Figure BDA00026665300100001220
步骤B95、从剩余样本集合LU中去掉NBv_s(w)中的样本，即LU＝LU-NBv_s(w)，w＝w+1，若
Figure BDA00026665300100001221
跳转至步骤B93；
步骤B10、根据矩阵U的聚类结果，将原始瞬态热响应样本分配类簇，即当矩阵U的第n行属于第j个类簇时，将原始瞬态热响应样本n分配到第j个类簇；输出聚类结果
Figure BDA0002666530010000122
步骤C、计算聚类结果中每个类别的聚类中心：
Figure BDA0002666530010000123
其中
Figure BDA0002666530010000124
为第
Figure BDA0002666530010000125
个聚类结果
Figure BDA0002666530010000126
的第k个代表瞬态热响应，
Figure BDA0002666530010000127
为第h个聚类结果中所包含的瞬态热响应总数；
基于聚类结果
Figure BDA0002666530010000128
从每个类别中选择能够充分描述当前PS形状和多样性的一个代表瞬态热响应
Figure BDA0002666530010000129
Figure BDA00026665300100001210
其中，
Figure BDA00026665300100001211
为第
Figure BDA00026665300100001212
个聚类结果
Figure BDA00026665300100001213
的聚类中心；这样一共得到W个能够充分描述当前PS形状和多样性的代表瞬态热响应，加上
Figure BDA00026665300100001214
解集整体瞬态热响应的中心，构成W+1个PS多方向预测集；
步骤S723、根据第m-1次和第m次环境的PS多方向预测集
Figure BDA00026665300100001215
和
Figure BDA00026665300100001216
其中，
Figure BDA00026665300100001217
按照步骤S721、步骤S722的方法获得，W'是
Figure BDA00026665300100001218
集中代表瞬态热响应的个数；
计算预测方向
Figure BDA00026665300100001219
Figure BDA0002666530010000131
其中，
Figure BDA0002666530010000132
是PS多方向预测集
Figure BDA0002666530010000133
中与
Figure BDA0002666530010000134
距离最近的瞬态热响应；
步骤S724、迭代次数g'＝0时，第m+1次环境下的近似前沿解集的初始化种群瞬态热响应个数为Np，其中，
Figure BDA0002666530010000135
个初始种群瞬态热响应在取值范围内随机生成，
Figure BDA0002666530010000136
个初始种群瞬态热响应由下式预测获得：
Figure BDA0002666530010000137
其中，hn为瞬态热响应
Figure BDA0002666530010000138
所属聚类结果
Figure BDA0002666530010000139
的序号，
Figure BDA00026665300100001310
是一个服从均值为0，方差为
Figure BDA00026665300100001311
的正态分布的随机数，方差
Figure BDA00026665300100001312
的计算公式为：
Figure BDA00026665300100001313
步骤S73、初始化相关参数
初始化迭代次数g'＝0，一组均匀分布的权重向量
Figure BDA00026665300100001314
其中
Figure BDA00026665300100001315
初始化参考点
Figure BDA00026665300100001316
是函数
Figure BDA00026665300100001317
对应的参考点；
Figure BDA00026665300100001318
最大迭代次数g'max；
初始化每个种群瞬态热响应的进化速度为
Figure BDA00026665300100001319
种群瞬态热响应的全局最优和局部最优满足
Figure BDA00026665300100001320
步骤S74、利用
Figure BDA00026665300100001321
构造基于惩罚项的边界交叉法下每个种群瞬态热响应的动态目标函数适应度值
Figure BDA00026665300100001322
Figure BDA00026665300100001323
其中
Figure BDA00026665300100001324
为预设参数，
Figure BDA00026665300100001325
和
Figure BDA00026665300100001326
的设定如下：
Figure BDA00026665300100001327
Figure BDA00026665300100001328
步骤S75、对n＝1,…,NP：按照粒子群算法更新速度
Figure BDA0002666530010000141
和种群瞬态热响应
Figure BDA0002666530010000142
按照多目标优化算法比较
Figure BDA0002666530010000143
更新全局最优
Figure BDA0002666530010000144
局部最优
Figure BDA0002666530010000145
和参考点
Figure BDA0002666530010000146
从
Figure BDA0002666530010000147
中保留支配
Figure BDA0002666530010000148
的解向量，移除所有被
Figure BDA0002666530010000149
支配的解向量，如果
Figure BDA00026665300100001410
中的向量都不支配
Figure BDA00026665300100001411
将
Figure BDA00026665300100001412
加入
Figure BDA00026665300100001413
同时n＝n+1，若n≤NP，则g'＝g'+1；
步骤S76、进化终止判断：若g'≤g'max，则重复步骤S75，若g'＞g'max，则得到第i'类温度瞬态热响应最终的前沿近似解集
Figure BDA00026665300100001414
步骤S77、基于加权隶属度方案从前沿近似解集
Figure BDA00026665300100001415
选出第i'类瞬态热响应的代表i'REP；
根据下列公式计算前沿近似解集
Figure BDA00026665300100001416
中第k个解对第l个目标函数的隶属度：
Figure BDA00026665300100001417
其中，Fl为第l个目标函数的值，Fl max、Fl min分别为对应目标函数的最大和最小值；
为目标函数设置权重λl,(l＝1,2,…,L)，计算前沿近似解集的隶属度加权值，取所得最大值对应的前沿解即为第i'类瞬态热响应的代表i'REP，公式如下：
Figure BDA00026665300100001418
其中，
Figure BDA00026665300100001419
为前沿近似解集
Figure BDA00026665300100001420
所包含的解集个数，L为目标函数个数，
Figure BDA00026665300100001421
为隶属度函数值；
所有L类的瞬态热响应代表按列放置，一列为T个时刻的像素值即温度值，构成一个T×L的矩阵Y。
优选的是，其中，所述步骤八计算瞬态热响应代表之间的斯皮尔曼相关系数的具体步骤包括：
以i'REP和j'REP，(i',j'＝1,2,…,L)表示任意两个瞬态热响应代表，将i'REP和j'REP对应的各个时刻的温度元素值i'REPt,(t＝1,2,…,T)和j'REPt,(t＝1,2,…,T)转换为其在各自瞬态热响应向量总体时刻温度值中的降序排名，将温度值最高的元素转换为1，温度值最低的元素转换为T，其余温度值元素按大小顺序转换为排名并记为Ra(i'REPt)和Ra(j'REPt)；根据公式：
Figure BDA0002666530010000151
计算两个瞬态热响应代表之间对应时刻元素值的差异Da；最后根据公式：
Figure BDA0002666530010000152
计算两个瞬态热响应代表i'REP和j'REP之间的斯皮尔曼相关系数；
设定相关性阈值θ，两两之间比较斯皮尔曼相关系数，保留相关性最小的瞬态热响应代表，若Rs(i'REPt,j'REPt)＜θ，则保留第i，j类瞬态热响应代表，否则剔除其中一类，得到L'类瞬态热响应代表；所有L'类的瞬态热响应代表按列放置，一列为T个时刻的像素值即温度值，构成一个T×L'的矩阵Y'；
将三维矩阵S中的每一帧从第一列开始，将后一列接在前一列的末尾，构成新的一列，得到T帧对应的T列像素值，然后，依据时间先后，将T列像素值依次放置，构成I×J行、T列二维图像矩阵O，用矩阵Y对二维矩阵O进行线性变换，即：
Figure BDA0002666530010000153
得到二维图像矩阵R，其中，
Figure BDA0002666530010000154
为L'×T矩阵，是矩阵Y'的伪逆矩阵，OT二维图像矩阵O的转置矩阵，得到的二维图像矩阵R为L'行、I×J列；
二维图像矩阵R的每一行，按J列依次进行截取，截取的J列按行依次放置，构成一张I×J二维图像，这样L'行，得到L'张I×J二维图像，这些图片都包含了缺陷区域，为方便缺陷轮廓提取，选择缺陷区域和非缺陷区域像素值差距最大的一张二维图像，并记为f(x,y)。
优选的是，其中，所述步骤九利用区域卷积神经网络R-CNN对步骤八中的二维图像f(x,y)进行图像中的缺陷检测和缺陷区域分割，从而实现缺陷部分的定位识别和特征提取。
本发明至少包括以下有益效果：
1、本发明采用多目标优化方法实现了差异性与相似性的综合考量，并准确的刻画缺陷轮廓，弥补了传统方法对于降维处理上的一些不足，比只基于差异性的算法提取缺陷特征更具代表性；
2、本发明采用多方向预测策略，结合基于样本邻域标准差的完全自适应的SC-SD谱聚类算法快速选择和引入多个代表瞬态热响应适当地描述PS(Pareto Set)的形状，记录每次环境PS的分布情况，以此来预测PS的新位置。在环境发生了变化后，用前两次环境的代表瞬态热响应来预测PS的新位置，在新位置生成若干新的初始种群瞬态热响应，由此加快了对环境变化的响应。基于样本邻域标准差的完全自适应的SC-SD谱聚类算法，在保证能够准确获得PS种群的形状信息和位置信息的前提下，自适应地确定样本的尺度参数，消除主观性，同时引入SD_K-medoids算法，进一步提升算法聚类准确度，实现PS新位置预测的快速性和有效性；
3、本发明采用基于惩罚项的边界交叉法进行多目标问题分解。使得在处理目标数超过两维的多目标优化问题时，相比切比雪夫方法得到的最优解分布更均匀，并在处理高维多目标优化问题时，基于惩罚项的边界交叉法要显著优于切比雪夫方法，更契合航天器的空间碎片复杂撞击损伤检测评估方面的要求。同时，由于惩罚项的引入，能够自由地选择进化算法获得最优解的收敛性与多样性间的平衡，以应对动态环境下动态多目标优化问题的不同要求。
本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明：
图1是本发明基于动态多目标优化的红外热图像缺陷特征提取方法一种具体实施方式的流程图；
图2是改进后的红外热图像数据大尺寸分块变步长的瞬态热响应搜索流程图；
图3是完全自适应的SC-SD谱聚类算法流程图；
图4是采用模糊C均值聚类对选取的瞬态热响应进行分类后的结果图；
图5是基于惩罚项的边界交叉法得到的缺陷1的前沿近似解集散点图，以及基于加权隶属度方案选择的缺陷1的瞬态热响应代表；
图6是基于惩罚项的边界交叉法得到的缺陷2的前沿近似解集散点图，以及基于加权隶属度方案选择的缺陷2的瞬态热响应代表；
图7是基于惩罚项的边界交叉法得到的背景区域的前沿近似解集散点图，以及基于加权隶属度方案选择的背景区域的瞬态热响应代表；
图8是缺陷1温度点的瞬态热响应曲线图；
图9是缺陷2温度点的瞬态热响应曲线图；
图10是背景区域温度点的瞬态热响应曲线图；
图11是基于本发明选取的对应缺陷1温度点的瞬态热响应曲线图；
图12是基于本发明选取的对应缺陷2温度点的瞬态热响应曲线图；
图13是基于本发明选取的对应背景区域温度点的瞬态热响应曲线图；
图14是基于本发明提取的缺陷1特征图；
图15是基于本发明提取的缺陷2特征图。
具体实施方式：
下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不配出一个或多个其它元件或其组合的存在或添加。
如图1-3所示：本发明的一种在轨航天器微小撞击损伤可视化检测评估方法，包括以下步骤：
步骤一、将红外热像仪获取的航天器撞击损伤试件的热图像序列用三维矩阵S表示，其中三维矩阵S的元素S(i,j,t)表示热图像序列的t帧热图像的第i行、第j列的像素值；
步骤二、从三维矩阵S的所有瞬态热响应中，选取上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)，其中，izz、jzz分别表示上升率最小的瞬态热响应对应的像素点所在行的行数、所在列的列数；
步骤三、基于最小上升率瞬态热响应确定分块大小尺寸；
步骤四、根据分块尺寸进行分块并确定块内搜索步长；
步骤五、分块分步长选取瞬态热响应；
步骤六、采用无监督聚类算法对选取的瞬态热响应进行分类；
步骤七、基于完全自适应SC-SD谱聚类算法实现动态预测并利用基于惩罚项的边界交叉法进行多目标优化选取每类瞬态热响应的代表，构成矩阵Y；
步骤八、基于斯皮尔曼相关系数进行实际缺陷类别数判别，并根据实际缺陷类别数矩阵变化得到二维图像；
步骤九、利用区域卷积神经网络R-CNN对二维图像进行特征提取，得到航天器撞击损伤试件的缺陷特征。
本发明基于动态多目标优化的红外热图像缺陷特征提取方法，通过对热图像序列变换步长选取出像素点的瞬态热响应，并采用FCM进行分类，得到各个像素点的瞬态热响应的所属类别，然后考虑每个类别像素点与同类像素点的像素值(温度值)相似性，同时考虑该像素点(温度点)与不同类别像素点(温度点)的差异性，构造相应的多目标函数，同时，在每次环境发生变化后，通过预测机制，采用结合了基于样本邻域标准差的完全自适应的SC-SD谱聚类算法，自适应地确定样本的尺度参数,以发现数据集样本的潜在分布规律,避免了人为给定参数带来的聚类结果不确定性，消除主观性，提高了聚类算法的运行效率，同时引入SD_K-medoids算法，进一步提升算法聚类准确性，使得聚类算法能够更有效地获得PS种群的形状信息和位置信息，保证了预测算法的准确性和有效性，实现更准确，更快速的种群跟踪，从而提升动态多目标算法整体的运行效率，为种群进化提供引导方向，帮助多目标优化算法对新变化做出快速响应。通过基于惩罚项的边界交叉法的多目标优化算法，使得进化算法获得的最优解分布更均匀，同时更方便地根据不同的多目标优化问题环境灵活地决定最优解间多样性和收敛性间地平衡。然后利用加权隶属度，从各个类别的前沿最优解集中选择各类别的瞬态热响应代表。再利用基于斯皮尔曼相关系数的缺陷类别数判别方法，判别实际的缺陷类别总数，从而更精确地获得热图像序列的降维结果，最后利用区域卷积神经网络进行特征提取，从而提取出红外热图像的缺陷特征。通过上述步骤，实现代表瞬态热相应(温度点)的精确选择，保证了缺陷特征提取的精准度，同时降低了动态环境下获取各个类别信息代表瞬态热相应的计算消耗。
在上述技术方案中，所述步骤三基于最小上升率瞬态热响应确定分块大小尺寸的具体步骤包括：
设置分块行阈值K_THVr，分块列阈值K_THVc；依次计算上升率最小的瞬态热响应对应的像素点所在列的其他温度点S(iN,jzz,t)与上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)的相关度
Figure BDA0002666530010000191
其中iN表示距离上升率最小的瞬态热响应对应的像素点最近的第N个点；一直寻找到第一个
Figure BDA0002666530010000192
小于分块行阈值K_THVr的像素点，统计像素点个数N；
依次计算上升率最小的瞬态热响应对应的像素点所在行的其他温度点S(izz,jM,t)与上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)的相关度
Figure BDA0002666530010000193
其中jM表示上升率最小的瞬态热响应对应的像素点最近的第M个点；一直寻找到第一个
Figure BDA0002666530010000194
小于分块列阈值K_THVc的像素点，统计像素点个数M。
在上述技术方案中，所述步骤四根据分块尺寸进行分块并确定块内搜索步长的具体步骤包括：
根据基于分块行、列阈值得到的像素点个数M、N，将三维矩阵依次分解为个数为K个，大小为N×M的子三维矩阵块kS(in,jm,t)，其中k表示第k个子三维矩阵块，in、jm、t分别表示第k个子三维矩阵块的第in行、第jm列、第t帧的像素值，n＝1,2,…,N，m＝1,2,…,M，t＝1,2,…T，T为三维矩阵S帧的总数量；
在第k个子三维数据块kS(in,jm,t)中，以子三维数据块中心点为圆心，向四周搜索，找到子三维数据块内中心点，记为kS(kiN/2,kjM/2,kt)，其中，kiN/2,kjM/2,kt分别表示第k个子三维数据块中块内极大值像素点所在行数、列数以及帧数；
设置第k个子三维数据块的块内行阈值R_THVk，依次计算子三维数据块内中心点所在帧、所在行的温度点kS(kin',kjM/2,kt)与子三维数据块内中心点kS(kiN/2,kjM/2,kt)的相关度
Figure BDA0002666530010000195
其中kin'表示在第k个子三维数据块中，子三维数据块内中心点最近的第n'个像素点；计算并统计
Figure BDA0002666530010000196
的温度点个数，记为kRSS，作为第k个子三维数据块的块内行步长；
设置第k个子三维数据块的块内列阈值C_THVk，依次计算子三维数据块内中心点所在帧、所在列的温度点kS(kiN/2,kjm',kt)与子三维数据块内中心点kS(kiN/2,kjM/2,kt)的相关度
Figure BDA0002666530010000201
其中kjm'表示在第k个子三维数据块中，距子三维数据块内中心点kS(kiN/2,kjM/2,kt)最近的第m'个像素点；计算并统计
Figure BDA0002666530010000202
的温度点个数，记为kCSS，作为第k个子三维数据块的块内列步长。
在上述技术方案中，所述步骤五分块分步长选取瞬态热响应的具体步骤包括：
步骤S51、根据步骤三统计的M、N个像素值对三维矩阵进行分块，得到K个个大小为N×M×T的子三维数据块，其中kS(in,jm,t)表示第k个子三维数据块内的第in行、第jm列像素的瞬态热响应，其中t＝1,2,…,T，T为原三维矩阵S帧的总数量；
步骤S52、对于每个子三维数据块内的像素点kS(in,jm,t)，设置阈值DD，初始化集合编号g＝1，初始化像素点位置in＝1,jm＝1，并将块内极大值kS(kizz,kjzz,ktzz)对应的瞬态热响应kS(kizz,kjzz,t)，t＝1,2,…,T，存储在集合X(g)中；然后计算子三维数据块内的像素点kS(in,jm,t)中的位于in行，jm列的瞬态热响应与集合X(g)间的相关度Rei,j，并判断：
如果Rei,j<DD，则g＝g+1，并将瞬态热响应kS(in,jm,t)作为一个新特征存储在集合X(g)中；否则，令in＝in+kRSS，继续计算下一个瞬态热响应kS(in,jm,t)，t＝1,2,…,T与集合X(g)的相关度；如果in＞N，则令in＝in-N，jm＝jm+kCSS，即变化到第jm+kCSS列进行计算，如果jm＞M,则第k个子三维矩阵的瞬态热响应选取完毕，k＝k+1；其中，N、M分别为第k个子三维矩阵数据块kS(:,:,:)的行数、列数
在上述技术方案中，所述步骤六采用无监督聚类算法对选取的瞬态热响应进行分类的具体方法为：将步骤五选取的所有K个数据块的所有集合X(g)即瞬态热响应采用模糊C均值聚类FCM算法分为L类，得到每个瞬态热响应所属的类别，具体包括以下步骤：
步骤S61、设置聚类数目L，初始化迭代次数c＝0，设定终止迭代条件阈值ε；
步骤S62、利用公式
Figure BDA0002666530010000211
计算隶属度矩阵U；
其中，i'＝1,2,…,L，c∈L，n'dk'＝||xk'-i'V||，n'＝i',j'，n'dk'表示第k'个像素点与第i'聚类中心i'V的欧氏距离，xk'表示第k'个像素点的坐标；τ为常数；i'uk'表示第k'个像素点隶属于第i'类的程度；
步骤S63、更新聚类中心
Figure BDA0002666530010000212
其中，
Figure BDA0002666530010000213
表示第k'个像素点的热响应值；
步骤S64、如果迭代次数到达最大值L或者前后两次聚类中心之差绝对值小于ε，则算法结束，并输出隶属度矩阵U和聚类中心V，再进入步骤S65；否则，令c＝c+1，返回步骤S62；
步骤S65、利用隶属度最大化准则对所有像素点去模糊化，得到每个像素点所属类别，即Mk'＝argi'max(i'uk')。
在上述技术方案中，所述步骤七基于动态多目标选取每类瞬态热响应的代表构成矩阵Y的具体步骤包括：
步骤S71、在第m+1次外部环境下，对第i'(i'＝1,…,L)个类瞬态热响应选代表时，定义多目标函数：
Figure BDA0002666530010000214
其中，
Figure BDA0002666530010000215
为第m+1次外部环境下第i'类瞬态热响应选出的一个瞬态热响应
Figure BDA0002666530010000216
的类内欧氏距离，表示为：
Figure BDA0002666530010000217
为第i'类瞬态热响应选出的一个瞬态热响应
Figure BDA0002666530010000218
的L-1个类间欧氏距离，由计算出的L-1个类间欧氏距离
Figure BDA0002666530010000219
重新编号组成，
Figure BDA00026665300100002110
表示为：
Figure BDA0002666530010000221
Figure BDA0002666530010000222
为瞬态热响应
Figure BDA0002666530010000223
在第t时刻的像素值即温度值，
Figure BDA0002666530010000224
为第i'类瞬态热响应聚类中心在第t时刻的像素值即温度值，
Figure BDA0002666530010000225
为第j'类瞬态热响应聚类中心在第t时刻的像素值即温度值；
步骤S72、第m-1次和第m次环境下得到的多目标函数近似前沿解集分别为
Figure BDA0002666530010000226
和
Figure BDA0002666530010000227
对应的种群瞬态热响应，解集分别为
Figure BDA0002666530010000228
和
Figure BDA0002666530010000229
其数目分别为
Figure BDA00026665300100002210
和
Figure BDA00026665300100002211
在环境变化后，根据第m-1次和第m次环境的历史信息，预测计算第m+1次环境下的近似前沿解集的初始化种群瞬态热响应，步骤如下：
步骤S721、
Figure BDA00026665300100002212
是从
Figure BDA00026665300100002213
解集中随机选择NE个瞬态热响应
Figure BDA00026665300100002214
构成的瞬态热响应集，n'＝1,2,..,NE，计算
Figure BDA00026665300100002215
集中代表瞬态热响应的个数W，用于获得第m+1次环境下多方向预测集：
Figure BDA00026665300100002216
Figure BDA00026665300100002217
Figure BDA00026665300100002218
其中，
Figure BDA00026665300100002219
是从
Figure BDA00026665300100002220
解集中随机选择NE个瞬态热响应
Figure BDA00026665300100002221
构成的瞬态热响应集，n'＝1,2,..,NE；
步骤S722、选择W+1个代表瞬态热响应，其具体步骤包括：
步骤A、代表瞬态热响应构成的PS多方向预测集由两部分组成：
一是
Figure BDA00026665300100002222
解集瞬态热响应的中心，记为
Figure BDA00026665300100002223
Figure BDA00026665300100002224
其中，
Figure BDA00026665300100002225
为解集
Figure BDA00026665300100002226
中第n个瞬态热响应；
二是基于完全自适应谱聚类算法得到的能够充分描述当前PS形状和多样性的W个代表瞬态热响应集
Figure BDA00026665300100002227
此时，集合
Figure BDA0002666530010000231
中代表瞬态热响应的个数
Figure BDA0002666530010000232
为W个；
步骤B、用完全自适应谱聚类算法SC_SD聚类算法将解集
Figure BDA0002666530010000233
中的瞬态热响应聚类成聚类集合
Figure BDA0002666530010000234
其具体步骤包括：
步骤B1、为消除瞬态热响应不同属性的不同量纲对聚类结果的影响，采用最大最小标准化方法对瞬态热响应点进行标准化：
Figure BDA0002666530010000235
其中
Figure BDA0002666530010000236
表示解集
Figure BDA0002666530010000237
中第n个瞬态热响应的第d个时刻的温度值取值d∈{1,2,…,T}，
Figure BDA0002666530010000238
表示解集
Figure BDA0002666530010000239
中所有瞬态热响应样本在第d个个时刻的温度值取值；
步骤B2、计算标准化后每个瞬态热响应样本的全局标准差，以瞬态热响应n为例，n＝1,2,…,NE，其全局标准差公式如下：
Figure BDA00026665300100002310
其中
Figure BDA00026665300100002311
为标准化后瞬态热响应样本n和标准化后瞬态热响应样本j之间的欧氏距离
Figure BDA00026665300100002312
NE为解集
Figure BDA00026665300100002313
中瞬态热响应样本总数；
步骤B3、以stdn为瞬态热响应样本n的领域半径，n＝1,2,…,NE，统计瞬态热响应样本n邻域内的瞬态热响应样本数Cn；利用瞬态热响应样本n在相应邻域的局部标准差获得瞬态热响应样本n对应的自适应局部尺度参数σn：
Figure BDA00026665300100002314
其中Cn为瞬态热响应样本n对应领域半径内瞬态热响应样本的数量；
步骤B4、计算基于自适应的局部尺度参数的亲和矩阵M，亲和矩阵M为由其第n行，j列的元素mnj,n,j∈{1,2,…,NE}构成的NE×NE的方阵，其元素mnj由以下公式确定：
Figure BDA00026665300100002315
步骤B5、计算亲和矩阵M的度矩阵D：
Figure BDA0002666530010000241
度矩阵D为
Figure BDA0002666530010000242
是由dn,n∈{1,2,…,NE}组成的大小为NE×NE对角矩阵，mnj为亲和矩阵的第n行、第j列元素；
步骤B6、计算拉普拉斯矩阵
Figure BDA0002666530010000243
步骤B7、计算拉普拉斯矩阵L的特征值，令λ1≥λ2≥…≥λW≥0为拉普拉斯矩阵L的前W个最大特征值，找到其对应的特征向量ξ1,ξ2,…,ξW，构造矩阵
Figure BDA0002666530010000244
其中ξn为列向量，n∈{1,2,…,W}，令vij表示矩阵V中的第i行、第j列个元素；
步骤B8、按行对矩阵V进行归一化，得到矩阵U，U的第i行，第j列元素值uij由以下公式确定：
Figure BDA0002666530010000245
其中，uij、vik分别表示矩阵U中的第i行、第j列个元素和矩阵V中的第i行、第k列个元素；
步骤B9、将矩阵U每行uN，N∈{1,2,…,NE}，当作一个样本，采用SD_K-medoids算法聚类，其具体步骤包括：
步骤B91、根据下面公式计算矩阵U中每个样本的密度：
Figure BDA0002666530010000247
其中uN,uj,ut分别表示矩阵U的第N，第j，第t行，选取样本密度值最小的样本作为初始聚类中心，将剩余样本作为集合LU，需聚类成的总簇数为W，设置聚类标记数w，置w＝1；
步骤B92、根据公式：
Figure BDA0002666530010000248
计算各个样本的方差VaN,N∈{1,2,…,NE}；
步骤B93、在剩余样本集LU中选择方差最小的样本，记为uv_s(w)，加入到初始聚类中心集合C中，即C＝C∪{uv_s(w)}；
步骤B94、根据公式：
Figure BDA0002666530010000251
计算样本uv_s(w)的标准差Sv_s(w)，以Sv_s(w)为样本uv_s(w)的领域半径rv_s(w)，即
Figure BDA0002666530010000252
计算样本uv_s(w)的领域NBv_s(w)，使其满足：
Figure BDA00026665300100002519
步骤B95、从剩余样本集合LU中去掉NBv_s(w)中的样本，即LU＝LU-NBv_s(w)，w＝w+1，若
Figure BDA00026665300100002520
跳转至步骤B93；
步骤B10、根据矩阵U的聚类结果，将原始瞬态热响应样本分配类簇，即当矩阵U的第n行属于第j个类簇时，将原始瞬态热响应样本n分配到第j个类簇；输出聚类结果
Figure BDA0002666530010000253
步骤C、计算聚类结果中每个类别的聚类中心：
Figure BDA0002666530010000254
其中
Figure BDA0002666530010000255
为第
Figure BDA0002666530010000256
个聚类结果
Figure BDA0002666530010000257
的第k个代表瞬态热响应，
Figure BDA0002666530010000258
为第h个聚类结果中所包含的瞬态热响应总数；
基于聚类结果
Figure BDA0002666530010000259
从每个类别中选择能够充分描述当前PS形状和多样性的一个代表瞬态热响应
Figure BDA00026665300100002510
Figure BDA00026665300100002511
其中，
Figure BDA00026665300100002512
为第
Figure BDA00026665300100002513
个聚类结果
Figure BDA00026665300100002514
的聚类中心；这样一共得到W个能够充分描述当前PS形状和多样性的代表瞬态热响应，加上
Figure BDA00026665300100002515
解集整体瞬态热响应的中心，构成W+1个PS多方向预测集；
步骤S723、根据第m-1次和第m次环境的PS多方向预测集
Figure BDA00026665300100002516
和
Figure BDA00026665300100002517
其中，
Figure BDA00026665300100002518
按照步骤S721、步骤S722的方法获得，W'是
Figure BDA0002666530010000261
集中代表瞬态热响应的个数；
计算预测方向
Figure BDA0002666530010000262
Figure BDA0002666530010000263
其中，
Figure BDA0002666530010000264
是PS多方向预测集
Figure BDA0002666530010000265
中与
Figure BDA0002666530010000266
距离最近的瞬态热响应；
步骤S724、迭代次数g'＝0时，第m+1次环境下的近似前沿解集的初始化种群瞬态热响应个数为Np，其中，
Figure BDA0002666530010000267
个初始种群瞬态热响应在取值范围内随机生成，
Figure BDA0002666530010000268
个初始种群瞬态热响应由下式预测获得：
Figure BDA0002666530010000269
其中，hn为瞬态热响应
Figure BDA00026665300100002610
所属聚类结果
Figure BDA00026665300100002611
的序号，
Figure BDA00026665300100002612
是一个服从均值为0，方差为
Figure BDA00026665300100002613
的正态分布的随机数，方差
Figure BDA00026665300100002614
的计算公式为：
Figure BDA00026665300100002615
步骤S73、初始化相关参数
初始化迭代次数g'＝0，一组均匀分布的权重向量
Figure BDA00026665300100002616
其中
Figure BDA00026665300100002617
初始化参考点
Figure BDA00026665300100002618
是函数
Figure BDA00026665300100002619
对应的参考点；
Figure BDA00026665300100002620
最大迭代次数g'max；
初始化每个种群瞬态热响应的进化速度为
Figure BDA00026665300100002621
种群瞬态热响应的全局最优和局部最优满足
Figure BDA00026665300100002622
步骤S74、利用
Figure BDA00026665300100002623
构造基于惩罚项的边界交叉法下每个种群瞬态热响应的动态目标函数适应度值
Figure BDA00026665300100002624
Figure BDA00026665300100002625
其中
Figure BDA00026665300100002626
为预设参数，
Figure BDA00026665300100002627
和
Figure BDA00026665300100002628
的设定如下：
Figure BDA00026665300100002629
Figure BDA0002666530010000271
步骤S75、对n＝1,…,NP：按照粒子群算法更新速度
Figure BDA0002666530010000272
和种群瞬态热响应
Figure BDA0002666530010000273
按照多目标优化算法比较
Figure BDA0002666530010000274
更新全局最优
Figure BDA0002666530010000275
局部最优
Figure BDA0002666530010000276
和参考点
Figure BDA0002666530010000277
从
Figure BDA0002666530010000278
中保留支配
Figure BDA0002666530010000279
的解向量，移除所有被
Figure BDA00026665300100002710
支配的解向量，如果
Figure BDA00026665300100002711
中的向量都不支配
Figure BDA00026665300100002712
将
Figure BDA00026665300100002713
加入
Figure BDA00026665300100002714
同时n＝n+1，若n≤NP，则g'＝g'+1；
步骤S76、进化终止判断：若g'≤g'max，则重复步骤S75，若g'＞g'max，则得到第i'类温度瞬态热响应最终的前沿近似解集
Figure BDA00026665300100002715
步骤S77、基于加权隶属度方案从前沿近似解集
Figure BDA00026665300100002716
选出第i'类瞬态热响应的代表i'REP；
根据下列公式计算前沿近似解集
Figure BDA00026665300100002717
中第k个解对第l个目标函数的隶属度：
Figure BDA00026665300100002718
其中，Fl为第l个目标函数的值，Fl max、Fl min分别为对应目标函数的最大和最小值；
为目标函数设置权重λl,(l＝1,2,…,L)，计算前沿近似解集的隶属度加权值，取所得最大值对应的前沿解即为第i'类瞬态热响应的代表i'REP，公式如下：
Figure BDA00026665300100002719
其中，
Figure BDA00026665300100002720
为前沿近似解集
Figure BDA00026665300100002721
所包含的解集个数，L为目标函数个数，
Figure BDA00026665300100002722
为隶属度函数值；
所有L类的瞬态热响应代表按列放置，一列为T个时刻的像素值即温度值，构成一个T×L的矩阵Y。
在上述技术方案中，所述步骤八计算瞬态热响应代表之间的斯皮尔曼相关系数的具体步骤包括：
以i'REP和j'REP，(i',j'＝1,2,…,L)表示任意两个瞬态热响应代表，将i'REP和j'REP对应的各个时刻的温度元素值i'REPt,(t＝1,2,…,T)和j'REPt,(t＝1,2,…,T)转换为其在各自瞬态热响应向量总体时刻温度值中的降序排名，将温度值最高的元素转换为1，温度值最低的元素转换为T，其余温度值元素按大小顺序转换为排名并记为Ra(i'REPt)和Ra(j'REPt)；根据公式：
Figure BDA0002666530010000281
计算两个瞬态热响应代表之间对应时刻元素值的差异Da；最后根据公式：
Figure BDA0002666530010000282
计算两个瞬态热响应代表i'REP和j'REP之间的斯皮尔曼相关系数；
设定相关性阈值θ，两两之间比较斯皮尔曼相关系数，保留相关性最小的瞬态热响应代表，若Rs(i'REPt,j'REPt)＜θ，则保留第i，j类瞬态热响应代表，否则剔除其中一类，得到L'类瞬态热响应代表；所有L'类的瞬态热响应代表按列放置，一列为T个时刻的像素值即温度值，构成一个T×L'的矩阵Y'；
将三维矩阵S中的每一帧从第一列开始，将后一列接在前一列的末尾，构成新的一列，得到T帧对应的T列像素值，然后依据时间先后，将T列像素值依次放置，构成I×J行、T列二维图像矩阵O，用矩阵Y对二维矩阵O进行线性变换，即：
Figure BDA0002666530010000283
得到二维图像矩阵R，其中，
Figure BDA0002666530010000284
为L'×T矩阵，是矩阵Y'的伪逆矩阵，OT二维图像矩阵O的转置矩阵，得到的二维图像矩阵R为L'行、I×J列；
二维图像矩阵R的每一行，按J列依次进行截取，截取的J列按行依次放置，构成一张I×J二维图像，这样L'行，得到L'张I×J二维图像，这些图片都包含了缺陷区域，为方便缺陷轮廓提取，选择缺陷区域和非缺陷区域像素值差距最大的一张二维图像，并记为f(x,y)。
在上述技术方案中，所述步骤九利用区域卷积神经网络R-CNN对步骤八中的二维图像f(x,y)进行图像中的缺陷检测和缺陷区域分割，从而实现缺陷部分的定位识别和特征提取。
实施例：
在本实施例中，在试件上有两种缺陷，即填充导热性差材料的缺陷1以及背部穿孔的缺陷2。
红外热图像数据大尺寸分块变步长的瞬态热响应搜索流程图如图2所示。
完全自适应的SC-SD谱聚类算法流程图如图3所示。
在本实例中，采用模糊C均值聚类对选取的瞬态热响应进行分类后的结果图如图4所示。
采集基于惩罚项的边界交叉法获得的缺陷1温度点、缺陷2温度点以及材料本身温度点的近似前沿解集为1AP、2AP以及3AP。同时采用加权隶属度方案选择的缺陷1温度点、缺陷2温度点以及材料本身温度点的代表瞬态热响应为AFV63、BFV13以及CFV25，如图5、图6、图7所示。
在试件的热像图序列中直接提取三个已知温度点，即材料本身温度点、缺陷1温度点以及缺陷2温度点的瞬态热响应曲线，分别记为TTRBackground、TTRDe1以及TTRDe2，如图8、图9、图10所示。
用本发明中动态多目标优化选瞬态热响应代表的方法，获得了三个瞬态热响应代表：AFV63、BFV13以及CFV25，它们分别对应缺陷1温度点、缺陷2温度点以及材料本身温度点，其曲线如图11、图12、图13所示。
由热响应曲线可知：缺陷1温度点峰值明显高于背景区域温度峰值，缺陷2温度点的峰值最高，大大高于背景区域峰值，背景区域温度点峰值最低。三个特征相比，缺陷1温度点吸热较多，缺陷2温度点吸收热量最多。
本发明的瞬态热响应曲线与直接从热像图序列提取对应的瞬态热响应曲线相关度如表1所示。
表1
本身温度点 缺陷1温度点 缺陷2温度点
本发明 0.997 0.998 0.998
从表1，可以看出，本发明方法选取的瞬态热响应曲线的相关性较好。
在本实施例中，提取的缺陷1和缺陷2的特征如图12、13所示。
这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (8)

1.一种在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，包括以下步骤：
步骤一、将红外热像仪获取的航天器撞击损伤试件的热图像序列用三维矩阵S表示，其中三维矩阵S的元素S(i,j,t)表示热图像序列的t帧热图像的第i行、第j列的像素值；
步骤二、从三维矩阵S的所有瞬态热响应中，选取上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)，其中，izz、jzz分别表示上升率最小的瞬态热响应对应的像素点所在行的行数、所在列的列数；
步骤三、基于最小上升率瞬态热响应确定分块大小尺寸；
步骤四、根据分块尺寸进行分块并确定块内搜索步长；
步骤五、分块分步长选取瞬态热响应；
步骤六、采用无监督聚类算法对选取的瞬态热响应进行分类；
步骤七、基于完全自适应SC-SD谱聚类算法实现动态预测并利用基于惩罚项的边界交叉法进行多目标优化选取每类瞬态热响应的代表，构成矩阵Y；
步骤八、基于斯皮尔曼相关系数进行实际缺陷类别数判别，并根据实际缺陷类别数矩阵变化得到二维图像；
步骤九、利用区域卷积神经网络R-CNN对二维图像进行特征提取，得到航天器撞击损伤试件的缺陷特征。
2.如权利要求1所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤三基于最小上升率瞬态热响应确定分块大小尺寸的具体步骤包括：
设置分块行阈值K_THVr，分块列阈值K_THVc；依次计算上升率最小的瞬态热响应对应的像素点所在列的其他温度点S(iN,jzz,t)与上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)的相关度
Figure FDA0003541315390000011
其中iN表示距离上升率最小的瞬态热响应对应的像素点最近的第N个点；一直寻找到第一个
Figure FDA0003541315390000012
小于分块行阈值K_THVr的像素点，统计像素点个数N；
依次计算上升率最小的瞬态热响应对应的像素点所在行的其他温度点S(izz,jM,t)与上升率最小的瞬态热响应对应的像素点S(izz,jzz,t)的相关度
Figure FDA0003541315390000021
其中jM表示上升率最小的瞬态热响应对应的像素点最近的第M个点；一直寻找到第一个
Figure FDA0003541315390000022
小于分块列阈值K_THVc的像素点，统计像素点个数M。
3.如权利要求1所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤四根据分块尺寸进行分块并确定块内搜索步长的具体步骤包括：
根据基于分块行、列阈值得到的像素点个数M、N，将三维矩阵依次分解为个数为K个，大小为N×M的子三维矩阵块kS(in,jm,t)，其中k表示第k个子三维矩阵块，in、jm、t分别表示第k个子三维矩阵块的第in行、第jm列、第t帧的像素值，n＝1,2,…,N，m＝1,2,…,M，t＝1,2,…T，T为三维矩阵S帧的总数量；
在第k个子三维数据块kS(in,jm,t)中，以子三维数据块中心点为圆心，向四周搜索，找到子三维数据块内中心点，记为kS(kiN/2,kjM/2,kt)，其中，kiN/2,kjM/2,kt分别表示第k个子三维数据块中块内极大值像素点所在行数、列数以及帧数；
设置第k个子三维数据块的块内行阈值R_THVk，依次计算子三维数据块内中心点所在帧、所在行的温度点kS(kin',kjM/2,kt)与子三维数据块内中心点kS(kiN/2,kjM/2,kt)的相关度
Figure FDA0003541315390000023
其中kin'表示在第k个子三维数据块中，子三维数据块内中心点最近的第n'个像素点；计算并统计
Figure FDA0003541315390000024
的温度点个数，记为kRSS，作为第k个子三维数据块的块内行步长；
设置第k个子三维数据块的块内列阈值C_THVk，依次计算子三维数据块内中心点所在帧、所在列的温度点kS(kiN/2,kjm',kt)与子三维数据块内中心点kS(kiN/2,kjM/2,kt)的相关度
Figure FDA0003541315390000025
其中kjm'表示在第k个子三维数据块中，距子三维数据块内中心点kS(kiN/2,kjM/2,kt)最近的第m'个像素点；计算并统计
Figure FDA0003541315390000026
的温度点个数，记为kCSS，作为第k个子三维数据块的块内列步长。
4.如权利要求1所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤五分块分步长选取瞬态热响应的具体步骤包括：
步骤S51、根据步骤三统计的M、N个像素值对三维矩阵进行分块，得到K个大小为N×M×T的子三维数据块，其中kS(in,jm,t)表示第k个子三维数据块内的第in行、第jm列像素的瞬态热响应，其中t＝1,2,…,T，T为原三维矩阵S帧的总数量；
步骤S52、对于每个子三维数据块内的像素点kS(in,jm,t)，设置阈值DD，初始化集合编号g＝1，初始化像素点位置in＝1,jm＝1，并将块内极大值kS(kizz,kjzz,ktzz)对应的瞬态热响应kS(kizz,kjzz,t)，t＝1,2,...,T，存储在集合X(g)中；然后计算子三维数据块内的像素点kS(in,jm,t)中的位于in行，jm列的瞬态热响应与集合X(g)间的相关度Rei,j，并判断：
如果Rei,j<DD，则g＝g+1，并将瞬态热响应kS(in,jm,t)作为一个新特征存储在集合X(g)中；否则，令in＝in+kRSS，继续计算下一个瞬态热响应kS(in,jm,t)，t＝1,2,...,T与集合X(g)的相关度；如果in＞N，则令in＝in-N，jm＝jm+kCSS，即变化到第jm+kCSS列进行计算，如果jm＞M,则第k个子三维矩阵的瞬态热响应选取完毕，k＝k+1；其中，N、M分别为第k个子三维矩阵数据块kS(:,:,:)的行数、列数。
5.如权利要求1所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤六采用无监督聚类算法对选取的瞬态热响应进行分类的具体方法为：将步骤五选取的所有K个数据块的所有集合X(g)即瞬态热响应采用模糊C均值聚类FCM算法分为L类，得到每个瞬态热响应所属的类别，具体包括以下步骤：
步骤S61、设置聚类数目L，初始化迭代次数c＝0，设定终止迭代条件阈值ε；
步骤S62、利用公式
Figure FDA0003541315390000031
计算隶属度矩阵U；
其中，i′＝1,2,...,L，c∈L，n'dk'＝||xk'-i'V||，n'＝i',j'，n'dk'表示第k'个像素点与第n′聚类中心i'V的欧氏距离，xk'表示第k'个像素点的坐标；τ为常数；i'uk'表示第k'个像素点隶属于第i'类的程度；
步骤S63、更新聚类中心i'V
Figure FDA0003541315390000041
其中，
Figure FDA0003541315390000042
表示第k'个像素点的热响应值；
步骤S64、如果迭代次数到达最大值L或者前后两次聚类中心之差绝对值小于ε，则算法结束，并输出隶属度矩阵U和聚类中心V，再进入步骤S65；否则，令c＝c+1，返回步骤S62；
步骤S65、利用隶属度最大化准则对所有像素点去模糊化，得到每个像素点所属类别，即Mk'＝argi'max(i'uk')。
6.如权利要求1所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤七基于动态多目标选取每类瞬态热响应的代表构成矩阵Y的具体步骤包括：
步骤S71、在第m+1次外部环境下，对第i'个类瞬态热响应选代表时，定义多目标函数：
Figure FDA0003541315390000043
其中，i′＝1,2,...,L，
Figure FDA0003541315390000044
为第m+1次外部环境下第i'类瞬态热响应选出的一个瞬态热响应
Figure FDA0003541315390000045
的类内欧氏距离，表示为：
Figure FDA0003541315390000046
Figure FDA0003541315390000047
为第i'类瞬态热响应选出的一个瞬态热响应
Figure FDA0003541315390000048
的L-1个类间欧氏距离，由计算出的L-1个类间欧氏距离
Figure FDA0003541315390000049
重新编号组成，
Figure FDA00035413153900000410
表示为：
Figure FDA00035413153900000411
Figure FDA00035413153900000412
为瞬态热响应
Figure FDA00035413153900000413
在第t时刻的像素值即温度值，
Figure FDA00035413153900000414
为第i'类瞬态热响应聚类中心在第t时刻的像素值即温度值，
Figure FDA00035413153900000415
为第j'类瞬态热响应聚类中心在第t时刻的像素值即温度值；
步骤S72、第m-1次和第m次环境下得到的多目标函数近似前沿解集分别为
Figure FDA0003541315390000051
和
Figure FDA0003541315390000052
对应的种群瞬态热响应，解集分别为
Figure FDA0003541315390000053
和
Figure FDA0003541315390000054
其数目分别为
Figure FDA0003541315390000055
和
Figure FDA0003541315390000056
在环境变化后，根据第m-1次和第m次环境的历史信息，预测计算第m+1次环境下的近似前沿解集的初始化种群瞬态热响应，步骤如下：
步骤S721、
Figure FDA0003541315390000057
是从
Figure FDA0003541315390000058
解集中随机选择NE个瞬态热响应
Figure FDA0003541315390000059
构成的瞬态热响应集，n'＝1,2,..,NE，计算
Figure FDA00035413153900000510
集中代表瞬态热响应的个数W，用于获得第m+1次环境下多方向预测集：
Figure FDA00035413153900000511
Figure FDA00035413153900000512
Figure FDA00035413153900000513
其中，
Figure FDA00035413153900000514
是从
Figure FDA00035413153900000515
解集中随机选择NE个瞬态热响应
Figure FDA00035413153900000516
构成的瞬态热响应集，n'＝1,2,..,NE；
步骤S722、选择W+1个代表瞬态热响应，其具体步骤包括：
步骤A、代表瞬态热响应构成的PS多方向预测集由两部分组成：
一是
Figure FDA00035413153900000517
解集瞬态热响应的中心，记为
Figure FDA00035413153900000518
Figure FDA00035413153900000519
其中，
Figure FDA00035413153900000520
为解集
Figure FDA00035413153900000521
中第n个瞬态热响应；
二是基于完全自适应谱聚类算法得到的能够充分描述当前PS形状和多样性的W个代表瞬态热响应集
Figure FDA00035413153900000522
此时，集合
Figure FDA00035413153900000523
中代表瞬态热响应的个数
Figure FDA00035413153900000524
为W个；
步骤B、用完全自适应谱聚类算法SC_SD聚类算法将解集
Figure FDA00035413153900000525
中的瞬态热响应聚类成聚类集合
Figure FDA00035413153900000526
其具体步骤包括：
步骤B1、为消除瞬态热响应不同属性的不同量纲对聚类结果的影响，采用最大最小标准化方法对瞬态热响应点进行标准化：
Figure FDA0003541315390000061
其中
Figure FDA0003541315390000062
NE表示解集
Figure FDA0003541315390000063
中第n个瞬态热响应的第d个时刻的温度值取值d∈{1,2,…,T}，
Figure FDA0003541315390000064
表示解集
Figure FDA0003541315390000065
中所有瞬态热响应样本在第d个时刻的温度值取值；
步骤B2、计算标准化后每个瞬态热响应样本的全局标准差，以瞬态热响应n为例，i′＝1,2,…,NE，其全局标准差公式如下：
Figure FDA0003541315390000066
其中
Figure FDA0003541315390000067
为标准化后瞬态热响应样本n和标准化后瞬态热响应样本j之间的欧氏距离
Figure FDA0003541315390000068
NE为解集
Figure FDA0003541315390000069
中瞬态热响应样本总数；
步骤B3、以stdn为瞬态热响应样本n的领域半径，统计瞬态热响应样本n邻域内的瞬态热响应样本数Cn；利用瞬态热响应样本n在相应邻域的局部标准差获得瞬态热响应样本n对应的自适应局部尺度参数σn：
Figure FDA00035413153900000610
其中Cn为瞬态热响应样本n对应领域半径内瞬态热响应样本的数量；
步骤B4、计算基于自适应的局部尺度参数的亲和矩阵M，亲和矩阵M为由其第n行，j列的元素mnj,n,j∈{1,2,…,NE}构成的NE×NE的方阵，其元素mnj由以下公式确定：
Figure FDA00035413153900000611
步骤B5、计算亲和矩阵M的度矩阵D：
Figure FDA00035413153900000612
度矩阵D为
Figure FDA00035413153900000613
是由dn,n∈{1,2,…,NE}组成的大小为NE×NE对角矩阵，mnj为亲和矩阵的第n行、第j列元素；
步骤B6、计算拉普拉斯矩阵
Figure FDA00035413153900000614
步骤B7、计算拉普拉斯矩阵L的特征值，令λ1≥λ2≥…≥λW≥0为拉普拉斯矩阵L的前W个最大特征值，找到其对应的特征向量ξ1,ξ2,…,ξW，构造矩阵
Figure FDA0003541315390000071
其中ξn为列向量，n∈{1,2,…,W}，令vij表示矩阵V中的第i行、第j列个元素；
步骤B8、按行对矩阵V进行归一化，得到矩阵U，U的第i行，第j列元素值uij由以下公式确定：
Figure FDA0003541315390000072
其中，uij、vik分别表示矩阵U中的第i行、第j列个元素和矩阵V中的第i行、第k列个元素；
步骤B9、将矩阵U每行uN,N∈{1,2,…,NE}，当作一个样本，采用SD_K-medoids算法聚类，其具体步骤包括：
步骤B91、根据下面公式计算矩阵U中每个样本的密度：
Figure FDA0003541315390000073
其中uN,uj,ut分别表示矩阵U的第N，第j，第t行，选取样本密度值最小的样本作为初始聚类中心，将剩余样本作为集合LU，需聚类成的总簇数为W，设置聚类标记数w，置w＝1；
步骤B92、根据公式：
Figure FDA0003541315390000074
计算各个样本的方差VaN,N∈{1,2,…,NE}；
步骤B93、在剩余样本集LU中选择方差最小的样本，记为uv_s(w)，加入到初始聚类中心集合C中，即C＝C∪{uv_s(w)}；
步骤B94、根据公式：
Figure FDA0003541315390000075
计算样本uv_s(w)的标准差Sv_s(w)，以Sv_s(w)为样本uv_s(w)的领域半径rv_s(w)，即
Figure FDA0003541315390000081
计算样本uv_s(w)的领域NBv_s(w)，使其满足：
Figure FDA00035413153900000818
步骤B95、从剩余样本集合LU中去掉NBv_s(w)中的样本，即LU＝LU-NBv_s(w)，w＝w+1，若
Figure FDA00035413153900000819
跳转至步骤B93；
步骤B10、根据矩阵U的聚类结果，将原始瞬态热响应样本分配类簇，即当矩阵U的第n行属于第j个类簇时，将原始瞬态热响应样本n分配到第j个类簇；输出聚类结果
Figure FDA0003541315390000082
步骤C、计算聚类结果中每个类别的聚类中心：
Figure FDA0003541315390000083
其中
Figure FDA0003541315390000084
为第h个聚类结果
Figure FDA0003541315390000085
的第k个代表瞬态热响应，
Figure FDA0003541315390000086
为第h个聚类结果中所包含的瞬态热响应总数；
基于聚类结果
Figure FDA0003541315390000087
从每个类别中选择能够充分描述当前PS形状和多样性的一个代表瞬态热响应
Figure FDA0003541315390000088
Figure FDA0003541315390000089
其中，
Figure FDA00035413153900000810
为第h个聚类结果
Figure FDA00035413153900000811
的聚类中心；这样一共得到W个能够充分描述当前PS形状和多样性的代表瞬态热响应，加上
Figure FDA00035413153900000812
解集整体瞬态热响应的中心，构成W+1个PS多方向预测集；
Figure FDA00035413153900000813
步骤S723、根据第m-1次和第m次环境的PS多方向预测集
Figure FDA00035413153900000814
和
Figure FDA00035413153900000815
其中，
Figure FDA00035413153900000816
按照步骤S721、步骤S722的方法获得，W'是
Figure FDA00035413153900000817
集中代表瞬态热响应的个数；
计算预测方向
Figure FDA0003541315390000091
Figure FDA0003541315390000092
其中，
Figure FDA0003541315390000093
是PS多方向预测集
Figure FDA0003541315390000094
中与
Figure FDA0003541315390000095
距离最近的瞬态热响应；
步骤S724、迭代次数g'＝0时，第m+1次环境下的近似前沿解集的初始化种群瞬态热响应个数为Np，其中，
Figure FDA0003541315390000096
个初始种群瞬态热响应在取值范围内随机生成，
Figure FDA0003541315390000097
个初始种群瞬态热响应由下式预测获得：
Figure FDA0003541315390000098
其中，hn为瞬态热响应
Figure FDA0003541315390000099
所属聚类结果
Figure FDA00035413153900000910
的序号，
Figure FDA00035413153900000911
是一个服从均值为0，方差为
Figure FDA00035413153900000912
的正态分布的随机数，方差
Figure FDA00035413153900000913
的计算公式为：
Figure FDA00035413153900000914
步骤S73、初始化相关参数
初始化迭代次数g'＝0，一组均匀分布的权重向量
Figure FDA00035413153900000915
其中：
Figure FDA00035413153900000916
初始化参考点
Figure FDA00035413153900000917
是函数
Figure FDA00035413153900000918
对应的参考点；
Figure FDA00035413153900000919
最大迭代次数g'max；
初始化每个种群瞬态热响应的进化速度为
Figure FDA00035413153900000920
种群瞬态热响应的全局最优和局部最优满足
Figure FDA00035413153900000921
步骤S74、利用
Figure FDA00035413153900000922
构造基于惩罚项的边界交叉法下每个种群瞬态热响应的动态目标函数适应度值
Figure FDA00035413153900000923
Figure FDA00035413153900000924
其中
Figure FDA00035413153900000925
为预设参数，
Figure FDA00035413153900000926
和
Figure FDA00035413153900000927
的设定如下：
Figure FDA00035413153900000928
Figure FDA0003541315390000101
步骤S75、对n＝1,K,NP：按照粒子群算法更新速度
Figure FDA0003541315390000102
和种群瞬态热响应
Figure FDA0003541315390000103
按照多目标优化算法比较
Figure FDA0003541315390000104
更新全局最优
Figure FDA0003541315390000105
局部最优
Figure FDA0003541315390000106
和参考点
Figure FDA0003541315390000107
从
Figure FDA0003541315390000108
中保留支配
Figure FDA0003541315390000109
的解向量，移除所有被
Figure FDA00035413153900001010
支配的解向量，如果
Figure FDA00035413153900001011
中的向量都不支配
Figure FDA00035413153900001012
将
Figure FDA00035413153900001013
加入
Figure FDA00035413153900001014
同时n＝n+1，若n≤NP，则g'＝g'+1；
步骤S76、进化终止判断：若g'≤g'max，则重复步骤S75，若g'＞g'max，则得到第i'类温度瞬态热响应最终的前沿近似解集
Figure FDA00035413153900001015
步骤S77、基于加权隶属度方案从前沿近似解集
Figure FDA00035413153900001016
选出第i'类瞬态热响应的代表i'REP；
根据下列公式计算前沿近似解集
Figure FDA00035413153900001017
中第k个解对第l个目标函数的隶属度：
Figure FDA00035413153900001018
其中，Fl为第l个目标函数的值，Fl max、Fl min分别为对应目标函数的最大和最小值；
为目标函数设置权重λl，l＝1,2,…,L，计算前沿近似解集的隶属度加权值，取所得最大值对应的前沿解即为第i'类瞬态热响应的代表i'REP，公式如下：
Figure FDA00035413153900001019
其中，
Figure FDA00035413153900001020
为前沿近似解集
Figure FDA00035413153900001021
所包含的解集个数，L为目标函数个数，
Figure FDA00035413153900001022
为隶属度函数值；
所有L类的瞬态热响应代表按列放置，一列为T个时刻的像素值即温度值，构成一个T×L的矩阵Y。
7.如权利要求1所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤八计算瞬态热响应代表之间的斯皮尔曼相关系数的具体步骤包括：
以i'REP和j'REP，i′,j′＝1,2…,L表示任意两个瞬态热响应代表，将i'REP和j'REP对应的各个时刻的温度元素值i′REPt,t＝1,2,L,T和j′REPt,t＝1,2,L,T转换为其在各自瞬态热响应向量总体时刻温度值中的降序排名，将温度值最高的元素转换为1，温度值最低的元素转换为T，其余温度值元素按大小顺序转换为排名并记为Ra(i'REPt)和Ra(j'REPt)；根据公式：
Figure FDA0003541315390000111
计算两个瞬态热响应代表之间对应时刻元素值的差异Da；最后根据公式：
Figure FDA0003541315390000112
计算两个瞬态热响应代表i'REP和j'REP之间的斯皮尔曼相关系数；
设定相关性阈值θ，两两之间比较斯皮尔曼相关系数，保留相关性最小的瞬态热响应代表，若Rs(i'REPt,j'REPt)＜θ，则保留第i，j类瞬态热响应代表，否则剔除其中一类，得到L'类瞬态热响应代表；所有L'类的瞬态热响应代表按列放置，一列为T个时刻的像素值即温度值，构成一个T×L'的矩阵Y'；
将三维矩阵S中的每一帧从第一列开始，将后一列接在前一列的末尾，构成新的一列，得到T帧对应的T列像素值，然后，依据时间先后，将T列像素值依次放置，构成I×J行、T列二维图像矩阵O，用矩阵Y对二维矩阵O进行线性变换，即：
Figure FDA0003541315390000113
得到二维图像矩阵R，其中，
Figure FDA0003541315390000114
为L'×T矩阵，是矩阵Y'的伪逆矩阵，OT二维图像矩阵O的转置矩阵，得到的二维图像矩阵R为L'行、I×J列；
二维图像矩阵R的每一行，按J列依次进行截取，截取的J列按行依次放置，构成一张I×J二维图像，这样L'行，得到L'张I×J二维图像，这些图片都包含了缺陷区域，为方便缺陷轮廓提取，选择缺陷区域和非缺陷区域像素值差距最大的一张二维图像，并记为f(x,y)。
8.如权利要求7所述的在轨航天器微小撞击损伤可视化检测评估方法，其特征在于，所述步骤九利用区域卷积神经网络R-CNN对步骤八中的二维图像f(x,y)进行图像中的缺陷检测和缺陷区域分割，从而实现缺陷部分的定位识别和特征提取。
CN202010920381.0A 2020-09-04 2020-09-04 一种在轨航天器微小撞击损伤可视化检测评估方法 Active CN112016627B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010920381.0A CN112016627B (zh) 2020-09-04 2020-09-04 一种在轨航天器微小撞击损伤可视化检测评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010920381.0A CN112016627B (zh) 2020-09-04 2020-09-04 一种在轨航天器微小撞击损伤可视化检测评估方法

Publications (2)

Publication Number Publication Date
CN112016627A CN112016627A (zh) 2020-12-01
CN112016627B true CN112016627B (zh) 2022-04-19

Family

ID=73515350

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010920381.0A Active CN112016627B (zh) 2020-09-04 2020-09-04 一种在轨航天器微小撞击损伤可视化检测评估方法

Country Status (1)

Country Link
CN (1) CN112016627B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112818822B (zh) * 2021-01-28 2022-05-06 中国空气动力研究与发展中心超高速空气动力研究所 一种航空航天复合材料损伤区域自动识别方法
CN112784847B (zh) * 2021-01-28 2022-03-04 中国空气动力研究与发展中心超高速空气动力研究所 一种超高速撞击损伤检测图像的分割识别方法
CN112944104B (zh) * 2021-03-03 2022-09-30 杭州申昊科技股份有限公司 一种检测缺陷的管道机器人及其控制方法与控制系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109767437A (zh) * 2019-01-09 2019-05-17 电子科技大学 基于k均值动态多目标的红外热图像缺陷特征提取方法
CN109767438A (zh) * 2019-01-09 2019-05-17 电子科技大学 一种基于动态多目标优化的红外热图像缺陷特征识别方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109767437A (zh) * 2019-01-09 2019-05-17 电子科技大学 基于k均值动态多目标的红外热图像缺陷特征提取方法
CN109767438A (zh) * 2019-01-09 2019-05-17 电子科技大学 一种基于动态多目标优化的红外热图像缺陷特征识别方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
分解多目标算法在航天器损伤红外检测中的应用;薛婷;《中国优秀硕士学位论文全文数据库 工程科技II辑》;20200715(第07期);第C031-193页 *

Also Published As

Publication number Publication date
CN112016627A (zh) 2020-12-01

Similar Documents

Publication Publication Date Title
CN112016628B (zh) 一种基于动态多目标预测的空间碎片撞击损伤判读方法
CN112016627B (zh) 一种在轨航天器微小撞击损伤可视化检测评估方法
CN112037211B (zh) 一种用于微小空间碎片撞击事件动态监测的损伤特征识别方法
US11475360B2 (en) System and method for relational time series learning with the aid of a digital computer
CN109767438B (zh) 一种基于动态多目标优化的红外热图像缺陷特征识别方法
Roustant et al. Group kernels for Gaussian process metamodels with categorical inputs
US20150185270A1 (en) Method for recognizing transformer partial discharge pattern based on singular value decomposition algorithm
CN109767437B (zh) 基于k均值动态多目标的红外热图像缺陷特征提取方法
CN102831474A (zh) 基于量子粒子群优化改进的模糊c-均值聚类方法
CN115131618B (zh) 基于因果推理的半监督图像分类方法
CN106228197A (zh) 一种基于自适应极限学习机的卫星图像云量识别方法
Yi et al. White blood cell image segmentation using on-line trained neural network
CN114595884A (zh) 一种遗传智能优化的神经网络风力发电设备温度预测方法
Wang et al. Particle swarm optimization for fuzzy c-means clustering
CN104463207B (zh) 知识自编码网络及其极化sar影像地物分类方法
CN115393631A (zh) 基于贝叶斯层图卷积神经网络的高光谱图像分类方法
CN113177576A (zh) 一种用于目标检测的多示例主动学习方法
Pourkamali-Anaraki et al. Neural networks and imbalanced learning for data-driven scientific computing with uncertainties
Portman et al. Sampling algorithms for validation of supervised learning models for Ising-like systems
CN110263808B (zh) 一种基于lstm网络和注意力机制的图像情感分类方法
CN105719006B (zh) 基于流特征的因果结构学习方法
CN111897735B (zh) 基于粒子群和多标准的模糊测试用例选择方法和装置
CN105678798A (zh) 一种结合局部空间信息的多目标模糊聚类图像分割方法
CN115640868A (zh) 一种新建风电场少数据风电功率短期预测方法
CN111626332A (zh) 一种基于图卷积极限学习机的快速半监督分类方法

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