一種基于LM算法的激光足印中心提取方法

王志文,李 松,羅 敏

(武漢大學(xué)電子信息學(xué)院,湖北 武漢 430072)

1 引 言

星載激光測高儀是一種可以實(shí)現(xiàn)對地表高精度三維觀測的天基測量設(shè)備,而其精細(xì)的真實(shí)性校驗(yàn)是保證激光測高儀數(shù)據(jù)產(chǎn)品在各個(gè)領(lǐng)域得到有效應(yīng)用的前提[1]。2003年,首個(gè)對地觀測的激光測高儀(Geoscience Laser Altimeter System,GLAS)于美國成功發(fā)射,隨后,相關(guān)研究人員提出了多種真實(shí)性校驗(yàn)的方式[2-4]。其中,L.A.Magruder等提出一種獨(dú)立于其他校驗(yàn)方式的地面探測系統(tǒng)[4],通過地面探測系統(tǒng)捕獲星載激光測高儀精確的激光足印位置,來實(shí)現(xiàn)對星載激光測高儀的真實(shí)性校驗(yàn),而高精度的激光足印中心提取方法是獲取激光足印精確位置的關(guān)鍵。地面探測系統(tǒng)捕獲到的激光足印,在傳輸?shù)倪^程中,會受到大氣湍流的影響,使得最終到達(dá)地面的激光足印的能量空間分布發(fā)生一定的變化,如偏轉(zhuǎn)、擴(kuò)展[5-6],因此,對于激光足印中心的提取精度,除了依賴于激光足印本身的高斯能量分布特性,還會受到大氣湍流噪聲的影響。

傳統(tǒng)的激光足印中心定位算法有灰度重心法[7](Gray Centroid Method,GCM)、圓或橢圓擬合法[8](Elliptic Fitting Method,EFM)、基于二維向量空間圓擬合法[9](2D Vector Space Circle Fitting,2DCF)等等,這些方法要么是基于整個(gè)激光足印的能量分布,如灰度重心法,要么是基于激光足印的幾何形狀,如圓擬合法?；诩す庾阌〉哪芰糠植贾行奶崛》椒▽τ诩す庾阌〉某跏寄芰糠植驾^為嚴(yán)格,當(dāng)噪聲比較小,提取的精度較好,但不適用于噪聲相對比較大的情形且沒有充分利用激光足印的空間分布信息；基于幾何形狀的提取方法,考慮了激光足印的幾何分布,但是當(dāng)邊緣或整個(gè)空間分布不是規(guī)則的圓或者橢圓,那么就會產(chǎn)生比較大的誤差。

本文提出一種基于LM優(yōu)化算法的激光足印中心提取方法,對于激光足印數(shù)據(jù)的擬合過程比較符合激光足印捕獲數(shù)據(jù)的數(shù)據(jù)格式以及激光足印的光強(qiáng)分布特點(diǎn),不會依賴于每個(gè)能級邊緣形狀的規(guī)則程度,能從全局上找到一個(gè)最優(yōu)的激光足印中心,提取精度高,并且抗干擾能力較強(qiáng)。

2 基于LM優(yōu)化算法的激光足印中心提取原理

星載激光測高儀發(fā)射出的激光經(jīng)過大氣傳輸后,被地面探測系統(tǒng)捕獲,然后通過相關(guān)算法提取可以得到激光足印的中心,從而實(shí)現(xiàn)星載激光測高儀的在軌校驗(yàn)。激光足印輪廓在大多數(shù)情況下服從橢圓高斯分布[10],可用下式表示:

(1)

其中,α=(x,y)為坐標(biāo)信息,β=(a,b,r1,r2,θ,E)為參數(shù)列表;a為x方向上激光足印中心坐標(biāo)；b為y方向上激光足印中心坐標(biāo)；r1為x方向上橢圓高斯激光足印的束腰半徑；r2為y方向上橢圓高斯激光足印的束腰半徑,當(dāng)r1=r2時(shí),足印的長軸和短軸相等,此時(shí)激光足印呈圓高斯分布；θ為對應(yīng)的橢圓高斯激光足印的方位角,取值范圍為:[0,π/4]；E表示橢圓高斯激光足印的強(qiáng)度信息。N(α)為相應(yīng)的大氣湍流噪聲。

根據(jù)采集到的橢圓高斯激光足印信號的特點(diǎn),采用一種基于LM算法的激光足印中心提取方法,首先,構(gòu)建最小二乘目標(biāo)函數(shù),確定擬合的橢圓高斯激光足印初始參數(shù):激光足印中心的位置信息、強(qiáng)度信息、長短軸的大小和方位角；然后計(jì)算最小二乘目標(biāo)函數(shù)的值,使用LM算法更新參數(shù),使每次更新后的目標(biāo)函數(shù)的值比上一次迭代的值小,否則,重新更新參數(shù),即按照目標(biāo)函數(shù)減小的方向進(jìn)行迭代；當(dāng)目標(biāo)函數(shù)的值小于預(yù)先設(shè)置的閾值,或者達(dá)到迭代次數(shù)時(shí),即停止迭代,此時(shí)擬合函數(shù)與輸入數(shù)據(jù)的殘差最小,而對應(yīng)的參數(shù)(a,b)就是待求的足印中心。整個(gè)處理的過程如圖1所示。

圖1 基于LM算法的激光足印中心提取方法Fig.1 Laser footprint center extraction method based on LM algorithm

2.1 地面探測系統(tǒng)及其相應(yīng)數(shù)據(jù)格式

針對我國自主研制的搭載于資源三號02星的激光測高儀的在軌校驗(yàn)實(shí)驗(yàn),文獻(xiàn)[1]中提出了相應(yīng)的地面足印探測系統(tǒng),其中,用于捕獲激光足印的能量探測器有8個(gè)量化能級。假設(shè)布設(shè)方式為30×30方陣,布設(shè)間距為4 m,激光足印半徑30 m,在激光足印過境的位置固定好探測器后會獲取到激光足印的8能級30×30的采樣數(shù)據(jù),并且包含足印的位置信息。

2.2 大氣湍流噪聲

星載激光測高儀所發(fā)射的激光探測信號在大氣中傳播的過程中,主要受到了大氣湍流的影響而發(fā)生畸變。Kolomogorov等建立了大氣湍流的理論模型后,由于澤尼克多項(xiàng)式低頻部分占主要部分,高頻成分較少,這與大氣湍流的理論模型表現(xiàn)一致,因而Noll等采用澤尼克多項(xiàng)式來描述大氣湍流[11]。選用澤尼克多項(xiàng)式對大氣湍流進(jìn)行描述的時(shí)候,每一項(xiàng)均有對應(yīng)的物理意義,如波前像斑的傾斜程度由第一、二項(xiàng)描述；離焦、像散由第三到第五項(xiàng)描述；第八、九項(xiàng)描述了像斑的擴(kuò)展、模糊程度，等等[12]。此處,采用40項(xiàng)澤尼克多項(xiàng)式來描述大氣湍流對激光足印信號產(chǎn)生的影響。

2.3 初始參數(shù)的確定

(1) 激光足印中心初始化

針對地面探測系統(tǒng)獲取的采樣數(shù)據(jù),對數(shù)據(jù)進(jìn)行相應(yīng)的處理。首先,尋找數(shù)據(jù)中的峰值區(qū)域,并將峰值區(qū)域?qū)?yīng)位置標(biāo)記為1,其余地方標(biāo)記為0。由于光斑半徑比較大,采樣數(shù)據(jù)為多個(gè)能級,峰值區(qū)域可能并不是一個(gè)單一的采樣點(diǎn),而是一個(gè)含有多個(gè)采樣點(diǎn)連通域。統(tǒng)計(jì)連通域的個(gè)數(shù)m,統(tǒng)計(jì)每個(gè)峰值區(qū)域的能量之和,將峰值區(qū)域能量最大的區(qū)域作為足印中心所在的區(qū)域。

足印中心初始化包括中心坐標(biāo)初始化和強(qiáng)度信息初始化。在上述所求峰值區(qū)域內(nèi)隨機(jī)選取一個(gè)采樣點(diǎn)的坐標(biāo)(as,bs)作為足印的初始中心,將其能量值Es作為中心位置的能量初始值。

(2)激光足印束腰半徑和方位角初始化

束腰半徑的大小直接反映了激光足印大小,選取合適的束腰半徑會使得算法效率更高。根據(jù)上述確定的激光足印中心位置,在x方向上,將峰值區(qū)域邊界點(diǎn)到足印初始中心(as,bs)的平均距離r1s表示激光足印的x方向上束腰半徑;在y方向上,將峰值區(qū)域邊界點(diǎn)到足印初始中心(as,bs)的平均距離r2s表示激光足印的y方向上束腰半徑。方位角θ是激光足印與選定的坐標(biāo)系的夾角,選定初始化參數(shù)θ=θs=0。

選取合適的初始參數(shù),是為了讓算法更快的達(dá)到最終的迭代結(jié)果,而初始參數(shù)與最終迭代結(jié)果的偏差不影響算法的定位精度。參數(shù)初始化后,式(1)可以表示為:

(2)

2.4 相關(guān)參數(shù)的更新與求解

(1)LM算法求參數(shù)最優(yōu)解與足印中心

Levenberg-Marquardt(LM)算法是一種常用的解決非線性最小二乘問題的方法[13-14],是一種介于高斯牛頓法和梯度下降法之間的優(yōu)化算法,同時(shí),LM算法綜合了兩者之間的優(yōu)勢,在引入阻尼因子之后,能夠調(diào)節(jié)每一次迭代的步長,提高了迭代的效率,使得目標(biāo)函數(shù)以更快的速度收斂。構(gòu)造如下式所示的最小二乘目標(biāo)函數(shù):

(3)

式中,α=(x,y)為坐標(biāo)信息,β=(as,bs,r1s,r2s,θs,Es)為待優(yōu)化輸出向量;Z(α)為測量值;f(α,β)為殘差方程,最終將求解激光足印中心的問題轉(zhuǎn)化為求解最優(yōu)參數(shù)β使得F(α,β)最小。對F(α,β)做二階泰勒展開,有:

(4)

式中,Jf為f(α,β)關(guān)于參數(shù)向量β的雅可比矩陣,Hf為f(α,β)的二階偏導(dǎo)數(shù)組成的關(guān)于參數(shù)向量β的Hessian矩陣?？梢詫⒔频慕Y(jié)果記為A(h)。引入阻尼因子μ后的LM迭代公式可由下式表示:

(5)

式中,μ為阻尼因子,當(dāng)μ≥0時(shí),LM算法的迭代方向是向著F(α,β)減小的方向進(jìn)行,也即是保證了迭代的正確方向；當(dāng)μ很大時(shí),LM算法轉(zhuǎn)化為梯度下降法,對應(yīng)的迭代公式為:

(6)

此時(shí)的算法具有梯度下降法的優(yōu)勢:下降量大、迭代迅速；當(dāng)μ>0且μ很小的時(shí)候,LM算法轉(zhuǎn)化為高斯牛頓法,對應(yīng)的迭代公式可以表示為:

(7)

此時(shí)LM算法具有高斯牛頓法的優(yōu)勢:具有二階收斂性的特點(diǎn)。而阻尼因子μ可以通過系數(shù)ρ來調(diào)節(jié):

(8)

通過系數(shù)ρ的大小對當(dāng)前迭代步驟進(jìn)行取舍,并對μ進(jìn)行調(diào)節(jié)。當(dāng)F(α,β)小于設(shè)定的閾值或者達(dá)到了設(shè)定的迭代次數(shù)即停止迭代。

通過上述的LM算法最終求出參數(shù)β的最優(yōu)解,由此可以得到激光足印中心的坐標(biāo)信息(a*,b*)和其他的相關(guān)參數(shù)。

(2) LM算法局部最優(yōu)解問題

LM算法求目標(biāo)函數(shù)的最小二乘解的過程中,如果初始參數(shù)選取不當(dāng),可能會使得最終求得的參數(shù)為局部最優(yōu)解。針對這種情況,有兩種處理策略:一是隨機(jī)選取多組初始參數(shù),比較最終的目標(biāo)函數(shù)的大小,最小的目標(biāo)函數(shù)對應(yīng)的參數(shù)β即為最優(yōu)解,那么對應(yīng)的(a*,b*)即為激光足印中心。峰值區(qū)域往往對應(yīng)著激光足印能量比較高的區(qū)域,因而第二種處理方法是為每個(gè)峰值區(qū)域選取多組初始參數(shù)進(jìn)行求解,最小的目標(biāo)函數(shù)對應(yīng)的參數(shù)β即為最優(yōu)解,那么對應(yīng)的(a*,b*)即為激光足印中心。

3 數(shù)據(jù)仿真與驗(yàn)證

3.1 大氣湍流仿真

由于地面探測器系統(tǒng)能夠在能級范圍內(nèi)捕獲到星載激光測高儀的激光足印,所以大氣湍流對于激光足印的影響的仿真數(shù)據(jù),只關(guān)注湍流分布。湍流的分布可以通過澤尼克多項(xiàng)式的項(xiàng)數(shù)或者澤尼克多項(xiàng)式的系數(shù)來改變[12]。此處,選取澤尼克項(xiàng)數(shù)為40項(xiàng),通過改變澤尼克多項(xiàng)式的系數(shù)來改變湍流分布,生成40項(xiàng)澤尼克多項(xiàng)式描述的大氣湍流隨機(jī)相位屏,如圖2所示:

圖2 前40項(xiàng)大氣湍流隨機(jī)相位屏Fig.2 The first 40 atmospheric turbulence random phase screens

3.2 激光足印仿真數(shù)據(jù)

根據(jù)地面探測系統(tǒng)捕獲激光足印的方式以及探測器8能級的特點(diǎn),生成大小為30×30的橢圓高斯激光足印仿真數(shù)據(jù),激光足印的強(qiáng)度從中心到邊緣呈現(xiàn)出高斯分布,相對強(qiáng)度從0～8,如圖3所示,激光足印中心坐標(biāo)為(16,16)。

圖3 理想橢圓高斯激光足印Fig.3 Ideal elliptical Gaussian laser footprint

將圖3所示的理想橢圓高斯激光足印信號通過40項(xiàng)澤尼克多項(xiàng)式描述的大氣湍流隨機(jī)相位屏,得到地面探測系統(tǒng)的輸入仿真數(shù)據(jù),如圖4所示,不斷改變澤尼克多項(xiàng)式的系數(shù),得到多組仿真數(shù)據(jù),用GCM、EFM、2DCF以及本文算法,對仿真數(shù)據(jù)進(jìn)行處理,提取激光足印中心,通過與實(shí)際中心坐標(biāo)(16,16)進(jìn)行對比,分析各個(gè)算法基于歐式距離的定位誤差,即各個(gè)算法所求的激光足印中心位置到激光足印實(shí)際的位置之間的歐式距離。

圖4 通過大氣湍流相位屏后的激光足印Fig.4 Laser footprint after passing through the atmospheric turbulence phase screen

3.3 結(jié)果分析

通過用GCM、EFM、2DCF以及本文算法對5000組上述激光足印仿真數(shù)據(jù)的處理,得到了如圖5所示的各算法基于歐式距離的定位誤差的分布情況和表1所示的各算法定位精度的均值和標(biāo)準(zhǔn)差。

圖5 各算法基于歐式距離的定位誤差分布Fig.5 The positioning error distribution of each algorithm based on Euclidean distance

由圖5可知,本文算法80 %以上的數(shù)據(jù)定位誤差在0.2個(gè)布設(shè)間距以下,只有不到5 %的數(shù)據(jù)定位誤差會超過0.4個(gè)布設(shè)間距,GCM方法的定位誤差在整個(gè)區(qū)間上分布比較均勻,EFM和2DCF方法的定位誤差在0.4個(gè)布設(shè)間距以上的占比接近40 %,在0.2個(gè)布設(shè)間距以內(nèi)的數(shù)據(jù)不足25 %。圖5表明了本文算法對于受到大氣湍流影響的激光足印信號中心的提取效果明顯優(yōu)于其他三種算法。

表1 各算法定位精度的均值和標(biāo)準(zhǔn)差Tab.1 Mean and standard deviation of positioning accuracy of each algorithm

由表1可以得到:本文算法對于激光足印中心定位的平均誤差為0.1529個(gè)布設(shè)間距,假設(shè)布設(shè)間距為4 m,那么此時(shí)的定位誤差為0.61 m,而GCM、EFM和2DCF三種算法的定位誤差分別為1.13 m、1.53 m和1.52 m,本文算法的定位精度明顯優(yōu)于其他三種算法,至少優(yōu)于0.5 m。本文算法的定位誤差的標(biāo)準(zhǔn)差為0.1016個(gè)布設(shè)間距,僅為其他三種算法的1/2左右,說明隨著大氣湍流隨機(jī)相位屏的不斷改變,本文算法的定位誤差沒有太大的起伏,而GCM、EFC和2DCF的定位誤差表現(xiàn)出不穩(wěn)定。

4 結(jié) 論

本文介紹了一種基于LM優(yōu)化算法的激光足印中心提取算法,采用了40項(xiàng)澤尼克多項(xiàng)式模擬大氣湍流隨機(jī)相位屏作為激光足印在大氣傳輸過程中的噪聲,利用仿真輸入數(shù)據(jù)研究了本文的算法。實(shí)驗(yàn)結(jié)果表明,本文算法比傳統(tǒng)的中心定位方法有著更高的定位精度并且相對穩(wěn)定,假設(shè)激光足印探測器的布設(shè)間距為4 m,本文算法比其他三種傳統(tǒng)算法的提取精度至少優(yōu)于0.5 m。根據(jù)本文算法的優(yōu)化結(jié)果,可以提取到更為精確的激光足印中心,因而能夠減少星載激光測高儀的測量誤差,也可以給星載激光測高儀的指向角誤差預(yù)留更大的空間或者可以適當(dāng)增大探測器的布設(shè)間距以減少成本。