基于OTSU二值化的航天器在轨环境预测方法

2020-04-29郜泽霖

计算机测量与控制 2020年4期

魏晋，郜泽霖

(吉林大学仪器科学与电气工程学院，长春 130000)

0 引言

航天器在飞行期间可能会受到由于起飞、入轨时序动作引起的噪声、振动、冲击等环境影响，使航天器出现偏离轨道的现象发生[1]。对于目前生活环境中，航空运输是非常重要的环节，无论是出行还是邮递东西，都少不了一些飞行工具，而如何能在最短时间内把飞机上的东西送达，尤为重要[2]。针对如何在最短路径内躲避开一些因为天气原因而不能通过的区域，是当下航天器在轨环境预测最为关键的问题。冲击环境大都以高频为主的宽频带、瞬时状态特点，很难利用分析方法精准预测航天器在轨环境。因此，有学者利用统计能量法和虚拟模态综合法预测在轨环境，虽然预测速度较快，但受到已知结构某些冲击环境影响，预测误差较大，采用上述数值类方法不太适用。基于此，提出了基于OTSU(最大类间方差法)二值化的航天器在轨环境预测方法，该方法需先获取航天器像素版图，应用OTSU算法对其进行二值化处理，将航天器的每一个模块都处理为一个小像素块，在其中设置一些不可以通过的点，模拟因为天气原因而使航天器无法通过的区域，应用BFS宽度优先搜索算法求解出航天器从一个地方到达另一个地方。该预测方法结合了OTSU算法处理图像与BFS算法寻找最短路径理论，可以求解像素点间最短路径，实现航天器在轨环境预测。

1 航天器像素版图获取

在航天器像素版图获取过程中，使用单目视觉光学测量法测量航天器在轨相对位置的姿态时，需先合理选择航天器特征标志[3]。在单目视觉光学测量时以圆周作为靶面像素测量特征，通过相机视野范围内得到的靶面的位置姿态信息，以靶面圆周直径进行测量不受空间圆位置姿态的影响，在相机的视野范围内都会找到任意一条与像面平行的直径，比采用其它测量方法要精准得多。由于这种特征直径测量法在透视投影中直径与像面平行，而且直径的长度与像长之间的比例关系不会发生变化，因此在透视投影中其像点的像面坐标与靶面圆心之间形成一种精确的比例关系[4]。基于这种以圆心的物、像空间坐标和特征直径的原长与其像长作为特征参数，来计算靶面的相对位置姿态参数的计算方式，具有简单精准、无误差的优点。再进行图像处理时应先准确找到特征直径和圆心的像的位置，才能完成精准测量。

利用摄像机获得的单目透视投影的图片中，空间圆周的透视投影像为椭圆形，具有特征直径的像的有形部分位于椭圆圆周上的两个点上，而圆心的像点位置却不在椭圆像的中心位置，在进行图像处理时无法快速准确地找到圆心和特征直径的像的位置[5]，因此设计了如图1所示的三圆靶进行点特征的测量和提取。

图1 三圆靶示意图

根据三圆靶示意图所示，当摄像头采集图片时，投影到摄像机二维平面上的坐标系为直角坐标系的形式，在等腰直角三角形的三个顶点上分别为三个圆的圆心，以它们之间相交的两直线段的交点作为圆心的标志特征[6]。利用这种特征直径的计算方法能够准确地求出这两条直线段在像面上的像的直线方程，就可以精确找到圆心的像点位置。在运用亚像素处理方法处理数字照片求出像面的椭圆方程，通过计算结果就可获得特征直径及其像长的准确比例。

为实现精准测量将三圆靶固定安装在目标航天器上，利用摄像机对靶进行拍照[7]。首先把目标坐标系的原点O以圆心O1为中心坐标圆点，把圆心O1和O2的连线方向设为x轴，O1和O3的连线方向则为y轴，对三个坐标轴的方向进行目标坐标系与目标航天器坐标系的相对位置姿态参数进行测量，然后在进行追踪航天器坐标系的相对位置和姿态参数的获取。以追踪航天器摄像机的光心作为相机坐标系的原点O，以相机主光轴向外方向为z轴，以x轴为平行于相机CCD面阵平行的方向，y轴平行于拍摄方向，求出追踪航天器的坐标系中的相对位置和姿态参数，获取像素板图。

2 基于OTSU二值化图像分割

为确保图像分割的精准性，对图像二值化分割进行改进研究，采用了一种在原有二维算法的基础上进行改进的OTSU的算法，并引入类内最小离散度与邻域平均灰度值，保证了域值分割的完整性[8]。两个一维OTSU算法是由一维直方图获取一个灰度阈值，再由相邻域一维直方图得到一个平均灰度阈值，不仅保障了图像本身的灰度信息，还考虑了其邻域像素点的信息，同时还起到了良好的去噪作用。两个一维OTSU算法解决了传统计算中存在的计算复杂、难度大、易误差的缺点，提高了算法的精准性和可靠性。

通过对二维直方图中(a，b)出现的频数nab，判断出像素灰度值a出现的频数和领域灰度值b出现的频数，分别为:

(1)

(2)

根据计算结果由概率分布Qab得到像素灰度值a和领域灰度值b的一维直方图分布，其概率分别为：

(3)

(4)

在公式(3)、(4)中，a=0,1,...,m-1；b=0,1,...,m-1。

一般情况下，在二维OTSU算法中，边缘点和噪声点较小，即噪声和边缘部分所对应的概率可以近似为0，也就是：

(5)

因此，可以得到像素灰度值a所对应的一维直方图中目标和背景所占的比例。

根据获取的相关信息，引入标准萤火虫算法分割图像，如图2所示。

图2 基于OTSU二值化图像分割流程

引入步长调整函数扩大了寻优区域，可以快速高效地搜寻到全局最优值[9]。读入待分割的图像，定义目标函数，设定最大迭代次数，并设置相关参数。初始化萤火虫位置，设定萤火虫位置目标函数，并计算每个萤火虫的绝对亮度值。根据新的位置计算新的亮度值，找到的最优解进行图像分割。

3 基于宽度优先搜索算法像素点间最短路径求取

在像素点最短途径求取过程中，采取优先搜索算法进行求取，通过计算机软件建立三角形网格模型。分别以三角形网格模型上的节点和边作为图，与航天器在轨的特征构成边的权值，来计算塑件流长比例[10]。根据边的厚度平均值的倒数作为边的权值，进而生成一个二值化图，然后采用宽度优先搜索算法，来计算图上任意两点之间的最短距离，从而得到一个最初的最短路径[11]。最后再对计算出的最短路径周围的三角形网格边进行细分，生成新的二值化子图，应是采取宽度优先搜索算法求取新的最短路径，经过反复迭代细分与计算过程，最终得到趋向真实的最短路径。

采用宽度优先搜索算法得到的近似最短路径实际上是三角形网格模型中首尾相连的边的集合，而并不是实际网格面上的最短路径[12]，最初近似最短路径如图3所示。

图3 最初近似最短路径

如图3所示，采用宽度优先搜索算法得到的近似路径，与虚线表示的实际路径存在较大距离，因此必须将最短路径周围所经过的顶点、三角形单元的边进行细分，再采用宽度优先搜索算法计算得到细分的顶点、边的最短路径，如此循环细分计算，直到得到符合标准的最短路径:

|H(m)-H(m+1)|<σ

(6)

由式(6)可知：H表示数组用于存放从初始原点到各个点上的二值化最短路径长度值；σ表示容许误差。

1)将二值化子图最短路径所经过节点加入到下一个子图之中；

2)对于最短路径上任意节点转移，若其相邻点满足二值化子图，则把点和边加入二值化子图，并做上相应的标记；

3)对于最短路径上任意节点，若其相邻点满足边界值为0，并将点和边加入二值化子图；

4)对此时子图边数进行记录；

5)把图中的边再对边进行细分，使其产生新的节点和边，并把这些新边和点加入到子图中，同时从子图中删去边，子图点集生成示意图如图4所示。

图4 子图点集生成示意图

图4中，点1和点2为二值化子图的最短路径上的点，且1、2、3属于原图上的节点，14、45、52、16、67、GC为原边上新生成的边，细实线为由不在同一条边上细分点相连生成的边；

6)计数K减一，若K不为零，则返回步骤5)；

7)对同一三角形网格上不同边上的点，需边加入子图。

根据上述流程，求取的基于宽度优先搜索算法像素点间最短路径如图5所示。

图5 最短路径生成

通过上述预测方法，可以生成最短路径，即两点之间最短路径，以此求取像素点间最短路径，完成航天器在轨环境预测。

4 结构动力学实验

在实际航天工程中，由于冲击环境限制条件较多，通常使用地面测量产品冲击环境模拟航天器在轨环境，以此验证基于OTSU二值化的航天器在轨环境预测方法研究的合理性。

4.1 航天器结构分析

根据航天器结构自身承受载荷功能，可将整个航天器结构进行主结构和此结构划分。其中主结构是整个航天器的主要承载力结构，将其与运载火箭对接，可将载荷传递到航天器之中。其结构如图6所示。

1.圆筒体;2.远地点发动机框;3.锥壳;4.包带分离框;5.天线支承框;6.外面板;7.内面板;8.燃料储箱支承框。图6 航天器主结构

次结构是主结构分支出来的航天器其余结构，包括仪器设备安装和卫星外壳，在考虑卫星发射载荷作用下，充分考虑在轨环境噪声作用，辅助航天器稳定运行。

航天器结构如图7所示。

图7 航天器结构

航天器由若干个不同系统组成，包括通信卫星转发器、通信天线、可见光照相机、电视摄像机、热控制系统、电源系统和姿态控制系统等。使用有限元模型设计的航天器结构，可以计算各种边界条件下的动力响应。针对实际对象建立的模型，往往会因为边界连接条件简化问题而产生在轨环境预测误差。为了将误差控制到允许范围内，必须进行结构动力学实验。

4.2 在轨力学环境测量平台设计

通过布置航天器内部的传感器，来监测在轨结构受到冲击响应。常规的冲击传感器容易受到振动信号的干扰，故此在安装冲击传感器的同时，需要配合高、低频振动传感器，以此辅助模数转换模块对振动信号与冲击信号的分类传送。在轨力学环境预测平台组成如图8所示。

图8 在轨力学环境预测平台

该平台是由传感器、电缆线和采集单元组成的，其中传感器负责采集冲击信号；采集单元对采集到的信号进行存储，并通过以太网传输到交换机处，供地面人员对采集到的数据进行分析；利用电缆进行各单元之间信号传输。

4.3 冲击环境预示

在轨力学环境测量平台共有30个传感器，150路模拟通道。表1给出了部分力学环境预测平台位置和平台测试点布局。

表1 测平台位置和平台测试点布局

鉴定级试验条件设置如表2所示。

表2 正弦振动试验条件

在试验过程中，设置了冲击环境，预测到的两点5次冲击响应谱如图9所示。

图9 两点5次冲击响应谱

预测平台-D1，对应平台-P1的冲击响应谱：在频率为100 Hz时，第4次的冲击响应谱为20 g，其余冲击响应谱基本为0；在频率为1 000 Hz时，5次冲击响应谱相差较小，都维持在600～800 g范围内波动；在频率为10 000 Hz时，5次冲击响应谱依次为1 100 g、990 g、900 g、990 g、600 g。

预测平台-D2，对应平台-P2的冲击响应谱：在频率为100 Hz时，5次冲击响应谱依次为25 g、20 g、45 g、50 g、90 g；在频率为1 000 Hz时，5次冲击响应谱依次为1 400 g、1 300 g、1 100 g、1 100 g、1 100 g；在频率为1 000 Hz时，5次冲击响应谱依次为1 400 g、1 400 g、1 400 g、1 300 g、1 200 g。