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

王志文,李 松,罗 敏

(武汉大学电子信息学院,湖北 武汉 430072)

1 引 言

星载激光测高仪是一种可以实现对地表高精度三维观测的天基测量设备,而其精细的真实性校验是保证激光测高仪数据产品在各个领域得到有效应用的前提[1]。2003年,首个对地观测的激光测高仪(Geoscience Laser Altimeter System,GLAS)于美国成功发射,随后,相关研究人员提出了多种真实性校验的方式[2-4]。其中,L.A.Magruder等提出一种独立于其他校验方式的地面探测系统[4],通过地面探测系统捕获星载激光测高仪精确的激光足印位置,来实现对星载激光测高仪的真实性校验,而高精度的激光足印中心提取方法是获取激光足印精确位置的关键。地面探测系统捕获到的激光足印,在传输的过程中,会受到大气湍流的影响,使得最终到达地面的激光足印的能量空间分布发生一定的变化,如偏转、扩展[5-6],因此,对于激光足印中心的提取精度,除了依赖于激光足印本身的高斯能量分布特性,还会受到大气湍流噪声的影响。

传统的激光足印中心定位算法有灰度重心法[7](Gray Centroid Method,GCM)、圆或椭圆拟合法[8](Elliptic Fitting Method,EFM)、基于二维向量空间圆拟合法[9](2D Vector Space Circle Fitting,2DCF)等等,这些方法要么是基于整个激光足印的能量分布,如灰度重心法,要么是基于激光足印的几何形状,如圆拟合法。基于激光足印的能量分布中心提取方法对于激光足印的初始能量分布较为严格,当噪声比较小,提取的精度较好,但不适用于噪声相对比较大的情形且没有充分利用激光足印的空间分布信息；基于几何形状的提取方法,考虑了激光足印的几何分布,但是当边缘或整个空间分布不是规则的圆或者椭圆,那么就会产生比较大的误差。

本文提出一种基于LM优化算法的激光足印中心提取方法,对于激光足印数据的拟合过程比较符合激光足印捕获数据的数据格式以及激光足印的光强分布特点,不会依赖于每个能级边缘形状的规则程度,能从全局上找到一个最优的激光足印中心,提取精度高,并且抗干扰能力较强。

2 基于LM优化算法的激光足印中心提取原理

星载激光测高仪发射出的激光经过大气传输后,被地面探测系统捕获,然后通过相关算法提取可以得到激光足印的中心,从而实现星载激光测高仪的在轨校验。激光足印轮廓在大多数情况下服从椭圆高斯分布[10],可用下式表示:

(1)

其中,α=(x,y)为坐标信息,β=(a,b,r1,r2,θ,E)为参数列表;a为x方向上激光足印中心坐标；b为y方向上激光足印中心坐标；r1为x方向上椭圆高斯激光足印的束腰半径；r2为y方向上椭圆高斯激光足印的束腰半径,当r1=r2时,足印的长轴和短轴相等,此时激光足印呈圆高斯分布；θ为对应的椭圆高斯激光足印的方位角,取值范围为:[0,π/4]；E表示椭圆高斯激光足印的强度信息。N(α)为相应的大气湍流噪声。

根据采集到的椭圆高斯激光足印信号的特点,采用一种基于LM算法的激光足印中心提取方法,首先,构建最小二乘目标函数,确定拟合的椭圆高斯激光足印初始参数:激光足印中心的位置信息、强度信息、长短轴的大小和方位角；然后计算最小二乘目标函数的值,使用LM算法更新参数,使每次更新后的目标函数的值比上一次迭代的值小,否则,重新更新参数,即按照目标函数减小的方向进行迭代；当目标函数的值小于预先设置的阈值,或者达到迭代次数时,即停止迭代,此时拟合函数与输入数据的残差最小,而对应的参数(a,b)就是待求的足印中心。整个处理的过程如图1所示。

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

2.1 地面探测系统及其相应数据格式

针对我国自主研制的搭载于资源三号02星的激光测高仪的在轨校验实验,文献[1]中提出了相应的地面足印探测系统,其中,用于捕获激光足印的能量探测器有8个量化能级。假设布设方式为30×30方阵,布设间距为4 m,激光足印半径30 m,在激光足印过境的位置固定好探测器后会获取到激光足印的8能级30×30的采样数据,并且包含足印的位置信息。

2.2 大气湍流噪声

星载激光测高仪所发射的激光探测信号在大气中传播的过程中,主要受到了大气湍流的影响而发生畸变。Kolomogorov等建立了大气湍流的理论模型后,由于泽尼克多项式低频部分占主要部分,高频成分较少,这与大气湍流的理论模型表现一致,因而Noll等采用泽尼克多项式来描述大气湍流[11]。选用泽尼克多项式对大气湍流进行描述的时候,每一项均有对应的物理意义,如波前像斑的倾斜程度由第一、二项描述；离焦、像散由第三到第五项描述；第八、九项描述了像斑的扩展、模糊程度，等等[12]。此处,采用40项泽尼克多项式来描述大气湍流对激光足印信号产生的影响。

2.3 初始参数的确定

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

针对地面探测系统获取的采样数据,对数据进行相应的处理。首先,寻找数据中的峰值区域,并将峰值区域对应位置标记为1,其余地方标记为0。由于光斑半径比较大,采样数据为多个能级,峰值区域可能并不是一个单一的采样点,而是一个含有多个采样点连通域。统计连通域的个数m,统计每个峰值区域的能量之和,将峰值区域能量最大的区域作为足印中心所在的区域。

足印中心初始化包括中心坐标初始化和强度信息初始化。在上述所求峰值区域内随机选取一个采样点的坐标(as,bs)作为足印的初始中心,将其能量值Es作为中心位置的能量初始值。

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

束腰半径的大小直接反映了激光足印大小,选取合适的束腰半径会使得算法效率更高。根据上述确定的激光足印中心位置,在x方向上,将峰值区域边界点到足印初始中心(as,bs)的平均距离r1s表示激光足印的x方向上束腰半径;在y方向上,将峰值区域边界点到足印初始中心(as,bs)的平均距离r2s表示激光足印的y方向上束腰半径。方位角θ是激光足印与选定的坐标系的夹角,选定初始化参数θ=θs=0。

选取合适的初始参数,是为了让算法更快的达到最终的迭代结果,而初始参数与最终迭代结果的偏差不影响算法的定位精度。参数初始化后,式(1)可以表示为:

(2)

2.4 相关参数的更新与求解

(1)LM算法求参数最优解与足印中心

Levenberg-Marquardt(LM)算法是一种常用的解决非线性最小二乘问题的方法[13-14],是一种介于高斯牛顿法和梯度下降法之间的优化算法,同时,LM算法综合了两者之间的优势,在引入阻尼因子之后,能够调节每一次迭代的步长,提高了迭代的效率,使得目标函数以更快的速度收敛。构造如下式所示的最小二乘目标函数:

(3)

式中,α=(x,y)为坐标信息,β=(as,bs,r1s,r2s,θs,Es)为待优化输出向量;Z(α)为测量值;f(α,β)为残差方程,最终将求解激光足印中心的问题转化为求解最优参数β使得F(α,β)最小。对F(α,β)做二阶泰勒展开,有:

(4)

式中,Jf为f(α,β)关于参数向量β的雅可比矩阵,Hf为f(α,β)的二阶偏导数组成的关于参数向量β的Hessian矩阵。可以将近似的结果记为A(h)。引入阻尼因子μ后的LM迭代公式可由下式表示:

(5)

式中,μ为阻尼因子,当μ≥0时,LM算法的迭代方向是向着F(α,β)减小的方向进行,也即是保证了迭代的正确方向；当μ很大时,LM算法转化为梯度下降法,对应的迭代公式为:

(6)

此时的算法具有梯度下降法的优势:下降量大、迭代迅速；当μ>0且μ很小的时候,LM算法转化为高斯牛顿法,对应的迭代公式可以表示为:

(7)

此时LM算法具有高斯牛顿法的优势:具有二阶收敛性的特点。而阻尼因子μ可以通过系数ρ来调节:

(8)

通过系数ρ的大小对当前迭代步骤进行取舍,并对μ进行调节。当F(α,β)小于设定的阈值或者达到了设定的迭代次数即停止迭代。

通过上述的LM算法最终求出参数β的最优解,由此可以得到激光足印中心的坐标信息(a*,b*)和其他的相关参数。

(2) LM算法局部最优解问题

LM算法求目标函数的最小二乘解的过程中,如果初始参数选取不当,可能会使得最终求得的参数为局部最优解。针对这种情况,有两种处理策略:一是随机选取多组初始参数,比较最终的目标函数的大小,最小的目标函数对应的参数β即为最优解,那么对应的(a*,b*)即为激光足印中心。峰值区域往往对应着激光足印能量比较高的区域,因而第二种处理方法是为每个峰值区域选取多组初始参数进行求解,最小的目标函数对应的参数β即为最优解,那么对应的(a*,b*)即为激光足印中心。

3 数据仿真与验证

3.1 大气湍流仿真

由于地面探测器系统能够在能级范围内捕获到星载激光测高仪的激光足印,所以大气湍流对于激光足印的影响的仿真数据,只关注湍流分布。湍流的分布可以通过泽尼克多项式的项数或者泽尼克多项式的系数来改变[12]。此处,选取泽尼克项数为40项,通过改变泽尼克多项式的系数来改变湍流分布,生成40项泽尼克多项式描述的大气湍流随机相位屏,如图2所示:

图2 前40项大气湍流随机相位屏Fig.2 The first 40 atmospheric turbulence random phase screens

3.2 激光足印仿真数据

根据地面探测系统捕获激光足印的方式以及探测器8能级的特点,生成大小为30×30的椭圆高斯激光足印仿真数据,激光足印的强度从中心到边缘呈现出高斯分布,相对强度从0～8,如图3所示,激光足印中心坐标为(16,16)。

图3 理想椭圆高斯激光足印Fig.3 Ideal elliptical Gaussian laser footprint

将图3所示的理想椭圆高斯激光足印信号通过40项泽尼克多项式描述的大气湍流随机相位屏,得到地面探测系统的输入仿真数据,如图4所示,不断改变泽尼克多项式的系数,得到多组仿真数据,用GCM、EFM、2DCF以及本文算法,对仿真数据进行处理,提取激光足印中心,通过与实际中心坐标(16,16)进行对比,分析各个算法基于欧式距离的定位误差,即各个算法所求的激光足印中心位置到激光足印实际的位置之间的欧式距离。

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

3.3 结果分析

通过用GCM、EFM、2DCF以及本文算法对5000组上述激光足印仿真数据的处理,得到了如图5所示的各算法基于欧式距离的定位误差的分布情况和表1所示的各算法定位精度的均值和标准差。

图5 各算法基于欧式距离的定位误差分布Fig.5 The positioning error distribution of each algorithm based on Euclidean distance

由图5可知,本文算法80 %以上的数据定位误差在0.2个布设间距以下,只有不到5 %的数据定位误差会超过0.4个布设间距,GCM方法的定位误差在整个区间上分布比较均匀,EFM和2DCF方法的定位误差在0.4个布设间距以上的占比接近40 %,在0.2个布设间距以内的数据不足25 %。图5表明了本文算法对于受到大气湍流影响的激光足印信号中心的提取效果明显优于其他三种算法。

表1 各算法定位精度的均值和标准差Tab.1 Mean and standard deviation of positioning accuracy of each algorithm

由表1可以得到:本文算法对于激光足印中心定位的平均误差为0.1529个布设间距,假设布设间距为4 m,那么此时的定位误差为0.61 m,而GCM、EFM和2DCF三种算法的定位误差分别为1.13 m、1.53 m和1.52 m,本文算法的定位精度明显优于其他三种算法,至少优于0.5 m。本文算法的定位误差的标准差为0.1016个布设间距,仅为其他三种算法的1/2左右,说明随着大气湍流随机相位屏的不断改变,本文算法的定位误差没有太大的起伏,而GCM、EFC和2DCF的定位误差表现出不稳定。

4 结 论

本文介绍了一种基于LM优化算法的激光足印中心提取算法,采用了40项泽尼克多项式模拟大气湍流随机相位屏作为激光足印在大气传输过程中的噪声,利用仿真输入数据研究了本文的算法。实验结果表明,本文算法比传统的中心定位方法有着更高的定位精度并且相对稳定,假设激光足印探测器的布设间距为4 m,本文算法比其他三种传统算法的提取精度至少优于0.5 m。根据本文算法的优化结果,可以提取到更为精确的激光足印中心,因而能够减少星载激光测高仪的测量误差,也可以给星载激光测高仪的指向角误差预留更大的空间或者可以适当增大探测器的布设间距以减少成本。