CN114913218A - 一种基于极小模网络的烟幕红外遮蔽面积测量方法 - Google Patents

一种基于极小模网络的烟幕红外遮蔽面积测量方法 Download PDF

Info

Publication number
CN114913218A
CN114913218A CN202210475105.7A CN202210475105A CN114913218A CN 114913218 A CN114913218 A CN 114913218A CN 202210475105 A CN202210475105 A CN 202210475105A CN 114913218 A CN114913218 A CN 114913218A
Authority
CN
China
Prior art keywords
smoke screen
area
infrared
point source
point
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
CN202210475105.7A
Other languages
English (en)
Other versions
CN114913218B (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.)
Unit 63891 Of Pla
Original Assignee
Unit 63891 Of Pla
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 Unit 63891 Of Pla filed Critical Unit 63891 Of Pla
Priority to CN202210475105.7A priority Critical patent/CN114913218B/zh
Publication of CN114913218A publication Critical patent/CN114913218A/zh
Application granted granted Critical
Publication of CN114913218B publication Critical patent/CN114913218B/zh
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/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/08—Projecting images onto non-planar surfaces, e.g. geodetic screens
    • 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/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Geometry (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开一种基于极小模网络的烟幕红外遮蔽面积测量方法，其包括以下步骤：对获取的烟幕下的红外原始图像进行滤波；对预处理后的图像进行分割；提取出的点源形状和位置，计算所有点源的重心位置坐标，作为每个点源的位置坐标；计算红外点源位置的遮蔽率作为插值节点的函数值；将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分三角单元；构造曲线网络建立求解方程组；构造平均平面估计各节点处一阶偏导数，计算出插值节点的一阶偏导数；根据每个三角单元的插值系数和插值公式逐点插值计算烟幕遮蔽率，最后根据所有点的烟幕遮蔽率统计烟幕遮蔽面积。本发明能够极大的减少插值算法的复杂度，提高插值精度，减少数据处理时间。

Description

一种基于极小模网络的烟幕红外遮蔽面积测量方法
技术领域
本发明属于烟幕干扰评估技术领域，尤其是涉及一种基于极小模网络的烟幕红外遮蔽面积测量方法。
背景技术
在军民融合、大数据分析等领域，随着现代科学技术的发展，新型可见光与红外波段侦察设备的广泛应用，使侦察与监视的水平和能力有了极大提高，对具有高效遮蔽能力的烟幕需求越来越迫切。烟幕遮蔽是一种快速、高效遮蔽目标的手段。因此，烟幕有效遮蔽性能的研究对于烟幕红外有效遮蔽面积评估具有重大意义。
对红外烟幕遮敝效果的测量主要是通过对烟幕透过率、遮蔽面积、形成时间和持续时间等的测量来定量描述的。目前，使用最多的是用凝视型红外成像热像仪的多点测量法。
现有技术中，在使用摄像法测量烟幕的红外遮蔽性能时，通常的方法是采用布设红外点源阵列，在红外点源阵列前方释放烟幕，并使用红外热像仪分别测量烟幕施放前、施放时红外点源阵列的红外辐射强度变化，得到点源位置处的遮蔽率(透过率)，再通过插值的方法得出烟幕遮蔽率的空间分布，最后利用烟幕遮蔽率阈值确定有效遮蔽面积。在布设红外点源阵列时，应使其在红外热像仪上的成像间距分布均匀，以便于插值。但在实际试验中，由于试验条件限制，很难使红外点源阵列的成像间距均匀分布，导致插值算法复杂、插值精度差，得出的有效遮蔽面积结果不精确，数据处理时间长，使得试验周期长、成本高、效率低。
发明内容
针对现有技术中在使用摄像法测量烟幕的红外遮蔽性能时，实际试验中存在红外点源阵列在红外热像仪上的成像间距分布不均匀的问题，本发明的目的是提供一种基于极小模网络的烟幕红外遮蔽面积测量方法，其利用样条函数能量极小原理构建插值函数进行插值，首先基于三角剖分构建曲线网络，在此基础上构造出双自变量三次自然样条插值函数，并通过迭代的方法计算插值节点处的偏导数，最后代入插值函数进行插值，得出插值结果。
为实现上述发明目的，本发明采用如下技术方案：
一种基于极小模网络的烟幕红外遮蔽面积测量方法，其包括以下步骤：
S1、图像预处理：采用数学形态学中的顶帽变换算法对获取的烟幕下的红外原始图像进行滤波，增强红外点源灰度值，抑制背景；
S2、图像分割：使用粗分割、精分割两步分割法对步骤S1中预处理后的图像进行分割；操作为：首先，使用全局阈值分割算法提取出每个红外点源的粗略位置；然后在每个点源位置处向外扩张一定距离，划分为局部待分割区域；最后，在每个局部待分割区域内使用最大类间方差法依次对每个点源进行分割，提取点源形状和位置；
S3、提取点源位置坐标：根据步骤S2中提取出的点源形状和位置，计算所有点源的重心位置坐标，作为每个点源的位置坐标；
S4、计算点源的烟幕遮蔽率或透过率：根据步骤S2中的图像分割结果中的目标区域和背景区域施放烟幕前、施放烟幕后的灰度变化，计算每个红外点源位置处的烟幕遮蔽率或透过率；
S5、三角剖分：根据步骤S2中提取出的点源位置坐标，将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分为基本三角形单元，三角形单元的每个顶点依次对应一个点源，所有点源形成的三角形网络构成插值曲面片；
S6、构造曲线网络建立求解方程组：在步骤S5中三角剖分的基础上，建立定义于三角剖分中所有边和节点函数值及一阶偏导数的插值曲线网络，将所有插值曲线网络联立建立一个求解方程组；
S7、迭代计算插值节点偏导数：构造平均平面估计各节点处一阶偏导数，将估计出的偏导数作为迭代法求解时的初始值，迭代求解方程组，计算得出各节点处一阶偏导数；
S8、坐标变换：将待插值点直角坐标(xi,yi)进行坐标变换，变换为局部重心坐标(b1,b2,b3)；
S9、计算遮蔽面积：将节点函数值和偏导数代入插值函数，遍历所有的待插值点，逐点插值，得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积；最后根据几何投影关系，计算出实际烟幕遮蔽面积。
进一步地，上述的步骤S1中，使用顶帽变换算法对图像预处理，顶帽变换算法的计算公式如下：
使用结构元素b对图像f(x,y)的灰度膨胀记为
Figure BDA0003625068440000031
定义为：
Figure BDA0003625068440000032
式中：Db是结构元素b的定义域；
使用结构元素b对图像f(x,y)的灰度腐蚀记为
Figure BDA0003625068440000033
定义为：
Figure BDA0003625068440000034
灰度腐蚀是一个局部最小值算子；
使用结构元素b对图像f(x,y)的开运算记为
Figure BDA0003625068440000036
，定义为：
Figure BDA0003625068440000035
即f先由b腐蚀，再由b膨胀；
使用结构元素b对图像f(x,y)的闭运算记为f·b，定义为：
Figure BDA0003625068440000041
即f先由b膨胀，再由b腐蚀；
灰度图像的顶帽变换定义为f减去其开运算结果，即：
Figure BDA0003625068440000042
顶帽变换算法用于暗背景上的亮目标，用来校正不均匀背景的影响。
进一步地，上述的步骤S2中，使用最大类间方差法对局部待分割区域进行分割，最大化类间方差
Figure BDA0003625068440000043
计算公式如下：
Figure BDA0003625068440000044
其中，
Figure BDA0003625068440000045
Figure BDA0003625068440000046
Figure BDA0003625068440000047
Figure BDA0003625068440000048
Figure BDA0003625068440000049
Figure BDA00036250684400000410
式中：nq是待分割区域Z中第q级灰度的像素数；N是待分割区域Z中的像素总数；L是待分割区域Z中的灰度级；
求出待分割区域所有可能灰度级的类间方差后，用下式求得分割阈值t：
Figure BDA00036250684400000411
分割出的目标区域表示为待分割区域Z中灰度值f(x,y)≥t的区域，即
T＝{(x,y)|f(x,y)≥t,(x,y)∈Z}。
进一步地，上述的步骤S3中，计算点源重心坐标的算法为：
Figure BDA0003625068440000051
其中，n为目标区域T的像素总数。
进一步地，上述的步骤S4中，点源的烟幕透过率计算公式如下：
Figure BDA0003625068440000052
点源的烟幕遮蔽率计算公式如下：
Figure BDA0003625068440000053
其中，
Figure BDA0003625068440000054
分别为施放烟幕前、施放烟幕后目标区域的平均灰度值，
Figure BDA0003625068440000055
Figure BDA0003625068440000056
分别为施放烟幕前、施放烟幕后背景区域的平均灰度值。
进一步地，上述的步骤S5中，划分基本三角形单元的方法，步骤为：
首先根据空间位置关系将所有点源按每列、每行点源连接划分为基本四边形单元；已知任意1个四边形能够沿对角顶点划分为2个三角形且有两种划分方法，分别按两种方法进行划分并计算2个三角形的最小内角，取最小内角最大的划分方法进行划分。
进一步地，上述的步骤S6中，构造曲线网络建立求解方程组的方法，步骤为：令S(x,y)为定义在三角形网络构成插值曲面片上的三次自然样条插值函数，将所有插值曲线网络联立，建立一个求解方程组；
Figure BDA0003625068440000057
式中，Ni是Ne中以节点Vi为端点的部分；Sx(Vi)、Sy(Vi)分别是节点Vi处插值函数x、y方向的一阶偏导数；S(Vi)＝zi为节点Vi处的遮蔽率或透过率；
若总红外点源数为n，即n个节点，每个节点均能够根据式(16)列出2个方程，共有2n个方程2n个未知数，这2n个方程相互耦合，必须联立形成方程组联合求解，用直接迭代法进行求解。
进一步地，上述的步骤S7中，迭代计算插值节点偏导数的方法，步骤为：
首先，构造平均平面对Sx(Vi)、Sy(Vi)进行初始化，平均平面表示为
Amix+Bmiy+Cmiz+Dmi＝0 (17)
其中
Figure BDA0003625068440000061
Figure BDA0003625068440000062
Figure BDA0003625068440000063
Ai＝(yi2-yi3)zi1+(yi3-yi1)zi2+(yi1-yi2)zi3 (21)
Bi＝(xi3-xi2)zi1+(xi1-xi3)zi2+(xi2-xi1)zi3 (22)
Ci＝(yi2-yi3)xi1+(yi3-yi1)xi2+(yi1-yi2)xi3 (23)
式中，Tijk表示三角形网络中以Vi、Vj、Vk为顶点的三角形单元；(xi1,yi1,zi1)、(xi1,yi1,zi1)、(xi1,yi1,zi1)为第i个三角单元的顶点坐标及函数值；M为以Vi为顶点的三角形个数；则有
Figure BDA0003625068440000064
Figure BDA0003625068440000065
由式(16)推导出Sx(Vi)、Sy(Vi)的表达式，用于直接迭代求解：
Figure BDA0003625068440000066
其中
Figure BDA0003625068440000071
Figure BDA0003625068440000072
Figure BDA0003625068440000073
将初始值
Figure BDA0003625068440000074
代入式(26)右边，求解得到
Figure BDA0003625068440000075
继续代入式(26)右边，求解得到
Figure BDA0003625068440000076
如此往复，直至计算结果满足精度要求；若令
Figure BDA0003625068440000077
其中，
Figure BDA0003625068440000078
为第k次迭代节点Vi处的一阶偏导数，则当||dk-dk-1||≤ε时停止迭代，ε为指定的迭代停止阈值，此时即认为求解精度满足要求，即
Figure BDA0003625068440000079
为满足精度的解。
进一步地，上述的步骤S8中，坐标变换的方法，步骤为：
△V1V2V3中，任意一点P(x,y)的位置由三角形的3个顶点确定的重心坐标系表示，即用下面三个比值来规定：
Figure BDA00036250684400000710
其中，b1、b2、b3为点P在重心坐标系中的坐标；S为△V1V2V3的面积；S1、S2、S3分别为△V1V2V4、△V4V2V3、△V1V4V3的面积；获知
b1+b2+b3＝1 (33)
Figure BDA00036250684400000711
Figure BDA0003625068440000081
根据(32)式～式(35)，得到重心坐标关系式
Figure BDA0003625068440000082
进一步地，上述的步骤S9中，在每个三角形单元中，插值函数为
Figure BDA0003625068440000083
Figure BDA0003625068440000084
Figure BDA0003625068440000085
其中，(i,j,k)＝{(1,2,3),(2,3,1),(3,1,2)}；
遍历所有的待插值点，逐点插值，即得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积；最后根据几何投影关系，计算出实际烟幕遮蔽面积；由烟幕遮蔽面积在红外热像仪上的投影面积，计算实际烟幕遮蔽面积的方法为
S＝NSp (42)
其中，S为实际烟幕遮蔽面积；N为图像中遮蔽率大于遮蔽率阈值的像素数；Sp为根据几何投影关系计算出的单个像素代表的实际面积。
由于采用如上所述的技术方案，本发明具有如下优越性：
该基于极小模网络的烟幕红外遮蔽面积测量方法，其利用样条函数能量极小原理构建插值函数进行插值，能够有效避免红外点源阵列分布不均匀造成数据点分布不均匀，进而导致插值算法复杂，插值精度差、计算速度慢等现象；能够极大的减少插值算法的复杂度，提高插值精度，减少数据处理时间，缩短试验周期，降低成本。
附图说明
图1是摄像法的典型试验布局示意图；
图2是理想情况下，红外点源阵列在红外热像仪上的投影分布示意图；
图3是本发明一实施例提供的试验中，红外点源阵列在红外热像仪上的实际投影分布图；
图4对图3进行三角剖分结果图；
图5是重心坐标变换示意图；
图6是本发明一实施例提供的一帧图像的最终插值结果及等位线图。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步详细说明。
如图1～6所示，基于极小模网络的烟幕红外遮蔽面积测量方法，其包括以下步骤：
S1、图像预处理：由于弥散作用，红外点源在红外热像仪上所成的像不是一个点，而是一个近似于圆形的小目标，因此需使用图像分割算法提取目标，分割之前需对图像进行预处理，以避免误分、漏分；
已知，形态学的基本运算有4个：膨胀(或扩张)、腐蚀(或侵蚀)、开启和闭合，基于上述基本运算能够推导和组合成各种数学形态学实用算法，用它们进行图像形状和结构的分析及处理，包括图像分割、特征抽取、边缘检测、图像滤波、图像增强和恢复等。
本发明采用数学形态学中的顶帽变换算法对获取的烟幕下的红外原始图像进行滤波，增强红外点源灰度值，抑制背景，以避免分割时出现误分、漏分；顶帽变换算法的计算公式如下：
使用结构元素b对图像f(x,y)的灰度膨胀记为
Figure BDA0003625068440000101
定义为：
Figure BDA0003625068440000102
式中：Db是结构元素b的定义域；
使用结构元素b对图像f(x,y)的灰度腐蚀记为
Figure BDA0003625068440000103
定义为：
Figure BDA0003625068440000104
灰度腐蚀是一个局部最小值算子；
使用结构元素b对图像f(x,y)的开运算记为
Figure BDA0003625068440000105
定义为：
Figure BDA0003625068440000106
即，f先由b腐蚀，再由b膨胀；
使用结构元素b对图像f(x,y)的闭运算记为f·b，定义为：
Figure BDA0003625068440000107
即，f先由b膨胀，再由b腐蚀；
灰度图像的顶帽变换定义为f减去其开运算结果，即：
Figure BDA0003625068440000108
顶帽变换算法用于暗背景上的亮目标，用来校正不均匀背景的影响；
上述式中，结构元素b的定义类似于邻域，区别的是具有各种不同的形状，如菱形、圆盘、线形、八边形、甚至2个以上的结构元素等，二值矩阵Db决定了哪些位置的元素包括在最大值运算中，可以将膨胀运算看作是一个局部最大值算子；本实施例中，结构元素b取半径为3的圆；
S2、图像分割：使用粗分割、精分割两步分割法对步骤S1中预处理后的图像进行分割；
步骤为：首先，使用全局阈值分割算法提取出每个红外点源的粗略位置；然后在每个点源位置处向外扩张一定距离d，划分为局部待分割区域，本实施例中，取距离d＝5；最后，在每个局部待分割区域内使用最大类间方差法依次对每个点源进行分割，提取点源形状和位置；
最大化类间方差
Figure BDA0003625068440000111
计算公式如下：
Figure BDA0003625068440000112
其中，
Figure BDA0003625068440000113
Figure BDA0003625068440000114
Figure BDA0003625068440000115
Figure BDA0003625068440000116
Figure BDA0003625068440000117
Figure BDA0003625068440000118
式中：nq是待分割区域Z中第q级灰度的像素数；N是待分割区域中的像素总数；求出待分割区域所有可能灰度级的类间方差后，用下式求得分割阈值t：
Figure BDA0003625068440000119
分割出的目标区域表示为待分割区域Z中灰度值f(x,y)≥t的区域，即
T＝{(x,y)|f(x,y)≥t,(x,y)∈Z}；
S3、提取点源位置坐标：根据步骤S2中提取出的点源形状和位置，计算所有点源的重心位置坐标，作为每个点源的位置坐标；
计算点源重心坐标的算法为：
Figure BDA00036250684400001110
其中，n为目标区域T的像素总数；
本实施例中，选取点源重心坐标作为点源位置坐标；
S4、计算点源的烟幕遮蔽率或透过率，根据试验要求确定：步骤S2中，每个局部待分割区域都能够根据最终分割结果划分为目标区域和背景区域两部分；施放烟幕前后，目标区域和背景区域的灰度变化不同，进而能够计算出每个目标位置处的烟幕遮蔽率或透过率；
点源的烟幕透过率计算公式如下：
Figure BDA0003625068440000121
点源的烟幕遮蔽率计算公式如下：
Figure BDA0003625068440000122
其中，
Figure BDA0003625068440000123
分别为施放烟幕前、施放烟幕后目标区域的平均灰度值，
Figure BDA0003625068440000124
分别为施放烟幕前、施放烟幕后背景区域的平均灰度值；
本实施例中，选择计算点源的烟幕遮蔽率；
S5、三角剖分：根据步骤S2中提取出的点源位置坐标，将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分为基本三角形单元，三角形单元的每个顶点依次对应一个点源，所有点源形成的三角形网络构成插值曲面片；
划分基本三角形单元的方法为：首先，根据空间位置关系将所有点源按每列、每行点源连接划分为基本四边形单元；已知任意1个四边形可沿对角顶点划分为2个三角形且有两种划分方法，根据三角剖分的相关理论，三角形单元越接近于等边三角形，三角形单元的质量越好，引起的误差越小，利用这一原则划分基本三角形单元；在划分完成基本四边形单元的基础上，已知任意1个四边形可沿对角顶点划分为2个三角形且有两种划分方法，在划分完成基本四边形单元的基础上，分别按两种方法进行划分并计算2个三角形的最小内角，取最小内角最大的划分方法进行划分(或取最长边与最短边之比最接近于1的划分方法进行划分，二者等价)；
S6、构造曲线网络建立求解方程组：在步骤S5中三角剖分的基础上，建立定义于三角剖分中所有边和节点函数值及一阶偏导数的插值曲线网络，将所有插值曲线网络联立建立一个求解方程组；建立求解方程组的方法，步骤为：
设f(x)∈C2(a,b)是任意插值函数，S(x)是三次自然样条插值函数，则下式成立：
Figure BDA0003625068440000131
式中，等号仅当f(x)＝S(x)时成立；
上式称为三次样条的极小模性质，又称为曲率极小性质，是由力学上的样条能量极小原理推出的；极小模性质表明，三次自然样条插值函数符合能量极小条件；
将三次自然样条函数的极小模性质扩展到三角剖分定义的网络上，即通过求解下式，找出构成曲线网络的各点的偏导数
Figure BDA0003625068440000132
式中，Ne＝{ij:i,j＝1,2,…,n,i≠j}为三角网络中以点Vi、Vj为端点的线段的组合；eij表示点Vi、Vj为端点的线段；dsij表示对应于线段eij的曲线弧长元素；若令
fij(t)＝F[(1-t)Vi+tVj]ij∈Ne,0≤t≤1 (13)
则有
Figure BDA0003625068440000133
式中，||eij||为eij的长度；
由式(12)～式(14)推导出
Figure BDA0003625068440000141
若令S(x,y)为定义在三角形网络构成插值曲面片上的三次自然样条插值函数，则推导出
Figure BDA0003625068440000142
式中，Ni是Ne中以节点Vi为端点的部分；Sx(Vi)、Sy(Vi)分别是节点Vi处插值函数x、y方向的一阶偏导数；S(Vi)＝zi为节点Vi处的遮蔽率或透过率；
本实施例中，总红外点源数为135，即135个节点，每个节点均能够根据式(11)列出2个方程，构成含有270个未知数的270个方程，这270个方程相互耦合，必须联合求解，用直接迭代法进行求解。
S7、迭代计算插值节点偏导数：构造平均平面估计各节点处一阶偏导数，将估计出的偏导数作为迭代法求解时的初始值，迭代求解方程组，计算得出各节点处一阶偏导数；迭代计算插值节点偏导数的方法为：
首先，构造平均平面对Sx(Vi)、Sy(Vi)进行初始化，平均平面表示为
Amix+Bmiy+Cmiz+Dmi＝0 (17)
其中
Figure BDA0003625068440000143
Figure BDA0003625068440000144
Figure BDA0003625068440000145
Ai＝(yi2-yi3)zi1+(yi3-yi1)zi2+(yi1-yi2)zi3 (21)
Bi＝(xi3-xi2)zi1+(xi1-xi3)zi2+(xi2-xi1)zi3 (22)
Ci＝(yi2-yi3)xi1+(yi3-yi1)xi2+(yi1-yi2)xi3 (23)
式中，Tijk表示三角形网络中以Vi、Vj、Vk为顶点的三角形单元；(xi1,yi1,zi1)、(xi1,yi1,zi1)、(xi1,yi1,zi1)为第i个三角单元的顶点坐标及函数值；M为以Vi为顶点的三角形个数；则有
Figure BDA0003625068440000151
Figure BDA0003625068440000152
由式(16)推导出Sx(Vi)、Sy(Vi)的表达式，用于直接迭代求解：
Figure BDA0003625068440000153
其中
Figure BDA0003625068440000154
Figure BDA0003625068440000155
Figure BDA0003625068440000156
将初始值
Figure BDA0003625068440000157
代入式(26)右边，求解得到
Figure BDA0003625068440000158
继续代入式(26)右边，求解得到
Figure BDA0003625068440000159
如此往复，直至计算结果满足精度要求；若令
Figure BDA00036250684400001510
其中，
Figure BDA00036250684400001511
为第k次迭代节点Vi处的一阶偏导数，则当||dk-dk-1||≤ε时停止迭代，ε为指定的迭代停止阈值，本实施例中，取ε＝2-52；此时即认为求解精度满足要求，即
Figure BDA0003625068440000161
为满足精度的解；
S8、坐标变换：将待插值点直角坐标(xi,yi)进行坐标变换，变换为局部重心坐标(b1,b2,b3)；坐标变换的方法为：
△V1V2V3中，任意一点P(x,y)的位置由三角形的3个顶点确定的重心坐标系表示，即用下面三个比值来规定：
Figure BDA0003625068440000162
其中，b1、b2、b3称为点P在重心坐标系中的坐标；S为△V1V2V3的面积；S1、S2、S3分别为△V1V2V4、△V4V2V3、△V1V4V3的面积；获知
b1+b2+b3＝1 (33)
Figure BDA0003625068440000163
Figure BDA0003625068440000164
根据(32)式～式(35)，得到重心坐标关系式
Figure BDA0003625068440000165
也写为
Figure BDA0003625068440000166
其逆变换为
Figure BDA0003625068440000171
S9、计算遮蔽面积：将节点函数值和偏导数代入插值函数，在每个三角形单元中，插值函数为
Figure BDA0003625068440000172
Figure BDA0003625068440000173
Figure BDA0003625068440000174
其中，(i,j,k)＝{(1,2,3),(2,3,1),(3,1,2)}；
遍历所有的待插值点，逐点插值，即得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积；最后根据几何投影关系，计算出实际烟幕遮蔽面积；由烟幕遮蔽面积在红外热像仪上的投影面积，计算实际烟幕遮蔽面积的方法为
S＝NSp (42)
其中，S为实际烟幕遮蔽面积；N为图像中遮蔽率大于遮蔽率阈值的像素数；Sp为根据几何投影关系计算出的单个像素代表的实际面积。
以上所述仅为本发明的较佳实施例，而非对本发明的限制，在不脱离本发明的精神和范围的情况下，凡依本发明申请专利范围所作的均等变化与修饰，皆应属本发明的专利保护范围之内。

Claims (10)

1.一种基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其包括以下步骤：
S1、图像预处理：采用数学形态学中的顶帽变换算法对获取的烟幕下的红外原始图像进行滤波，增强红外点源灰度值，抑制背景；
S2、图像分割：使用粗分割、精分割两步分割法对步骤S1中预处理后的图像进行分割；操作为：首先，使用全局阈值分割算法提取出每个红外点源的粗略位置；然后在每个点源位置处向外扩张一定距离，划分为局部待分割区域；最后，在每个局部待分割区域内使用最大类间方差法依次对每个点源进行分割，提取点源形状和位置；
S3、提取点源位置坐标：根据步骤S2中提取出的点源形状和位置，计算所有点源的重心位置坐标，作为每个点源的位置坐标；
S4、计算点源的烟幕遮蔽率或透过率：根据步骤S2中的图像分割结果中的目标区域和背景区域施放烟幕前、施放烟幕后的灰度变化，计算每个红外点源位置处的烟幕遮蔽率或透过率；
S5、三角剖分：根据步骤S2中提取出的点源位置坐标，将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分为基本三角形单元，三角形单元的每个顶点依次对应一个点源，所有点源形成的三角形网络构成插值曲面片；
S6、构造曲线网络建立求解方程组：在步骤S5中三角剖分的基础上，建立定义于三角剖分中所有边和节点函数值及一阶偏导数的插值曲线网络，将所有插值曲线网络联立建立一个求解方程组；
S7、迭代计算插值节点偏导数：构造平均平面估计各节点处一阶偏导数，将估计出的偏导数作为迭代法求解时的初始值，迭代求解方程组，计算得出各节点处一阶偏导数；
S8、坐标变换：将待插值点直角坐标(xi,yi)进行坐标变换，变换为局部重心坐标(b1,b2,b3)；
S9、计算遮蔽面积：将节点函数值和偏导数代入插值函数，遍历所有的待插值点，逐点插值，得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积；最后根据几何投影关系，计算出实际烟幕遮蔽面积。
2.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S1中，使用顶帽变换算法对图像预处理，顶帽变换算法的计算公式如下：
使用结构元素b对图像f(x,y)的灰度膨胀记为
Figure FDA0003625068430000021
定义为：
Figure FDA0003625068430000022
式中：Db是结构元素b的定义域；
使用结构元素b对图像f(x,y)的灰度腐蚀记为
Figure FDA0003625068430000023
定义为：
Figure FDA0003625068430000024
灰度腐蚀是一个局部最小值算子；
使用结构元素b对图像f(x,y)的开运算记为
Figure FDA0003625068430000025
定义为：
Figure FDA0003625068430000026
即f先由b腐蚀，再由b膨胀；
使用结构元素b对图像f(x,y)的闭运算记为f·b，定义为：
Figure FDA0003625068430000027
即f先由b膨胀，再由b腐蚀；
灰度图像的顶帽变换定义为f减去其开运算结果，即：
Figure FDA0003625068430000031
3.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S2中，使用最大类间方差法对局部待分割区域进行分割，最大化类间方差
Figure FDA0003625068430000032
计算公式如下：
Figure FDA0003625068430000033
其中，
Figure FDA0003625068430000034
Figure FDA0003625068430000035
Figure FDA0003625068430000036
Figure FDA0003625068430000037
Figure FDA0003625068430000038
Figure FDA0003625068430000039
式中：nq是待分割区域Z中第q级灰度的像素数；N是待分割区域Z中的像素总数；L是待分割区域Z中的灰度级；
求出待分割区域所有可能灰度级的类间方差后，用下式求得分割阈值t：
Figure FDA00036250684300000310
分割出的目标区域表示为待分割区域Z中灰度值f(x,y)≥t的区域，即
T＝{(x,y)|f(x,y)≥t,(x,y)∈Z}。
4.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S3中，计算点源重心坐标的算法为：
Figure FDA0003625068430000041
其中，n为目标区域T的像素总数。
5.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S4中，点源的烟幕透过率计算公式如下：
Figure FDA0003625068430000042
点源的烟幕遮蔽率计算公式如下：
Figure FDA0003625068430000043
其中，
Figure FDA0003625068430000044
分别为施放烟幕前、施放烟幕后目标区域的平均灰度值，
Figure FDA0003625068430000045
Figure FDA0003625068430000046
分别为施放烟幕前、施放烟幕后背景区域的平均灰度值。
6.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S5中，划分基本三角形单元的方法，步骤为：
首先根据空间位置关系将所有点源按每列、每行点源连接划分为基本四边形单元；已知任意1个四边形能够沿对角顶点划分为2个三角形且有两种划分方法，分别按两种方法进行划分并计算2个三角形的最小内角，取最小内角最大的划分方法进行划分。
7.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S6中，构造曲线网络建立求解方程组的方法，步骤为：令S(x,y)为定义在三角形网络构成插值曲面片上的三次自然样条插值函数，将所有插值曲线网络联立，建立一个求解方程组：
Figure FDA0003625068430000051
式中，Ni是Ne中以节点Vi为端点的部分；Sx(Vi)、Sy(Vi)分别是节点Vi处插值函数x、y方向的一阶偏导数；S(Vi)＝zi为节点Vi处的遮蔽率或透过率；
若总红外点源数为n，即n个节点，每个节点均能够根据式(16)列出2个方程，共有2n个方程2n个未知数，这2n个方程相互耦合，必须联立形成方程组联合求解，用直接迭代法进行求解。
8.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S7中，迭代计算插值节点偏导数的方法，步骤为：
首先，构造平均平面对Sx(Vi)、Sy(Vi)进行初始化，平均平面表示为
Amix+Bmiy+Cmiz+Dmi＝0 (17)
其中
Figure FDA0003625068430000052
Figure FDA0003625068430000053
Figure FDA0003625068430000054
Ai＝(yi2-yi3)zi1+(yi3-yi1)zi2+(yi1-yi2)zi3 (21)
Bi＝(xi3-xi2)zi1+(xi1-xi3)zi2+(xi2-xi1)zi3 (22)
Ci＝(yi2-yi3)xi1+(yi3-yi1)xi2+(yi1-yi2)xi3 (23)
式中，Tijk表示三角形网络中以Vi、Vj、Vk为顶点的三角形单元；(xi1,yi1,zi1)、(xi1,yi1,zi1)、(xi1,yi1,zi1)为第i个三角单元的顶点坐标及函数值；M为以Vi为顶点的三角形个数；则有
Figure FDA0003625068430000061
Figure FDA0003625068430000062
由式(16)推导出Sx(Vi)、Sy(Vi)的表达式，用于直接迭代求解：
Figure FDA0003625068430000063
其中
Figure FDA0003625068430000064
Figure FDA0003625068430000065
Figure FDA0003625068430000066
将初始值
Figure FDA0003625068430000067
代入式(26)右边，求解得到
Figure FDA0003625068430000068
继续代入式(26)右边，求解得到
Figure FDA0003625068430000069
如此往复，直至计算结果满足精度要求；若令
Figure FDA00036250684300000610
其中，
Figure FDA00036250684300000611
为第k次迭代节点Vi处的一阶偏导数，则当||dk-dk-1||≤ε时停止迭代，ε为指定的迭代停止阈值，此时即认为求解精度满足要求，即
Figure FDA00036250684300000612
为满足精度的解。
9.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S8中，坐标变换的方法，步骤为：
△V1V2V3中，任意一点P(x,y)的位置由三角形的3个顶点确定的重心坐标系表示，即用下面三个比值来规定：
Figure FDA0003625068430000071
其中，b1、b2、b3为点P在重心坐标系中的坐标；S为△V1V2V3的面积；S1、S2、S3分别为△V1V2V4、△V4V2V3、△V1V4V3的面积；获知
b1+b2+b3＝1 (33)
Figure FDA0003625068430000072
Figure FDA0003625068430000073
根据(32)式～式(35)，得到重心坐标关系式
Figure FDA0003625068430000074
10.根据权利要求1所述的基于极小模网络的烟幕红外遮蔽面积测量方法，其特征是：其步骤S9中，在每个三角形单元中，插值函数为
Figure FDA0003625068430000075
Figure FDA0003625068430000076
Figure FDA0003625068430000077
其中，(i,j,k)＝{(1,2,3),(2,3,1),(3,1,2)}；
遍历所有的待插值点，逐点插值，即得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积；最后根据几何投影关系，计算出实际烟幕遮蔽面积；由烟幕遮蔽面积在红外热像仪上的投影面积，计算实际烟幕遮蔽面积的方法为
S＝NSp (42)
其中，S为实际烟幕遮蔽面积；N为图像中遮蔽率大于遮蔽率阈值的像素数；Sp为根据几何投影关系计算出的单个像素代表的实际面积。
CN202210475105.7A 2022-04-29 2022-04-29 一种基于极小模网络的烟幕红外遮蔽面积测量方法 Active CN114913218B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210475105.7A CN114913218B (zh) 2022-04-29 2022-04-29 一种基于极小模网络的烟幕红外遮蔽面积测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210475105.7A CN114913218B (zh) 2022-04-29 2022-04-29 一种基于极小模网络的烟幕红外遮蔽面积测量方法

Publications (2)

Publication Number Publication Date
CN114913218A true CN114913218A (zh) 2022-08-16
CN114913218B CN114913218B (zh) 2024-04-30

Family

ID=82765069

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210475105.7A Active CN114913218B (zh) 2022-04-29 2022-04-29 一种基于极小模网络的烟幕红外遮蔽面积测量方法

Country Status (1)

Country Link
CN (1) CN114913218B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114926408A (zh) * 2022-04-29 2022-08-19 中国人民解放军63891部队 烟幕对红外遮蔽效果的测量方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012106435A (ja) * 2010-11-18 2012-06-07 National Printing Bureau 透過潜像画像を有する印刷媒体
US20190005171A1 (en) * 2017-04-26 2019-01-03 Central South University Method for constructing finite element interpolation function
CN111242912A (zh) * 2020-01-08 2020-06-05 北京电子工程总体研究所 一种烟幕干扰效能获取方法
CN113554694A (zh) * 2021-09-22 2021-10-26 中国人民解放军国防科技大学 一种烟幕释放过程中红外有效遮蔽面积的获取方法和系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012106435A (ja) * 2010-11-18 2012-06-07 National Printing Bureau 透過潜像画像を有する印刷媒体
US20190005171A1 (en) * 2017-04-26 2019-01-03 Central South University Method for constructing finite element interpolation function
CN111242912A (zh) * 2020-01-08 2020-06-05 北京电子工程总体研究所 一种烟幕干扰效能获取方法
CN113554694A (zh) * 2021-09-22 2021-10-26 中国人民解放军国防科技大学 一种烟幕释放过程中红外有效遮蔽面积的获取方法和系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
孙杜娟;胡以华;王磊;徐凯;: "一种新的红外烟幕遮蔽效果时域特性评估方法", 航天电子对抗, no. 01, 28 February 2010 (2010-02-28) *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114926408A (zh) * 2022-04-29 2022-08-19 中国人民解放军63891部队 烟幕对红外遮蔽效果的测量方法
CN114926408B (zh) * 2022-04-29 2024-05-28 中国人民解放军63891部队 烟幕对红外遮蔽效果的测量方法

Also Published As

Publication number Publication date
CN114913218B (zh) 2024-04-30

Similar Documents

Publication Publication Date Title
WO2021098083A1 (zh) 基于显著特征的多光谱相机动态立体标定算法
CN110941999B (zh) 一种人群计数系统中自适应计算高斯核大小的方法
CN105976330A (zh) 一种嵌入式雾天实时视频稳像方法
CN107462897B (zh) 基于激光雷达的三维建图的方法
CN104182968B (zh) 宽基线多阵列光学探测系统模糊动目标分割方法
CN110796691B (zh) 一种基于形状上下文和hog特征的异源图像配准方法
CN112652020B (zh) 一种基于AdaLAM算法的视觉SLAM方法
CN110910456B (zh) 基于Harris角点互信息匹配的立体相机动态标定方法
CN108022261B (zh) 一种基于改进光流场模型的非刚性图像配准方法
CN111914695B (zh) 一种基于机器视觉的涌潮监测方法
CN111210396A (zh) 一种结合天空光偏振模型的多光谱偏振图像去雾方法
Liu et al. Generic distortion model for metrology under optical microscopes
CN115170619B (zh) 一种基于密集光流法的云遮挡预测方法
CN114913218B (zh) 一种基于极小模网络的烟幕红外遮蔽面积测量方法
CN110580715B (zh) 一种基于照度约束和格网变形的图像对齐方法
CN110390338B (zh) 一种基于非线性引导滤波与比率梯度的sar高精度匹配方法
CN114926524B (zh) 一种用于提高烟幕红外有效遮蔽面积测量精度的方法
Stent et al. Precise deterministic change detection for smooth surfaces
Ljubičić et al. SSIMS-Flow: Image velocimetry workbench for open-channel flow rate estimation
CN112837353A (zh) 一种基于多阶特征点线匹配的异源影像匹配方法
CN110910457B (zh) 基于角点特征的多光谱立体相机外参计算方法
CN116863357A (zh) 一种无人机遥感堤坝影像定标与智能分割变化检测方法
CN116958416A (zh) 一种三维建模方法、装置、系统以及存储介质
CN116402904A (zh) 一种基于激光雷达间和单目相机的联合标定方法
CN108010076B (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