[发明专利]基于星图匹配的星空图像畸变检测与估计方法

权利要求书

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

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

xn=Σ(x,y)∈ΩnI(x,y)·xΣ(x,y)∈ΩnI(x,y),yn=Σ(x,y)∈ΩnI(x,y)·yΣ(x,y)∈ΩnI(x,y)---(1)]]>

式中，(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的邻接矩阵，其中每个元素表示中两点之间在图像平面的欧氏距离，对该矩阵各行进行升序排序得到特征矩阵：

FMatD=0d1,2...d1,I0d2,2...d2,I............0dI,2...dI,I---(2)]]>

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

IMatD=1p1,2...p1,I2p2,2...p1,I............IpI,2...p1,I---(3)]]>

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

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

c[i1,j1]=0,i1=0orj1=0c[i1-1,j1-1]+1,i1,j1>0and||dx,i1-dy,j1||≤ϵf(dx,i1,dy,j1)max(c[i1,j1-1],c[i1-1,j1]),i1,j1>0and||dx,i1-dy,j1||>ϵf(dx,i1,dy,j1)---(4)]]>

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

f(dx,i1,dy,j1)=1Lexp{||Dx-C||+||DIMatD[x,i1]-C||+||Ry-C||+||RIMatR[y,j1]-C||4σ2}---(5)]]>

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

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

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

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

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

DxmDym=T(RxmRym)=abcd·RxmRym+ef---(6)]]>

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

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

Ax＝b          (7)

A=Rx1Ry11000000Rx1Ry11..................RxMRxM1000000RxMRxM1---(8)]]>

x＝[a b e c d f]^T          (9)

b＝[D_x1 D_y1...D_xM D_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)计算投影误差ek=Σ<Rq,Dq>∈{Pq}con||Dq-Tcon(Rq)||;]]>

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

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

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

DxnDyn=GRxn′Ryn′=Gx(Rxn′)Gy(Ryn′)=L(r)Rxn′-CxRyn′-Cy+CxCy---(12)]]>

其中，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=(Rxn′-Cx)2+(Ryn′-Cy)2]]>

其中，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之后均有：

Jk,n(Rn;Pk)=∂Gx(Rnx;P)∂P[k,n]∂Gy(Rny;P)∂P[k,n]---(16)]]>

其中

∂Gx(Rx;P)∂P[k,n]=[∂Gx(Rx;k1)∂k1[k,n],∂Gx(Rx;k2)∂k2[k,n],∂Gx(Rx;Cx)∂Cx[k,n],∂Gx(Rx;Cy)∂Cy[k,n]]---(17)]]>而类似地表示；针对G(·)中的G_x(·)，具体地有G_x(·)针对不同参数求偏导：

∂Gx(Rx;P)∂k1[k,n]=rk,n2Rxn′]]>

∂Gx(Rx;P)∂k2[k,n]=rk,n4Rxn′]]>    (18)

∂Gx(Rx;P)∂Cx[k,n]=-k1,krk,n2-2k1,k(Rxn-Cx,k)2-k2,krk,n4-4k2,krn,k2(Rnx-Cx,k)2]]>

∂Gx(Rx;P)∂Cy[k,n]=-2k1,k(Rxn-Cx,k)(Ryn-Cy,k)-4k2,kr2(Rnx-Cx,k)(Rny-Cy,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次迭代：

(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)计算误差ϵk=1NΣn=1N||Gk(Rn′)-Dn||,]]>如果ε_k<ε_com,则T^(·)=Tk(·),]]>G^(·)=Gk(·),]]>停止循环，输出结果，否则，进行下一步；

(f)更新{Rj*}j=1J,]]>有{Rj*}j=1J=Tk-1(Gk(Tk({Rj}j=1J))),]]>进入下一次循环；

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