CN114926524A - 一种用于提高烟幕红外有效遮蔽面积测量精度的方法 - Google Patents

一种用于提高烟幕红外有效遮蔽面积测量精度的方法 Download PDF

Info

Publication number
CN114926524A
CN114926524A CN202210473572.6A CN202210473572A CN114926524A CN 114926524 A CN114926524 A CN 114926524A CN 202210473572 A CN202210473572 A CN 202210473572A CN 114926524 A CN114926524 A CN 114926524A
Authority
CN
China
Prior art keywords
smoke screen
area
infrared
point source
coordinates
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
CN202210473572.6A
Other languages
English (en)
Other versions
CN114926524B (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 CN202210473572.6A priority Critical patent/CN114926524B/zh
Publication of CN114926524A publication Critical patent/CN114926524A/zh
Application granted granted Critical
Publication of CN114926524B publication Critical patent/CN114926524B/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
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G06T17/205—Re-meshing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/20—Image enhancement or restoration using local operators
    • G06T5/30—Erosion or dilatation, e.g. thinning
    • 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/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/155—Segmentation; Edge detection involving morphological operators
    • 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

Landscapes

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

Abstract

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

Description

一种用于提高烟幕红外有效遮蔽面积测量精度的方法
技术领域
本发明属于烟幕干扰评估技术领域，尤其是涉及一种用于提高烟幕红外有效遮蔽面积测量精度的方法。
背景技术
在军民融合、大数据分析等领域，随着现代科学技术的发展，新型可见光与红外波段侦察设备的广泛应用，使侦察与监视的水平和能力有了极大提高，对具有高效遮蔽能力的烟幕需求越来越迫切。烟幕遮蔽是一种快速、高效遮蔽目标的手段。因此，烟幕有效遮蔽性能的研究对于烟幕红外有效遮蔽面积评估具有重大意义。
对红外烟幕遮敝效果的测量主要是通过对烟幕透过率、遮蔽面积、形成时间和持续时间等的测量来定量描述的。目前，使用最多的是用凝视型红外成像热像仪的多点测量法。
现有技术中，在使用摄像法测量烟幕的红外遮蔽性能时，通常的方法是采用布设红外点源阵列，在红外点源阵列前方释放烟幕，并使用红外热像仪分别测量烟幕施放前、施放时红外点源阵列的红外辐射强度变化，得到点源位置处的遮蔽率(透过率)，再通过插值的方法得出烟幕遮蔽率的空间分布，最后利用烟幕遮蔽率阈值确定有效遮蔽面积。在布设红外点源阵列时，应使其在红外热像仪上的成像间距分布均匀，以便于插值。但在实际试验中，由于试验条件限制，很难使红外点源阵列的成像间距均匀分布，导致插值算法复杂、插值精度差，得出的有效遮蔽面积结果不精确，数据处理时间长，使得试验周期长、成本高、效率低。
发明内容
针对现有技术中在使用摄像法测量烟幕的红外遮蔽性能时，实际试验中存在红外点源阵列在红外热像仪上的成像间距分布不均匀的问题，本发明的目的是提供一种用于提高烟幕红外有效遮蔽面积测量精度的方法，其利用重心坐标变换与形状无关的特性，将原本散乱的插值节点转换为易插值的重心坐标节点，从而降低插值的算法复杂度，提高精度和计算速度。
为实现上述发明目的，本发明采用如下技术方案：
一种用于提高烟幕红外有效遮蔽面积测量精度的方法，其包括以下步骤：
S1、图像预处理：采用数学形态学中的顶帽变换算法对获取的烟幕下的红外原始图像进行滤波，增强红外点源灰度值，抑制背景；
S2、图像分割：使用粗分割、精分割两步分割法对步骤S1中预处理后的图像进行分割；操作为：首先，使用全局阈值分割算法提取出每个红外点源的粗略位置；然后在每个点源位置处向外扩张一定距离，划分为局部待分割区域；最后，在每个局部待分割区域内使用最大类间方差法依次对每个点源进行分割，提取点源形状和位置；
S3、提取点源位置坐标：根据步骤S2中提取出的点源形状和位置，计算所有点源的形心或重心位置坐标，作为每个点源的位置坐标；
S4、计算点源的烟幕遮蔽率或透过率：根据步骤S2中的图像分割结果中的目标区域和背景区域施放烟幕前、施放烟幕后的灰度变化，计算每个红外点源位置处的烟幕遮蔽率或透过率；
S5、三角剖分：根据步骤S2中提取出的点源位置坐标，将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分为基本三角形单元，三角形单元的每个顶点依次对应一个点源，所有点源形成的三角形网络构成插值曲面片；
S6、坐标变换：将待插值点直角坐标(xi,yi)进行坐标变换，变换为局部重心坐标(b1,b2,b3,b4)；
S7、计算插值系数：根据步骤S4中计算出的红外点源位置处烟幕遮蔽率或透过率，估计每个红外点源位置处烟幕遮蔽率或透过率的偏导数值，并计算插值系数；最后，利用插值公式插值计算坐标(b1,b2,b3,b4)处的烟幕遮蔽率，坐标(b1,b2,b3,b4)处的烟幕遮蔽率即为坐标(xi,yi)处的烟幕遮蔽率；
S8、计算遮蔽面积：遍历所有的待插值点，对所有待插值点执行步骤S6、S7，得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积，最后，根据几何投影关系，计算出实际烟幕遮蔽面积。
进一步地，上述的步骤S1中，使用顶帽变换算法对图像预处理，顶帽变换算法的计算公式如下：
使用结构元素b对图像f(x,y)的灰度膨胀记为
Figure BDA0003624153580000031
定义为：
Figure BDA0003624153580000032
式中：Db是结构元素b的定义域；
使用结构元素b对图像f(x,y)的灰度腐蚀记为
Figure BDA0003624153580000033
定义为：
Figure BDA0003624153580000034
灰度腐蚀是一个局部最小值算子；
使用结构元素b对图像f(x,y)的开运算记为
Figure BDA0003624153580000035
定义为：
Figure BDA0003624153580000036
即f先由b腐蚀，再由b膨胀；
使用结构元素b对图像f(x,y)的闭运算记为f·b，定义为：
Figure BDA0003624153580000037
即f先由b膨胀，再由b腐蚀；
灰度图像的顶帽变换定义为f减去其开运算结果，即：
Figure BDA0003624153580000041
进一步地，上述的步骤S2中，使用最大类间方差法对局部待分割区域进行分割，最大化类间方差
Figure BDA0003624153580000042
计算公式如下：
Figure BDA0003624153580000043
其中，
Figure BDA0003624153580000044
Figure BDA0003624153580000045
Figure BDA0003624153580000046
Figure BDA0003624153580000047
Figure BDA0003624153580000048
Figure BDA0003624153580000049
式中：nq是待分割区域Z中第q级灰度的像素数；N是待分割区域Z中的像素总数；L是待分割区域Z中的灰度级；
求出待分割区域所有可能灰度级的类间方差后，用下式求得分割阈值t：
Figure BDA00036241535800000410
分割出的目标区域表示为待分割区域Z中灰度值f(x,y)≥t的区域，即
T＝{(x,y)|f(x,y)≥t,(x,y)∈Z}。
进一步地，上述的步骤S3中，计算点源形心坐标的算法为：
Figure BDA00036241535800000411
计算点源重心坐标的算法为：
Figure BDA0003624153580000051
其中，n为目标区域T的像素总数。
进一步地，上述的步骤S4中，点源的烟幕透过率计算公式如下：
Figure BDA0003624153580000052
点源的烟幕遮蔽率计算公式如下：
Figure BDA0003624153580000053
其中，
Figure BDA0003624153580000054
分别为施放烟幕前、施放烟幕后目标区域的平均灰度值，
Figure BDA0003624153580000055
Figure BDA0003624153580000056
分别为施放烟幕前、施放烟幕后背景区域的平均灰度值。
进一步地，上述的步骤S5中，划分基本三角形单元的方法，步骤为：
首先根据空间位置关系将所有点源按每列、每行点源连接划分为基本四边形单元；已知任意1个四边形能够沿对角顶点划分为2个三角形且有两种划分方法，分别按两种方法进行划分并计算2个三角形的最小内角，取最小内角最大的划分方法进行划分。
进一步地，上述的步骤S6中，坐标变换的方法，步骤为：
△V1V2V3中，任意一点P(x,y)的位置由三角形的3个顶点确定的重心坐标系表示，即用下面三个比值来规定：
Figure BDA0003624153580000057
其中，b1、b2、b3为点P在重心坐标系中的坐标；S为△V1V2V3的面积；S1、S2、S3分别为△V1V2V4、△V4V2V3、△V1V4V3的面积；获知
b1+b2+b3＝1 (13)
Figure BDA0003624153580000061
Figure BDA0003624153580000062
根据式(12)～式(15)，得到重心坐标关系式
Figure BDA0003624153580000063
也写为
Figure BDA0003624153580000064
其逆变换为
Figure BDA0003624153580000065
若在每个三角单元中取一点V4(x4,y4)，其中，x4＝(x1+x2+x3)/3，y4＝(y1+y2+y3)/3，进一步将每个三角单元划分为3个子三角单元，则三坐标(b1,b2,b3)能够由四坐标(b1,b2,b3,b4)表示，即在每个子三角单元Ti中，有bi＝0，其余3个坐标由坐标变换公式(17)求得，有
Figure BDA0003624153580000066
其中，b′i为在3个子三角单元中用式(17)求得的坐标，bi为四坐标表示重心坐标，Ti为顶点Vi对应的子三角单元编号。
进一步地，上述的步骤S7中，偏导数值Fx(xi,yi)、Fy(xi,yi)的估计方法为：假设F(x,y)为二次函数，即
F(x,y)＝ai1(x-xi)2+ai2(x-xi)(y-yi)+ai3(y-yi)2+ai4(x-xi)+ai5(y-yi)+ai6 (20)
则
Fx(xi,yi)＝ai4 (21)
Fy(xi,yi)＝ai5 (22)
系数a,b,c,d,e,f由加权最小二乘法求得，即求解使下列方程最小化问题；
Figure BDA0003624153580000071
其中，Di称为节点(xi,yi)的影响域；若令ri(x,y)表示任一点(x,y)到节点(xi,yi)的距离，则有Di＝{(x,y)|ri(x,y)≤Ri}，Ri称为节点(xi,yi)的影响半径；权值wij由下式定义
Figure BDA0003624153580000072
令
Figure BDA0003624153580000073
Figure BDA0003624153580000074
Figure BDA0003624153580000075
则有
Figure BDA0003624153580000076
遍历所有的三角单元节点，即能够由式(20)～式(28)估算出所有节点的偏导数值；
插值公式为
Figure BDA0003624153580000081
其中，n＝3，i,j,k,l＝0,1,...,n，λ:(i,j,k,l)，且|λ|＝i+j+k+l，在各控制点处有坐标(b1,b2,b3,b4)＝(i/n,j/n,k/n,l/n)，且
Figure BDA0003624153580000082
Figure BDA0003624153580000083
令eij表示三角单元上从顶点Vi到Vj的方向矢量，若已知三角单元顶点(xi,yi)处的函数值F(xi,yi)和偏导数值Fx(xi,yi)、Fy(xi,yi)，则有方向导数
Figure BDA0003624153580000084
其中||eij||为矢量eij的欧式长度；
插值公式(29)中，系数cijkl由下列公式求得：
c3000＝F(x1,y1) (33)
c0300＝F(x2,y2) (34)
c0030＝F(x3,y3) (35)
Figure BDA0003624153580000085
Figure BDA0003624153580000086
Figure BDA0003624153580000087
Figure BDA0003624153580000088
Figure BDA0003624153580000091
Figure BDA0003624153580000092
Figure BDA0003624153580000093
Figure BDA0003624153580000094
Figure BDA0003624153580000095
Figure BDA0003624153580000096
Figure BDA0003624153580000097
Figure BDA0003624153580000098
Figure BDA0003624153580000099
Figure BDA00036241535800000910
Figure BDA00036241535800000911
Figure BDA00036241535800000912
Figure BDA00036241535800000913
Figure BDA00036241535800000914
Figure BDA00036241535800000915
其中，F(xi,yi)即为三角单元顶点Vi处的烟幕遮蔽率或透过率。
进一步地，上述的步骤S8中，由烟幕遮蔽面积在红外热像仪上的投影面积计算实际烟幕遮蔽面积的方法为
S＝NSp (55)
其中，S为实际烟幕遮蔽面积；N为图像中遮蔽率大于遮蔽率阈值的像素数；Sp为根据几何投影关系计算出的单个像素代表的实际面积。
由于采用如上所述的技术方案，本发明具有如下优越性：
该用于提高烟幕红外有效遮蔽面积测量精度的方法，其利用坐标变换的方法，将不均匀的成像阵列网格转换为结构化的均匀网格后再进行插值，能够有效避免红外点源阵列分布不均匀造成数据点分布不均匀，进而导致插值算法复杂，插值精度差、计算速度慢等现象；能够极大的减少插值算法的复杂度，提高插值精度，减少数据处理时间，缩短试验周期，降低成本。
附图说明
图1是红外源阵列沿地面水平布设时，摄像法的典型试验布局示意图；
图2是理想情况下，红外点源阵列在红外热像仪上的投影分布示意图；
图3是本发明一实施例提供的试验中，红外点源阵列在红外热像仪上的实际投影分布图；
图4对图3进行三角剖分结果图；
图5是重心坐标变换示意图；
图6是三角单元继续划分为3个子单元示意图及其各控制点位置及编号；
图7是本发明一实施例提供的一帧图像的最终插值结果及等位线图。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步详细说明。
如图1～7所示，一种用于提高烟幕红外有效遮蔽面积测量精度的方法，其包括以下步骤：
S1、图像预处理：由于弥散作用，红外点源在红外热像仪上所成的像不是一个点，而是一个近似于圆形的小目标，因此需使用图像分割算法提取目标，分割之前需对图像进行预处理，以避免误分、漏分；
已知，形态学的基本运算有4个：膨胀(或扩张)、腐蚀(或侵蚀)、开启和闭合，基于上述基本运算能够推导和组合成各种数学形态学实用算法，用它们进行图像形状和结构的分析及处理，包括图像分割、特征抽取、边缘检测、图像滤波、图像增强和恢复等；
本发明采用数学形态学中的顶帽变换算法对获取的烟幕下的红外原始图像进行滤波，增强红外点源灰度值，抑制背景，以避免分割时出现误分、漏分；顶帽变换算法的计算公式如下：
使用结构元素b对图像f(x,y)的灰度膨胀记为
Figure BDA0003624153580000111
定义为：
Figure BDA0003624153580000112
式中，Db是结构元素b的定义域；
使用结构元素b对图像f(x,y)的灰度腐蚀记为
Figure BDA0003624153580000113
定义为：
Figure BDA0003624153580000114
灰度腐蚀是一个局部最小值算子；
使用结构元素b对图像f(x,y)的开运算记为
Figure BDA0003624153580000115
定义为：
Figure BDA0003624153580000116
即，f先由b腐蚀，再由b膨胀；
使用结构元素b对图像f(x,y)的闭运算记为f·b，定义为：
Figure BDA0003624153580000117
即，f先由b膨胀，再由b腐蚀；
灰度图像的顶帽变换定义为f减去其开运算结果，即：
Figure BDA0003624153580000118
顶帽变换算法用于暗背景上的亮目标，用来校正不均匀背景的影响；
上述式中，结构元素b的定义类似于邻域，区别的是具有各种不同的形状，如菱形、圆盘、线形、八边形、甚至2个以上的结构元素等，二值矩阵Db决定了哪些位置的元素包括在最大值运算中，可以将膨胀运算看作是一个局部最大值算子；本实施例中，结构元素b取半径为3的圆；
S2、图像分割：使用粗分割、精分割两步分割法对步骤S1中预处理后的图像进行分割；
步骤为：首先，使用全局阈值分割算法提取出每个红外点源的粗略位置；然后在每个点源位置处向外扩张一定距离d，划分为局部待分割区域，本实施例中，取距离d＝5；最后，在每个局部待分割区域内使用最大类间方差法依次对每个点源进行分割，提取点源形状和位置；
最大化类间方差
Figure BDA0003624153580000121
计算公式如下：
Figure BDA0003624153580000122
其中，
Figure BDA0003624153580000123
Figure BDA0003624153580000124
Figure BDA0003624153580000125
Figure BDA0003624153580000126
Figure BDA0003624153580000127
Figure BDA0003624153580000128
式中，nq是待分割区域Z中第q级灰度的像素数；N是待分割区域中的像素总数；L是待分割区域Z中的灰度级；
求出待分割区域所有可能灰度级的类间方差后，用下式求得分割阈值t：
Figure BDA0003624153580000129
分割出的目标区域表示为待分割区域Z中灰度值f(x,y)≥t的区域，即
T＝{(x,y)|f(x,y)≥t,(x,y)∈Z}；
S3、提取点源位置坐标：根据步骤S2中提取出的点源形状和位置，计算所有点源的形心或重心位置坐标，作为每个点源的位置坐标；
计算点源形心坐标的算法为：
Figure BDA0003624153580000131
计算点源重心坐标的算法为：
Figure BDA0003624153580000132
式中，n为目标区域T的像素总数；
本实施例中，选取点源重心坐标作为点源位置坐标；
S4、计算点源的烟幕遮蔽率或透过率，根据试验要求确定：步骤S2中，每个局部待分割区域都能够根据最终分割结果划分为目标区域和背景区域两部分；施放烟幕前后，目标区域和背景区域的灰度变化不同，进而能够计算出每个目标位置处的烟幕遮蔽率或透过率；
点源的烟幕透过率计算公式如下：
Figure BDA0003624153580000133
点源的烟幕遮蔽率计算公式如下：
Figure BDA0003624153580000134
其中，
Figure BDA0003624153580000135
分别为施放烟幕前、施放烟幕后目标区域的平均灰度值，
Figure BDA0003624153580000136
分别为施放烟幕前、施放烟幕后背景区域的平均灰度值；
本实施例中，选择计算点源的烟幕遮蔽率；
S5、三角剖分：根据步骤S2中提取出的点源位置坐标，将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分为基本三角形单元，三角形单元的每个顶点依次对应一个点源，所有点源形成的三角形网络构成插值曲面片；
划分基本三角形单元的方法为：首先，根据空间位置关系将所有点源按每列、每行点源连接划分为基本四边形单元；已知任意1个四边形可沿对角顶点划分为2个三角形且有两种划分方法，根据三角剖分的相关理论，三角形单元越接近于等边三角形，三角形单元的质量越好，引起的误差越小，利用这一原则划分基本三角形单元；在划分完成基本四边形单元的基础上，已知任意1个四边形可沿对角顶点划分为2个三角形且有两种划分方法，在划分完成基本四边形单元的基础上，分别按两种方法进行划分并计算2个三角形的最小内角，取最小内角最大的划分方法进行划分(或取最长边与最短边之比最接近于1的划分方法进行划分，二者等价)；
S6、坐标变换：将待插值点直角坐标(xi,yi)进行坐标变换，变换为局部重心坐标(b1,b2,b3,b4)；坐标变换的方法为：
△V1V2V3中，任意一点P(x,y)的位置由三角形的3个顶点确定的重心坐标系表示，即用下面三个比值来规定：
Figure BDA0003624153580000141
其中，b1、b2、b3称为点P在重心坐标系中的坐标；S为△V1V2V3的面积；S1、S2、S3分别为△V1V2V4、△V4V2V3、△V1V4V3的面积；获知
b1+b2+b3＝1 (13)
Figure BDA0003624153580000142
Figure BDA0003624153580000143
根据(12)式～式(15)，得到重心坐标关系式
Figure BDA0003624153580000151
也写为
Figure BDA0003624153580000152
其逆变换为
Figure BDA0003624153580000153
若在每个三角单元中取一点V4(x4,y4)，其中，x4＝(x1+x2+x3)/3，y4＝(y1+y2+y3)/3，进一步将每个三角单元划分为3个子三角单元，则三坐标(b1,b2,b3)能够由四坐标(b1,b2,b3,b4)表示，即在每个子三角单元Ti中，有bi＝0，其余3个坐标由坐标变换公式(17)求得，有
Figure BDA0003624153580000154
其中，b′i为在3个子三角单元中用式(17)求得的坐标，bi为四坐标表示重心坐标，Ti为顶点Vi对应的子三角单元编号；
S7、计算插值系数：根据步骤S4中计算出的红外点源位置处烟幕遮蔽率，估计每个红外点源位置处烟幕遮蔽率的偏导数值，并计算插值系数，最后利用插值公式插值计算坐标(b1,b2,b3,b4)处的烟幕遮蔽率，坐标(b1,b2,b3,b4)处的烟幕遮蔽率即为坐标(xi,yi)处的烟幕遮蔽率；
偏导数值Fx(xi,yi)、Fy(xi,yi)的估计方法为：假设F(x,y)为二次函数，即
F(x,y)＝ai1(x-xi)2+ai2(x-xi)(y-yi)+ai3(y-yi)2+ai4(x-xi)+ai5(y-yi)+ai6 (20)
那么
Fx(xi,yi)＝ai4 (21)
Fy(xi,yi)＝ai5 (22)
系数ai1,ai2,ai3,ai4,ai5,ai6可由加权最小二乘法求得，即求解使下列方程最小化问题。
Figure BDA0003624153580000161
其中，Di称为节点(xi,yi)的影响域，若令ri(x,y)表示任一点(x,y)到节点(xi,yi)的距离，则有Di＝{(x,y)|ri(x,y)≤Ri}，Ri称为节点(xi,yi)的影响半径，优选地，取距节点(xi,yi)最近的k个节点中距离的最大值，在本实施例中，取距当前节点(xi,yi)最近的8个节点中距离的最大值；权值wij由下式定义
Figure BDA0003624153580000162
令
Figure BDA0003624153580000163
Figure BDA0003624153580000164
Figure BDA0003624153580000165
则有
Figure BDA0003624153580000166
遍历所有的节点，即能够由式(20)～式(28)估算出所有节点的偏导数值；
插值方法为用下列插值公式进行插值，
Figure BDA0003624153580000171
其中，n＝3，i,j,k,l＝0,1,...,n，λ:(i,j,k,l)，且|λ|＝i+j+k+l，在各控制点处有坐标(b1,b2,b3,b4)＝(i/n,j/n,k/n,l/n)，且
Figure BDA0003624153580000172
Figure BDA0003624153580000173
令eij表示三角单元上从顶点Vi到Vj的方向矢量，若已知三角单元顶点(xi,yi)处的函数值F(xi,yi)和偏导数值Fx(xi,yi)、Fy(xi,yi)，则有方向导数
Figure BDA0003624153580000174
其中||eij||为矢量eij的欧式长度；
插值公式(29)中，系数cijkl由下列公式求得：
c3000＝F(x1,y1) (33)
c0300＝F(x2,y2) (34)
c0030＝F(x3,y3) (35)
Figure BDA0003624153580000175
Figure BDA0003624153580000176
Figure BDA0003624153580000177
Figure BDA0003624153580000178
Figure BDA0003624153580000179
Figure BDA0003624153580000181
Figure BDA0003624153580000182
Figure BDA0003624153580000183
Figure BDA0003624153580000184
Figure BDA0003624153580000185
Figure BDA0003624153580000186
Figure BDA0003624153580000187
Figure BDA0003624153580000188
Figure BDA0003624153580000189
Figure BDA00036241535800001810
Figure BDA00036241535800001811
Figure BDA00036241535800001812
Figure BDA00036241535800001813
Figure BDA00036241535800001814
其中，F(xi,yi)即为三角单元顶点Vi处的烟幕遮蔽率或透过率。
S8、计算遮蔽面积：遍历所有的待插值点，对所有待插值点执行步骤S6、S7，得到整个待计算区域的烟幕遮蔽率，根据烟幕遮蔽率阈值，可在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积，最后根据几何投影关系即可计算出实际烟幕遮蔽面积；
由烟幕遮蔽面积在红外热像仪上的投影面积计算实际烟幕遮蔽面积的方法为
S＝NSp (55)
其中，S为实际烟幕遮蔽面积，N为图像中遮蔽率大于遮蔽率阈值的像素数，Sp为根据几何投影关系计算出的单个像素代表的实际面积。
以上所述仅为本发明的较佳实施例，而非对本发明的限制，在不脱离本发明的精神和范围的情况下，凡依本发明申请专利范围所作的均等变化与修饰，皆应属本发明的专利保护范围之内。

Claims (9)

1.一种用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其其包括以下步骤：
S1、图像预处理：采用数学形态学中的顶帽变换算法对获取的烟幕下的红外原始图像进行滤波，增强红外点源灰度值，抑制背景；
S2、图像分割：使用粗分割、精分割两步分割法对步骤S1中预处理后的图像进行分割；操作为：首先，使用全局阈值分割算法提取出每个红外点源的粗略位置；然后在每个点源位置处向外扩张一定距离，划分为局部待分割区域；最后，在每个局部待分割区域内使用最大类间方差法依次对每个点源进行分割，提取点源形状和位置；
S3、提取点源位置坐标：根据步骤S2中提取出的点源形状和位置，计算所有点源的形心或重心位置坐标，作为每个点源的位置坐标；
S4、计算点源的烟幕遮蔽率或透过率：根据步骤S2中的图像分割结果中的目标区域和背景区域施放烟幕前、施放烟幕后的灰度变化，计算每个红外点源位置处的烟幕遮蔽率或透过率；
S5、三角剖分：根据步骤S2中提取出的点源位置坐标，将红外点源位置坐标作为三角单元的顶点按最小内角最大化准则划分为基本三角形单元，三角形单元的每个顶点依次对应一个点源，所有点源形成的三角形网络构成插值曲面片；
S6、坐标变换：将待插值点直角坐标(xi,yi)进行坐标变换，变换为局部重心坐标(b1,b2,b3,b4)；
S7、计算插值系数：根据步骤S4中计算出的红外点源位置处烟幕遮蔽率或透过率，估计每个红外点源位置处烟幕遮蔽率或透过率的偏导数值，并计算插值系数；最后，利用插值公式插值计算坐标(b1,b2,b3,b4)处的烟幕遮蔽率，坐标(b1,b2,b3,b4)处的烟幕遮蔽率即为坐标(xi,yi)处的烟幕遮蔽率；
S8、计算遮蔽面积：遍历所有的待插值点，对所有待插值点执行步骤S6、S7，得到整个待计算区域的烟幕遮蔽率；根据烟幕遮蔽率阈值，在图像中分割出烟幕遮蔽率大于遮蔽率阈值的区域，统计该区域的像素数，即为烟幕遮蔽面积在红外热像仪上的投影面积，最后，根据几何投影关系，计算出实际烟幕遮蔽面积。
2.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S1中，使用顶帽变换算法对图像预处理，顶帽变换算法的计算公式如下：
使用结构元素b对图像f(x,y)的灰度膨胀记为
Figure FDA0003624153570000021
定义为：
Figure FDA0003624153570000022
式中：Db是结构元素b的定义域；
使用结构元素b对图像f(x,y)的灰度腐蚀记为
Figure FDA0003624153570000026
定义为：
Figure FDA0003624153570000023
灰度腐蚀是一个局部最小值算子；
使用结构元素b对图像f(x,y)的开运算记为
Figure FDA0003624153570000027
定义为：
Figure FDA0003624153570000024
即f先由b腐蚀，再由b膨胀；
使用结构元素b对图像f(x,y)的闭运算记为f·b，定义为：
Figure FDA0003624153570000025
即f先由b膨胀，再由b腐蚀；
灰度图像的顶帽变换定义为f减去其开运算结果，即：
Figure FDA0003624153570000028
3.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S2中，使用最大类间方差法对局部待分割区域进行分割，最大化类间方差
Figure FDA0003624153570000031
计算公式如下：
Figure FDA0003624153570000032
其中，
Figure FDA0003624153570000033
Figure FDA0003624153570000034
Figure FDA0003624153570000035
Figure FDA0003624153570000036
Figure FDA0003624153570000037
Figure FDA0003624153570000038
式中：nq是待分割区域Z中第q级灰度的像素数；N是待分割区域Z中的像素总数；L是待分割区域Z中的灰度级；
求出待分割区域所有可能灰度级的类间方差后，用下式求得分割阈值t：
Figure FDA0003624153570000039
分割出的目标区域表示为待分割区域Z中灰度值f(x,y)≥t的区域，即
T＝{(x,y)|f(x,y)≥t,(x,y)∈Z}。
4.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S3中，计算点源形心坐标的算法为：
Figure FDA0003624153570000041
计算点源重心坐标的算法为：
Figure FDA0003624153570000042
其中，n为目标区域T的像素总数。
5.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S4中，点源的烟幕透过率计算公式如下：
Figure FDA0003624153570000043
点源的烟幕遮蔽率计算公式如下：
Figure FDA0003624153570000044
其中，
Figure FDA0003624153570000045
分别为施放烟幕前、施放烟幕后目标区域的平均灰度值，
Figure FDA0003624153570000046
Figure FDA0003624153570000047
分别为施放烟幕前、施放烟幕后背景区域的平均灰度值。
6.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S5中，划分基本三角形单元的方法，步骤为：
首先根据空间位置关系将所有点源按每列、每行点源连接划分为基本四边形单元；已知任意1个四边形能够沿对角顶点划分为2个三角形且有两种划分方法，分别按两种方法进行划分并计算2个三角形的最小内角，取最小内角最大的划分方法进行划分。
7.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S6中，坐标变换的方法，步骤为：
△V1V2V3中，任意一点P(x,y)的位置由三角形的3个顶点确定的重心坐标系表示，即用下面三个比值来规定：
Figure FDA0003624153570000051
其中，b1、b2、b3为点P在重心坐标系中的坐标；S为△V1V2V3的面积；S1、S2、S3分别为△V1V2V4、△V4V2V3、△V1V4V3的面积；获知
b1+b2+b3＝1 (13)
Figure FDA0003624153570000052
Figure FDA0003624153570000053
根据式(12)～式(15)，得到重心坐标关系式
Figure FDA0003624153570000054
也写为
Figure FDA0003624153570000055
其逆变换为
Figure FDA0003624153570000056
若在每个三角单元中取一点V4(x4,y4)，其中，x4＝(x1+x2+x3)/3，y4＝(y1+y2+y3)/3，进一步将每个三角单元划分为3个子三角单元，则三坐标(b1,b2,b3)能够由四坐标(b1,b2,b3,b4)表示，即在每个子三角单元Ti中，有bi＝0，其余3个坐标由坐标变换公式(17)求得，有
Figure FDA0003624153570000061
其中，b′i为在3个子三角单元中用式(17)求得的坐标，bi为四坐标表示重心坐标，Ti为顶点Vi对应的子三角单元编号。
8.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S7中，偏导数值Fx(xi,yi)、Fy(xi,yi)的估计方法为：假设F(x,y)为二次函数，即
F(x,y)＝ai1(x-xi)2+ai2(x-xi)(y-yi)+ai3(y-yi)2+ai4(x-xi)+ai5(y-yi)+ai6 (20)
则
Fx(xi,yi)＝ai4 (21)
Fy(xi,yi)＝ai5 (22)
系数a,b,c,d,e,f由加权最小二乘法求得，即求解使下列方程最小化问题；
Figure FDA0003624153570000062
其中，Di称为节点(xi,yi)的影响域；若令ri(x,y)表示任一点(x,y)到节点(xi,yi)的距离，则有Di＝{(x,y)|ri(x,y)≤Ri}，Ri称为节点(xi,yi)的影响半径；权值wij由下式定义
Figure FDA0003624153570000063
令
Figure FDA0003624153570000064
Figure FDA0003624153570000071
Figure FDA0003624153570000072
则有
Figure FDA0003624153570000073
遍历所有的三角单元节点，即能够由式(20)～式(28)估算出所有节点的偏导数值；
插值公式为
Figure FDA0003624153570000074
其中，n＝3，i,j,k,l＝0,1,...,n，λ:(i,j,k,l)，且|λ|＝i+j+k+l，在各控制点处有坐标(b1,b2,b3,b4)＝(i/n,j/n,k/n,l/n)，且
Figure FDA0003624153570000075
Figure FDA0003624153570000076
令eij表示三角单元上从顶点Vi到Vj的方向矢量，若已知三角单元顶点(xi,yi)处的函数值F(xi,yi)和偏导数值Fx(xi,yi)、Fy(xi,yi)，则有方向导数
Figure FDA0003624153570000077
其中||eij||为矢量eij的欧式长度；
插值公式(29)中，系数cijkl由下列公式求得：
c3000＝F(x1,y1) (33)
c0300＝F(x2,y2) (34)
c0030＝F(x3,y3) (35)
Figure FDA0003624153570000081
Figure FDA0003624153570000082
Figure FDA0003624153570000083
Figure FDA0003624153570000084
Figure FDA0003624153570000085
Figure FDA0003624153570000086
Figure FDA0003624153570000087
Figure FDA0003624153570000088
Figure FDA0003624153570000089
Figure FDA00036241535700000810
Figure FDA00036241535700000811
Figure FDA00036241535700000812
Figure FDA00036241535700000813
Figure FDA00036241535700000814
Figure FDA00036241535700000815
Figure FDA00036241535700000816
Figure FDA00036241535700000817
Figure FDA0003624153570000091
Figure FDA0003624153570000092
其中，F(xi,yi)即为三角单元顶点Vi处的烟幕遮蔽率或透过率。
9.根据权利要求1所述的用于提高烟幕红外有效遮蔽面积测量精度的方法，其特征是：其步骤S8中，由烟幕遮蔽面积在红外热像仪上的投影面积计算实际烟幕遮蔽面积的方法为
S＝NSp (55)
其中，S为实际烟幕遮蔽面积；N为图像中遮蔽率大于遮蔽率阈值的像素数；Sp为根据几何投影关系计算出的单个像素代表的实际面积。
CN202210473572.6A 2022-04-29 2022-04-29 一种用于提高烟幕红外有效遮蔽面积测量精度的方法 Active CN114926524B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210473572.6A CN114926524B (zh) 2022-04-29 2022-04-29 一种用于提高烟幕红外有效遮蔽面积测量精度的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210473572.6A CN114926524B (zh) 2022-04-29 2022-04-29 一种用于提高烟幕红外有效遮蔽面积测量精度的方法

Publications (2)

Publication Number Publication Date
CN114926524A true CN114926524A (zh) 2022-08-19
CN114926524B CN114926524B (zh) 2024-04-26

Family

ID=82807049

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210473572.6A Active CN114926524B (zh) 2022-04-29 2022-04-29 一种用于提高烟幕红外有效遮蔽面积测量精度的方法

Country Status (1)

Country Link
CN (1) CN114926524B (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 (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200043228A1 (en) * 2018-08-01 2020-02-06 Nvidia Corporation System-generated stable barycentric coordinates and direct plane equation access
CN111242912A (zh) * 2020-01-08 2020-06-05 北京电子工程总体研究所 一种烟幕干扰效能获取方法
CN113554694A (zh) * 2021-09-22 2021-10-26 中国人民解放军国防科技大学 一种烟幕释放过程中红外有效遮蔽面积的获取方法和系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200043228A1 (en) * 2018-08-01 2020-02-06 Nvidia Corporation System-generated stable barycentric coordinates and direct plane equation access
CN111242912A (zh) * 2020-01-08 2020-06-05 北京电子工程总体研究所 一种烟幕干扰效能获取方法
CN113554694A (zh) * 2021-09-22 2021-10-26 中国人民解放军国防科技大学 一种烟幕释放过程中红外有效遮蔽面积的获取方法和系统

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
吕俊伟;卢毅;张鹏;谢勇;: "一种用热像仪获得烟幕透过率的方法", 红外技术, no. 10, 28 October 2006 (2006-10-28) *
李华;魏文俭;康华超;张宝东;王敏;: "基于多分辨率分析的烟幕遮蔽矢量提取方法分析", 光电技术应用, no. 02, 15 April 2007 (2007-04-15) *
李志国;白云塔;张思将;: "基于热像仪的红外烟幕遮蔽率测试方法研究", 红外, no. 11, 10 November 2008 (2008-11-10) *
窦茂森;高文静;王继光;耿春萍;: "红外烟幕遮蔽效果测量与评估", 光电技术应用, no. 05, 30 October 2006 (2006-10-30) *
高志扬;祝利;韩书键;: "烟幕弹对毫米波末制导雷达的有效遮蔽区域研究", 弹箭与制导学报, no. 06, 15 December 2015 (2015-12-15) *

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
CN114926524B (zh) 2024-04-26

Similar Documents

Publication Publication Date Title
CN109146930B (zh) 一种电力机房设备红外与可见光图像配准方法
CN106548462B (zh) 基于薄板样条插值的非线性sar图像几何校正方法
US11783457B2 (en) Multispectral camera dynamic stereo calibration algorithm based on saliency features
CN106886748B (zh) 一种基于tld的适用于无人机的变尺度目标跟踪方法
CN106887016B (zh) 一种gf-4卫星序列图像自动相对配准方法
CN111914695B (zh) 一种基于机器视觉的涌潮监测方法
WO2019010932A1 (zh) 一种利于模糊核估计的图像区域选择方法和系统
CN108257153B (zh) 一种基于方向梯度统计特征的目标跟踪方法
CN114926524A (zh) 一种用于提高烟幕红外有效遮蔽面积测量精度的方法
CN111062895B (zh) 一种基于多视场分割的显微图像复原方法
Ioli et al. UAV photogrammetry for metric evaluation of concrete bridge cracks
CN106846250B (zh) 一种基于多尺度滤波的超分辨率重建方法
CN106097292B (zh) 一种用于遥感影像处理的时空邻域高斯加权中值滤波方法
CN114913218B (zh) 一种基于极小模网络的烟幕红外遮蔽面积测量方法
CN107369163B (zh) 一种基于最佳熵双阈值分割的快速sar图像目标检测方法
Ljubičić et al. SSIMS-Flow: Image velocimetry workbench for open-channel flow rate estimation
CN108830273A (zh) 基于图像局部对比度的能见度测量方法
CN112862866A (zh) 基于麻雀搜索算法的图像配准方法及系统、计算设备
CN114926408B (zh) 烟幕对红外遮蔽效果的测量方法
US20230186594A1 (en) Information processing apparatus, information processing method, and non-transitory computer readable medium
CN113030968B (zh) 基于csar模式提取dem的方法、装置及存储介质
CN115082359A (zh) 基于海岸线数据的同步轨道光学卫星几何精校正方法
CN115049549A (zh) 一种基于稳健估计的红外图像条状噪声去除方法
CN115131555A (zh) 基于超像素分割的叠掩阴影检测方法及装置
CN110223250B (zh) 基于单应变换的sar几何校正方法

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