CN103440659B - 基于星图匹配的星空图像畸变检测与估计方法 - Google Patents

Abstract

本发明公开了一种基于星图匹配的星空图像畸变检测与估计方法，用于解决现有检测方法在图像遭受畸变影响的情况下检测精度低的技术问题。技术方案是利用已知的星表、成像系统参数和成像系统指向，对图像中星点与星表中星点进行匹配，进而对图像中可能存在的相机的旋转、平移和畸变进行检测与估计。由于考虑了畸变对星点匹配的影响、指向不精确和相机镜头平移旋转造成的图像中星点的偏移旋转、平移旋转与畸变综合作用造成的估计困难，利用基于畸变模型改进的最长LCS星点配准方法对星空图像与星表中星点进行配准，对匹配点对集合中可能存在的错误匹配外点进行剔除并估计相对的旋转平移关系，进而估计星图与星标之间的畸变。提高了检测精度。

Description

基于星图匹配的星空图像畸变检测与估计方法

技术领域

本发明涉及一种星空图像畸变检测与估计方法，特别是涉及一种基于星图匹配的星空图像畸变检测与估计方法。

背景技术

文献“基于星图匹配的小视场空间图像畸变校正[J].北京航空航天大学学报,2008,09:1016-1019.”提出了一种基于星图匹配的小视场空间图像畸变校正方法。该方法基于Hausdorff距离进行多帧空间图像中星点的匹配，并在通过误差项引入针对指向偏差的补偿。该方法通过匹配点之间的偏差对全局进行校正。但是由于该方法主要针对小视场图像，视场中星点较少，因此使用Hausdorff距离进行适用的星点匹配，并且简单地使用星点之间的偏差对全局进行建模。但是对于一般情况下的星空图像，视场内星点较多，尤其是在图像遭受畸变影响的情况下，Hausdorff距离难以进行准确的匹配，并且仅仅使用误差的插值难以对全局图像畸变进行准确描述。因而该方法无法直接应用于对含有几何畸变和指向偏差的一般星空图像畸变的检测和估计。

发明内容

为了克服现有检测方法在图像遭受畸变影响的情况下检测精度低的不足，本发明提供一种基于星图匹配的星空图像畸变检测与估计方法。该方法利用已知的星表、成像系统参数和成像系统指向，对图像中星点与星表中星点进行匹配，进而对图像中可能存在的相机的旋转、平移和畸变进行检测与估计。该方法充分考虑了畸变对星点匹配的影响、指向不精确和相机镜头平移旋转造成的图像中星点的偏移旋转、平移旋转与畸变综合作用造成的估计困难，利用基于畸变模型改进的最长公共子序列(LCS)星点配准方法对星空图像与星表中星点进行配准，对匹配点对集合中可能存在的错误匹配外点进行剔除并估计相对的旋转平移关系，在考虑旋转平移的情况下估计星图与星标之间的畸变。在以上过程中完成对星图中旋转平移和畸变的检测，同时通过联合迭代求解对旋转平移和畸变模型进行估计。

本发明解决其技术问题所采用的技术方案是：一种基于星图匹配的星空图像畸变检测与估计方法，其特点是包括以下步骤：

步骤一、对星空图像中的星点进行提取，并计算星点在图像上坐标。该过程包括，利用阈值对星空图像进行二值化分割，即根据像素值是否大于阈值而进行分割；对二值图像进行联通区域提取，得到星点光斑区域集合Ω_n；利用提取到的光斑区域，从原始图像中计算星点的亚像素坐标(x_n,y_n)，对于每个光斑，计算方法为

$x_n=\frac{\underset{(x,y)\∈{\mathrm{\Ω}}_n}{\Σ}I(x,y)\·x}{\underset{(x,y)\∈{\mathrm{\Ω}}_n}{\Σ}I(x,y)},y_n=\frac{\underset{(x,y)\∈{\mathrm{\Ω}}_n}{\Σ}I(x,y)\·y}{\underset{(x,y)\∈{\mathrm{\Ω}}_n}{\Σ}I(x,y)}---\left(1\right)$

式中，(x,y)∈Ω_n表示属于光斑Ω_n的像素坐标，I(x,y)为相应的像素值。对图像中的星点进行计算得到星点集合I为集合中星点数目。

步骤二、根据拍摄时的输入指向、视场大小、相机参数和拍摄时间，通过星空坐标转换得到对应天区的星体在CCD成像平面上的理想成像位置作为参考星点集合其中，J为集合中星点数目。

步骤三、像畸变检测与估计框架为在迭代中完成星图匹配、旋转平移估计和畸变估计，具体过程包括：

(a).使用基于LCS的星图配准方法根据星图畸变模型进行改进，对和进行匹配，若匹配失败，说明图像与输入的指向之间存在较大偏差或者图像畸变过于严重，完成星图畸变检测；若匹配成功，得到可能存在匹配错误的点对集合，进行进一步的检测与估计。

(b).使用一致随机采样对匹配外点进行剔除，并同时估计图像星点与理想星表星点之间的近似旋转平移关系模型，若剔除外点的匹配点对集合中点对个数小于求解畸变模型所需点对个数，则说明图像与输入的指向之间存在较大偏差或者图像畸变过于严重，完成星图畸变检测；否则，利用剔除外点后的点对集合与估计得到的旋转平移模型进行进一步的畸变检测与估计。

(c)利用无外点的、消除空间旋转平移差异的匹配集合估计图像中的几何畸变。

(d)迭代步骤(a)-(c)估计得到和之间的旋转平移模型和畸变模型，通过估计得的旋转平移和畸变模型得到图像的畸变程度的量化评价。

其中，步骤(a)-(c)的迭代中，包括基于畸变模型的LCS畸变星图配准、基于RANSAC的外点剔除与旋转平移模型估计、迭代估计几何畸变。

第一部分，基于畸变模型的LCS畸变星图配准。假设从星空图像中提取得到的星点集合为而通过给定参数得到的市场内星点在CCD平面的投影的星点集合为对于得到一个大小为I×I的邻接矩阵，其中每个元素表示中两点之间在图像平面的欧氏距离，对该矩阵各行进行升序排序得到特征矩阵：

${\mathrm{FMat}}_D=\left[\begin{array}{cccc}0& d_{1,2}& ...& d_{1,I}\\ 0& d_{2,2}& ...& d_{2,I}\\ ...& ...& ...& ...\\ 0& d_{I,2}& ...& d_{I,I}\end{array}\right]---\left(2\right)$

并且同时将FMat_D中元素在各行中的原始序号记录在矩阵中得到：

${\mathrm{IMat}}_D=\left[\begin{array}{cccc}1& p_{1,2}& ...& p_{1,I}\\ 2& p_{2,2}& ...& p_{1,1}\\ ...& ...& ...& ...\\ I& p_{I,2}& ...& p_{1,I}\end{array}\right]---\left(3\right)$

因此，元素p_i,j表示FMat₀中的元素d_i,j元素为中D_i与之间的距离。同样，对于星点集合得到相应的矩阵FMat_R和IMat_R。其中FMat_D和FMat_R中每行为相应集合中各个星点的特征。

在匹配过程中为比较星点D_x与R_y之间的相似性，需要计算序列<0,d_x,i,...,d_x,I>_D和序列<0,d_y,j,...,d_x,J>_R之间的公共子序列的长度，其中x和y为星点序号。基于一般动态规划求解LCS的算法框架，设定其中的公共子序列更新公式为：

$c\&lsqb;i_1,j_1\&rsqb;=\left\{\begin{array}{cc}0,& i_1=0or{j}_1=0\\ c\&lsqb;i_1-1,j_1-1\&rsqb;+1,& i_1,j_1>0and\left|\right|d_{x,i_1}-d_{y,j_1}\left|\right|\&le;\mathrm{\&epsiv;}f(d_{x,i_1}-d_{y,j_1})\\ \max (c\&lsqb;i_1,j_1-1\&rsqb;,c\&lsqb;i_1-1,j_1\&rsqb;),& i_1,j_1>0and\left|\right|d_{x,i_1}-d_{y,j_1}\left|\right|>\mathrm{\&epsiv;}f(d_{x,i_1}-d_{y,j_1})\end{array}\right.---\left(4\right)$

式中，i₁和j₁为动态规划算法表中的序号，ε为距离阈值，为[0,1]之间的函数，成立时，表示在畸变情况下而函数f(·)为：

$f(d_{x,i_1},d_{y,j_1})=\frac{1}{L}\exp \left\{\frac{\left|\right|C_x-C\left|\right|+\left|\right|D_{{\mathrm{IMat}}_D\&lsqb;x,i_1\&rsqb;}-C\left|\right|+\left|\right|R_y-C\left|\right|+\left|\right|R_{{\mathrm{IMat}}_R\&lsqb;y,j_1\&rsqb;}-C\left|\right|}{4{\mathrm{\&sigma;}}^2}\right\}---\left(5\right)$

式中，C为图像畸变中心，在迭代中更新，||·||表示欧式距离，σ²为方差，L为归一化因子。

得到与中任两个星点之间的LCS相似性度量。对于中星点，依次在中搜索LCS最近邻，如果与最近邻之间的关系满足阈值要求，则作为与相应D_i匹配的R_i，从而得到点对，否则舍弃对相应D_i的匹配。

经过星图匹配，得到匹配点对的集合其中P_m＝<R_m,D_m>。中包含足够多的正确配对点，但仍包含匹配错误点对。

第二部分，基于RANSAC的外点剔除与旋转平移模型估计。中依然存在匹配错误点对，但其中正确匹配的点对之间存在稳定的旋转平移关系，而不符合该模型的匹配错误点对为外点。使用RANSAC估计其旋转平移关系，同时剔除其中匹配错误的外点。

和之间的旋转平移关系描述为二维空间中的旋转平移关系。由到的转换关系T(·)描述为：

$\left[\begin{array}{c}{D}_{xm}\\ D_{ym}\end{array}\right]=T\left(\left[\begin{array}{c}{R}_{xm}\\ R_{ym}\end{array}\right]\right)=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\&CenterDot;\left[\begin{array}{c}{R}_{xm}\\ R_{ym}\end{array}\right]+\left[\begin{array}{c}e\\ f\end{array}\right]---\left(6\right)$

其中，D_xm表示第m个点的x方向坐标，D_ym、R_xm和R_ym同理。

在估计模型T(·)过程中，对于给定的N对匹配点对，由于存在测量误差和畸变未知引起的误差，点对中点坐标存在噪声，使用足够多的点对求取最优估计。根据式(6)建立方程

Ax＝b(7)

$A=\left[\begin{array}{cccccc}{R}_{x1}& R_{y1}& 1& 0& 0& 0\\ 0& 0& 0& R_{x1}& R_{y1}& 1\\ ...& ...& ...& ...& ...& ...\\ R_{xM}& R_{xM}& 1& 0& 0& 0\\ 0& 0& 0& R_{xM}& R_{xM}& 1\end{array}\right]---\left(8\right)$

x＝[abecdf]^T(9)

b＝[D_x1D_y1...D_xMD_yM]^T(10)

方程(7)有封闭形式解，因此旋转平移模型参数求解为

x＝(A^TA)^-1A^TB(11)

RANSAC算法通过在每次迭代中对包含外点的点对集合随机采样估计模型参数，得到符合该模型的候选点对集，并计算误差。最终从所有迭代中取出拟合误差最小的模型作为估计结果，相应的候选点对集合为剔除外电的集合。在所有迭代中选取最佳。

RANSAC的算法输入为：含有外点的点对集合每次采样最少点数目Q，迭代次数K，判断一个点对是否符合当前模型的误差阈值T，判断估计得到的模型是否足够好的候选集合大小D。算法流程为：迭代K次，对于第k次迭代，1≤k≤K，

(a)在中随机选取Q个点对得到候选点集{P_q}_con，此时，1≤q≤Q；

(b)使用式(11)估计参数x_k，得到投影模型T_k(·)；

(c)对于除{P_q}_con中点对外的其他点对<R_m,D_m>：计算投影误差e_k,m＝||D_n-T_k(R_n)||，如果e_k,m≤T，则将<R_m,D_m>加入候选点对集合{P_q}_con，每次更新使{P_q}_con变大，使q的取值范围上限变大；

(d)如果{P_q}_con元素数目大于D，则使用{P_q}_con估计式(11)中参数x_con，得到投影模型T_con(·)；

(e)计算投影误差 $e_k=\underset{<R_q,D_q>\&Element;{\left\{P_q\right\}}_{con}}{\&Sigma;}\left|\right|D_q-T_{con}\left(R_q\right)\left|\right|;$

对于每次迭代，求取e_k最小的T_con(·)作为输出模型参数，同时相应的{P_n}_con作为剔除外点的点集。

使用RANSAC剔除中外点得到集合估计得到相应旋转平移参数x，即和之间的转换包含的旋转平移函数T(·)。

第三部分，迭代估计几何畸变。在得到不包含外点的匹配点对集合和旋转平移函数T(·)之后，中的经过T(·)变换有与相应的组合得到点对集合因此，中与之间仅包含集合畸变的变换，被用来估计几何畸变。使用径向畸变来描述中的几何畸变。由到的几何畸变为

$\left[\begin{array}{c}{D}_{xn}\\ D_{yn}\end{array}\right]=G\left(\left[\begin{array}{c}{R}_{xn}^{\&prime;}\\ R_{yn}^{\&prime;}\end{array}\right]\right)=\left[\begin{array}{c}{G}_x\left(R_{xn}^{\&prime;}\right)\\ G_y\left(R_{yn}^{\&prime;}\right)\end{array}\right]=L\left(r\right)\left[\begin{array}{c}{R}_{xn}^{\&prime;}-C_x\\ R_{yn}^{\&prime;}-C_y\end{array}\right]+\left[\begin{array}{c}{C}_x\\ C_y\end{array}\right]---\left(12\right)$

其中，D_xn表示第n个点的x方向坐标，D_yn、R′_xn和R′_yn同理，C_x为畸变中心x方向坐标，C_y为畸变中心y方向坐标，G(·)表示畸变函数，表示到畸变中心C的距离，L(r)为失真函数，具有性质为L(0)＝1，使用其在0处的泰勒展开表示，仅使用偶次项并舍弃高次项，因此L(r)表示为。

L(r)＝1+k₁r²+k₂r⁴(13)

模型仅由畸变中心C＝[C_x,C_y]^T和畸变程度k₁、k₂确定，为方便描述，将参数表示为P＝[k₁,k₂,C_x,C_y]。

与y方向无差别地，x方向的畸变模型G_x(·)表示为：

$\begin{array}{c}{D}_{xn}=G_x(R_{xn}^{\&prime;};P)=(R_{xn}^{\&prime;}-C_x)+k_1{r}^2(R_{xn}^{\&prime;}-C_x)+k_2{r}^4(R_{xn}^{\&prime;}-C_x)\\ r=\sqrt{{(R_{xn}^{\&prime;}-C_x)}^2+{(R_{yn}^{\&prime;}-C_y)}^2}\end{array}---\left(14\right)$

其中，n为点对序号。同理，G_y(·)也由此表示。

使用牛顿迭代对该非线性模型进行求解。给定初始值参数初始值P₀，在迭代中求得最优解。迭代中参数更新方式为：

P_k+1＝P_k+Δ_k(15)

其中，k为迭代次数序号，且只有当k出现在下标与[]中时表示迭代次数，Δ_k为方程J_kΔ_k＝ε_k的最小二乘解，为取值在P_k上的雅各比矩阵，而ε_k＝D-G(R；P_k)为当前估计的误差。在每次迭代中对于每个值代入J之后均有：

$J_{k,n}(R_n;P_k)=\left[\begin{array}{c}\frac{\&part;G_x(R_{nx};P)}{\&part;P}\&lsqb;k,n\&rsqb;\\ \frac{\&part;G_y(R_{ny};P)}{\&part;P}\&lsqb;k,n\&rsqb;\end{array}\right]---\left(16\right)$

其中

$\frac{\&part;G_x(R_x;P)}{\&part;P}\&lsqb;k,n\&rsqb;=\&lsqb;\frac{\&part;G_x(R_x;k_1)}{\&part;k_1}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;k_2)}{\&part;k_2}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;C_x)}{\&part;C_x}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;C_y)}{\&part;C_y}\&lsqb;k,n\&rsqb;\&rsqb;---\left(17\right)$

而类似地表示。针对G(·)中的G_x(·)，具体地有G_x(·)针对不同参数求偏导：

$\begin{array}{c}\frac{\&part;G_x(R_x;P)}{\&part;k_1}\&lsqb;k,n\&rsqb;={r_{k,n}}^2{R}_{xn}^{\&prime;}\\ \frac{\&part;G_x(R_x;P)}{\&part;k_2}\&lsqb;k,n\&rsqb;={r_{k,n}}^4{R}_{xn}^{\&prime;}\\ \frac{\&part;G_x(R_x;P)}{\&part;C_x}\&lsqb;k,n\&rsqb;=-k_{1,k}{r_{k,n}}^2-2k_{1,k}{(R_{xn}-C_{x,k})}^2-k_{2,k}{r_{k,n}}^4-4k_{2,k}{r_{n,k}}^2{(R_{nx}-C_{x,k})}^2\\ \frac{\&part;G_x(R_x;P)}{\&part;C_y}\&lsqb;k,n\&rsqb;=-2k_{1,k}(R_{xn}-C_{x,k}){(R_{yn}-C_{y,k})}^2-4k_{2,k}{r}^2(R_{nx}-C_{x,k})(R_{ny}-C_{y,k})\end{array}---\left(18\right)$

其中，下标中的k均表示第k次迭代中相应的值，下标中的n表示相应的量为代入第n点对之后的相应值，对于G_y(·)得到相同的求偏导结果。

基于以上，在每次迭代中求解方程J_iΔ_i＝ε_i，在利用到所有点对的情况下有：

J_i＝[J_i(R₁；P_i),J_i(R₂；P_i),...,J_i(R_N；P_i)]^T(19)

ε_i＝[(D₁-G(R₁；P_i))^T,(D₂-G(R₂；P_i))^T,...,(D_N-G(R_N；P_i))^T]^T(20)

其中，(D₁-G(R₁；P_i))＝[D_x1-G_x(R_x1；P_i),D_y1-G_y(R_y1；P_i)]^T。于是，解Δ_i为

Δ_i＝J⁺ε_i＝(J^TJ)^-1J^Tε_i(21)

其中，J⁺＝(J^TJ)^-1J^T。

根据以上，迭代优化求解畸变模型的具体算法表示为，算法输入为输入：匹配点对集合模型G(·)，参数初始值P₀，迭代最大次数K，收敛误差ε_con。算法迭代过程中，如果迭代次数小于最大迭代次数k<K，则进行如下迭代，对于第k次迭代：

(a)利用式(20)代入当前参数值P_k，计算ε_k，如果||ε_k||<ε_con，停止迭代，输出结果，否则继续下一步；

(b)将和P_i代入式(19)，计算J_k；

(c)利用式(21)计算得到Δ_k；

(d)更新参数P_k+1＝P_k+Δ_k；

停止迭代，如果停止迭代时k的值为K，否则即为满足步骤(a)中迭代停止条件时的迭代步骤中的输出结果。

通过综合以上三部分，利用联合迭代框架检测和估计旋转平移和畸变模型。整体迭代算法框架输入为：星表参考星点集合星空图像星点集合星图匹配、旋转平移估计、畸变估计所需参数，最大迭代次数K，迭代停止误差ε_con。迭代算法首先初始化，令 为算法中用于传递每一步畸变估计结果的点集，令k从1到K进行循环迭代：

(a)使用改进的LCS匹配算法对和进行点对匹配得到匹配点对集合如果为空集，则说明和之间变换程度过大，检测到成像系统指向偏差过大或者图像畸变程度过大，退出迭代，输出结果，终止程序；

(b)将作为输入使用RANSAC方法估计旋转平移模型T_k(·)，同时得到剔除匹配外点的点对集合若估计失败或过小不足以进行后续运算，则说明和之间畸变程度过大，检测到成像系统指向偏差过大或者图像畸变程度过大，退出迭代，输出结果，终止程序；

(c)利用中点对的序号关系，得到对应的点对集合使用T_k(·)对{R_n}进行投影得到{P′_n}＝{<R′_n,D_n>}＝{<T(R_n),D_n>}；

(d)将{P′_n}代入，迭代估计得到畸变模型G_k(·)；如果估计失败，说明图像畸变程度过大，，退出迭代，输出结果，终止程序；

(e)计算误差 ${\mathrm{\&epsiv;}}_k=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\right|G_k\left(R_n^{\&prime;}\right)-D_n\left|\right|,$ 如果ε_k<ε_con，则 $\hat{T}(\&CenterDot;)=T_k(\&CenterDot;),\hat{G}(\&CenterDot;)=G_k(\&CenterDot;),$ 停止循环，输出结果，否则，进行下一步；

(f)更新 ${\left\{R_j^{*}\right\}}_{j=1}^J,$ 有 ${\left\{R_j^{*}\right\}}_{j=1}^J=T_k^{-1}\left(G_k\right(T_k\left({\left\{R_j\right\}}_{j=1}^J\right)\left)\right),$ 进入下一次循环。

如果最终循环进行了K次，则输出结果根据输出的和对图像与星表相比的旋转平移和畸变作出量化评价。

本发明的有益效果是：该方法利用已知的星表、成像系统参数和成像系统指向，对图像中星点与星表中星点进行匹配，进而对图像中可能存在的相机的旋转、平移和畸变进行检测与估计。该方法充分考虑了畸变对星点匹配的影响、指向不精确和相机镜头平移旋转造成的图像中星点的偏移旋转、平移旋转与畸变综合作用造成的估计困难，利用基于畸变模型改进的最长公共子序列(LCS)星点配准方法对星空图像与星表中星点进行配准，对匹配点对集合中可能存在的错误匹配外点进行剔除并估计相对的旋转平移关系，在考虑旋转平移的情况下估计星图与星标之间的畸变。在以上过程中完成对星图中旋转平移和畸变的检测，同时通过联合迭代求解对旋转平移和畸变模型进行估计。

下面结合具体实施方式对本发明作详细说明。

具体实施方式

本发明基于星图匹配的星空图像畸变检测与估计方法具体步骤如下：

1.图像星点坐标提取。

对星空图像中的星点进行提取，并计算星点在图像上坐标。该过程包括，利用阈值对星空图像进行二值化分割，即根据像素值是否大于阈值ε_bw而进行分割，图像(x,y)处像素点灰度值为I(x,y)，图像中像素灰度值值I(x,y)≥ε_bw的分割为星点光斑区域，而I(x,y)<ε_bw的像素点分割为背景区域。通过对星图数据进行鲁棒统计得到分割阈值

ε_bw＝μ+δσ

其中μ为图像灰度中值，σ为图像平均绝对偏差，δ为阈值参数，本实现中取为δ＝6。

在二值化图像中对二值图像进行联通区域提取，得到星点光斑区域集合{C_n}；利用提取到的光斑区域，从原始图像中计算星点的亚像素坐标(x_n,y_n)，对于每个光斑，计算方法为

$x_n=\frac{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)\&CenterDot;x}{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)},y_n=\frac{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)\&CenterDot;y}{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)}---\left(1\right)$

其中，(x,y)∈Ω_n表示属于光斑Ω_n的像素坐标，I(x,y)为相应的像素值。对图像中的星点进行计算得到星点集合其中，I为集合中星点数目。

2.星表星点提取。

由于已知图像拍摄时输入的摄像机指向、市场大小和拍摄时间，可以在星图中进行搜索得到相应天区应该出现的星体。在对星表中星体坐标进行修正后，利用摄像机参数，可以得到星表中星点在CCD平面上的理想位置，即参考星点集合其中，J为集合中星点数目。

3.迭代估计相对于的旋转平移变换和畸变。

将到的变换中的旋转平移和畸变成分进行分解，进行交替迭代估计。通过固定一部分来估计另一部分。在估计的迭代过程中同时对畸变进行检测。迭代过程可以描述为：

a)使用基于LCS的星图配准方法根据星图畸变模型进行了改进，对中和进行匹配，若匹配失败，说明图像与输入的指向之间存在较大偏差或者图像畸变过于严重，完成星图畸变检测；若匹配成功，得到可能存在匹配错误的点对集合，可以进行进一步细致的检测与估计。

b)使用一致随机采样(RANSAC)的方法匹配对匹配外点进行剔除，并同时估计星空图像与理想星表图像之间的近似旋转平移关系模型，若剔除外点的匹配点对集合中点对个数小于求解畸变模型所需点对个数，则说明图像与输入的指向之间存在较大偏差或者图像畸变过于严重，完成星图畸变检测；否则，利用剔除外点后的点对集合与估计得到的旋转平移模型进行进一步的畸变检测与估计。

c)利用无外点的匹配集合估计图像中的几何畸变。

d)迭代(1)-(3)估计得到和之间的旋转平移模型和畸变模型，通过估计得的旋转平移和畸变模型得到图像的畸变程度的量化评价。

其中a)-b)的迭代中，包括基于畸变模型的LCS畸变星图配准、基于RANSAC的外点剔除与旋转平移模型估计、迭代估计几何畸变。以下对迭代中的各个模块和迭代的具体实现进行详细描述。

(1)基于畸变模型的LCS畸变星图配准。

假设从星空图像中提取得到的星点集合为而通过给定参数得到的市场内星点在CCD平面的投影的星点集合为对于我们可以得到一个大小为I×I的邻接矩阵，其中每个元素表示中两点之间在图像平面的欧氏距离，对该矩阵各行进行升序排序可以得到特征矩阵：

${\mathrm{FMat}}_D=\left[\begin{array}{cccc}0& d_{1,2}& ...& d_{1,I}\\ 0& d_{2,2}& ...& d_{2,I}\\ ...& ...& ...& ...\\ 0& d_{I,2}& ...& d_{I,I}\end{array}\right]---\left(2\right)$

并且同时将FMat_D中元素在各行中的原始序号记录在矩阵中得到：

${\mathrm{IMat}}_D=\left[\begin{array}{cccc}1& p_{1,2}& ...& p_{1,I}\\ 2& p_{2,2}& ...& p_{1,1}\\ ...& ...& ...& ...\\ I& p_{I,2}& ...& p_{1,I}\end{array}\right]---\left(3\right)$

因此元素p_i,j表示FMat_D中的元素d_i,j元素为中D_i与之间的距离。同样，对于星点集合也可以得到相应的矩阵FMat_R和IMat_R。其中FMat_D和FMat_R中每行为相应集合中各个星点的特征。

在匹配过程中为比较星点D_x与R_y之间的相似性，需要计算序列<0,d_x,i,...,d_x,I>_D和序列<0,d_y,j,...,d_x,J>_R之间的公共子序列的长度，其中x和y为星点序号。基于一般动态规划求解LCS的算法框架，我们设定其中的公共子序列更新公式为：

$c\&lsqb;i_1,j_1\&rsqb;=\left\{\begin{array}{cc}0,& i_1=0or{j}_1=0\\ c\&lsqb;i_1-1,j_1-1\&rsqb;+1,& i_1,j_1>0and\left|\right|d_{x,i_1}-d_{y,j_1}\left|\right|\&le;\mathrm{\&epsiv;}f(d_{x,i_1}-d_{y,j_1})\\ \max (c\&lsqb;i_1,j_1-1\&rsqb;,c\&lsqb;i_1-1,j_1\&rsqb;),& i_1,j_1>0and\left|\right|d_{x,i_1}-d_{y,j_1}\left|\right|>\mathrm{\&epsiv;}f(d_{x,i_1}-d_{y,j_1})\end{array}\right.---\left(4\right)$

其中，i₁和j₁为一般动态规划算法表中的序号，ε为距离阈值，为[0,1]之间的函数，成立时，表示在畸变情况下而函数f(·)为：

$f(d_{x,i_1},d_{y,j_1})=\frac{1}{L}\exp \left\{\frac{\left|\right|C_x-C\left|\right|+\left|\right|D_{{\mathrm{IMat}}_D\&lsqb;x,i_1\&rsqb;}-C\left|\right|+\left|\right|R_y-C\left|\right|+\left|\right|R_{{\mathrm{IMat}}_R\&lsqb;y,j_1\&rsqb;}-C\left|\right|}{4{\mathrm{\&sigma;}}^2}\right\}---\left(5\right)$

其中，C为图像畸变中心，本实施例初始假设为图像中心，在迭代中更新，||·||表示欧式距离，σ²为方差，本文中取0.5，L为归一化因子。

经过计算，得到与中任两个星点之间的LCS相似性度量。对于中星点，依次在中搜索LCS最近邻，如果与最近邻之间的关系满足阈值要求，则作为与相应D_i匹配的R_i，从而得到点对，否则舍弃对相应D_i的匹配。本实施例设定阈值为与之间最长公共子序列的长度的1/2。

经过星图匹配，得到匹配点对的集合其中P_m＝<R_m,D_m>。中包含足够多的正确配对点，但仍包含匹配错误点对。

(2)基于RANSAC的外点剔除与旋转平移模型估计。

中依然存在匹配错误点对，但其中正确匹配的点对之间存在稳定的旋转平移关系，而不符合该模型的匹配错误点对为外点。使用RANSAC估计其旋转平移关系，同时剔除其中匹配错误的外点。

和之间的旋转平移关系可以描述为二维空间中的旋转平移关系。由到的转换关系T(·)可以描述为：

$\left[\begin{array}{c}{D}_{xm}\\ D_{ym}\end{array}\right]=T\left(\left[\begin{array}{c}{R}_{xm}\\ R_{ym}\end{array}\right]\right)=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\&CenterDot;\left[\begin{array}{c}{R}_{xm}\\ R_{ym}\end{array}\right]+\left[\begin{array}{c}e\\ f\end{array}\right]---\left(6\right)$

其中，D_xm表示第m个点的x方向坐标，D_ym、R_xm和R_ym同理。对于给定的N对匹配点对，根据式(6)建立方程

Ax＝b(7)

$A=\left[\begin{array}{cccccc}{R}_{x1}& R_{y1}& 1& 0& 0& 0\\ 0& 0& 0& R_{x1}& R_{y1}& 1\\ ...& ...& ...& ...& ...& ...\\ R_{xM}& R_{xM}& 1& 0& 0& 0\\ 0& 0& 0& R_{xM}& R_{xM}& 1\end{array}\right]---\left(8\right)$

x＝[abecdf]^T(9)

b＝[D_x1D_y1...D_xMD_yM]^T(10)

方程(7)有封闭形式解，因此旋转平移模型参数可以求解为

x＝(A^TA)^-1A^TB(11)

RANSAC算法通过在每次迭代中对包含外点的点对集合随机采样估计模型参数，得到符合该模型的候选点对集，并计算误差。最终从所有迭代中取出拟合误差最小的模型作为估计结果，相应的候选点对集合为剔除外电的集合。在所有迭代中选取最佳。

RANSAC的算法输入为：含有外点的点对集合每次采样最少点数目Q，迭代次数K，判断一个点对是否符合当前模型的误差阈值T，判断估计得到的模型是否足够好的候选集合大小D。算法流程为：迭代K次，对于第k次迭代，1≤k≤K，

1)在中随机选取Q个点对得到候选点集{P_q}_con，此时，1≤q≤Q；

2)使用式(11)估计参数x_k，得到投影模型T_k(·)；

3)对于除{P_q}_con中点对外的其他点对<R_m,D_m>：计算投影误差e_k,m＝||D_n-T_k(R_n)||，如果e_k,m≤T，则将<R_m,D_m>加入候选点对集合{P_q}_con，每次更新使{P_q}_con变大，使q的取值范围上限变大；

4)如果{P_q}_con元素数目大于D，则使用{P_q}_con估计式(11)中参数x_con，得到投影模型T_con(·)；

5)计算投影误差 $e_k=\underset{<R_q,D_q>\&Element;{\left\{p_q\right\}}_{con}}{\&Sigma;}\left|\right|D_q-T_{con}\left(R_q\right)\left|\right|;$

对于每次迭代，求取e_k最小的T_con(·)作为输出模型参数，同时相应的{P_n}_con作为剔除外点的点集。本实施例参数选取为，N≥6，K＝300，T＝1，D与星点数目、图像尺度和畸变程度相关，当星点数目约为50时取为10。

使用RANSAC剔除中外点得到集合估计得到相应旋转平移参数x，即和之间的转换包含的旋转平移函数T(·)。

(3)迭代估计几何畸变。

在得到不包含外点的匹配点对集合和旋转平移函数T(·)之后，中的经过T(·)变换有与相应的组合得到点对集合因此，中与之间仅包含集合畸变的变换，可以被用来估计几何畸变。

通常情况下，仅使用径向畸变来描述{P′_n}中的几何畸变。由{R′_n}到{D_n}的几何畸变可以为

$\left[\begin{array}{c}{D}_{xn}\\ D_{yn}\end{array}\right]=G\left(\left[\begin{array}{c}{R}_{xn}^{\&prime;}\\ R_{yn}^{\&prime;}\end{array}\right]\right)=\left[\begin{array}{c}{G}_x\left(R_{xn}^{\&prime;}\right)\\ G_y\left(R_{yn}^{\&prime;}\right)\end{array}\right]=L\left(r\right)\left[\begin{array}{c}{R}_{xn}^{\&prime;}-C_x\\ R_{yn}^{\&prime;}-C_y\end{array}\right]+\left[\begin{array}{c}{C}_x\\ C_y\end{array}\right]---\left(12\right)$

其中，D_xn表示第n个点的x方向坐标，D_yn、R′_xn和R′_yn同理，C_x为畸变中心x方向坐标，C_y为畸变中心y方向坐标，G(·)表示畸变函数，表示[R′_xn,R′_yn]^T到畸变中心C的距离，一般情况下，畸变中心不一定与图像几何中心重叠，但通常接近集合中心。L(r)为失真函数，具有性质为L(0)＝1，通常使用其在0处的泰勒展开表示，仅使用偶次项进行表示，并舍弃高次项，因此L(r)表示为

L(r)＝1+k₁r²+k₂r⁴(13)

模型仅由畸变中心C＝[C_x,C_y]^T和畸变程度k₁、k₂确定，为方便描述，将参数表示为P＝[k₁,k₂,C_x,C_y]。

与y方向无差别地，x方向的畸变模型G_x(·)可以表示为：

$\begin{array}{c}{D}_{xn}=G_x(R_{xn}^{\&prime;};P)=(R_{xn}^{\&prime;}-C_x)+k_1{r}^2(R_{xn}^{\&prime;}-C_x)+k_2{r}^4(R_{xn}^{\&prime;}-C_x)\\ r=\sqrt{{(R_{xn}^{\&prime;}-C_x)}^2+{(R_{yn}^{\&prime;}-C_y)}^2}\end{array}---\left(14\right)$

其中，n为点对序号。同理，G_y(·)也可以由此表示。

使用牛顿迭代对该非线性模型进行求解。牛顿迭代是一个迭代优化算法，给定初始值参数初始值P₀，在迭代中求得最优解。迭代中参数更新方式为：

P_k+1＝P_k+Δ_k(15)

其中，k为迭代次数序号，且只有当k出现在下标与[]中时表示迭代次数，Δ_k为方程J_kΔ_k＝ε_k的最小二乘解，为取值在P_k上的雅各比矩阵，而ε_k＝D-G(R；P_k)为当前估计的误差。在每次迭代中对于每个值代入J之后均有：

$J_{k,n}(R_n;P_k)=\left[\begin{array}{c}\frac{\&part;G_x(R_{nx};P)}{\&part;P}\&lsqb;k,n\&rsqb;\\ \frac{\&part;G_y(R_{ny};P)}{\&part;P}\&lsqb;k,n\&rsqb;\end{array}\right]---\left(16\right)$

其中

$\frac{\&part;G_x(R_x;P)}{\&part;P}\&lsqb;k,n\&rsqb;=\&lsqb;\frac{\&part;G_x(R_x;k_1)}{\&part;k_1}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;k_2)}{\&part;k_2}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;C_x)}{\&part;C_x}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;C_y)}{\&part;C_y}\&lsqb;k,n\&rsqb;\&rsqb;---\left(17\right)$

而可以类似地表示。针对G(·)中的G_x(·)，具体地有G_x(·)针对不同参数求偏导：

$\begin{array}{c}\frac{\&part;G_x(R_x;P)}{\&part;k_1}\&lsqb;k,n\&rsqb;={r_{k,n}}^2{R}_{xn}^{\&prime;}\\ \frac{\&part;G_x(R_x;P)}{\&part;k_2}\&lsqb;k,n\&rsqb;={r_{k,n}}^4{R}_{xn}^{\&prime;}\\ \frac{\&part;G_x(R_x;P)}{\&part;C_x}\&lsqb;k,n\&rsqb;=-k_{1,k}{r_{k,n}}^2-2k_{1,k}{(R_{xn}-C_{x,k})}^2-k_{2,k}{r_{k,n}}^4-4k_{2,k}{r_{n,k}}^2{(R_{nx}-C_{x,k})}^2\\ \frac{\&part;G_x(R_x;P)}{\&part;C_y}\&lsqb;k,n\&rsqb;=-2k_{1,k}(R_{xn}-C_{x,k}){(R_{yn}-C_{y,k})}^2-4k_{2,k}{r}^2(R_{nx}-C_{x,k})(R_{ny}-C_{y,k})\end{array}---\left(18\right)$

其中，下标出现的k均表示第k次迭代中相应的值，下标中的n表示相应的量为代入第n点对之后的相应值，对于G_y(·)可以得到相同的求偏导结果。

基于以上，在每次迭代中求解方程J_iΔ_i＝ε_i，在利用到所有点对的情况下有：

J_i＝[J_i(R₁；P_i),J_i(R₂；P_i),...,J_i(R_N；P_i)]^T(19)

ε_i＝[(D₁-G(R₁；P_i))^T,(D₂-G(R₂；P_i))^T,...,(D_N-G(R_N；P_i))^T]^T(20)

其中(D₁-G(R₁；P_i))＝[D_x1-G_x(R_x1；P_i),D_y1-G_y(R_y1；P_i)]^T。于是，解Δ_i为

Δ_i＝J⁺ε_i＝(J^TJ)^-1J^Tε_i(21)

根据以上，迭代优化求解畸变模型的具体算法表示为，算法输入为输入：匹配点对集合模型G(·)，参数初始值P₀，迭代最大次数K，收敛误差ε_con。算法迭代过程中，如果迭代次数小于最大迭代次数k<K，则进行如下迭代，对于第k次迭代：

1)利用式(20)代入当前参数值P_k，计算ε_k，如果||ε_k||<ε_con，停止迭代，输出结果，否则继续下一步；

2)将和P_i代入式(19)，计算J_k；

3)利用式(21)计算得到Δ_k；

4)更新参数P_k+1＝P_k+Δ_k；

停止迭代，如果停止迭代时k的值为K，否则即为满足1)中迭代停止条件时的迭代步骤中的输出结果。对于以上参数，本实施例设置P₀中的k₁＝k₂＝7×10^-7，[C_x,C_y]初始值为图像集合中心，迭代最大次数为K＝20，ε_con与具体数据相关，可以默认设置为ε_con＝0.5。

(4)联合优化迭代框架。

通过综合以上的星点匹配、RANSAC估计旋转平移和畸变估计方法，利用联合迭代框架检测和估计旋转平移和畸变模型。整体迭代算法框架输入为星表参考星点集合星空图像星点集合星图匹配、旋转平移估计、畸变估计所需参数，最大迭代次数K，迭代停止误差ε_con。迭代算法首先初始化，令 为算法中用于传递每一步畸变估计结果的点集。令k从1到K进行循环迭代：

1)使用改进的LCS匹配算法对和进行点对匹配得到匹配点对集合如果为空集，则说明和之间变换程度过大，检测到成像系统指向偏差过大或者图像畸变程度过大，退出迭代，输出结果，终止程序；

2)将作为输入使用RANSAC方法估计旋转平移模型T_k(·)，同时得到剔除匹配外点的点对集合若估计失败或过小不足以进行后续运算，则说明和之间畸变程度过大，检测到成像系统指向偏差过大或者图像畸变程度过大，退出迭代，输出结果，终止程序；

3)利用中点对的序号关系，得到对应的点对集合使用T_k(·)对{R_n}进行投影得到{P′_n}＝{<R′_n,D_n>}＝{<T(R_n),D_n>}；

4)将{P′_n}代入，迭代估计得到畸变模型G_k(·)；如果估计失败，说明图像畸变程度过大，退出迭代，输出结果，终止程序；

5)计算误差 ${\mathrm{\&epsiv;}}_k=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\right|G_k\left(R_n^{\&prime;}\right)-D_n\left|\right|,$ 如果ε_k<ε_con，则 $\hat{T}(\&CenterDot;)=T_k(\&CenterDot;),\hat{G}(\&CenterDot;)=G_k(\&CenterDot;),$ 停止循环，输出结果，否则，进行下一步；

6)更新有同时更新匹配中使用的图像畸变中心为[C_x,C_y]，进入下一次循环。

如果最终循环进行了K次，则输出结果根据输出的和对图像与星表相比的旋转平移和畸变作出量化评价。参数选择与畸变程度相关，本实施例参数选择为，K＝10，ε_con＝0.005。

Claims (1)

1.一种基于星图匹配的星空图像畸变检测与估计方法，其特征在于包括以下步骤：

步骤一、对星空图像中的星点进行提取，并计算星点在图像上坐标；该过程包括，利用阈值对星空图像进行二值化分割，即根据像素值是否大于阈值而进行分割；对二值图像进行联通区域提取，得到星点光斑区域集合Ω_n；利用提取到的光斑区域，从原始图像中计算星点的亚像素坐标(x_n,y_n)，对于每个光斑，计算方法为

$x_n=\frac{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)\&CenterDot;x}{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)},y_n=\frac{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)\&CenterDot;y}{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_n}{\&Sigma;}I(x,y)}---\left(1\right)$ 式中，(x,y)∈Ω_n表示属于光斑Ω_n的像素坐标，I(x,y)为相应的像素值；对图像中的星点进行计算得到星点集合I为集合中星点数目；

步骤二、根据拍摄时的输入指向、视场大小、相机参数和拍摄时间，通过星空坐标转换得到对应天区的星体在CCD成像平面上的理想成像位置作为参考星点集合其中，J为集合中星点数目；

步骤三、像畸变检测与估计框架为在迭代中完成星图匹配、旋转平移估计和畸变估计，具体过程包括：

(a).使用基于LCS的星图配准方法根据星图畸变模型进行改进，对和进行匹配，若匹配失败，说明图像与输入的指向之间存在较大偏差或者图像畸变过于严重，完成星图畸变检测；若匹配成功，得到可能存在匹配错误的点对集合，进行进一步的检测与估计；

(b).使用一致随机采样对匹配外点进行剔除，并同时估计图像星点与理想星表星点之间的旋转平移模型，若剔除外点的匹配点对集合中点对个数小于求解畸变模型所需点对个数，则说明图像与输入的指向之间存在较大偏差或者图像畸变过于严重，完成星图畸变检测；否则，利用剔除外点后的点对集合与估计得到的旋转平移模型进行进一步的畸变检测与估计；

(c)利用无外点的、消除空间旋转平移差异的匹配集合估计图像中的几何畸变；

(d)迭代步骤(a)-(c)估计得到和之间的旋转平移模型和畸变模型，通过估计得的旋转平移和畸变模型得到图像的畸变程度的量化评价；

其中，步骤(a)-(c)的迭代中，包括基于畸变模型的LCS畸变星图配准、基于RANSAC的外点剔除与旋转平移模型估计、迭代估计几何畸变；

第一部分，基于畸变模型的LCS畸变星图配准；假设从星空图像中提取得到的星点集合为而通过给定参数得到的市场内星点在CCD平面的投影的星点集合为对于得到一个大小为I×I的邻接矩阵，其中每个元素表示中两点之间在图像平面的欧氏距离，对该矩阵各行进行升序排序得到特征矩阵：

${\mathrm{FMat}}_D=\left[\begin{array}{cccc}0& d_{1,2}& ...& d_{1,I}\\ 0& d_{2,2}& ...& d_{2,I}\\ ...& ...& ...& ...\\ 0& d_{I,2}& ...& d_{I,I}\end{array}\right]---\left(2\right)$

并且同时将FMat_D中元素在各行中的原始序号记录在矩阵中得到：

${\mathrm{IMat}}_D=\left[\begin{array}{cccc}1& p_{1,2}& ...& p_{1,I}\\ 2& p_{2,2}& ...& p_{1,1}\\ ...& ...& ...& ...\\ I& p_{I,2}& ...& p_{1,I}\end{array}\right]---\left(3\right)$

因此，元素p_i,j表示FMat_D中的元素为中D_i与之间的距离；同样，对于星点集合得到相应的矩阵FMat_R和IMat_R；其中FMat_D和FMat_R中每行为相应集合中各个星点的特征；

在匹配过程中为比较星点D_x与R_y之间的相似性，需要计算序列<0,d_x,i,...,d_x,I>_D和序列<0,d_y,j,...,d_x,J>_R之间的公共子序列的长度，其中x和y为星点序号；基于一般动态规划求解LCS的算法框架，设定其中的公共子序列更新公式为：

$c\&lsqb;i_1,j_1\&rsqb;=\left\{\begin{array}{cc}0,& i_1=0or{j}_1=0\\ c\&lsqb;i_1-1,j_1-1\&rsqb;+1,& i_1,j_1>0and\left|\right|d_{x,i_1}-d_{y,j_1}\left|\right|\&le;\mathrm{\&epsiv;}f(d_{x,i_1}-d_{y,j_1})\\ \max (c\&lsqb;i_1,j_1-1\&rsqb;,c\&lsqb;i_1-1,j_1\&rsqb;),& i_1,j_1>0and\left|\right|d_{x,i_1}-d_{y,j_1}\left|\right|>\mathrm{\&epsiv;}f(d_{x,i_1}-d_{y,j_1})\end{array}\right.---\left(4\right)$

式中，i₁和j₁为动态规划算法表中的序号，ε为距离阈值，为[0,1]之间的函数，成立时，表示在畸变情况下而函数f(·)为：

$f(d_{x,i_1},d_{y,j_1})=\frac{1}{L}\exp \left\{\frac{\left|\right|D_x-C\left|\right|+\left|\right|D_{{\mathrm{IMat}}_D\&lsqb;x,i_1\&rsqb;}-C\left|\right|+\left|\right|R_y-C\left|\right|+\left|\right|R_{{\mathrm{IMat}}_R\&lsqb;y,j_1\&rsqb;}-C\left|\right|}{4{\mathrm{\&sigma;}}^2}\right\}---\left(5\right)$

式中，C为图像畸变中心，在迭代中更新，||·||表示欧式距离，σ²为方差，L为归一化因子；

得到与中任两个星点之间的LCS相似性度量；对于中星点，依次在中搜索LCS最近邻，如果与最近邻之间的关系满足阈值要求，则作为与相应D_i匹配的R_i，从而得到点对，否则舍弃对相应D_i的匹配；

经过星图匹配，得到匹配点对的集合其中P_m＝<R_m,D_m>；中包含足够多的正确配对点，但仍包含匹配错误点对；

第二部分，基于RANSAC的外点剔除与旋转平移模型估计；中依然存在匹配错误点对，但其中正确匹配的点对之间存在稳定的旋转平移关系，而不符合该模型的匹配错误点对为外点；使用RANSAC估计其旋转平移关系，同时剔除其中匹配错误的外点；

和之间的旋转平移关系描述为二维空间中的旋转平移关系；由到的转换关系T(·)描述为：

$\left[\begin{array}{c}{D}_{xm}\\ D_{ym}\end{array}\right]=T\left(\left[\begin{array}{c}{R}_{xm}\\ R_{ym}\end{array}\right]\right)=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\&CenterDot;\left[\begin{array}{c}{R}_{xm}\\ R_{ym}\end{array}\right]+\left[\begin{array}{c}e\\ f\end{array}\right]---\left(6\right)$

其中，D_xm表示第m个点的x方向坐标，D_ym、R_xm和R_ym同理；

在估计模型T(·)过程中，对于给定的N对匹配点对，由于存在测量误差和畸变未知引起的误差，点对中点坐标存在噪声，使用足够多的点对求取最优估计；根据式(6)建立方程

Ax＝b(7)

$A=\left[\begin{array}{cccccc}{R}_{x1}& R_{y1}& 1& 0& 0& 0\\ 0& 0& 0& R_{x1}& R_{y1}& 1\\ ...& ...& ...& ...& ...& ...\\ R_{xM}& R_{xM}& 1& 0& 0& 0\\ 0& 0& 0& R_{xM}& R_{xM}& 1\end{array}\right]---\left(8\right)$

x＝[abecdf]^T(9)

b＝[D_x1D_y1...D_xMD_yM]^T(10)

方程(7)有封闭形式解，因此旋转平移模型参数求解为

x＝(A^TA)^-1A^TB(11)

RANSAC算法通过在每次迭代中对包含外点的点对集合随机采样估计模型参数，得到符合该模型的候选点对集，并计算误差；最终从所有迭代中取出拟合误差最小的模型作为估计结果，相应的候选点对集合为剔除外电的集合；在所有迭代中选取最佳；

RANSAC的算法输入为：含有外点的点对集合每次采样最少点数目Q，迭代次数K，判断一个点对是否符合当前模型的误差阈值T，判断估计得到的模型是否足够好的候选集合大小D；算法流程为：迭代K次，对于第k次迭代，1≤k≤K，

(e)在中随机选取Q个点对得到候选点集{P_q}_con，此时，1≤q≤Q；

(f)使用式(11)估计参数x_k，得到投影模型T_k(·)；

(g)对于除{P_q}_con中点对外的其他点对<R_m,D_m>：计算投影误差e_k,m＝||D_n-T_k(R_n)||，如果e_k,m≤T，则将<R_m,D_m>加入候选点对集合{P_q}_con，每次更新使{P_q}_con变大，使q的取值范围上限变大；

(h)如果{P_q}_con元素数目大于D，则使用{P_q}_con估计式(11)中参数x_con，得到投影模型T_con(·)；

(i)计算投影误差 $e_k=\underset{<R_q,D_q>\&Element;{\left\{P_q\right\}}_{con}}{\&Sigma;}\left|\right|D_q-T_{con}\left(R_q\right)\left|\right|;$

对于每次迭代，求取e_k最小的T_con(·)作为输出模型参数，同时相应的{P_n}_con作为剔除外点的点集；

使用RANSAC剔除中外点得到集合估计得到相应旋转平移参数x，即和之间的转换包含的旋转平移函数T(·)；

第三部分，迭代估计几何畸变；在得到不包含外点的匹配点对集合和旋转平移函数T(·)之后，中的经过T(·)变换有与相应的组合得到点对集合因此，中与之间仅包含集合畸变的变换，被用来估计几何畸变；使用径向畸变来描述中的几何畸变；由到的几何畸变为

$\left[\begin{array}{c}{D}_{xn}\\ D_{yn}\end{array}\right]=G\left(\left[\begin{array}{c}{R}_{xn}^{\&prime;}\\ R_{yn}^{\&prime;}\end{array}\right]\right)=\left[\begin{array}{c}{G}_x\left(R_{xn}^{\&prime;}\right)\\ G_y\left(R_{yn}^{\&prime;}\right)\end{array}\right]=L\left(r\right)\left[\begin{array}{c}{R}_{xn}^{\&prime;}-C_x\\ R_{yn}^{\&prime;}-C_y\end{array}\right]+\left[\begin{array}{c}{C}_x\\ C_y\end{array}\right]---\left(12\right)$

其中，D_xn表示第n个点的x方向坐标，D_yn、R′_xn和R′_yn同理，C_x为畸变中心x方向坐标，C_y为畸变中心y方向坐标，G(·)表示畸变函数，表示[R′_xn,R′_yn]^T到畸变中心C的距离，L(r)为失真函数，具有性质为L(0)＝1，使用其在0处的泰勒展开表示，仅使用偶次项并舍弃高次项，因此L(r)表示为；

L(r)＝1+k₁r²+k₂r⁴(13)

模型仅由畸变中心C＝[C_x,C_y]^T和畸变程度k₁、k₂确定，为方便描述，将参数表示为P＝[k₁,k₂,C_x,C_y]；

与y方向无差别地，x方向的畸变模型G_x(·)表示为：

D_xn＝G_x(R′_xn；P)＝(R′_xn-C_x)+k₁r²(R′_xn-C_x)+k₂r⁴(R′_xn-C_x)

(14)

$r=\sqrt{{(R_{xn}^{\&prime;}-C_x)}^2+{(R_{yn}^{\&prime;}-C_y)}^2}$

其中，n为点对序号；同理，G_y(·)也由此表示；

使用牛顿迭代对畸变模型G_x(·)进行求解；给定初始值参数初始值P₀，在迭代中求得最优解；迭代中参数更新方式为：

P_k+1＝P_k+△_k(15)

其中，k为迭代次数序号，且只有当k出现在下标与[]中时表示迭代次数，△_k为方程J_k△_k＝ε_k的最小二乘解，为取值在P_k上的雅各比矩阵，而ε_k＝D-G(R；P_k)为当前估计的误差；在每次迭代中对于每个值代入J之后均有：

$J_{k,n}(R_n;P_k)=\left[\begin{array}{c}\frac{\&part;G_x(R_{nx};P)}{\&part;P}\&lsqb;k,n\&rsqb;\\ \frac{\&part;G_y(R_{ny};P)}{\&part;P}\&lsqb;k,n\&rsqb;\end{array}\right]---\left(16\right)$

其中

$\frac{\&part;G_x(R_x;P)}{\&part;P}\&lsqb;k,n\&rsqb;=\&lsqb;\frac{\&part;G_x(R_x;k_1)}{\&part;k_1}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;k_2)}{\&part;k_2}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;k_x)}{\&part;C_x}\&lsqb;k,n\&rsqb;,\frac{\&part;G_x(R_x;k_y)}{\&part;C_y}\&lsqb;k,n\&rsqb;\&rsqb;---\left(17\right)$

而类似地表示；针对G(·)中的G_x(·)，具体地有G_x(·)针对不同参数求偏导：

$\frac{\&part;G_x(R_x;P)}{\&part;k_1}\&lsqb;k,n\&rsqb;={r_{k,n}}^2{R}_{xn}^{\&prime;}$

$\frac{\&part;G_x(R_x;P)}{\&part;k_2}\&lsqb;k,n\&rsqb;={r_{k,n}}^4{R}_{xn}^{\&prime;}$

(18)

$\frac{\&part;G_x(R_x;P)}{\&part;C_x}\&lsqb;k,n\&rsqb;=-k_{1,k}{r_{k,n}}^2-2k_{1,k}{(R_{xn}-C_{x,k})}^2-k_{2,k}{r_{k,n}}^4-4k_{2,k}{r_{n,k}}^2{(R_{nx}-C_{x,k})}^2$

$\frac{\&part;G_x(R_x;P)}{\&part;C_y}\&lsqb;k,n\&rsqb;=-2k_{1,k}(R_{xn}-C_{x,k})(R_{yn}-C_{y,k})-4k_{2,k}{r}^2(R_{nx}-C_{x,k})(R_{ny}-C_{y,k})$

其中，下标中的k均表示第k次迭代中相应的值，下标中的n表示相应的量为代入第n点对之后的相应值，对于G_y(·)得到相同的求偏导结果；

基于以上，在每次迭代中求解方程J_i△_i＝ε_i，在利用到所有点对的情况下有：

J_i＝[J_i(R₁；P_i),J_i(R₂；P_i),...,J_i(R_N；P_i)]^T(19)

ε_i＝[(D₁-G(R₁；P_i))^T,(D₂-G(R₂；P_i))^T,...,(D_N-G(R_N；P_i))^T]^T(20)

其中，(D₁-G(R₁；P_i))＝[D_x1-G_x(R_x1；P_i),D_y1-G_y(R_y1；P_i)]^T；于是，解△_i为

△_i＝J⁺ε_i＝(J^TJ)^-1J^Tε_i(21)其中，J⁺＝(J^TJ)^-1J^T；

根据以上，迭代优化求解畸变模型的具体算法表示为，算法输入为输入：匹配点对集合模型G(·)，参数初始值P₀，迭代最大次数K，收敛误差ε_con；算法迭代过程中，如果迭代次数小于最大迭代次数k<K，则进行如下迭代，对于第k次迭代：

(j)利用式(20)代入当前参数值P_k，计算ε_k，如果||ε_k||<ε_con，停止迭代，输出结果，否则继续下一步；

(k)将和P_i代入式(19)，计算J_k；

(l)利用式(21)计算得到△_k；

(m)更新参数P_k+1＝P_k+△_k；

停止迭代，如果停止迭代时k的值为K，否则即为满足步骤(j)中迭代停止条件时的迭代步骤中的输出结果；

通过综合以上三部分，利用联合迭代框架检测和估计旋转平移和畸变模型；整体迭代算法框架输入为：星表参考星点集合星空图像星点集合星图匹配、旋转平移估计、畸变估计所需参数，最大迭代次数K，迭代停止误差ε_con；迭代算法首先初始化，令为算法中用于传递每一步畸变估计结果的点集，令k从1到K进行循环迭代：

(n)使用改进的LCS匹配算法对和进行点对匹配得到匹配点对集合如果为空集，则说明和之间变换程度过大，检测到成像系统指向偏差过大或者图像畸变程度过大，退出迭代，输出结果，终止程序；

(o)将作为输入使用RANSAC方法估计旋转平移模型T_k(·)，同时得到剔除匹配外点的点对集合若估计失败或过小不足以进行后续运算，则说明和之间畸变程度过大，检测到成像系统指向偏差过大或者图像畸变程度过大，退出迭代，输出结果，终止程序；

(p)利用中点对的序号关系，得到对应的点对集合使用T_k(·)对{R_n}进行投影得到{P′_n}＝{<R′_n,D_n>}＝{<T(R_n),D_n>}；

(q)将{P′_n}代入，迭代估计得到畸变模型G_k(·)；如果估计失败，说明图像畸变程度过大，退出迭代，输出结果，终止程序；

(r)计算误差 ${\mathrm{\&epsiv;}}_k=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}\left|\right|G_k\left(R_n^{\&prime;}\right)-D_n\left|\right|,$ 如果ε_k＜ε_con，则 $\hat{T}(\&CenterDot;)=T_k(\&CenterDot;),\hat{G}(\&CenterDot;)=G_k(\&CenterDot;),$ 停止循环，输出结果，否则，进行下一步；

(s)更新有 ${\left\{R_j^{*}\right\}}_{j=1}^J=T_k^{-1}\left(G_k\right(T_k\left({\left\{R_j\right\}}_{j=1}^J\right)\left)\right),$ 进入下一次循环；

如果最终循环进行了K次，则输出结果根据输出的T和对图像与星表相比的旋转平移和畸变作出量化评价。