航天器太阳能板振动单演多尺度相位测量*

郑少杰，薛斌党，张国琪

0 引 言

航天器通常配备多块太阳能帆板、可伸缩大型天线等柔性部件.柔性部件在航天器轨道机动或进、出入太阳照射区域时，将发生严重形变和振动，导致航天器主体振动，影响航天器的任务及其运行安全.因此，研究航天器挠性振动在轨测量，并在此基础上进行振动抑制、姿态镇定和姿态机动过程控制，对于保障航天器安全运行和航天器任务的顺利完成非常重要[1-2].同时，可以利用柔性部件微动信息，感知非合作航天器的运行态势及姿态机动，有利于实现太空态势感知.

太阳能板是航天器的重要能量来源，在轨监测和辨识其动力学特性及模态参数具有重要的应用价值.早在1989年，喷气推进实验室(JPL)就对某大型空间柔性航天器平台进行了在轨模态参数辨识试验[3].1994年,哈伯望远镜共进行了两项试验以获得外干扰的频率特性和航天器动力学参数[4].NASA的Langley研究中心[5]在1993年成功将时域辨识方法应用在“自由号”空间站的模态参数辨识试验中.从1996年到1997年，研究者们又在和平号空间站开展了一系列模态参数在轨辨识试验[6].从1995年到1999年，日本在工程试验卫星ETS-Ⅵ上[7]采用姿态控制发动机进行激励，开展了一系列在轨模态参数辨识试验.

光学测量是基于专业摄像设备和数字图像处理算法的测量技术，被广泛应用于物体几何尺寸和空间位置的测量，已经在航天器交会对接和姿态测量中得到应用[8-11].另外，NASA在1996年就通过航天飞机搭载摄影测量设备对和平号空间站的太阳能帆板进行了模态测试[12].由于太空成像环境的特殊性以及太阳能板的弹性机械结构，使得常用的基于图像灰度特征点匹配的形变测量方法[13-15]处理航天器太阳能帆板振动的在轨测量时，精度受影响较大.已有研究结果表明，在频域中，图像的相位信息包含了场景的主要结构信息[16]，同时基于图像相位的处理方法对噪声和光照变化的鲁棒性优于空域中基于图像灰度的处理方法,因此广泛应用于解决图像配准、图像形变测量等问题[17-22].基于Riesz变换构造的单演信号是二维解析信号[23]，能够表述图像的局部特征，可以将图像分解为局部幅度、局部相位和局部方向，分别表示图像的能量、结构和几何信息.其中，局部相位信息能够有效地表示信号的局部结构特征，因此适用于微动分析、图像特征提取等[24-26].单演尺度空间是Felsberg基于三维拉普拉斯方程推导出的非高斯线性尺度空间[27]，由泊松尺度空间和其对应的共轭泊松尺度空间构成.单演尺度空间是由单演信号演变而来，具有多尺度多可控方向的优点，因此能够分析多尺度多方向上图像序列局部幅值和局部相位[28]，已应用于目标识别[29]和图像分析[30].

针对航天器太阳能板振动测量问题，本文提出了一种基于单演多尺度相位信息的光学成像量方法，用局部相位不变替换传统光流算法中的局部亮度不变，通过局部相位计算太阳能板振动信息.最后，通过仿真实验和地面实验验证了方法的有效性.

1 方法描述

1.1 太阳能板振动的光流表达

太阳能板在三维空间内发生形变或振动时，其在相机成像表面的对应点也会发生改变.如图1所示：

图1 相机成像等效针孔模型Fig.1 The pin-hole model of out-plane camera imaging

图2为计算两帧太阳能板振动图像光流场的流程.步骤如下：a)在单演尺度空间中表示相邻两帧图像；b)在对应低空间频率的粗糙尺度下提取图像局部相位和局部方向特征；c)用图像间的局部相位差计算位移相邻两帧图像之间的变形量；d)基于粗略估计的形变量调整浮动图像，继续计算初步变形后的浮动图像与参考图像的局部相位差；e)调整单演尺度空间滤波器的尺度因子，在更精细的空间频率下重复局部相位差实现多尺度逐级迭代计算.

图2 基于单演多尺度相位的光流计算Fig.2 Monogenic multiscale phase based optical flow computing

1.2 图像的单演尺度空间表达

单演尺度空间通过3个2D球状正交滤波器(spherical quadrature filters,SQF)计算图像的带有尺度特性的局部幅度，局部相位和局部方向特征[27].球状滤波器包含一个偶对称和旋转不变的带通滤波器Be(ω)和两个奇对称带通滤波器Bo1(ω)、Bo2(ω)(这里指频域内滤波器形式).

图3 不同尺度带通频域滤波器示意图Fig.3 Schematic diagrams of different scale bandpass frequency domain filters

奇对称滤波器由偶对称滤波器的Riesz变换得到，Riesz变换的频域表达为：

(1)

其中，大写字母B代表傅里叶变换，ω为频域因子，Be(ω)表达式为

Be(ω)=e-|ω|s1-e-|ω|s2

(2)

其中，s1,s2是尺度参数，带通滤波后的波长由尺度参数来表示

(3)

基于3个滤波器响函数，可以得到在尺度参数λ0下原信号I的单演相位φ(x;λ0),方向θ(x;λ0)和幅度A(x;λ0),其中x=(x,y) 为空间坐标.为了表示方便，省去下标λ0:

φ(x)=arg(p(x)+iq(x)),

(4)

以及单演信号MSI(λ0)：

MSI=p(x)+iq(x)

(5)

其中,p(x)=(I*be)(x)，q1(x)=(I*bo1)(x)，q2(x)=(I*bo2)(x)，q(x)=[q1(x)q2(x)]T，be、bo1、bo2为球状滤波器的空域形式，*是二维卷积符号.局部相位和局部方向成相位矢量的形式为：

(6)

根据局部相位的空间变化可以计算出局部频率，即相位在局部方向上的方向导数，为了避免相位信息的缠绕性，采用间接计算的方式：

(7)

其中▽=[∂x∂y]T.通过式(6)可以根据3个SQF的响应函数及对应的4个偏导数∂xp,∂yp,∂xq,∂yq来估计局部频率.

1.3 单演多尺度相位光流计算

设相邻两帧图像为I(x,t)和I(x,t+1),对应单演相位矢量为r(x,t)和r(x,t+1),两帧图像之间的位移场为d(x)=[d1(x)d2(x)]T.相位光流法是基于相位不变性的假设，用相位矢量表示为r(x,t+1)=r(x-d(x),t).微小位移的假设下，相位矢量的一阶泰勒级数展开为：

r(x-d(x),t)≈r(x,t)-J(x,t)·d(x)

(8)

J是r的雅克比行列式.假设中心：x0=[x0y0],大小为w的窗口中所有像素以相同的位移d0移动，则线性系统的方程可以写作：

wd0=-w,

(9)

其中，rik=∂kri,rt=r(x,t+1)-r(x,t)表示r关于时间的导数.由公式(7)可以发现，J的特征值为单演频率f，特征向量为n=[cosθsinθ]T,所以J也可以表示为：

(10)

为了提升光流计算精度，本文借鉴多尺度金字塔估计大位移光流的方法，并结合高斯牛顿的迭代方法优化：在空间频率较低的尺度估计一个粗糙的位移，再将这个粗糙的位移代入更加精细的尺度进一步计算残余位移.这种由粗到精的迭代方式能够有效地提升光流估计的精度.在实际算法实现过程中，由粗到精的尺度划分是通过调整SQF带通滤波器的中心频率实现的，先利用低频率滤波得到估计一个粗糙的位移，然后将这个位移代入更加精细的尺度(对应中心频率高一些)进行下一步计算.图4给出了单演多尺度相位光流估计计算的流程图.

图4 单演多尺度相位光流估计计算流程图Fig.4 Monogenic multi-scale phase based optical flow estimation flow chart

2 微动测量实验

我们用仿真微动数据和地面实物实验数据验证方法的有效性.所有实验都是在PC机上完成.

2.1 仿真微动测量

为了定量评估方法，我们仿真了一段微动视频.实验所用4组图像如图5所示，针对不同光照强度下，以及添加高斯白噪声后的图像(图5(a)、(b)、(c)光照强度依次增强，图5(d)为(c)添加高斯白噪声后的图像)施加一定强度的位移场进行变形，根据变形前后的图像计算光流，选取光流场和像素位移真值的终点绝对误差评价算法的精度，统计结果见表1(单位是像素).另外，我们与基于图像灰度特征的传统光流法[31]进行了对比.仿真数据实验结果表明，本文提出的方法具有亚像素级的精度，对光照和噪声更加鲁棒.

图5 仿真微动太阳能板图像序列Fig.5 Panel images with synthetic motion

表1 4组仿真图像光流计算的绝对误差(μ±σ)Tab.1 Endpoint error (μ±σ) on 4 simulated images

位移(0.5/pixel)传统光流法[31]本文方法图2(a)0.011±0.0050.0053±0.0032图2(b)0.023±0.0210.0031±0.0020图2(c)0.056±0.0140.0060±0.0040图2(d)0.127±0.0360.0063±0.0056

另外，我们针对图5(b)的光照环境下，利用正弦函数y=0.5sin(ωt)仿真了不同振动频率下的3组微动视频，然后将本文方法的测量结果与真值作对比，误差(单位：Hz)分别为0.0013、0.0021、0.0016.

2.2 地面实验太阳能板微动测量

地面太阳能板微动实验中，在合作模式下，我们模拟了光照变化条件下单相机拍摄的太阳能板微动图像序列.图6所示为两帧不同光照下地面实物太阳能板微动图像.基于单演多尺度相位光流法测量出图6中黄色原点微动信息，如图7所示.

图6 地面实验太阳能板微动图像Fig.6 Images of micro-vibration of solar panels on ground

另外，我们可以将测量出的微动信息进行放大，便于观察分析和后续的模态分析及姿态机动感知.

图7是太阳能板微动视频中某点跟踪的坐标变换，7(a)是原始视频的测量结果；7(b)是微动放大3倍数后的测量结果.图8为太阳能板局部区域微动测量及微动放大后测量结果(单位是像素).

图7 太阳能板某点微动测量结果Fig.7 Trend chart of tracking point coordinate

图8 太阳能局部区域微动测量结果Fig.8 Local region micro-vibration field vision measurement

为了定量评估测量振动信息的准确性，比较了图6中圆点所示测点单相机测量的帆板振动振型、振动频率结果和激光测距仪测量的结果.在待测点和相机的空间位置已知情况下，可以将计算出的像素位移解算后与激光测距仪的测量结果做比较，结果如图9所示.我们统计了前两个振动周期的振型拟合平均误差，为4.5 mm.

从实验结果可以看出，我们的单演多尺度相位光流测量结果的振型与激光测距仪测量结果的振型基本吻合，因此能够有效测量太阳能板的振动振型.另外，我们又分别在不同激励下进行了3组实验，单相机测量的振动频率结果和激光测距仪测量结果绝对偏差的平均值(单位：Hz),分别为：0.0020,0.0043和0.0036.因此，该单演多尺度相位光流方法不仅能够提取振型信息，还能够提取振动频率信息，有了振型和振动频率，便可以用于振动的模型分析.

图9 本文方法与激光测距仪比较结果Fig.9 Comparison of the results of this method with laser range finder

4 结 论

本文提出了一种单演多尺度相位光流微动测量方法，通过单演多尺度空间中粗糙尺度局部相位到精细尺度局部相位迭代计算相位差，保证了测量精度，同时图像相位对成像光照变化和噪声影响鲁棒性较好，因此适合于空间飞行器柔性部件的微动测量.本文将相位局部不变假设来代替传统光流法中的灰义不变假设，可以有效计算出亚像素级的激动.另外，可以利用这些亚像素级的位移均实现微动放大.本文提出的航天器太阳能板微动测量方法可以应用于挠性航天器姿态控制和振动抑制，还可以通过微动测量及微动模态分析，感知非合作航天器姿态动机.