Partial Recursive Segmentation Algorithm Based on Otsu and Image Entropy
-
摘要: 在图像分割中,门限的选取至关重要,本文在Otsu门限选择准则[1]的基础上,研究了一种基于图像熵的局部递归的门限选择及分割算法,该算法在图像背景不均或图像不是简单的单峰、双峰图像的情况下可以有效的分割,分割后图像的细节更加丰富,有利于分割后的特征提取.对几种不同的目标进行了实验,获得了较好的实验结果.
-
0. 引言
由于红外探测器探测距离限制，红外体制的导引头多用于末端防空领域或与其他模式组合成复合制导体制完成打击任务。空中红外截获是复合制导体制中的关键技术，而目标的落入概率和检测识别概率是空中截获技术中决定着打击成功率的重要技术指标[1-2]。对于低目标指示精度下的空中红外截获，框架扫描搜索技术可以有效提高目标的落入概率。而传统的信息处理技术只能对探测器视场内的目标进行探测和识别，对于处于探测视场外的目标不具备识别能力，故本文设计了一种实时扫描和视场拼接算法，首先根据红外导引头特点选择合适的框架扫描控制方法，在扫描过程中得到红外图像信息，通过将惯性信息与图像特征偏移信息相融合的图像拼接方法，构造了虚拟的大角度探测视场，为后期信息处理提供了基础。
1. 位标器两轴框架稳定扫描搜索算法
框架扫描搜索是增加目标落入概率的重要实现途径[3]，扫描方式有定点扫描、圆形扫描和匀速矩形扫描。定点扫描需要机构快速的运动到设定的位置点并在停留的位置点进行探测，这对控制系统的快速性要求极高，同时在机构快速运动的过程中探测到的图像会出现模糊的现象，故此定点扫描不适用于小型化位标器机构；圆形扫描由于不同位置点的扫描速度不同，也会造成扫描外圈图像模糊的现象。所以从工程实现角度考虑，匀速矩形扫描既可保证视场各点的扫描速度相同，也降低了对控制系统的性能要求，有利于提高图像拼接效果。
匀速矩形扫描方式下，探测器在惯性空间内按照固定的速度进行匀速运动，如图 1所示。
![]()
图 1 匀速矩形扫描的视场移动方式Figure 1. Field of view movement method for uniform rectangular scanning在图像拼接应用的领域，根据应用场景的区别，会存在特定的先验知识，对图像拼接的复杂程度进行降维，在设定拼接算法的过程中，图像的变化模型的简化可以通过增加先验知识来实现，扫描方式就是一种重要的先验知识。扫描方式的简化可以有效降低图像拼接算法的复杂度。综合以上分析，采用匀速矩形扫描作为框架的扫描搜索方案。框架扫描算法如下所示：
$$\begin{aligned} & \theta(t)=\arctan \left(y_{\mathrm{p}}(t) / \sqrt{x_{\mathrm{p}}^2(t)+z_{\mathrm{p}}^2(t)}\right)+\alpha(t) \\ & \varphi(t)=\arctan \left(z_{\mathrm{p}}(t) / x_{\mathrm{p}}(t)\right)+\beta(t)\end{aligned}$$ (1) 其中：
$$ \left\{ \begin{array}{l} \alpha (t) = \left( \begin{array}{l} 8(t - k);\, \, \;\quad \, k \leqslant t < k + 0.25\, \hfill \\ 2;\, \quad \;\, \;\, \;\, \quad \;\;k + 0.25 \leqslant t < k + 0.5 \hfill \\ 2 - 8(t - k);\, \;k + 0.5 \leqslant t < k + 0.75 \hfill \\ - 2;\, \, \;\, \quad \, \;\, \, \, \, \;\, \, k + 0.75 \leqslant t < k + 1 \hfill \\ \end{array} \right) \hfill \\ \beta (t) = \left( \begin{array}{l} 2;\, \quad \quad \quad \quad k \leqslant t < k + 0.25 \hfill \\ 2 - 8(t - k);\, \;k + 0.25 \leqslant t < k + 0.5 \hfill \\ - 2;\, \quad \;\quad \quad \, k + 0.5 \leqslant t < k + 0.75\, \hfill \\ 8(t - k);\, \, \;\quad k + 0.75 \leqslant t < k + 1 \hfill \\ \end{array} \right) \hfill \\ \end{array} \right. $$ 式中：θ(t)、ϕ(t)为框架运动的输入角度；α(t)、β(t)为随时间变化的矩形扫描补偿角度；xp(t)、yp(t)、zp(t)为惯导指向的目标三坐标位置，t为时间；k为递增整数。按上述公式，导引头框架可以完成矩形扫描。对于远距目标而言，框架的扫描在惯性空间内进行，以惯性系下目指角度进行扫描和跟踪[4]。
在扫描方式确定后，框架运动扫描周期是影响后续成像效果的重要指标，扫描速度过快，会造成图像模糊，扫描周期过大，扩展红外视场的实时性就会降低。根据自身红外探测系统特点，以探测器能量积累过程的探测器运动角度小于1个像素为原则，设计框架扫描速度为24°/s，完成一个周期3°×3°扫描需要0.5 s，故扩展视场的更新周期为0.5 s；故对框架伺服跟踪系统的进动性能要求较高，需要保证系统有较高的快速性。
两轴框架控制回路由位置环、速率环和电流环组成，在速率回路中引入MEMS陀螺传感器输出，有效提高了伺服系统的解耦能力，对弹体的震动和摆动有较高的隔离作用，为图像拼接算法的应用创造了较好的力学环境。通过增加校正环节的形式提高系统的动态性能。校正环节传递函数形式如下：
$$ \begin{array}{l} {\omega _1}(s) = {k_1}{{({T_{\text{1}}}S + {\text{1}})} \mathord{\left/ {\vphantom {{({T_{\text{1}}}S + {\text{1}})} {({T_{\text{2}}}S + {\text{1}}}}} \right. } {({T_{\text{2}}}S + {\text{1}}}}) \hfill \\ {\omega _{\text{2}}}(s) = {k_{\text{2}}}{{({T_{\text{3}}}S + {\text{1}})} \mathord{\left/ {\vphantom {{({T_{\text{3}}}S + {\text{1}})} {({T_{\text{4}}}S + {\text{1}}}}} \right. } {({T_{\text{4}}}S + {\text{1}}}}) \hfill \\ \end{array} $$ (2) 式中：S为传递函数中的复频率。
经过校正后的系统的稳定性和快速性有明显改善，具备较为稳定的矩形扫描功能，扫描状态下码盘输出如图 3所示。
2. 姿态补偿图像快速拼接算法
在框架扫描完成后，框架扫描的一个周期内就会生成连续的红外图像信息。图像拼接的主要研究内容就是确定相邻帧红外图像间的相对位置关系，通过设定融合算法处理具有相同位置映射的重复像素。两张相邻红外图像的变化模型主要取决于图像的变换方式，旋转变换、平移变换、等距变换、仿射变换、相似变换等是常用的红外图像变换方式。对于两轴框架式的扫描方式，主要存在的红外帧间图像变换方式为平移变换。平移变换指相邻两帧红外图像之间只存在俯仰和方位两自由度的位移，即一张红外图像可以看作另一张图像进行了一个二维矢量的平移变化，平移变化表征公式如下所示：
$$ \left( {\begin{array}{*{20}{c}} {x'} \\ {y'} \\ {{P_{{\text{gray}}}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 1&0&{\Delta {l_x}} \\ 0&1&{\Delta {l_y}} \\ 0&0&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} x \\ y \\ {{P_{{\text{gray}}}}} \end{array}} \right) $$ (3) 式中：x、y、Pgray分别为图像像素二维坐标和图像灰度值；Δlx、Δly为帧间图像的偏移量。更加规整的表达形式如下：
$$x^{\prime}=\boldsymbol{H}_{\mathrm{R}} X=\left(\begin{array}{cc}E & t \\ \mathbf{0}^{\mathrm{T}} & 1\end{array}\right) X$$ (4) 式中：0T代表0矩阵的转置，如此图像拼接问题可以转化为求取HR矩阵的问题。
对于两轴框架式位标器来说，由于探测器旁安装有MEMS器件，故探测器的角度信息可以通过惯性导航解算，所以：
$$ {\boldsymbol{H}_{\text{R}}} = \left( {\begin{array}{*{20}{c}} 1&0&{\int {_t^{t + \Delta T}{\omega _y}{\text{d}}t} } \\ 0&1&{\int {_t^{t + \Delta T}{\omega _x}{\text{d}}t} } \\ 0&0&1 \end{array}} \right) $$ (5) 式中：ωx、ωy为探测惯性组件陀螺输出。
根据需求探测视场的大小设定扩大的拼接图像矩阵，在计算出帧转换矩阵后，将偏移量叠加到原始图像中，采用加权平均或融合的方法得出拼接图像[6]，图 4为用于拼接的原始图像，图 5为使用惯性姿态信息拼接得到的图像。
从图 5可以看出，采用惯性姿态补偿计算偏移量的方法可以有效快速地拼接出视场较大的图像，但是存在明显的拼接瑕疵，如图 5左上侧红框标注所示，柱子边缘存在明显扭曲，代表拼接过程惯性姿态信息存在瞬时偏差导致拼接效果变差；同时如图 5中间红框所示，拼接整周期的最后一幅图像的电脑位置和第一幅图像的电脑位置严重错位，这表明惯性积累误差导致图像拼接效果变差。
3. 角点检测和互信息的图像偏移量提取算法
图像自身特征信息匹配也是图像偏移矩阵计算的重要途径。从整体（区域）和局部信息来看，图像特征信息可以分为互信息和特征点信息，通过相邻幅图像的互信息和特征点信息匹配可以得到红外图像间的偏移矩阵。互信息是基于区域的方法采用某种相似度量来比较图像间预定义尺寸下像素的相似性，当两幅图像精准匹配后，代表求得的平移矩阵最优，两幅红外图像具有最大的互信息值，因此基于此方法测量图像间的互信息，并通过微调平移矩阵，搜索其最大值。两幅图像的互信息为两幅图像间共享信息量的度量，两幅图像无关，则一幅图像完全不能给出另外一幅图像的任何信息，反之，若两幅图像相等，则互信息会达到极大值，都等于两幅图像各自的信息熵[7]。两个随机变量X和Y的互信息可以定义为：
$$ I\left( {X;Y} \right) = \mathop \sum \limits_{} p\left( {x, y} \right){\text{log}}\left( {\frac{{p\left( {x, y} \right)}}{{p\left( x \right)p\left( y \right)}}} \right) $$ (6) 式中：p(x, y)是随机变量X和Y的联合概率分布函数；p(x)和p(y)分别是随机变量X和Y的边缘概率分布函数。
互信息的计算需要对一定尺寸内的原始图像进行统计学运算，若图像自身较大，互信息极值的计算需要遍历较多的平移尺度，故计算量较大，不适宜工程上应用。所以需要图像的局部特征来初步确定偏移矩阵的范围，在偏移矩阵的可能范围内求取互信息的极值可以有效减小图像拼接的计算量。
对于红外图像来说，通过二维红外图像反应出的物体的形状信息主要集中在图像的高曲率点的汇集处。因此采用两幅红外图像的曲率特征匹配是计算偏移矩阵的常用方法。Harris角点检测算法于1998年由Harris和Stephens提出[8]，基于自相关矩阵，描述采样点附近邻域内的梯度分布：
$$ \begin{array}{l} M = \sigma _{\text{D}}^{\text{2}}g \hfill \\ \left\{ {{\sigma _{\text{t}}} \cdot \left[ {\begin{array}{*{20}{c}} {I_x^2(x, y, {\sigma _\text{D}})}&{I_x(x, y, {\sigma _\text{D}})I_y(x, y, {\sigma _\text{D}})} \\ {I_x(x, y, {\sigma _\text{D}})I_y(x, y, {\sigma _\text{D}})}&{I_y^2(x, y, {\sigma _\text{D}})} \end{array}} \right]} \right\} \hfill \\ \end{array} $$ (7) 式中：
$$ {I}_{x}(x, y, {\sigma }_{\text{D}}）=\frac{\partial }{\partial x}g({\sigma }_{\text{D}})I(x, y) \text{；} $$ $$ g(x, y, {\sigma }_{\text{D}}）=\frac{1}{2\pi {\sigma }^{2}}{\text{e}}^{-\frac{{\text{x}}^{2}+{\text{y}}^{2}}{2{\sigma }^{2}}} $$ 红外图像微分的机选应用尺度为σD的高斯核计算；σt的高斯窗在采样点的一个邻域中对其微分平滑处理表征了自相关矩阵的特征值。自相关矩阵的特征值代表了邻域中对该微分的平滑处理，角点在图像信号多个方向均存在显著变化的位置[9-10]，提取原始图像的角点图如图 6所示。
通过求取红外图像的Harris角点算子，通过图像特征点的变化量，初步识别帧间红外图像偏移量，由于红外图像自身分辨率不高，且角点的计算存在误差，所以在误差的范围内，通过寻找互信息的极大值的方法来寻求偏移矩阵的最优解[11]。
虽然通过帧间红外图像局部特征和互信息匹配的方法在一般情况下可以求取出较为精准的帧间图像偏移矩阵，合成品质较高的瞬时探测视场，但是对于低特征的红外图像存在图像特征匹配存在偏差的现象，即特征点信息和互信息匹配存在多组奇异解，此时就会造成转移矩阵计算跳动。同时，由于后续视场合成依赖于前序图像的拼接效果，故帧间的错误偏移量会在后续图片拼接中积累，引起误差积累，造成拼接效果变差而无法供后续信息处理系统使用的现象，如图 7所示。从不同偏移的帧间方差图可以得出互信息匹配矩阵图，如图 8所示。
故此，单纯使用图像信息进行帧间图像拼接，依赖于图像特征点的孤立性，为消除此问题，可以使用更为复杂的图像特征点求取方式，比如ORB（oriented FAST and rotated BRIEF）特征匹配算法，或最小核值相似区SUSAN（smallest univalue segment assimilating nucleus）角点检测匹配算法。ORB特征点检测基于BRIEF（特征描述子算法）特征点描述和FAST（feature from accelerated segment test）特征检测算法，在二者基础上，融合了二者的简洁高效的特性，具备了旋转不变特性和尺度不变换特性[12]。而SUSAN角点检测方法，应用了一种形态学法来检测特征点，而不是计算量大且对红外图像噪声敏感的求取局部梯度的计算方法[13]。虽然通过对特征点检测算法的优化，可以有效应对某些特殊情况，但是仅从图像自身维度，很难提出一种普遍使用所有工况的特征点提取算法，如图 9所示，采用互信息提取出的偏移关系因为图像复杂点的存在而引入孤立点误差，从而导致拼接效果较差，如图 10所示。
![]()
图 10 通过图像特征匹配方法拼接处的原始图像Figure 10. The original image at the stitching point using image feature matching method故本文提出引入速度维度的传感器信息进入特征点匹配算法，来提高对复杂工况的普适性。
4. 基于惯性器件的图像偏移量卡尔曼滤波补偿修正
在导引头位标器中，可以通过与探测器固连的MEMS器件测量到红外探测器在惯性空间内俯仰和方位的运动角速度。引入红外探测器速度信号对红外帧间图像偏差量进行滤波修正来消除图像特征量匹配跳点带来的影响。图像的偏移量体现了探测器在惯性系下的速度信息，将角速度积分可以得到在惯性系下的角度信息[14]。
$$ \begin{array}{l} {\varphi _x} = \int {_t^{t + \Delta T}{\omega _x}{\text{d}}t} \hfill \\ {\varphi _y} = \int {_t^{t + \Delta T}{\omega _y}{\text{d}}t} \hfill \\ \end{array} $$ (8) 式中：φx、φy为图像的偏移角度，可以作为卡尔曼滤波器的一阶信号；而陀螺输出ωx、ωy作为滤波的二阶信号。
建立框架运动的状态方程和观测方程。卡尔曼滤波理论是在高斯白噪声下对随机过程的最优估计，假定在框架扫描运动的过程中，惯性系下探测器运动的角度变化为线性变化，用线性卡尔曼滤波器可以得到估计误差最小的全局估计结果，设定系统的状态变量$X=(\dot{\boldsymbol{q}} \ddot{\boldsymbol{q}} \dddot{\boldsymbol{q}})^{\mathrm{T}}$，设定系统采样时间为T0，由运动学公式可以得出：
$$ \begin{array}{l} \dot{\boldsymbol{q}}(k + 1) = \dot{\boldsymbol{q}}(k) + {T_{\text{0}}}\ddot{\boldsymbol{q}}(k) + \frac{{\text{1}}}{{\text{2}}}{T_{\text{0}}}^{\text{2}}\dddot{\boldsymbol{q}}(k) \hfill \\ \ddot{\boldsymbol{q}}(k + 1) = \ddot{\boldsymbol{q}}(k) + {T_{\text{0}}}\dddot{\boldsymbol{q}}(k) + {T_{\text{0}}}a(k) \hfill \\ \dddot{\boldsymbol{q}}(k + 1) = \dddot{\boldsymbol{q}}(k) \hfill \\ a(k) = u(k) + \omega (k) \hfill \\ \end{array} $$ (9) 式中：$ \dot{\boldsymbol{q}}\; $、$ \ddot{\boldsymbol{q}} $、$ \dddot{\boldsymbol{q}} $为探测器运动过程的惯性系下角度、惯性系下探测器运动的角速度和角加速度；T0为系统解算周期时间；u(k)为电机输入力矩；ω(k)为干扰力矩引起的随机加速度；
加速度可以用驱动框架运动的电机力矩u(k)和框架干扰力矩引起的随机加速度ω(k)得到，ω(k)为白噪声，由此可以得出采样时刻的角速度和角加速度信息。进而得出系统状态转移矩阵如下：
$$\begin{array}{l}\left[\begin{array}{c}\hat{\dot{q}}(k) \\ \hat{\dddot{q}}(k) \\ \hat{\ddot{q}}(k)\end{array}\right]=\left[\begin{array}{ccc}1 & T_0 & \frac{1}{2} T_0^2 \\ 0 & 1 & T_0 \\ 0 & 0 & 1\end{array}\right]\left[\begin{array}{c}\hat{\dot{q}}(k-1) \\ \hat{\ddot{q}}(k-1) \\ \hat{\dddot{q}}(k-1)\end{array}\right]+ \\ \left[\begin{array}{l}0 \\ 1 \\ 0\end{array}\right] u_{k-1}+\left[\begin{array}{l}0 \\ 1 \\ 0\end{array}\right] \omega_{k-1} \end{array}$$ (10) $${{\boldsymbol{Z}}_{k}}=\left[ \begin{array}{*{35}{l}} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} \right]\left[ \begin{matrix} \hat{\dot{q}}(k) \\ \hat{\ddot{q}}(k) \\ \hat{\dddot{q}}(k) \\ \end{matrix} \right]+{{v}_{k}}$$ (11) 建立探测器运动状态方程：
$$X_k=\boldsymbol{A} X_{k-1}+\boldsymbol{B} u_{k-1}+\boldsymbol{C} \omega_{k-1}$$ 建立角度、角速度信号观测方程：
$$Z_k=\boldsymbol{H} X_k+v_k$$ 式中：vk设定为白噪声；
$$ \boldsymbol{A} = \left[ {\begin{array}{*{20}{c}} 1&{{T_{\text{0}}}}&{\frac{{\text{1}}}{{\text{2}}}{T_{\text{0}}}^{\text{2}}} \\ {\text{0}}&{\text{1}}&{{T_{\text{0}}}} \\ {\text{0}}&{\text{0}}&{\text{1}} \end{array}} \right]\text{，}\boldsymbol{B} = \left[ {\begin{array}{*{20}{c}} {\text{0}} \\ {\text{1}} \\ {\text{0}} \end{array}} \right]\text{，}\boldsymbol{C} = \left[ {\begin{array}{*{20}{c}} {\text{0}} \\ {\text{1}} \\ {\text{0}} \end{array}} \right]\text{，}\boldsymbol{H} = \left[ {\begin{array}{*{20}{c}} 1 \\ {0\, } \end{array}\begin{array}{*{20}{c}} {\, 0} \\ 1 \end{array}\, \begin{array}{*{20}{c}} {\, 0} \\ {\, 0} \end{array}} \right]。 $$ 在得出框架运动模型后，红外探测器运动帧间位移的求解问题就变为在存在二阶角速度信号的前提下，红外探测器运动角度的最优估计。
采用卡尔曼滤波融合的方式来对运动角度进行最优估计，依次建立时间更新方程和状态更新方程。
时间更新方程：
$$\begin{aligned} & \hat{X}_{\bar{k}}=\boldsymbol{A} \hat{X_{\bar{k}-1}}+\boldsymbol{B} u_{k-1} \\ & P_{\bar{k}}=\boldsymbol{A} P_{\bar{k}-1} \boldsymbol{A}^{\mathrm{T}}+\boldsymbol{Q}\end{aligned}$$ (12) 状态更新方程：
$$\begin{aligned} & K_k=\frac{P_{\overline{\mathrm{k}}} \boldsymbol{H}^{\mathrm{T}}}{H P_{\bar{K}} \boldsymbol{H}^{\mathrm{T}}+R} \\ & \hat{X}_{\bar{k}}=\hat{X_{\bar{k}-1}}+K_k\left(Z_k-\boldsymbol{H} \hat{X}_{\bar{k}}\right) \\ & \boldsymbol{P}_k=\left(I-K_k \boldsymbol{H}\right) P_{\bar{k}} \end{aligned}$$ (13) $$ L_{\mathrm{kal} X / Y}=\hat{X}_{\bar{k}}(0) $$ 式中：LkalX/Y为合成的视线角速度；Xk为状态向量；$\mathop {{X_{\bar k}}}\limits^ \wedge (0)$为状态向量的第一个元素；Pk为协方差矩阵；Q为系统噪声协方差矩阵；R为测量噪声协方差矩阵；通过框架运动的最优卡尔曼滤波估计，对帧间红外图像的位移量进行滤波，将滤波后的位移量进行图像匹配，可以有效应对低红外特征或帧间红外特征匹配存在奇异点的问题。
5. 仿真结果分析
原始红外图像（如图 11所示）经过角点检测算法，求取出角点特征图像，通过设定分割方法，得出图像的Harris特征点，如图 12所示；对特征点信息进行二值化，如图 13所示；求取特征点坐标，并对帧间的特征点变化进行匹配，得出初步的帧间偏移量，如图 14所示。
在初始帧间偏移量的基础上，构造卡尔曼修正模型，将探测器中的陀螺速度信息引入修正回路，对通过角点识别出的帧间位移量进行融合滤波，可以有效去除低红外特征图像的偏移量计算奇异点，拟合后的卡尔曼修正图像偏移量如图 14所示，卡尔曼修正后的图像偏移量有效消除了原图像特征算法偏移量中段存在跳变的问题，可以有效避免互信息奇异性带来的拼接错误。
将滤波修正后的偏移量作为帧间图像变化，实时对探测视场进行合成，在一个扫描周期后，采用新算法的拼接图像如图 15所示，相对于惯性姿态拼接算法（图 5），墙壁的扭曲和显示器的错位效果得到了明显的改善；相对于图像特征拼接算法（图 10），中间由于特征点复杂存在奇异解而造成拼接错误的问题得到了有效的规避，图像质量明显提高。
6. 结论
使用MEMS角速度信息构建速率环的框架扫描算法可以有效隔离弹体的震动；通过惯性信息将图像特征量进行融合滤波的拼接方式，可以有效对抗复杂场景对图像特征信息带来的干扰，使得图像拼接效果有了极大的提高。该算法可以有效提高拼接视场的成像质量，扩大红外视场便于目标识别，为低指向精度红外交班提供了一种可行的技术途径。该算法对于扩展红外武器的作战距离有着重要的意义。
-
期刊类型引用(3)
1. 张慧，韩新宁，韩惠丽，常莉红. 基于引导滤波二尺度分解的红外与可见光图像融合. 红外技术. 2023(12): 1216-1222 . 
本站查看
2. 杨磊，李爱华，李宜平，黄刚，李金哲，伍煌彬. 输配电线路零值绝缘子红外检测技术研究. 电工技术. 2022(13): 79-82 . 百度学术
3. 潘翀，沈鹏飞，张忠，王博，朱如桂，张颖. 基于无人机巡检图像的绝缘子串实时定位研究. 电瓷避雷器. 2020(01): 234-240 . 百度学术
其他类型引用(5)
下载:
计量
- 文章访问数: 94
- HTML全文浏览量: 31
- PDF下载量: 5
- 被引次数: 8
