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

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

Info

Publication number
CN103440659B
CN103440659B CN201310390034.1A CN201310390034A CN103440659B CN 103440659 B CN103440659 B CN 103440659B CN 201310390034 A CN201310390034 A CN 201310390034A CN 103440659 B CN103440659 B CN 103440659B
Authority
CN
China
Prior art keywords
distortion
point
iteration
asterism
model
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201310390034.1A
Other languages
English (en)
Other versions
CN103440659A (zh
Inventor
张艳宁
巩东
孙瑾秋
李海森
丁王斌
黄建余
惠建江
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN201310390034.1A priority Critical patent/CN103440659B/zh
Publication of CN103440659A publication Critical patent/CN103440659A/zh
Application granted granted Critical
Publication of CN103440659B publication Critical patent/CN103440659B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

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

Description

基于星图匹配的星空图像畸变检测与估计方法
技术领域
本发明涉及一种星空图像畸变检测与估计方法,特别是涉及一种基于星图匹配的星空图像畸变检测与估计方法。
背景技术
文献“基于星图匹配的小视场空间图像畸变校正[J].北京航空航天大学学报,2008,09:1016-1019.”提出了一种基于星图匹配的小视场空间图像畸变校正方法。该方法基于Hausdorff距离进行多帧空间图像中星点的匹配,并在通过误差项引入针对指向偏差的补偿。该方法通过匹配点之间的偏差对全局进行校正。但是由于该方法主要针对小视场图像,视场中星点较少,因此使用Hausdorff距离进行适用的星点匹配,并且简单地使用星点之间的偏差对全局进行建模。但是对于一般情况下的星空图像,视场内星点较多,尤其是在图像遭受畸变影响的情况下,Hausdorff距离难以进行准确的匹配,并且仅仅使用误差的插值难以对全局图像畸变进行准确描述。因而该方法无法直接应用于对含有几何畸变和指向偏差的一般星空图像畸变的检测和估计。
发明内容
为了克服现有检测方法在图像遭受畸变影响的情况下检测精度低的不足,本发明提供一种基于星图匹配的星空图像畸变检测与估计方法。该方法利用已知的星表、成像系统参数和成像系统指向,对图像中星点与星表中星点进行匹配,进而对图像中可能存在的相机的旋转、平移和畸变进行检测与估计。该方法充分考虑了畸变对星点匹配的影响、指向不精确和相机镜头平移旋转造成的图像中星点的偏移旋转、平移旋转与畸变综合作用造成的估计困难,利用基于畸变模型改进的最长公共子序列(LCS)星点配准方法对星空图像与星表中星点进行配准,对匹配点对集合中可能存在的错误匹配外点进行剔除并估计相对的旋转平移关系,在考虑旋转平移的情况下估计星图与星标之间的畸变。在以上过程中完成对星图中旋转平移和畸变的检测,同时通过联合迭代求解对旋转平移和畸变模型进行估计。
本发明解决其技术问题所采用的技术方案是:一种基于星图匹配的星空图像畸变检测与估计方法,其特点是包括以下步骤:
步骤一、对星空图像中的星点进行提取,并计算星点在图像上坐标。该过程包括,利用阈值对星空图像进行二值化分割,即根据像素值是否大于阈值而进行分割;对二值图像进行联通区域提取,得到星点光斑区域集合Ωn;利用提取到的光斑区域,从原始图像中计算星点的亚像素坐标(xn,yn),对于每个光斑,计算方法为
x n = Σ ( x , y ) ∈ Ω n I ( x , y ) · x Σ ( x , y ) ∈ Ω n I ( x , y ) , y n = Σ ( x , y ) ∈ Ω n I ( x , y ) · y Σ ( x , y ) ∈ Ω n I ( 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的邻接矩阵,其中每个元素表示中两点之间在图像平面的欧氏距离,对该矩阵各行进行升序排序得到特征矩阵:
FMat D = 0 d 1 , 2 ... d 1 , I 0 d 2 , 2 ... d 2 , I ... ... ... ... 0 d I , 2 ... d I , I - - - ( 2 )
并且同时将FMatD中元素在各行中的原始序号记录在矩阵中得到:
IMat D = 1 p 1 , 2 ... p 1 , I 2 p 2 , 2 ... p 1 , 1 ... ... ... ... I p I , 2 ... p 1 , I - - - ( 3 )
因此,元素pi,j表示FMat0中的元素di,j元素为中Di之间的距离。同样,对于星点集合得到相应的矩阵FMatR和IMatR。其中FMatD和FMatR中每行为相应集合中各个星点的特征。
在匹配过程中为比较星点Dx与Ry之间的相似性,需要计算序列<0,dx,i,...,dx,I>D和序列<0,dy,j,...,dx,J>R之间的公共子序列的长度,其中x和y为星点序号。基于一般动态规划求解LCS的算法框架,设定其中的公共子序列更新公式为:
c &lsqb; i 1 , j 1 &rsqb; = 0 , i 1 = 0 o r j 1 = 0 c &lsqb; i 1 - 1 , j 1 - 1 &rsqb; + 1 , i 1 , j 1 > 0 a n d | | d x , i 1 - d y , j 1 | | &le; &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 > 0 a n d | | d x , i 1 - d y , j 1 | | > &epsiv; f ( d x , i 1 - d y , j 1 ) - - - ( 4 )
式中,i1和j1为动态规划算法表中的序号,ε为距离阈值,为[0,1]之间的函数,成立时,表示在畸变情况下而函数f(·)为:
f ( d x , i 1 , d y , j 1 ) = 1 L exp { | | C x - C | | + | | D IMat D &lsqb; x , i 1 &rsqb; - C | | + | | R y - C | | + | | R IMat R &lsqb; y , j 1 &rsqb; - C | | 4 &sigma; 2 } - - - ( 5 )
式中,C为图像畸变中心,在迭代中更新,||·||表示欧式距离,σ2为方差,L为归一化因子。
得到中任两个星点之间的LCS相似性度量。对于中星点,依次在中搜索LCS最近邻,如果与最近邻之间的关系满足阈值要求,则作为与相应Di匹配的Ri,从而得到点对,否则舍弃对相应Di的匹配。
经过星图匹配,得到匹配点对的集合其中Pm=<Rm,Dm>。中包含足够多的正确配对点,但仍包含匹配错误点对。
第二部分,基于RANSAC的外点剔除与旋转平移模型估计。中依然存在匹配错误点对,但其中正确匹配的点对之间存在稳定的旋转平移关系,而不符合该模型的匹配错误点对为外点。使用RANSAC估计其旋转平移关系,同时剔除其中匹配错误的外点。
之间的旋转平移关系描述为二维空间中的旋转平移关系。由的转换关系T(·)描述为:
D x m D y m = T ( R x m R y m ) = a b c d &CenterDot; R x m R y m + e f - - - ( 6 )
其中,Dxm表示第m个点的x方向坐标,Dym、Rxm和Rym同理。
在估计模型T(·)过程中,对于给定的N对匹配点对,由于存在测量误差和畸变未知引起的误差,点对中点坐标存在噪声,使用足够多的点对求取最优估计。根据式(6)建立方程
Ax=b(7)
A = R x 1 R y 1 1 0 0 0 0 0 0 R x 1 R y 1 1 ... ... ... ... ... ... R x M R x M 1 0 0 0 0 0 0 R x M R x M 1 - - - ( 8 )
x=[abecdf]T(9)
b=[Dx1Dy1...DxMDyM]T(10)
方程(7)有封闭形式解,因此旋转平移模型参数求解为
x=(ATA)-1ATB(11)
RANSAC算法通过在每次迭代中对包含外点的点对集合随机采样估计模型参数,得到符合该模型的候选点对集,并计算误差。最终从所有迭代中取出拟合误差最小的模型作为估计结果,相应的候选点对集合为剔除外电的集合。在所有迭代中选取最佳。
RANSAC的算法输入为:含有外点的点对集合每次采样最少点数目Q,迭代次数K,判断一个点对是否符合当前模型的误差阈值T,判断估计得到的模型是否足够好的候选集合大小D。算法流程为:迭代K次,对于第k次迭代,1≤k≤K,
(a)在中随机选取Q个点对得到候选点集{Pq}con,此时,1≤q≤Q;
(b)使用式(11)估计参数xk,得到投影模型Tk(·);
(c)对于除{Pq}con中点对外的其他点对<Rm,Dm>:计算投影误差ek,m=||Dn-Tk(Rn)||,如果ek,m≤T,则将<Rm,Dm>加入候选点对集合{Pq}con,每次更新使{Pq}con变大,使q的取值范围上限变大;
(d)如果{Pq}con元素数目大于D,则使用{Pq}con估计式(11)中参数xcon,得到投影模型Tcon(·);
(e)计算投影误差 e k = &Sigma; < R q , D q > &Element; { P q } c o n | | D q - T c o n ( R q ) | | ;
对于每次迭代,求取ek最小的Tcon(·)作为输出模型参数,同时相应的{Pn}con作为剔除外点的点集。
使用RANSAC剔除中外点得到集合估计得到相应旋转平移参数x,即之间的转换包含的旋转平移函数T(·)。
第三部分,迭代估计几何畸变。在得到不包含外点的匹配点对集合和旋转平移函数T(·)之后,中的经过T(·)变换有与相应的组合得到点对集合因此,之间仅包含集合畸变的变换,被用来估计几何畸变。使用径向畸变来描述中的几何畸变。由的几何畸变为
D x n D y n = G ( R x n &prime; R y n &prime; ) = G x ( R x n &prime; ) G y ( R y n &prime; ) = L ( r ) R x n &prime; - C x R y n &prime; - C y + C x C y - - - ( 12 )
其中,Dxn表示第n个点的x方向坐标,Dyn、R′xn和R′yn同理,Cx为畸变中心x方向坐标,Cy为畸变中心y方向坐标,G(·)表示畸变函数,表示到畸变中心C的距离,L(r)为失真函数,具有性质为L(0)=1,使用其在0处的泰勒展开表示,仅使用偶次项并舍弃高次项,因此L(r)表示为。
L(r)=1+k1r2+k2r4(13)
模型仅由畸变中心C=[Cx,Cy]T和畸变程度k1、k2确定,为方便描述,将参数表示为P=[k1,k2,Cx,Cy]。
与y方向无差别地,x方向的畸变模型Gx(·)表示为:
D x n = G x ( R x n &prime; ; P ) = ( R x n &prime; - C x ) + k 1 r 2 ( R x n &prime; - C x ) + k 2 r 4 ( R x n &prime; - C x ) r = ( R x n &prime; - C x ) 2 + ( R y n &prime; - C y ) 2 - - - ( 14 )
其中,n为点对序号。同理,Gy(·)也由此表示。
使用牛顿迭代对该非线性模型进行求解。给定初始值参数初始值P0,在迭代中求得最优解。迭代中参数更新方式为:
Pk+1=Pkk(15)
其中,k为迭代次数序号,且只有当k出现在下标与[]中时表示迭代次数,Δk为方程JkΔk=εk的最小二乘解,取值在Pk上的雅各比矩阵,而εk=D-G(R;Pk)为当前估计的误差。在每次迭代中对于每个值代入J之后均有:
J k , n ( R n ; P k ) = &part; G x ( R n x ; P ) &part; P &lsqb; k , n &rsqb; &part; G y ( R n y ; P ) &part; P &lsqb; k , n &rsqb; - - - ( 16 )
其中
&part; G x ( R x ; P ) &part; P &lsqb; k , n &rsqb; = &lsqb; &part; G x ( R x ; k 1 ) &part; k 1 &lsqb; k , n &rsqb; , &part; G x ( R x ; k 2 ) &part; k 2 &lsqb; k , n &rsqb; , &part; G x ( R x ; C x ) &part; C x &lsqb; k , n &rsqb; , &part; G x ( R x ; C y ) &part; C y &lsqb; k , n &rsqb; &rsqb; - - - ( 17 )
类似地表示。针对G(·)中的Gx(·),具体地有Gx(·)针对不同参数求偏导:
&part; G x ( R x ; P ) &part; k 1 &lsqb; k , n &rsqb; = r k , n 2 R x n &prime; &part; G x ( R x ; P ) &part; k 2 &lsqb; k , n &rsqb; = r k , n 4 R x n &prime; &part; G x ( R x ; P ) &part; C x &lsqb; k , n &rsqb; = - k 1 , k r k , n 2 - 2 k 1 , k ( R x n - C x , k ) 2 - k 2 , k r k , n 4 - 4 k 2 , k r n , k 2 ( R n x - C x , k ) 2 &part; G x ( R x ; P ) &part; C y &lsqb; k , n &rsqb; = - 2 k 1 , k ( R x n - C x , k ) ( R y n - C y , k ) 2 - 4 k 2 , k r 2 ( R n x - C x , k ) ( R n y - C y , k ) - - - ( 18 )
其中,下标中的k均表示第k次迭代中相应的值,下标中的n表示相应的量为代入第n点对之后的相应值,对于Gy(·)得到相同的求偏导结果。
基于以上,在每次迭代中求解方程JiΔi=εi,在利用到所有点对的情况下有:
Ji=[Ji(R1;Pi),Ji(R2;Pi),...,Ji(RN;Pi)]T(19)
εi=[(D1-G(R1;Pi))T,(D2-G(R2;Pi))T,...,(DN-G(RN;Pi))T]T(20)
其中,(D1-G(R1;Pi))=[Dx1-Gx(Rx1;Pi),Dy1-Gy(Ry1;Pi)]T。于是,解Δi
Δi=J+εi=(JTJ)-1JTεi(21)
其中,J+=(JTJ)-1JT
根据以上,迭代优化求解畸变模型的具体算法表示为,算法输入为输入:匹配点对集合模型G(·),参数初始值P0,迭代最大次数K,收敛误差εcon。算法迭代过程中,如果迭代次数小于最大迭代次数k<K,则进行如下迭代,对于第k次迭代:
(a)利用式(20)代入当前参数值Pk,计算εk,如果||εk||<εcon停止迭代,输出结果,否则继续下一步;
(b)将和Pi代入式(19),计算Jk
(c)利用式(21)计算得到Δk
(d)更新参数Pk+1=Pkk
停止迭代,如果停止迭代时k的值为K,否则即为满足步骤(a)中迭代停止条件时的迭代步骤中的输出结果。
通过综合以上三部分,利用联合迭代框架检测和估计旋转平移和畸变模型。整体迭代算法框架输入为:星表参考星点集合星空图像星点集合星图匹配、旋转平移估计、畸变估计所需参数,最大迭代次数K,迭代停止误差εcon。迭代算法首先初始化,令 为算法中用于传递每一步畸变估计结果的点集,令k从1到K进行循环迭代:
(a)使用改进的LCS匹配算法对进行点对匹配得到匹配点对集合如果为空集,则说明之间变换程度过大,检测到成像系统指向偏差过大或者图像畸变程度过大,退出迭代,输出结果,终止程序;
(b)将作为输入使用RANSAC方法估计旋转平移模型Tk(·),同时得到剔除匹配外点的点对集合若估计失败或过小不足以进行后续运算,则说明之间畸变程度过大,检测到成像系统指向偏差过大或者图像畸变程度过大,退出迭代,输出结果,终止程序;
(c)利用中点对的序号关系,得到对应的点对集合使用Tk(·)对{Rn}进行投影得到{P′n}={<R′n,Dn>}={<T(Rn),Dn>};
(d)将{P′n}代入,迭代估计得到畸变模型Gk(·);如果估计失败,说明图像畸变程度过大,,退出迭代,输出结果,终止程序;
(e)计算误差 &epsiv; k = 1 N &Sigma; n = 1 N | | G k ( R n &prime; ) - D n | | , 如果εkcon,则 T ^ ( &CenterDot; ) = T k ( &CenterDot; ) , G ^ ( &CenterDot; ) = G k ( &CenterDot; ) , 停止循环,输出结果,否则,进行下一步;
(f)更新 { R j * } j = 1 J , { R j * } j = 1 J = T k - 1 ( G k ( T k ( { R j } j = 1 J ) ) ) , 进入下一次循环。
如果最终循环进行了K次,则输出结果根据输出的对图像与星表相比的旋转平移和畸变作出量化评价。
本发明的有益效果是:该方法利用已知的星表、成像系统参数和成像系统指向,对图像中星点与星表中星点进行匹配,进而对图像中可能存在的相机的旋转、平移和畸变进行检测与估计。该方法充分考虑了畸变对星点匹配的影响、指向不精确和相机镜头平移旋转造成的图像中星点的偏移旋转、平移旋转与畸变综合作用造成的估计困难,利用基于畸变模型改进的最长公共子序列(LCS)星点配准方法对星空图像与星表中星点进行配准,对匹配点对集合中可能存在的错误匹配外点进行剔除并估计相对的旋转平移关系,在考虑旋转平移的情况下估计星图与星标之间的畸变。在以上过程中完成对星图中旋转平移和畸变的检测,同时通过联合迭代求解对旋转平移和畸变模型进行估计。
下面结合具体实施方式对本发明作详细说明。
具体实施方式
本发明基于星图匹配的星空图像畸变检测与估计方法具体步骤如下:
1.图像星点坐标提取。
对星空图像中的星点进行提取,并计算星点在图像上坐标。该过程包括,利用阈值对星空图像进行二值化分割,即根据像素值是否大于阈值εbw而进行分割,图像(x,y)处像素点灰度值为I(x,y),图像中像素灰度值值I(x,y)≥εbw的分割为星点光斑区域,而I(x,y)<εbw的像素点分割为背景区域。通过对星图数据进行鲁棒统计得到分割阈值
εbw=μ+δσ
其中μ为图像灰度中值,σ为图像平均绝对偏差,δ为阈值参数,本实现中取为δ=6。
在二值化图像中对二值图像进行联通区域提取,得到星点光斑区域集合{Cn};利用提取到的光斑区域,从原始图像中计算星点的亚像素坐标(xn,yn),对于每个光斑,计算方法为
x n = &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) &CenterDot; x &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) , y n = &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) &CenterDot; y &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) - - - ( 1 )
其中,(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的邻接矩阵,其中每个元素表示中两点之间在图像平面的欧氏距离,对该矩阵各行进行升序排序可以得到特征矩阵:
FMat D = 0 d 1 , 2 ... d 1 , I 0 d 2 , 2 ... d 2 , I ... ... ... ... 0 d I , 2 ... d I , I - - - ( 2 )
并且同时将FMatD中元素在各行中的原始序号记录在矩阵中得到:
IMat D = 1 p 1 , 2 ... p 1 , I 2 p 2 , 2 ... p 1 , 1 ... ... ... ... I p I , 2 ... p 1 , I - - - ( 3 )
因此元素pi,j表示FMatD中的元素di,j元素为中Di之间的距离。同样,对于星点集合也可以得到相应的矩阵FMatR和IMatR。其中FMatD和FMatR中每行为相应集合中各个星点的特征。
在匹配过程中为比较星点Dx与Ry之间的相似性,需要计算序列<0,dx,i,...,dx,I>D和序列<0,dy,j,...,dx,J>R之间的公共子序列的长度,其中x和y为星点序号。基于一般动态规划求解LCS的算法框架,我们设定其中的公共子序列更新公式为:
c &lsqb; i 1 , j 1 &rsqb; = 0 , i 1 = 0 o r j 1 = 0 c &lsqb; i 1 - 1 , j 1 - 1 &rsqb; + 1 , i 1 , j 1 > 0 a n d | | d x , i 1 - d y , j 1 | | &le; &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 > 0 a n d | | d x , i 1 - d y , j 1 | | > &epsiv; f ( d x , i 1 - d y , j 1 ) - - - ( 4 )
其中,i1和j1为一般动态规划算法表中的序号,ε为距离阈值,为[0,1]之间的函数,成立时,表示在畸变情况下而函数f(·)为:
f ( d x , i 1 , d y , j 1 ) = 1 L exp { | | C x - C | | + | | D IMat D &lsqb; x , i 1 &rsqb; - C | | + | | R y - C | | + | | R IMat R &lsqb; y , j 1 &rsqb; - C | | 4 &sigma; 2 } - - - ( 5 )
其中,C为图像畸变中心,本实施例初始假设为图像中心,在迭代中更新,||·||表示欧式距离,σ2为方差,本文中取0.5,L为归一化因子。
经过计算,得到中任两个星点之间的LCS相似性度量。对于中星点,依次在中搜索LCS最近邻,如果与最近邻之间的关系满足阈值要求,则作为与相应Di匹配的Ri,从而得到点对,否则舍弃对相应Di的匹配。本实施例设定阈值为之间最长公共子序列的长度的1/2。
经过星图匹配,得到匹配点对的集合其中Pm=<Rm,Dm>。中包含足够多的正确配对点,但仍包含匹配错误点对。
(2)基于RANSAC的外点剔除与旋转平移模型估计。
中依然存在匹配错误点对,但其中正确匹配的点对之间存在稳定的旋转平移关系,而不符合该模型的匹配错误点对为外点。使用RANSAC估计其旋转平移关系,同时剔除其中匹配错误的外点。
之间的旋转平移关系可以描述为二维空间中的旋转平移关系。由的转换关系T(·)可以描述为:
D x m D y m = T ( R x m R y m ) = a b c d &CenterDot; R x m R y m + e f - - - ( 6 )
其中,Dxm表示第m个点的x方向坐标,Dym、Rxm和Rym同理。对于给定的N对匹配点对,根据式(6)建立方程
Ax=b(7)
A = R x 1 R y 1 1 0 0 0 0 0 0 R x 1 R y 1 1 ... ... ... ... ... ... R x M R x M 1 0 0 0 0 0 0 R x M R x M 1 - - - ( 8 )
x=[abecdf]T(9)
b=[Dx1Dy1...DxMDyM]T(10)
方程(7)有封闭形式解,因此旋转平移模型参数可以求解为
x=(ATA)-1ATB(11)
RANSAC算法通过在每次迭代中对包含外点的点对集合随机采样估计模型参数,得到符合该模型的候选点对集,并计算误差。最终从所有迭代中取出拟合误差最小的模型作为估计结果,相应的候选点对集合为剔除外电的集合。在所有迭代中选取最佳。
RANSAC的算法输入为:含有外点的点对集合每次采样最少点数目Q,迭代次数K,判断一个点对是否符合当前模型的误差阈值T,判断估计得到的模型是否足够好的候选集合大小D。算法流程为:迭代K次,对于第k次迭代,1≤k≤K,
1)在中随机选取Q个点对得到候选点集{Pq}con,此时,1≤q≤Q;
2)使用式(11)估计参数xk,得到投影模型Tk(·);
3)对于除{Pq}con中点对外的其他点对<Rm,Dm>:计算投影误差ek,m=||Dn-Tk(Rn)||,如果ek,m≤T,则将<Rm,Dm>加入候选点对集合{Pq}con,每次更新使{Pq}con变大,使q的取值范围上限变大;
4)如果{Pq}con元素数目大于D,则使用{Pq}con估计式(11)中参数xcon,得到投影模型Tcon(·);
5)计算投影误差 e k = &Sigma; < R q , D q > &Element; { p q } c o n | | D q - T c o n ( R q ) | | ;
对于每次迭代,求取ek最小的Tcon(·)作为输出模型参数,同时相应的{Pn}con作为剔除外点的点集。本实施例参数选取为,N≥6,K=300,T=1,D与星点数目、图像尺度和畸变程度相关,当星点数目约为50时取为10。
使用RANSAC剔除中外点得到集合估计得到相应旋转平移参数x,即之间的转换包含的旋转平移函数T(·)。
(3)迭代估计几何畸变。
在得到不包含外点的匹配点对集合和旋转平移函数T(·)之后,中的经过T(·)变换有与相应的组合得到点对集合因此,之间仅包含集合畸变的变换,可以被用来估计几何畸变。
通常情况下,仅使用径向畸变来描述{P′n}中的几何畸变。由{R′n}到{Dn}的几何畸变可以为
D x n D y n = G ( R x n &prime; R y n &prime; ) = G x ( R x n &prime; ) G y ( R y n &prime; ) = L ( r ) R x n &prime; - C x R y n &prime; - C y + C x C y - - - ( 12 )
其中,Dxn表示第n个点的x方向坐标,Dyn、R′xn和R′yn同理,Cx为畸变中心x方向坐标,Cy为畸变中心y方向坐标,G(·)表示畸变函数,表示[R′xn,R′yn]T到畸变中心C的距离,一般情况下,畸变中心不一定与图像几何中心重叠,但通常接近集合中心。L(r)为失真函数,具有性质为L(0)=1,通常使用其在0处的泰勒展开表示,仅使用偶次项进行表示,并舍弃高次项,因此L(r)表示为
L(r)=1+k1r2+k2r4(13)
模型仅由畸变中心C=[Cx,Cy]T和畸变程度k1、k2确定,为方便描述,将参数表示为P=[k1,k2,Cx,Cy]。
与y方向无差别地,x方向的畸变模型Gx(·)可以表示为:
D x n = G x ( R x n &prime; ; P ) = ( R x n &prime; - C x ) + k 1 r 2 ( R x n &prime; - C x ) + k 2 r 4 ( R x n &prime; - C x ) r = ( R x n &prime; - C x ) 2 + ( R y n &prime; - C y ) 2 - - - ( 14 )
其中,n为点对序号。同理,Gy(·)也可以由此表示。
使用牛顿迭代对该非线性模型进行求解。牛顿迭代是一个迭代优化算法,给定初始值参数初始值P0,在迭代中求得最优解。迭代中参数更新方式为:
Pk+1=Pkk(15)
其中,k为迭代次数序号,且只有当k出现在下标与[]中时表示迭代次数,Δk为方程JkΔk=εk的最小二乘解,取值在Pk上的雅各比矩阵,而εk=D-G(R;Pk)为当前估计的误差。在每次迭代中对于每个值代入J之后均有:
J k , n ( R n ; P k ) = &part; G x ( R n x ; P ) &part; P &lsqb; k , n &rsqb; &part; G y ( R n y ; P ) &part; P &lsqb; k , n &rsqb; - - - ( 16 )
其中
&part; G x ( R x ; P ) &part; P &lsqb; k , n &rsqb; = &lsqb; &part; G x ( R x ; k 1 ) &part; k 1 &lsqb; k , n &rsqb; , &part; G x ( R x ; k 2 ) &part; k 2 &lsqb; k , n &rsqb; , &part; G x ( R x ; C x ) &part; C x &lsqb; k , n &rsqb; , &part; G x ( R x ; C y ) &part; C y &lsqb; k , n &rsqb; &rsqb; - - - ( 17 )
可以类似地表示。针对G(·)中的Gx(·),具体地有Gx(·)针对不同参数求偏导:
&part; G x ( R x ; P ) &part; k 1 &lsqb; k , n &rsqb; = r k , n 2 R x n &prime; &part; G x ( R x ; P ) &part; k 2 &lsqb; k , n &rsqb; = r k , n 4 R x n &prime; &part; G x ( R x ; P ) &part; C x &lsqb; k , n &rsqb; = - k 1 , k r k , n 2 - 2 k 1 , k ( R x n - C x , k ) 2 - k 2 , k r k , n 4 - 4 k 2 , k r n , k 2 ( R n x - C x , k ) 2 &part; G x ( R x ; P ) &part; C y &lsqb; k , n &rsqb; = - 2 k 1 , k ( R x n - C x , k ) ( R y n - C y , k ) 2 - 4 k 2 , k r 2 ( R n x - C x , k ) ( R n y - C y , k ) - - - ( 18 )
其中,下标出现的k均表示第k次迭代中相应的值,下标中的n表示相应的量为代入第n点对之后的相应值,对于Gy(·)可以得到相同的求偏导结果。
基于以上,在每次迭代中求解方程JiΔi=εi,在利用到所有点对的情况下有:
Ji=[Ji(R1;Pi),Ji(R2;Pi),...,Ji(RN;Pi)]T(19)
εi=[(D1-G(R1;Pi))T,(D2-G(R2;Pi))T,...,(DN-G(RN;Pi))T]T(20)
其中(D1-G(R1;Pi))=[Dx1-Gx(Rx1;Pi),Dy1-Gy(Ry1;Pi)]T。于是,解Δi
Δi=J+εi=(JTJ)-1JTεi(21)
根据以上,迭代优化求解畸变模型的具体算法表示为,算法输入为输入:匹配点对集合模型G(·),参数初始值P0,迭代最大次数K,收敛误差εcon。算法迭代过程中,如果迭代次数小于最大迭代次数k<K,则进行如下迭代,对于第k次迭代:
1)利用式(20)代入当前参数值Pk,计算εk,如果||εk||<εcon停止迭代,输出结果,否则继续下一步;
2)将和Pi代入式(19),计算Jk
3)利用式(21)计算得到Δk
4)更新参数Pk+1=Pkk
停止迭代,如果停止迭代时k的值为K,否则即为满足1)中迭代停止条件时的迭代步骤中的输出结果。对于以上参数,本实施例设置P0中的k1=k2=7×10-7,[Cx,Cy]初始值为图像集合中心,迭代最大次数为K=20,εcon与具体数据相关,可以默认设置为εcon=0.5。
(4)联合优化迭代框架。
通过综合以上的星点匹配、RANSAC估计旋转平移和畸变估计方法,利用联合迭代框架检测和估计旋转平移和畸变模型。整体迭代算法框架输入为星表参考星点集合星空图像星点集合星图匹配、旋转平移估计、畸变估计所需参数,最大迭代次数K,迭代停止误差εcon。迭代算法首先初始化,令 为算法中用于传递每一步畸变估计结果的点集。令k从1到K进行循环迭代:
1)使用改进的LCS匹配算法对进行点对匹配得到匹配点对集合如果为空集,则说明之间变换程度过大,检测到成像系统指向偏差过大或者图像畸变程度过大,退出迭代,输出结果,终止程序;
2)将作为输入使用RANSAC方法估计旋转平移模型Tk(·),同时得到剔除匹配外点的点对集合若估计失败或过小不足以进行后续运算,则说明之间畸变程度过大,检测到成像系统指向偏差过大或者图像畸变程度过大,退出迭代,输出结果,终止程序;
3)利用中点对的序号关系,得到对应的点对集合使用Tk(·)对{Rn}进行投影得到{P′n}={<R′n,Dn>}={<T(Rn),Dn>};
4)将{P′n}代入,迭代估计得到畸变模型Gk(·);如果估计失败,说明图像畸变程度过大,退出迭代,输出结果,终止程序;
5)计算误差 &epsiv; k = 1 N &Sigma; n = 1 N | | G k ( R n &prime; ) - D n | | , 如果εkcon,则 T ^ ( &CenterDot; ) = T k ( &CenterDot; ) , G ^ ( &CenterDot; ) = G k ( &CenterDot; ) , 停止循环,输出结果,否则,进行下一步;
6)更新同时更新匹配中使用的图像畸变中心为[Cx,Cy],进入下一次循环。
如果最终循环进行了K次,则输出结果根据输出的对图像与星表相比的旋转平移和畸变作出量化评价。参数选择与畸变程度相关,本实施例参数选择为,K=10,εcon=0.005。

Claims (1)

1.一种基于星图匹配的星空图像畸变检测与估计方法,其特征在于包括以下步骤:
步骤一、对星空图像中的星点进行提取,并计算星点在图像上坐标;该过程包括,利用阈值对星空图像进行二值化分割,即根据像素值是否大于阈值而进行分割;对二值图像进行联通区域提取,得到星点光斑区域集合Ωn;利用提取到的光斑区域,从原始图像中计算星点的亚像素坐标(xn,yn),对于每个光斑,计算方法为
x n = &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) &CenterDot; x &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) , y n = &Sigma; ( x , y ) &Element; &Omega; n I ( x , y ) &CenterDot; y &Sigma; ( x , y ) &Element; &Omega; n I ( 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的邻接矩阵,其中每个元素表示中两点之间在图像平面的欧氏距离,对该矩阵各行进行升序排序得到特征矩阵:
FMat D = 0 d 1 , 2 ... d 1 , I 0 d 2 , 2 ... d 2 , I ... ... ... ... 0 d I , 2 ... d I , I - - - ( 2 )
并且同时将FMatD中元素在各行中的原始序号记录在矩阵中得到:
IMat D = 1 p 1 , 2 ... p 1 , I 2 p 2 , 2 ... p 1 , 1 ... ... ... ... I p I , 2 ... p 1 , I - - - ( 3 )
因此,元素pi,j表示FMatD中的元素中Di之间的距离;同样,对于星点集合得到相应的矩阵FMatR和IMatR;其中FMatD和FMatR中每行为相应集合中各个星点的特征;
在匹配过程中为比较星点Dx与Ry之间的相似性,需要计算序列<0,dx,i,...,dx,I>D和序列<0,dy,j,...,dx,J>R之间的公共子序列的长度,其中x和y为星点序号;基于一般动态规划求解LCS的算法框架,设定其中的公共子序列更新公式为:
c &lsqb; i 1 , j 1 &rsqb; = 0 , i 1 = 0 o r j 1 = 0 c &lsqb; i 1 - 1 , j 1 - 1 &rsqb; + 1 , i 1 , j 1 > 0 a n d | | d x , i 1 - d y , j 1 | | &le; &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 > 0 a n d | | d x , i 1 - d y , j 1 | | > &epsiv; f ( d x , i 1 - d y , j 1 ) - - - ( 4 )
式中,i1和j1为动态规划算法表中的序号,ε为距离阈值,为[0,1]之间的函数,成立时,表示在畸变情况下而函数f(·)为:
f ( d x , i 1 , d y , j 1 ) = 1 L exp { | | D x - C | | + | | D IMat D &lsqb; x , i 1 &rsqb; - C | | + | | R y - C | | + | | R IMat R &lsqb; y , j 1 &rsqb; - C | | 4 &sigma; 2 } - - - ( 5 )
式中,C为图像畸变中心,在迭代中更新,||·||表示欧式距离,σ2为方差,L为归一化因子;
得到中任两个星点之间的LCS相似性度量;对于中星点,依次在中搜索LCS最近邻,如果与最近邻之间的关系满足阈值要求,则作为与相应Di匹配的Ri,从而得到点对,否则舍弃对相应Di的匹配;
经过星图匹配,得到匹配点对的集合其中Pm=<Rm,Dm>;中包含足够多的正确配对点,但仍包含匹配错误点对;
第二部分,基于RANSAC的外点剔除与旋转平移模型估计;中依然存在匹配错误点对,但其中正确匹配的点对之间存在稳定的旋转平移关系,而不符合该模型的匹配错误点对为外点;使用RANSAC估计其旋转平移关系,同时剔除其中匹配错误的外点;
之间的旋转平移关系描述为二维空间中的旋转平移关系;由的转换关系T(·)描述为:
D x m D y m = T ( R x m R y m ) = a b c d &CenterDot; R x m R y m + e f - - - ( 6 )
其中,Dxm表示第m个点的x方向坐标,Dym、Rxm和Rym同理;
在估计模型T(·)过程中,对于给定的N对匹配点对,由于存在测量误差和畸变未知引起的误差,点对中点坐标存在噪声,使用足够多的点对求取最优估计;根据式(6)建立方程
Ax=b(7)
A = R x 1 R y 1 1 0 0 0 0 0 0 R x 1 R y 1 1 ... ... ... ... ... ... R x M R x M 1 0 0 0 0 0 0 R x M R x M 1 - - - ( 8 )
x=[abecdf]T(9)
b=[Dx1Dy1...DxMDyM]T(10)
方程(7)有封闭形式解,因此旋转平移模型参数求解为
x=(ATA)-1ATB(11)
RANSAC算法通过在每次迭代中对包含外点的点对集合随机采样估计模型参数,得到符合该模型的候选点对集,并计算误差;最终从所有迭代中取出拟合误差最小的模型作为估计结果,相应的候选点对集合为剔除外电的集合;在所有迭代中选取最佳;
RANSAC的算法输入为:含有外点的点对集合每次采样最少点数目Q,迭代次数K,判断一个点对是否符合当前模型的误差阈值T,判断估计得到的模型是否足够好的候选集合大小D;算法流程为:迭代K次,对于第k次迭代,1≤k≤K,
(e)在中随机选取Q个点对得到候选点集{Pq}con,此时,1≤q≤Q;
(f)使用式(11)估计参数xk,得到投影模型Tk(·);
(g)对于除{Pq}con中点对外的其他点对<Rm,Dm>:计算投影误差ek,m=||Dn-Tk(Rn)||,如果ek,m≤T,则将<Rm,Dm>加入候选点对集合{Pq}con,每次更新使{Pq}con变大,使q的取值范围上限变大;
(h)如果{Pq}con元素数目大于D,则使用{Pq}con估计式(11)中参数xcon,得到投影模型Tcon(·);
(i)计算投影误差 e k = &Sigma; < R q , D q > &Element; { P q } c o n | | D q - T c o n ( R q ) | | ;
对于每次迭代,求取ek最小的Tcon(·)作为输出模型参数,同时相应的{Pn}con作为剔除外点的点集;
使用RANSAC剔除中外点得到集合估计得到相应旋转平移参数x,即之间的转换包含的旋转平移函数T(·);
第三部分,迭代估计几何畸变;在得到不包含外点的匹配点对集合和旋转平移函数T(·)之后,中的经过T(·)变换有与相应的组合得到点对集合因此,之间仅包含集合畸变的变换,被用来估计几何畸变;使用径向畸变来描述中的几何畸变;由的几何畸变为
D x n D y n = G ( R x n &prime; R y n &prime; ) = G x ( R x n &prime; ) G y ( R y n &prime; ) = L ( r ) R x n &prime; - C x R y n &prime; - C y + C x C y - - - ( 12 )
其中,Dxn表示第n个点的x方向坐标,Dyn、R′xn和R′yn同理,Cx为畸变中心x方向坐标,Cy为畸变中心y方向坐标,G(·)表示畸变函数,表示[R′xn,R′yn]T到畸变中心C的距离,L(r)为失真函数,具有性质为L(0)=1,使用其在0处的泰勒展开表示,仅使用偶次项并舍弃高次项,因此L(r)表示为;
L(r)=1+k1r2+k2r4(13)
模型仅由畸变中心C=[Cx,Cy]T和畸变程度k1、k2确定,为方便描述,将参数表示为P=[k1,k2,Cx,Cy];
与y方向无差别地,x方向的畸变模型Gx(·)表示为:
Dxn=Gx(R′xn;P)=(R′xn-Cx)+k1r2(R′xn-Cx)+k2r4(R′xn-Cx)
(14)
r = ( R x n &prime; - C x ) 2 + ( R y n &prime; - C y ) 2
其中,n为点对序号;同理,Gy(·)也由此表示;
使用牛顿迭代对畸变模型Gx(·)进行求解;给定初始值参数初始值P0,在迭代中求得最优解;迭代中参数更新方式为:
Pk+1=Pk+△k(15)
其中,k为迭代次数序号,且只有当k出现在下标与[]中时表示迭代次数,△k为方程Jkk=εk的最小二乘解,取值在Pk上的雅各比矩阵,而εk=D-G(R;Pk)为当前估计的误差;在每次迭代中对于每个值代入J之后均有:
J k , n ( R n ; P k ) = &part; G x ( R n x ; P ) &part; P &lsqb; k , n &rsqb; &part; G y ( R n y ; P ) &part; P &lsqb; k , n &rsqb; - - - ( 16 )
其中
&part; G x ( R x ; P ) &part; P &lsqb; k , n &rsqb; = &lsqb; &part; G x ( R x ; k 1 ) &part; k 1 &lsqb; k , n &rsqb; , &part; G x ( R x ; k 2 ) &part; k 2 &lsqb; k , n &rsqb; , &part; G x ( R x ; k x ) &part; C x &lsqb; k , n &rsqb; , &part; G x ( R x ; k y ) &part; C y &lsqb; k , n &rsqb; &rsqb; - - - ( 17 )
类似地表示;针对G(·)中的Gx(·),具体地有Gx(·)针对不同参数求偏导:
&part; G x ( R x ; P ) &part; k 1 &lsqb; k , n &rsqb; = r k , n 2 R x n &prime;
&part; G x ( R x ; P ) &part; k 2 &lsqb; k , n &rsqb; = r k , n 4 R x n &prime;
(18)
&part; G x ( R x ; P ) &part; C x &lsqb; k , n &rsqb; = - k 1 , k r k , n 2 - 2 k 1 , k ( R x n - C x , k ) 2 - k 2 , k r k , n 4 - 4 k 2 , k r n , k 2 ( R n x - C x , k ) 2
&part; G x ( R x ; P ) &part; C y &lsqb; k , n &rsqb; = - 2 k 1 , k ( R x n - C x , k ) ( R y n - C y , k ) - 4 k 2 , k r 2 ( R n x - C x , k ) ( R n y - C y , k )
其中,下标中的k均表示第k次迭代中相应的值,下标中的n表示相应的量为代入第n点对之后的相应值,对于Gy(·)得到相同的求偏导结果;
基于以上,在每次迭代中求解方程Jii=εi,在利用到所有点对的情况下有:
Ji=[Ji(R1;Pi),Ji(R2;Pi),...,Ji(RN;Pi)]T(19)
εi=[(D1-G(R1;Pi))T,(D2-G(R2;Pi))T,...,(DN-G(RN;Pi))T]T(20)
其中,(D1-G(R1;Pi))=[Dx1-Gx(Rx1;Pi),Dy1-Gy(Ry1;Pi)]T;于是,解△i
i=J+εi=(JTJ)-1JTεi(21)其中,J+=(JTJ)-1JT
根据以上,迭代优化求解畸变模型的具体算法表示为,算法输入为输入:匹配点对集合模型G(·),参数初始值P0,迭代最大次数K,收敛误差εcon;算法迭代过程中,如果迭代次数小于最大迭代次数k<K,则进行如下迭代,对于第k次迭代:
(j)利用式(20)代入当前参数值Pk,计算εk,如果||εk||<εcon停止迭代,输出结果,否则继续下一步;
(k)将和Pi代入式(19),计算Jk
(l)利用式(21)计算得到△k
(m)更新参数Pk+1=Pk+△k
停止迭代,如果停止迭代时k的值为K,否则即为满足步骤(j)中迭代停止条件时的迭代步骤中的输出结果;
通过综合以上三部分,利用联合迭代框架检测和估计旋转平移和畸变模型;整体迭代算法框架输入为:星表参考星点集合星空图像星点集合星图匹配、旋转平移估计、畸变估计所需参数,最大迭代次数K,迭代停止误差εcon;迭代算法首先初始化,令为算法中用于传递每一步畸变估计结果的点集,令k从1到K进行循环迭代:
(n)使用改进的LCS匹配算法对进行点对匹配得到匹配点对集合如果为空集,则说明之间变换程度过大,检测到成像系统指向偏差过大或者图像畸变程度过大,退出迭代,输出结果,终止程序;
(o)将作为输入使用RANSAC方法估计旋转平移模型Tk(·),同时得到剔除匹配外点的点对集合若估计失败或过小不足以进行后续运算,则说明之间畸变程度过大,检测到成像系统指向偏差过大或者图像畸变程度过大,退出迭代,输出结果,终止程序;
(p)利用中点对的序号关系,得到对应的点对集合使用Tk(·)对{Rn}进行投影得到{P′n}={<R′n,Dn>}={<T(Rn),Dn>};
(q)将{P′n}代入,迭代估计得到畸变模型Gk(·);如果估计失败,说明图像畸变程度过大,退出迭代,输出结果,终止程序;
(r)计算误差 &epsiv; k = 1 N &Sigma; n = 1 N | | G k ( R n &prime; ) - D n | | , 如果εk<εcon,则 T ^ ( &CenterDot; ) = T k ( &CenterDot; ) , G ^ ( &CenterDot; ) = G k ( &CenterDot; ) , 停止循环,输出结果,否则,进行下一步;
(s)更新 { R j * } j = 1 J = T k - 1 ( G k ( T k ( { R j } j = 1 J ) ) ) , 进入下一次循环;
如果最终循环进行了K次,则输出结果根据输出的T对图像与星表相比的旋转平移和畸变作出量化评价。
CN201310390034.1A 2013-08-30 2013-08-30 基于星图匹配的星空图像畸变检测与估计方法 Active CN103440659B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310390034.1A CN103440659B (zh) 2013-08-30 2013-08-30 基于星图匹配的星空图像畸变检测与估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310390034.1A CN103440659B (zh) 2013-08-30 2013-08-30 基于星图匹配的星空图像畸变检测与估计方法

Publications (2)

Publication Number Publication Date
CN103440659A CN103440659A (zh) 2013-12-11
CN103440659B true CN103440659B (zh) 2016-04-13

Family

ID=49694352

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310390034.1A Active CN103440659B (zh) 2013-08-30 2013-08-30 基于星图匹配的星空图像畸变检测与估计方法

Country Status (1)

Country Link
CN (1) CN103440659B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105631870B (zh) * 2015-12-25 2018-08-24 北京理工大学 基于NoC架构的星图配准与目标轨迹提取方法及装置
CN107341824B (zh) * 2017-06-12 2020-07-28 西安电子科技大学 一种图像配准的综合评价指标生成方法
CN107506772A (zh) * 2017-08-28 2017-12-22 西北工业大学 基于自适应排序的快速鲁棒模型估计方法
CN107609547B (zh) * 2017-09-06 2021-02-19 其峰科技有限公司 星体快速识别方法、装置及望远镜
CN109917974B (zh) * 2019-03-20 2022-03-22 安徽慧视金瞳科技有限公司 一种交互式投影系统非线性点坐标映射方法
CN110930329B (zh) * 2019-11-20 2023-04-21 维沃移动通信有限公司 星空图像处理方法及装置
CN113034394B (zh) * 2021-03-25 2022-09-06 中国科学院紫金山天文台 一种基于恒星星表的望远镜畸变校正方法
CN114820738B (zh) * 2022-06-30 2022-09-23 中国人民解放军国防科技大学 一种星图精确配准方法、装置、计算机设备和存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101246590A (zh) * 2008-03-03 2008-08-20 北京航空航天大学 星载相机空间畸变图像几何校正方法
WO2011003735A1 (de) * 2009-07-08 2011-01-13 Robert Bosch Gmbh Verzeichnungskorrektur von videosystemen
WO2012063467A1 (ja) * 2010-11-11 2012-05-18 パナソニック株式会社 画像処理装置、画像処理方法、プログラム及び撮影装置
CN103069433A (zh) * 2010-08-26 2013-04-24 索尼公司 具有图像对准机构的图像处理系统及其操作方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101246590A (zh) * 2008-03-03 2008-08-20 北京航空航天大学 星载相机空间畸变图像几何校正方法
WO2011003735A1 (de) * 2009-07-08 2011-01-13 Robert Bosch Gmbh Verzeichnungskorrektur von videosystemen
CN103069433A (zh) * 2010-08-26 2013-04-24 索尼公司 具有图像对准机构的图像处理系统及其操作方法
WO2012063467A1 (ja) * 2010-11-11 2012-05-18 パナソニック株式会社 画像処理装置、画像処理方法、プログラム及び撮影装置

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A method of optimization for the distorted model of star map based on improved genetic algorithm;Quan Wei 等;《Aerospace Science and Technology》;20100618;103-107 *
一种帧转移型CCD传感器Smear效应消除方法;高建伟 等;《中国体视学与图像分析》;20101031;第15卷(第4期);394-398 *
一种星空背景的运动补偿方法;孙瑾秋 等;《机械科学与技术》;20081231;第27卷(第12期);1563-1566 *
星图中基于小波变换的CCD传感器Smear现象消除方法;姚睿 等;《光子学报》;20110331;第40卷(第3期);413-418 *

Also Published As

Publication number Publication date
CN103440659A (zh) 2013-12-11

Similar Documents

Publication Publication Date Title
CN103440659B (zh) 基于星图匹配的星空图像畸变检测与估计方法
CN107564062B (zh) 位姿异常检测方法及装置
CN104484648B (zh) 基于轮廓识别的机器人可变视角障碍物检测方法
Mnih et al. Learning to label aerial images from noisy data
CN106408591B (zh) 一种抗遮挡的目标跟踪方法
CN107563438A (zh) 一种快速鲁棒的多模态遥感影像匹配方法和系统
CN103854283A (zh) 一种基于在线学习的移动增强现实跟踪注册方法
JP6397379B2 (ja) 変化領域検出装置、方法、及びプログラム
CN103048331B (zh) 基于柔性模板配准的印刷缺陷检测方法
CN107424161B (zh) 一种由粗至精的室内场景图像布局估计方法
CN104881029B (zh) 基于一点ransac和fast算法的移动机器人导航方法
CN106485751B (zh) 应用于基桩检测中的无人机摄影成像及数据处理方法及系统
CN113361542B (zh) 一种基于深度学习的局部特征提取方法
CN111145227B (zh) 一种地下隧道空间多视点云的可迭代整体配准方法
CN110197505B (zh) 基于深度网络及语义信息的遥感图像双目立体匹配方法
CN107229920B (zh) 基于整合深度典型时间规整及相关修正的行为识别方法
CN101398934A (zh) 对图像中的对象进行定位的方法和系统
CN105205818A (zh) 一种电气设备红外图像和可见光图像配准的方法
CN101692284A (zh) 基于量子免疫克隆算法的三维人体运动跟踪方法
Dal Poz et al. Object-space road extraction in rural areas using stereoscopic aerial images
Wagstaff et al. Self-supervised scale recovery for monocular depth and egomotion estimation
Jung et al. A line-based progressive refinement of 3D rooftop models using airborne LiDAR data with single view imagery
CN110764504A (zh) 一种用于变电站电缆沟道巡检的机器人导航方法及系统
CN104217209A (zh) 一种遥感图像自动配准点错误匹配消除方法
Zhang et al. Feature matching for multi-epoch historical aerial images

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant