一種地球紫外中心快速提取方法

張翰墨，盧 山，楊文博，呂春光

(1．上海市空間智能控制技術(shù)重點實驗室，上海201109;2．上海航天控制技術(shù)研究所，上海201109;3．南京大學(xué)國際地球系統(tǒng)科學(xué)研究所，南京210023)

0 引言

航天器的自主導(dǎo)航技術(shù)已經(jīng)成為各國研究的熱點，其抗干擾性和高可靠性是研究的目標(biāo)。地球紫外輻射特性具有可全天時觀測的特征，在地球邊緣處的“亮環(huán)”不受地貌和大氣的影響，可觀測到較穩(wěn)定且明顯的臨邊特征，這使得紫外敏感器具有提供高穩(wěn)定性導(dǎo)航信息的優(yōu)勢。通常，紫外敏感器通過對地球成像并對地心矢量估計來提供衛(wèi)星的導(dǎo)航信息［1－3］。紫外敏感器的研制和導(dǎo)航應(yīng)用研究始于20世紀(jì)90年代，隨著對星體紫外輻射特性研究的深入與積累，導(dǎo)航應(yīng)用研究也從中低軌發(fā)展為高軌或者深空探測的巡航段任務(wù)［4］。國內(nèi)研究的紫外敏感器應(yīng)用于我國的探月工程，該應(yīng)用首先對月球的紫外敏感器成像，然后估計月心矢量從而實現(xiàn)導(dǎo)航。該應(yīng)用中采用其8個環(huán)形不連續(xù)子視場拼接得到的大視場成像敏感器對月球成像并整合得到月球的邊緣信息，利用月球圓盤和黑暗太空背景之間的亮度反差估計得到200 km軌道高度的月心矢量及導(dǎo)航信息［2］。在此類大視場紫外敏感器基礎(chǔ)上，學(xué)者們開展了諸如坐標(biāo)系誤差［5］、姿態(tài)角影響［6］等相關(guān)的算法以補償導(dǎo)航精度。

面向高軌任務(wù)的地球紫外敏感器研究方向中，公開的研究成果集中在導(dǎo)航算法上。研究結(jié)果表明導(dǎo)航誤差的主要影響來自于地心方向的測量誤差，而不是地心距的測量誤差［4］，并且考慮到軌道高度，如何從地球紫外特征圖像中準(zhǔn)確估計地心方向是紫外導(dǎo)航敏感器研究領(lǐng)域的重點。目前行星形心的估計方法主要分為重心法［7］和擬合法［8－9］。重心法的有效性取決于目標(biāo)能量分布的均勻性、目標(biāo)邊緣的完整性，以及邊緣特征提取的準(zhǔn)確性?？紤]到地球紫外輻射特征的不均勻性，擬合法更適用于地球中心的提取。擬合法主要利用的是提取行星邊緣的特征，而梯度算子可以有效增強圖像的邊緣信息，更穩(wěn)定地描述局部圖像［10］。邊緣提取算法需要同時考慮計算機計算能力和內(nèi)存空間對其有效性的限制。本文以對地球完整一次成像的高軌敏感器為基礎(chǔ)，提出了一種基于梯度統(tǒng)計特征的邊緣提取算子，在有效檢測出地球紫外臨邊"亮環(huán)"邊緣后通過最小二乘擬合算法解出地平圓盤中心，實驗表明本文提出的算法能夠得到亞像素級精度的地心位置。

1 地球紫外圖像特征

地球的紫外輻射能量值具有穩(wěn)定性，其能量值主要與大氣臭氧分子濃度的垂直分布和吸收特性相關(guān)。

考慮到地球大氣層分子隨著海拔高度越來越少，而紫外譜段的能量主要被臭氧吸收，在280 nm～300 nm紫外波段范圍內(nèi)(見圖1)，大氣層40 km以下，大氣透過能量可近似忽略，而到海拔高度55 km時，大氣透過率可達(dá)到90%以上。此譜段的紫外輻射能量在地球表面55 km以上衰減非常小，因此可觀測到較明顯的地球臨邊特征;并且地球臨邊輻射能量由于不受到地球表面和大氣對流層、大氣平流層的影響，具有較高的穩(wěn)定性，是用于確定地心位置的有利依據(jù)。

圖1 大氣模型中臭氧值隨高度變化曲線Fig．1 Oxygen change curves with increasing altitudes in atmospheric model

由于目前國際上沒有公開的在軌的高軌紫外圖像數(shù)據(jù)，研究大多采用軟件模擬生成高軌的紫外圖像。綜合考慮軌道高度、觀測時間、觀測點得到的模擬數(shù)據(jù)中，地球的紫外輻射能量在地球臨邊呈現(xiàn)出“亮環(huán)”的現(xiàn)象。雖然由于同一天中觀測時間的不同，導(dǎo)致輻射亮度分布各異，能量較大的區(qū)域均分布在地球邊緣。

2 地球邊緣提取

由于受到軌道特性、大氣海拔高度、觀測時間等影響，地球紫外輻射特征在逐采樣點間的特性各異，并且地球邊緣“亮環(huán)”在圖像表達(dá)上呈現(xiàn)非均勻性，故采用梯度特征可以更魯棒地描述地球“亮環(huán)”特征。

2．1梯度圖計算

傳統(tǒng)的邊緣提取算法主要包括微分檢測算子和變換域檢測算子。其中微分檢測算子利用的是圖像邊緣的空間信息，基于邊緣點鄰域內(nèi)灰度信息在一階導(dǎo)數(shù)或者二階導(dǎo)數(shù)處的過零點特性實現(xiàn)邊緣提取，原理清晰且計算簡便。另一類變換域檢測算子主要是利用圖像邊緣結(jié)構(gòu)信息，計算量較大，且往往需要投票機制，占用較多的資源。結(jié)合星上計算機的存儲空間特性，微分檢測算子更加具有實用性。

選擇采用計算簡單且有效檢測的Sobel算子。通過X方向和Y方向的模板(見圖2)計算得到，具有低內(nèi)存要求的特點。

式中:G(x，y)表示圖像坐標(biāo)為(x，y)的像素的梯度值，GX(x，y)表示該點X方向梯度分量，GY(x，y)表示該點Y方向梯度分量。

圖2 Sobel算子模板Fig．2 Sobel operator templates

2．2梯度統(tǒng)計量計算

考慮到視場的因素，星空背景會使得邊緣提取中包含非地球邊緣點，進(jìn)而影響地平圓盤中心方向提取精度。因此需要采用有效的濾波算法以降低誤邊緣點的數(shù)量。

然而考慮到自主導(dǎo)航的實時性和星載計算機有限運算空間特性，傳統(tǒng)的濾波算法在取得較好去噪效果的同時，犧牲了計算效率或者內(nèi)存空間。

本文提出一種基于梯度統(tǒng)計的邊緣提取算子，在計算所得梯度圖的基礎(chǔ)上，通過計算局部二值模式(Local Binary Pattern，LBP)特征統(tǒng)計量來達(dá)到有效的濾波效果，實現(xiàn)邊緣提取。并且由于算法的簡單清晰，具有實時性和低開銷的優(yōu)點，易于星載計算機的移植。

1)LBP算子

LBP特征由Ojala等提出［11］，能有效描述局部圖像紋理的結(jié)構(gòu)信息。采用LBP特征構(gòu)建的直方圖來進(jìn)行匹配，在人臉識別、圖像分類等領(lǐng)域具有廣泛的應(yīng)用。

LBP計算公式為:

八鄰域LBP的基本模式有36個［12］，更大的鄰域計算半徑R和更多的模板采樣數(shù)P包含更多的模板算子(見圖3)。

雖然地平圓盤的紫外臨邊“亮環(huán)”特征結(jié)構(gòu)較為明顯，然而成像過程中CCD離散采樣和噪聲的影響，使得紫外輻射“亮環(huán)”梯度圖中表現(xiàn)出不均勻鋸齒狀的“毛刺”(見圖4(b))。在邊緣圖上明顯可見的“凸”型邊緣點會降低中心點的提取精度。

圖3 LBP計算模板示意圖Fig．3 Binary LBP template

圖4 邊緣“毛刺”Fig．4 Blur on edges

由于視角等因素，部分星空會進(jìn)入視場，使得梯度圖(見圖5)中會出現(xiàn)噪聲點，如圖5(b)中出現(xiàn)在地球大氣層邊緣的星點會直接影響地心提取精度。

圖5 部分原始Sobel算子邊緣梯度圖Fig．5 Partial gradient map calculated by Sobel operators

2)改進(jìn)的LBP統(tǒng)計量

本文結(jié)合邊緣點鄰域中像素點梯度方向的一致性和LBP算子的旋轉(zhuǎn)不變性，改進(jìn)的LBP統(tǒng)計量描述了地平圓盤邊界點的局部對稱結(jié)構(gòu)，可以有效濾除“凸”型邊緣點和離散噪聲點。選取LBP模板中的對稱結(jié)構(gòu)來描述地平圓盤邊界點。

分析2013年秋分零點的地球紫外輻射特性，以0．5 km間距采樣，計算得到大氣層90 km高的輻射量分布曲線(見圖6)?？梢杂^察到輻射量在15 km高度左右達(dá)到最大后急速下降，故梯度方向在邊緣處方向非突變，而存在一個漸變的過程。

圖6 2013年秋分零點某點輻射量曲線Fig．6 Atmosphere UV energy curve at AM 0 20130922

選取8個梯度方向作為梯度方向(見圖7)，計算每個鄰域內(nèi)在此8個方向所代表區(qū)域內(nèi)的梯度方向值，在LBP描述上利用X方向和Y方向的連續(xù)結(jié)構(gòu)的特點，實現(xiàn)邊緣點提取。典型邊緣點統(tǒng)計量包括:10000000，01000000等八個單方向特征量，以及10001000等連續(xù)性4個特征量。選擇5×5的鄰域尺寸進(jìn)行計算以達(dá)到濾波的效果。

圖7 LBP鄰域八方向模型Fig．7 Binary LBP template in 8 directions

2．3地平圓盤中心擬合

由于本文討論的是面向高軌的自主導(dǎo)航，在建立地平圓盤中心擬合模型時采用不考慮地球扁率的橢圓模型。當(dāng)敏感器工作姿態(tài)角為零時，成像的方程為一個正圓，而通常情況下敏感器工作下，需要通過三個姿態(tài)角俯仰角，滾轉(zhuǎn)角，偏航角將敏感器坐標(biāo)系轉(zhuǎn)換到軌道坐標(biāo)系，此時的成像為一個橢圓。目前研究多采用圓模型和橢圓模型進(jìn)行擬合［13］，本文采用橢圓模型。

定義誤差函數(shù)ES:

式中:

n為提取邊緣點的數(shù)量，z為實際焦距，k為視角半徑。

基于梯度統(tǒng)計算子提取的邊緣點坐標(biāo)(xi，yi)，通過最小化誤差函數(shù)ES，估算出地平圓盤的中心位置

式中:f為紫外敏感器焦距。

3 實驗結(jié)果

采用MOTRAN軟件對春分、夏至、秋分、冬至四天的地球紫外輻射特性進(jìn)行模擬。生成的圖像按照本文提出的方法進(jìn)行處理。

實驗比較了Sobel算子提取結(jié)果、SURF算子提取結(jié)果和本文提出改進(jìn)算法的邊緣點提取及圓心擬合結(jié)果。圖8(a)(c)(e)為較典型的紫外圖像，包括秋分傍晚、夏至凌晨、春分中午的偽彩色模擬數(shù)據(jù)。圖8(b)(d)(f)是對應(yīng)地球紫外圖像的地心提取結(jié)果，其中“×”表示提取得到的邊緣點，“·”表示通過擬合模型得到的橢圓結(jié)果。

模擬圖像大小為1046×746像素，模擬地球中心點位置為(522．5，372．5)。模擬觀測時間為春分、夏至、秋分、冬至這四節(jié)氣中每天的4個典型時間作為實驗圖像，考慮到紫外敏感器成像過程中的噪聲等因素的影響，對提取特征點數(shù)據(jù)采用隨機噪聲干擾，表1給出16幅模擬圖像在隨機噪聲干擾下的中心坐標(biāo)擬合結(jié)果。實驗結(jié)果表明地心提取精度誤差均在1個像素內(nèi)，并且處理時間在20 ms內(nèi)，滿足自主導(dǎo)航實時性的需求。本文算法僅需遍歷圖像像素點一遍，簡單快捷，內(nèi)存存儲的開銷僅為一幅圖像的大小，適合星上計算機移植。

表1 紫外圖像地平圓盤中心提取精度Table 1 Earth’s center detection results on ultraviolet image

圖8 紫外模擬數(shù)據(jù)及地球球心提取結(jié)果Fig．8 Simulated ultraviolet images and the earth extraction result

4 結(jié)論

本文實現(xiàn)了面向高軌自主導(dǎo)航的紫外特征地平圓盤中心提取。該方法采用基于梯度統(tǒng)計特征提取地球邊緣紫外輻射“亮環(huán)”，并通過最小二乘擬合估計橢圓模型參數(shù)，進(jìn)而求解得到地心位置。實驗表明，本文提出的方法可以實現(xiàn)亞像素級精度的地心提取，快速、低內(nèi)存開銷，具有較強的魯棒性。

［1］ 王裙，熊凱，張斌，等．一種基于紫外姿態(tài)敏感器和星間相對測量的星座組合導(dǎo)航方法［J］．航天控制，2013，31(5):83－89．［Wang Qun，Xiong Kai，Zhang Bin，et al．An integrated navigation method for the constellation based on ultraviolet sensor and inter-satellite relative measurements［J］．Aerospace Control，2013，31(5):83－89．］

［2］ 宋琛，張蓬蓬，張劍波，等．基于紫外敏感器的衛(wèi)星自主導(dǎo)航［J］．計算機仿真，2010，27(11):14－17．［Song Chen，Zhang Peng-peng，Zhang Jian-bo，et al．Satellite autonomous navigation based on ultraviolet sensors［J］．Computer Simulation，2010，27(11):14－17．］

［3］ 孫俊，張世杰，李葆華．利用地球紫外和恒星可見光的衛(wèi)星自主導(dǎo)航［J］．光學(xué)精密工程，2013，21(5):1192－1198．［Sun Jun，Zhang Shi-jie，Li Bao-hua．Autonomous navigation based on star light and ultraviolet Earth sensors［J］．Optics and Precision Engineering，2013，21(5):1192－1198．］

［4］ 喬國棟，李鐵壽，王大軼．基于紫外敏感器的地月轉(zhuǎn)移軌道慢旋探測器自主導(dǎo)航算法［J］．宇航學(xué)報，2009，30(2):492－496． ［Qiao Guo-dong，Li Tei-shou，Wang Da-yi．Autonomous navigation algorithm for slow-spin probe in earthlunar transfer orbit based on ultraviolet sensors［J］．Journal of Astronautics，2009，30(2):492－496．］

［5］ 魏春嶺，張斌，張春青．一種姿態(tài)機動輔助下的天文導(dǎo)航系統(tǒng)偏差自校準(zhǔn)方法［J］．宇航學(xué)報，2010，31(1):93－97．［Wei Chun-ling，Zhang Bin，Zhang Chun-qing．An attitude maneuvering aided self-calibration algorithm for celestial autonomous navigation system［J］．Journal of Astronautics，2010，31(1):93－97．］

［6］ 李茂登，荊武興，黃翔宇．考慮月球扁率修正的月球衛(wèi)星自主導(dǎo)航［J］．宇航學(xué)報，2012，33(7):896－902．［Li Maodeng，Jing Wu-xing，Huang Xiang-yu．Autonomous navigation for Lunar satellite with Lunar oblateness correction［J］．Journal of Astronautics，2012，33(7):896－902．］

［7］ 楊鵬，謝立，劉濟林，等．基于Zernike矩的高精度太陽圖像質(zhì)心提取算法［J］．宇航學(xué)報，2011，32(9):1963－1970．［Yang Peng，Xie Li，Liu Ji-lin，et al．Zernike moment based high-accuracy Sun image centroid algorithm［J］．Journal of Astronautics，2011，32(9):1963－1970．］

［8］ 陳闊，馮華君，徐之海，等．亞像素精度的行星中心定位算法［J］．光學(xué)精密工程，2013，21(7):1181－1190．［Chen Kuo，F(xiàn)eng Hua-jun，Xu Zhi-hai，et al．Sub-pixel location algorithm for planetary center measurement［J］．Optics and Precision Engineering，2013，21(7):1181－1190．］

［9］ 張瑩瑾，秦世引，胡曉惠．航天器自主交會對接的高精度目標(biāo)定位方法［J］．宇航學(xué)報，2014，35(12):1430－1436．［Zhang Ying-jin，Qin Shi-yin，Hu Xiao-hui．A high accuracy target positioning method for autonomous rendezvous and docking between spacecrafts［J］．Journal of Astronautics，2014，35(12):1430－1436．］

［10］Ahmed F．Gradient directional pattern:a robust feature descriptor for facial expression recognition［J］．Electronics Letters，2012，48(19):1203－1204．

［11］Ojala T，Valkealahti K，Oja E，et al．Texture discrimination with multidimensional distributions of signed gray level differences［J］．Pattern Recognition，2001，34(3):727－739．

［12］Grubler A C．Multiresolution Gray-scale and rotation invariant texture classification with local binary patterns［J］．IEEE Transactions on PAMI，2001，24(7):971–987．

［13］ 劉軍，韓潮．基于圓球與橢球的地球敏感器地心矢量算法分析［J］．航天控制，2008，26(2):86－91．［Liu Jun，Han Chao．Arithmetic analysis of geocentric vector based on Earth sensor for spheral model and oblate model［J］．Aerospace Control，2008，26(2):86－91．］