鄭少杰，薛斌黨，張國(guó)琪

0 引 言

航天器通常配備多塊太陽(yáng)能帆板、可伸縮大型天線等柔性部件.柔性部件在航天器軌道機(jī)動(dòng)或進(jìn)、出入太陽(yáng)照射區(qū)域時(shí)，將發(fā)生嚴(yán)重形變和振動(dòng)，導(dǎo)致航天器主體振動(dòng)，影響航天器的任務(wù)及其運(yùn)行安全.因此，研究航天器撓性振動(dòng)在軌測(cè)量，并在此基礎(chǔ)上進(jìn)行振動(dòng)抑制、姿態(tài)鎮(zhèn)定和姿態(tài)機(jī)動(dòng)過(guò)程控制，對(duì)于保障航天器安全運(yùn)行和航天器任務(wù)的順利完成非常重要[1-2].同時(shí)，可以利用柔性部件微動(dòng)信息，感知非合作航天器的運(yùn)行態(tài)勢(shì)及姿態(tài)機(jī)動(dòng)，有利于實(shí)現(xiàn)太空態(tài)勢(shì)感知.

太陽(yáng)能板是航天器的重要能量來(lái)源，在軌監(jiān)測(cè)和辨識(shí)其動(dòng)力學(xué)特性及模態(tài)參數(shù)具有重要的應(yīng)用價(jià)值.早在1989年，噴氣推進(jìn)實(shí)驗(yàn)室(JPL)就對(duì)某大型空間柔性航天器平臺(tái)進(jìn)行了在軌模態(tài)參數(shù)辨識(shí)試驗(yàn)[3].1994年,哈伯望遠(yuǎn)鏡共進(jìn)行了兩項(xiàng)試驗(yàn)以獲得外干擾的頻率特性和航天器動(dòng)力學(xué)參數(shù)[4].NASA的Langley研究中心[5]在1993年成功將時(shí)域辨識(shí)方法應(yīng)用在“自由號(hào)”空間站的模態(tài)參數(shù)辨識(shí)試驗(yàn)中.從1996年到1997年，研究者們又在和平號(hào)空間站開(kāi)展了一系列模態(tài)參數(shù)在軌辨識(shí)試驗(yàn)[6].從1995年到1999年，日本在工程試驗(yàn)衛(wèi)星ETS-Ⅵ上[7]采用姿態(tài)控制發(fā)動(dòng)機(jī)進(jìn)行激勵(lì)，開(kāi)展了一系列在軌模態(tài)參數(shù)辨識(shí)試驗(yàn).

光學(xué)測(cè)量是基于專業(yè)攝像設(shè)備和數(shù)字圖像處理算法的測(cè)量技術(shù)，被廣泛應(yīng)用于物體幾何尺寸和空間位置的測(cè)量，已經(jīng)在航天器交會(huì)對(duì)接和姿態(tài)測(cè)量中得到應(yīng)用[8-11].另外，NASA在1996年就通過(guò)航天飛機(jī)搭載攝影測(cè)量設(shè)備對(duì)和平號(hào)空間站的太陽(yáng)能帆板進(jìn)行了模態(tài)測(cè)試[12].由于太空成像環(huán)境的特殊性以及太陽(yáng)能板的彈性機(jī)械結(jié)構(gòu)，使得常用的基于圖像灰度特征點(diǎn)匹配的形變測(cè)量方法[13-15]處理航天器太陽(yáng)能帆板振動(dòng)的在軌測(cè)量時(shí)，精度受影響較大.已有研究結(jié)果表明，在頻域中，圖像的相位信息包含了場(chǎng)景的主要結(jié)構(gòu)信息[16]，同時(shí)基于圖像相位的處理方法對(duì)噪聲和光照變化的魯棒性優(yōu)于空域中基于圖像灰度的處理方法,因此廣泛應(yīng)用于解決圖像配準(zhǔn)、圖像形變測(cè)量等問(wèn)題[17-22].基于Riesz變換構(gòu)造的單演信號(hào)是二維解析信號(hào)[23]，能夠表述圖像的局部特征，可以將圖像分解為局部幅度、局部相位和局部方向，分別表示圖像的能量、結(jié)構(gòu)和幾何信息.其中，局部相位信息能夠有效地表示信號(hào)的局部結(jié)構(gòu)特征，因此適用于微動(dòng)分析、圖像特征提取等[24-26].單演尺度空間是Felsberg基于三維拉普拉斯方程推導(dǎo)出的非高斯線性尺度空間[27]，由泊松尺度空間和其對(duì)應(yīng)的共軛泊松尺度空間構(gòu)成.單演尺度空間是由單演信號(hào)演變而來(lái)，具有多尺度多可控方向的優(yōu)點(diǎn)，因此能夠分析多尺度多方向上圖像序列局部幅值和局部相位[28]，已應(yīng)用于目標(biāo)識(shí)別[29]和圖像分析[30].

針對(duì)航天器太陽(yáng)能板振動(dòng)測(cè)量問(wèn)題，本文提出了一種基于單演多尺度相位信息的光學(xué)成像量方法，用局部相位不變替換傳統(tǒng)光流算法中的局部亮度不變，通過(guò)局部相位計(jì)算太陽(yáng)能板振動(dòng)信息.最后，通過(guò)仿真實(shí)驗(yàn)和地面實(shí)驗(yàn)驗(yàn)證了方法的有效性.

1 方法描述

1.1 太陽(yáng)能板振動(dòng)的光流表達(dá)

太陽(yáng)能板在三維空間內(nèi)發(fā)生形變或振動(dòng)時(shí)，其在相機(jī)成像表面的對(duì)應(yīng)點(diǎn)也會(huì)發(fā)生改變.如圖1所示：

圖1 相機(jī)成像等效針孔模型Fig.1 The pin-hole model of out-plane camera imaging

圖2為計(jì)算兩幀太陽(yáng)能板振動(dòng)圖像光流場(chǎng)的流程.步驟如下：a)在單演尺度空間中表示相鄰兩幀圖像；b)在對(duì)應(yīng)低空間頻率的粗糙尺度下提取圖像局部相位和局部方向特征；c)用圖像間的局部相位差計(jì)算位移相鄰兩幀圖像之間的變形量；d)基于粗略估計(jì)的形變量調(diào)整浮動(dòng)圖像，繼續(xù)計(jì)算初步變形后的浮動(dòng)圖像與參考圖像的局部相位差；e)調(diào)整單演尺度空間濾波器的尺度因子，在更精細(xì)的空間頻率下重復(fù)局部相位差實(shí)現(xiàn)多尺度逐級(jí)迭代計(jì)算.

圖2 基于單演多尺度相位的光流計(jì)算Fig.2 Monogenic multiscale phase based optical flow computing

1.2 圖像的單演尺度空間表達(dá)

單演尺度空間通過(guò)3個(gè)2D球狀正交濾波器(spherical quadrature filters,SQF)計(jì)算圖像的帶有尺度特性的局部幅度，局部相位和局部方向特征[27].球狀濾波器包含一個(gè)偶對(duì)稱和旋轉(zhuǎn)不變的帶通濾波器Be(ω)和兩個(gè)奇對(duì)稱帶通濾波器Bo1(ω)、Bo2(ω)(這里指頻域內(nèi)濾波器形式).

圖3 不同尺度帶通頻域?yàn)V波器示意圖Fig.3 Schematic diagrams of different scale bandpass frequency domain filters

奇對(duì)稱濾波器由偶對(duì)稱濾波器的Riesz變換得到，Riesz變換的頻域表達(dá)為：

(1)

其中，大寫(xiě)字母B代表傅里葉變換，ω為頻域因子，Be(ω)表達(dá)式為

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

(2)

其中，s1,s2是尺度參數(shù)，帶通濾波后的波長(zhǎng)由尺度參數(shù)來(lái)表示

(3)

基于3個(gè)濾波器響函數(shù)，可以得到在尺度參數(shù)λ0下原信號(hào)I的單演相位φ(x;λ0),方向θ(x;λ0)和幅度A(x;λ0),其中x=(x,y) 為空間坐標(biāo).為了表示方便，省去下標(biāo)λ0:

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

(4)

以及單演信號(hào)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為球狀濾波器的空域形式，*是二維卷積符號(hào).局部相位和局部方向成相位矢量的形式為：

(6)

根據(jù)局部相位的空間變化可以計(jì)算出局部頻率，即相位在局部方向上的方向?qū)?shù)，為了避免相位信息的纏繞性，采用間接計(jì)算的方式：

(7)

其中▽=[?x?y]T.通過(guò)式(6)可以根據(jù)3個(gè)SQF的響應(yīng)函數(shù)及對(duì)應(yīng)的4個(gè)偏導(dǎo)數(shù)?xp,?yp,?xq,?yq來(lái)估計(jì)局部頻率.

1.3 單演多尺度相位光流計(jì)算

設(shè)相鄰兩幀圖像為I(x,t)和I(x,t+1),對(duì)應(yīng)單演相位矢量為r(x,t)和r(x,t+1),兩幀圖像之間的位移場(chǎng)為d(x)=[d1(x)d2(x)]T.相位光流法是基于相位不變性的假設(shè)，用相位矢量表示為r(x,t+1)=r(x-d(x),t).微小位移的假設(shè)下，相位矢量的一階泰勒級(jí)數(shù)展開(kāi)為：

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

(8)

J是r的雅克比行列式.假設(shè)中心：x0=[x0y0],大小為w的窗口中所有像素以相同的位移d0移動(dòng)，則線性系統(tǒng)的方程可以寫(xiě)作：

wd0=-w,

(9)

其中，rik=?kri,rt=r(x,t+1)-r(x,t)表示r關(guān)于時(shí)間的導(dǎo)數(shù).由公式(7)可以發(fā)現(xiàn)，J的特征值為單演頻率f，特征向量為n=[cosθsinθ]T,所以J也可以表示為：

(10)

為了提升光流計(jì)算精度，本文借鑒多尺度金字塔估計(jì)大位移光流的方法，并結(jié)合高斯牛頓的迭代方法優(yōu)化：在空間頻率較低的尺度估計(jì)一個(gè)粗糙的位移，再將這個(gè)粗糙的位移代入更加精細(xì)的尺度進(jìn)一步計(jì)算殘余位移.這種由粗到精的迭代方式能夠有效地提升光流估計(jì)的精度.在實(shí)際算法實(shí)現(xiàn)過(guò)程中，由粗到精的尺度劃分是通過(guò)調(diào)整SQF帶通濾波器的中心頻率實(shí)現(xiàn)的，先利用低頻率濾波得到估計(jì)一個(gè)粗糙的位移，然后將這個(gè)位移代入更加精細(xì)的尺度(對(duì)應(yīng)中心頻率高一些)進(jìn)行下一步計(jì)算.圖4給出了單演多尺度相位光流估計(jì)計(jì)算的流程圖.

圖4 單演多尺度相位光流估計(jì)計(jì)算流程圖Fig.4 Monogenic multi-scale phase based optical flow estimation flow chart

2 微動(dòng)測(cè)量實(shí)驗(yàn)

我們用仿真微動(dòng)數(shù)據(jù)和地面實(shí)物實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證方法的有效性.所有實(shí)驗(yàn)都是在PC機(jī)上完成.

2.1 仿真微動(dòng)測(cè)量

為了定量評(píng)估方法，我們仿真了一段微動(dòng)視頻.實(shí)驗(yàn)所用4組圖像如圖5所示，針對(duì)不同光照強(qiáng)度下，以及添加高斯白噪聲后的圖像(圖5(a)、(b)、(c)光照強(qiáng)度依次增強(qiáng)，圖5(d)為(c)添加高斯白噪聲后的圖像)施加一定強(qiáng)度的位移場(chǎng)進(jìn)行變形，根據(jù)變形前后的圖像計(jì)算光流，選取光流場(chǎng)和像素位移真值的終點(diǎn)絕對(duì)誤差評(píng)價(jià)算法的精度，統(tǒng)計(jì)結(jié)果見(jiàn)表1(單位是像素).另外，我們與基于圖像灰度特征的傳統(tǒng)光流法[31]進(jìn)行了對(duì)比.仿真數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明，本文提出的方法具有亞像素級(jí)的精度，對(duì)光照和噪聲更加魯棒.

圖5 仿真微動(dòng)太陽(yáng)能板圖像序列Fig.5 Panel images with synthetic motion

表1 4組仿真圖像光流計(jì)算的絕對(duì)誤差(μ±σ)Tab.1 Endpoint error (μ±σ) on 4 simulated images

位移(0.5/pixel)傳統(tǒng)光流法[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

另外，我們針對(duì)圖5(b)的光照環(huán)境下，利用正弦函數(shù)y=0.5sin(ωt)仿真了不同振動(dòng)頻率下的3組微動(dòng)視頻，然后將本文方法的測(cè)量結(jié)果與真值作對(duì)比，誤差(單位：Hz)分別為0.0013、0.0021、0.0016.

2.2 地面實(shí)驗(yàn)太陽(yáng)能板微動(dòng)測(cè)量

地面太陽(yáng)能板微動(dòng)實(shí)驗(yàn)中，在合作模式下，我們模擬了光照變化條件下單相機(jī)拍攝的太陽(yáng)能板微動(dòng)圖像序列.圖6所示為兩幀不同光照下地面實(shí)物太陽(yáng)能板微動(dòng)圖像.基于單演多尺度相位光流法測(cè)量出圖6中黃色原點(diǎn)微動(dòng)信息，如圖7所示.

圖6 地面實(shí)驗(yàn)太陽(yáng)能板微動(dòng)圖像Fig.6 Images of micro-vibration of solar panels on ground

另外，我們可以將測(cè)量出的微動(dòng)信息進(jìn)行放大，便于觀察分析和后續(xù)的模態(tài)分析及姿態(tài)機(jī)動(dòng)感知.

圖7是太陽(yáng)能板微動(dòng)視頻中某點(diǎn)跟蹤的坐標(biāo)變換，7(a)是原始視頻的測(cè)量結(jié)果；7(b)是微動(dòng)放大3倍數(shù)后的測(cè)量結(jié)果.圖8為太陽(yáng)能板局部區(qū)域微動(dòng)測(cè)量及微動(dòng)放大后測(cè)量結(jié)果(單位是像素).

圖7 太陽(yáng)能板某點(diǎn)微動(dòng)測(cè)量結(jié)果Fig.7 Trend chart of tracking point coordinate

圖8 太陽(yáng)能局部區(qū)域微動(dòng)測(cè)量結(jié)果Fig.8 Local region micro-vibration field vision measurement

為了定量評(píng)估測(cè)量振動(dòng)信息的準(zhǔn)確性，比較了圖6中圓點(diǎn)所示測(cè)點(diǎn)單相機(jī)測(cè)量的帆板振動(dòng)振型、振動(dòng)頻率結(jié)果和激光測(cè)距儀測(cè)量的結(jié)果.在待測(cè)點(diǎn)和相機(jī)的空間位置已知情況下，可以將計(jì)算出的像素位移解算后與激光測(cè)距儀的測(cè)量結(jié)果做比較，結(jié)果如圖9所示.我們統(tǒng)計(jì)了前兩個(gè)振動(dòng)周期的振型擬合平均誤差，為4.5 mm.

從實(shí)驗(yàn)結(jié)果可以看出，我們的單演多尺度相位光流測(cè)量結(jié)果的振型與激光測(cè)距儀測(cè)量結(jié)果的振型基本吻合，因此能夠有效測(cè)量太陽(yáng)能板的振動(dòng)振型.另外，我們又分別在不同激勵(lì)下進(jìn)行了3組實(shí)驗(yàn)，單相機(jī)測(cè)量的振動(dòng)頻率結(jié)果和激光測(cè)距儀測(cè)量結(jié)果絕對(duì)偏差的平均值(單位：Hz),分別為：0.0020,0.0043和0.0036.因此，該單演多尺度相位光流方法不僅能夠提取振型信息，還能夠提取振動(dòng)頻率信息，有了振型和振動(dòng)頻率，便可以用于振動(dòng)的模型分析.

圖9 本文方法與激光測(cè)距儀比較結(jié)果Fig.9 Comparison of the results of this method with laser range finder

4 結(jié) 論

本文提出了一種單演多尺度相位光流微動(dòng)測(cè)量方法，通過(guò)單演多尺度空間中粗糙尺度局部相位到精細(xì)尺度局部相位迭代計(jì)算相位差，保證了測(cè)量精度，同時(shí)圖像相位對(duì)成像光照變化和噪聲影響魯棒性較好，因此適合于空間飛行器柔性部件的微動(dòng)測(cè)量.本文將相位局部不變假設(shè)來(lái)代替?zhèn)鹘y(tǒng)光流法中的灰義不變假設(shè)，可以有效計(jì)算出亞像素級(jí)的激動(dòng).另外，可以利用這些亞像素級(jí)的位移均實(shí)現(xiàn)微動(dòng)放大.本文提出的航天器太陽(yáng)能板微動(dòng)測(cè)量方法可以應(yīng)用于撓性航天器姿態(tài)控制和振動(dòng)抑制，還可以通過(guò)微動(dòng)測(cè)量及微動(dòng)模態(tài)分析，感知非合作航天器姿態(tài)動(dòng)機(jī).