CN101246590A - 星载相机空间畸变图像几何校正方法 - Google Patents

Abstract

一种基于星图匹配的星载相机空间图像畸变校正方法，将图像上的恒星点作为参考点，实现畸变图像的校正，其步骤为：(1)图像分割，提取图像中恒星光斑部分；(2)通过重心法计算各光斑的重心坐标，作为实际坐标；(3)利用相机在地面标定时得到的畸变参数进行预校正；(4)利用单帧或多帧Hausdorff距离星图匹配方法进行星图匹配，计算各实际坐标对应的理论坐标；(5)利用实际坐标与理论坐标对应关系，计算图像的畸变参数；(6)通过三阶多项式拟合及插值法得到校正后图像。另外，本发明针对相机光轴指向系统误差对校正的影响，给出一种自适应补偿方法。本发明方法运算量小，图像校正效果明显，并能有效抑制图像噪声及输入系统误差给校正所带来的影响。

Description

星载相机空间畸变图像几何校正方法

技术领域

本发明主要涉及到图像处理领域中一种畸变图像几何校正方法，针对星载相机所拍摄的空间图像的畸变情况。

背景技术

由于相机实际镜头组成像过程难以达到理想透镜成像的效果，在像方垂直主轴平面不同点的切向放大率一般不完全相同，造成相机畸变的产生。一般认为，相机的畸变主要为径向畸变和切向畸变，其中径向畸变由于镜头各透镜曲率存在误差，而切向畸变则归咎于各镜头的光学主轴没有重合。

星载CCD相机虽然在地面可以进行精确的标定，但由于其焦距长、分辨率高、作用距离远，通光口径较大，光学零件尺寸比较大等特点，容易受到空间环境和相机内部环境的影响。卫星内部发热以及从外部吸收的热辐射会使光学元件发生变形。另外，卫星发射时的惯性力和冲击振动也会给光学元件带来影响。所以，星载相机所拍摄的空间图像会产生畸变，并且畸变模型、参数仍有可能发生变化。若不考虑这些因素，会给后续的图像恢复、目标定位等图像处理工作带来不利的影响。所以对于星载相机所拍摄的图像必须首先进行校正处理。

目前常用的畸变参数的测定方法，如相机非线性标定方法通过建立相机小孔透视模型，利用反复迭代实现相机的标定和畸变校正；等效曲面法将畸变的等效曲面设置为球面，通过求取球面的中心和半径实现图像的径向畸变校正；神经网络法在无需计算相机内外参数的场合利用多层前馈神经元网络实现图像非线性畸变的补偿。以上方法均为采用将相机固定在高精度的转台上，通过拍摄与相机垂直放置的标准模板的方法确定，其精度也比较高。但空间相机在太空中不可能有可供拍摄的模版，所以以上方法不适用。根据曲率变化检校畸变差的方法通过检测图像中的共线“直线”的曲率，确定图像的畸变程度。但要求被拍摄目标存在丰富的直线。由于空间图像背景单一，可观测物体少，所以该方法亦不适用。基于主动视觉的相机自标定技术是另一个研究的热点，其优点是仅利用摄像机在运动过程中周围环境的图像与图像之间的对应关系对摄像机进行标定。但由于卫星的运动受其轨道限制，不可能任意更改，同时星载相机所拍摄的图像中甚至不存在固定的可供拍摄的三维物体，所以基于主动视觉或kruppa方程求解的自标定技术对于星载相机也不适用。

因此，针对星载相机拍摄图像的特点，开发一种能够直接利用星载相机传回的图像对图像本身进行畸变参数计算和畸变校正的方法变得尤为重要。

发明内容

本发明的目的在于：克服现有技术的不足，提供一种星载相机空间畸变图像几何校正方法，该方法将恒星点作为特征点，具有计算量小，校正精度高，校正效果明显的优点，并能有效抑制图像噪声及输入系统误差给校正所带来的影响。

本发明的技术解决方案：一种基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于步骤如下：

(1)读取的图像进行分割处理，提取图像中恒星光斑部分；

(2)通过计算各光斑的重心坐标，作为实际坐标；

(3)利用相机在地面标定时得到的畸变参数，对步骤(2)得到的实际坐标进行预校正，使之更接近无畸变情况下的坐标，从而提高星图匹配的成功率；

(4)利用单帧或多帧Hausdorff距离星图匹配方法将星载相机所拍摄的空间图像进行星图匹配，计算各实际坐标对应的理论坐标；

(5)利用实际坐标与理论坐标对应关系，计算图像的畸变参数；

(6)利用插值法得到校正后图像。

本发明与现有技术相比的优点在于：本发明由于利用相机在太空中拍摄的图像中的恒星信息进行畸变参数计算，所以不再使用传统的模板拍摄方式。所计算的相机参数亦是相机在工作环境和时间下的参数。相比于地面测算参数后再发射进入太空的相机参数测算与工作时空分开的方法，本方法避免了卫星发射升空时振动和太空环境对相机参数的影响，因而校正精度更高。另外，在相机在太空长时间工作后，其参数可能发生漂移，本方法可以方便相机参数的重新测定，从而也避免了参数漂移的影响。

附图说明

图1为本发明方法的整体流程图；

图2为视场角为2°的星空卫星图像；

图3为6幅有畸变的卫星图像；

图4为6幅经过畸变校正的卫星图像；

图5为光轴指向各有1′误差时图像X、Y轴的偏差，其中(a)为光轴指向误差导致X轴的偏差，(b)为光轴指向误差导致Y轴的偏差；

图6为未考虑误差补偿校正图像；

图7为考虑误差补偿校正图像；

图8为校正前后误差对比仿真图像；

图9(a)、(b)、(c)分别为多帧图像中第一帧的无畸变、存在畸变以及校正后的图像，为了便于观察，图像中分别加入了与自身畸变程度相同的网格，并且星点用方格标出。

具体实施方式

如图1所示，本发明的步骤如下：

1.通过图像分割提取恒星光斑

首先读入足够帧的星空背景图像，保证图像中有10颗以上恒星，其次通过图像分割将恒星光斑部分分离出来。在本发明中利用改进型大津阈值的方法进行图像分割，其具体流程为：

(1)令T_min为图像灰度均值，即令图像大小为M×N，图像上(i，j)点的灰度值为G(i，j)，则 $T_{\min }=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}G(i,j)}{M\×N}.$

(2)在灰度级从T_min～255区间中利用大津阈值选取最优阈值T，步骤为：

a)令p(k)为灰度级k的频率，即 $p\left(k\right)=\frac{\underset{G(i,j)=k}{\mathrm{\Σ}}1}{M\×N},$ k∈[0，255]。

b)若t为分割阈值，将图像分割为两部分：{G(i，j)|G(i，j)≤t}和{G(i，j)|G(i，j)＞t}。ω₀和ω₁分别为目标和背景部分比例。

${\mathrm{\&omega;}}_0=\underset{t<i\&le;255}{\mathrm{\&Sigma;}}p\left(i\right)$

${\mathrm{\&omega;}}_1=\underset{0\&le;i\&le;t}{\mathrm{\&Sigma;}}p\left(i\right)$

μ₀和μ₁为目标和背景部分的灰度均值。

${\mathrm{\&mu;}}_0=\frac{\underset{t<i\&le;255}{\mathrm{\&Sigma;}}\mathrm{ip}\left(i\right)}{{\mathrm{\&omega;}}_1\left(t\right)}$

${\mathrm{\&mu;}}_1=\frac{\underset{0\&le;i\&le;t}{\mathrm{\&Sigma;}}\mathrm{ip}\left(i\right)}{{\mathrm{\&omega;}}_0\left(t\right)}$

总均值               μ＝ω₀(t)μ₀(t)+ω₁(t)μ₁(t)

a)最佳阈值T可由式(1)计算得到：

$T=\mathrm{Arg}\underset{T_{\min }\&le;t\&le;255}{\mathrm{Max}}[{\mathrm{\&omega;}}_0\left(t\right){({\mathrm{\&mu;}}_0\left(t\right)-\mathrm{\&mu;})}^2+{\mathrm{\&omega;}}_1\left(t\right){({\mathrm{\&mu;}}_1\left(t\right)-\mathrm{\&mu;})}^2]---\left(1\right)$

(3)令μ₀′和μ₁′为利用最佳阈值分割后图像目标和背景的灰度均值，σ₀和σ₁为目标和背景灰度标准差。考察式(2)是否成立，若成立则T即为所求的最佳阈值，否则令T_min＝T_min+1并返回步骤(2)。

μ₀′-μ₁′＞a(σ₀+σ₁)其中取a＝3～5                      (2)

实验表明，该方法不但可以准确地将图像中星点部分分割出来，同时在一定程度上能够避免由于整张图像光照不均匀造成的分割错误。

使检测出的认为是恒星部分的像素灰度保持不变，背景灰度全部置为零，即完成图像分割的工作。

2.利用重心法提取恒星实际坐标

利用八连通域标记的方法计算图中每颗星的面积及位置。首先从内存中申请与图像矩阵G同样大小的标记矩阵S，用以标记连通域。设通过上一步改进大津阈值方法确定的阈值为T，计算各星点重心的步骤为：

(1)初始化，矩阵S＝[0]，sign＝1用于记录连通域当前标号；SAME[N][2]＝[0]

用以记录占用多个标号的同一连通域；其中N根据图像情况可取某较大值；n＝0配合SAME矩阵使用。从图像矩阵左上角开始逐行扫描。

(2)若某像素点灰度值G(i，j)大于阈值T，进行如下操作：

a)if该点为图像左上角点，即i＝0；j＝0，then S(i，j)＝sign；endif

b)else if该点为第一行中的点，即i＝0，then：

i.if S(i，j-1)＝0 then S(i，j)＝sign sign++；endif

ii.if S(i，j-1)≠0 then S(i，j)＝S(i，j-1)；endif

endif

c)else if该点为第一列中的点，即j＝0，then：

iii.if S(i-1，j)＝0 and S(i-1，j+1)＝0 then S(i，j)＝sign sign++；

endif

iv.if S(i-1，j)＝0 and S(i-1，j+1)≠0 then S(i，j)＝S(i-1，j+1)；

endif

v.(if S(i-1，j)≠0 and S(i-1，j+1)＝0)

or(if S(i-1，j)≠0 and S(i-1，j+1)≠0 and S(i-1，j)＝S(i-1，j+1))then

            S(i，j)＝S(i-1，j)；endif

               vi.if S(i-1，j)≠0 and S(i-1，j+1)≠0 and S(i-1，j)≠S(i-1，j+1)

               then

                      (i).S(i，j)＝S(i-1，j)

                      (ii).SAME[n][0]＝S(i-1，j)；SAME[n][1]＝S(i-1，j+1)；

                      (iii).n++；

                   endif

                endif

         d)   else

            vii.if S(i，j-1)＝S(i-1，j-1)＝S(i-1，j)＝S(i-1，j+1)＝0

            then S(i，j)＝sign sign++；endif

            viii.if S(i，j-1)、S(i-1，j-1)、S(i-1，j)、S(i-1，j+1)中有且仅有一个

            S(p，q)不为0，then S(i，j)＝S(p，q)；endif

            ix.if S(i，j-1)、S(i-1，j-1)、S(i-1，j)、S(i-1，j+1)中有多个不为

        0，但相等，如S(p，q)＝S(p′，q′)≠0 then S(i，j)＝S(p，q)；endif

            x.if S(i，j-1)、S(i-1，j-1)、S(i-1，j)、S(i-1，j+1)中有多个不为

        0，且不相等，如S(p，q)≠S(p′，q′)≠0 then

                (i).S(i，j)＝S(p，q)；

                (ii).SAME[n][0]＝S(p，q)；SAME[n][1]＝S(p′，q′)；

                (iii).n++；

             endif

(3)if(i，j)不是图像最后一行最后一列像素点，then goto步骤(2)；endif

(4)for k＝0：n-1；

遍历S，if S(i，j)＝SAME[k][1]then S(i，j)＝SAME[k][0]；

end

(5)若矩阵S中某两点S(p，q)＝S(p′，q′)，则在图像G中相对应的两点G(p，q)、G(p′，q′)为同一星点光斑中的两像素点。通过矩阵S求G中各光斑的区域Ω及像素面积大小。为了进一步去除孤立噪声点和卫星等其它目标，可以设定恒星光斑面积大小双阈值，当连通成分不在此范围内时，便认为不是恒星目标而舍弃。

(6)将保留的恒星区域利用式(3)的重心法计算其在图像中的精确坐标，即每颗恒星的实际坐标。式(3)

为所求星点的重心坐标。

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{x}=\frac{{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{xG}(x,y)\mathrm{dx}}{{\&Integral;}_{\mathrm{\&Omega;}}G(x,y)\mathrm{dx}}\\ \stackrel{\&OverBar;}{y}=\frac{{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{yG}(x,y)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;}}G(x,y)\mathrm{dy}}\end{array}\right.---\left(3\right)$

3.对实际坐标进行预校正

由于HD方法对图像的畸变较为敏感，为了保证星图匹配的成功率，在进行匹配之前要对上一步计算得到的图像中各星点的实际坐标利用相机在地面测定时的畸变模型和参数进行畸变校正，称之为预校正，使之更接近无畸变的坐标。同时记录预校正后的坐标与实际坐标的对应关系。

设通过相机在地面标定的结果，图像某点存在畸变坐标(x_d，y_d)与准确坐标(x_u，y_u)之间的关系为：

$\left\{\begin{array}{c}{x}_u=x_d+k_1{x}_d(x_d^2+y_d^2)+p_1{x}_d({3x}_d^2+y_d^2)+2p_2{x}_d{y}_d\\ y_u=y_d+k_2{y}_d(x_d^2+y_d^2)+2p_1{x}_d{y}_d+p_2\mathrm{yd}(x_d^2+3y_d^2)\end{array}\right.---\left(4\right)$

式(4)中各参数已测知。

则通过步骤2所得到的各星点的实际坐标(X_d，Y_d)可利用式(5)进行预校正，得到预校正后坐标

$\left\{\begin{array}{c}{\tilde{X}}_u=X_d+k_1{X}_d(X_d^2+Y_d^2)+p_1{X}_d({3X}_d^2+Y_d^2)+2p_2{X}_d{Y}_d\\ {\tilde{Y}}_u=Y_d+k_2{Y}_d(X_d^2+Y_d^2)+2p_1{X}_d{Y}_d+p_2\mathrm{Yd}(X_d^2+3Y_d^2)\end{array}\right.---\left(5\right)$

4.单帧或多帧HD星图匹配

基于Hausdorff距离(HD)星图匹配方法是一种较为成熟的航天器自主导航技术，其基本原理是采用相对坐标的有向最小距离方法对星敏感器图像中的点和星表中的星点进行匹配，从而将二者建立一一对应的关系。HD星图匹配方法最突出的特点是抗噪能力强，即使图像中星点被遮掩或退化，该方法仍能够得到较好的匹配结果。

恒星星表是记载恒星各种参数(如位置、运动、星等、光谱型等)的表册。根据相机参数(光轴指向、旋转角、视场大小等)搜索星表，找到理论上应该落入视场的恒星，再利用式(6)将恒星从天球坐标转换为图像坐标，即得到理想图像控制点坐标，称之为理论坐标。式(6)中(X，Y)是恒星图像坐标；(N_x，N_y)是图像横纵方向分辨率；(FOV_x×FOV_y)是相机视场角；(α₀，δ₀)是相机主轴指向；(α，δ)是恒星天球坐标。

考虑相机旋转的情况，可将通过式(6)得到的X、Y坐标值代入式(7)中，便可以计算出相机沿视轴旋转ψ角度时星点的理论坐标。

$\left\{\begin{array}{c}X=\frac{N_x}{2}\frac{1}{\tan ({\mathrm{FOV}}_x/2)}\frac{\cos \mathrm{\&delta;}\sin \left({\mathrm{\&alpha;}-\mathrm{\&alpha;}}_0\right)}{\sin \mathrm{\&delta;}\sin {\mathrm{\&delta;}}_0+\cos \mathrm{\&delta;}\cos {\mathrm{\&delta;}}_0\cos (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)}\\ Y=\frac{N_y}{2}\frac{1}{\tan ({\mathrm{FOV}}_y/2)}\frac{\sin \mathrm{\&delta;}\cos {\mathrm{\&delta;}}_0-\cos {\mathrm{\&delta;}\sin \mathrm{\&delta;}}_0\cos (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)}{\sin \mathrm{\&delta;}\sin {\mathrm{\&delta;}}_0+\cos \mathrm{\&delta;}\cos {\mathrm{\&delta;}}_0\cos (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)}\end{array}\right.---\left(6\right)$

$\left\{\begin{array}{c}{X}^{\&prime;}=X\cos \mathrm{\&psi;}+Y\sin \mathrm{\&psi;}\\ Y^{\&prime;}=-X\sin \mathrm{\&psi;}+Y\cos \mathrm{\&psi;}\end{array}\right.---\left(7\right)$

当相机视场角FOV比较小，其所拍摄的空间图像中能观测到的星点非常少。若相机能敏感到7等星，当FOV_x和FOV_y均为2°时，每帧图像一般仅能拍摄到一两颗恒星。在这种情况下，一般的星图匹配方法无法使用。在星图匹配时不能简单地将各星点的理论坐标与距离最近星点的实际坐标联系起来作为匹配结果。因为图像存在畸变和可能错检测为星点的噪声、卫星等干扰，以及输入参数的误差，简单的利用单张图像中少量星点进行星图匹配很可能会造成匹配错误，而利用多帧图像进行星图匹配，可以很好地解决这一难题。

HD经常利用点的距离描述两幅图的相似程度，这里通过点间距描述两个点的相似程度，假设A＝{A₁，A₂，…，A_p}和B＝{B₁，B₂，…，B_q}分别是理想坐标数组和实际坐标数组。首先，对A中点A_i用其到A中其它点的距离表征。式(8)中da_ik为A_i到点A_k的欧氏距离。

A_i＝{da_ik|k＝1，2，…，p，k≠i}(i＝1，2，…，p)                    (8)

同样可得数组B中点的表征方法。再按式(9)、式(10)产生一个p×q的矩阵H表示A中点A_i和B中点B_j的HD。最后，搜索H找出HD较小的值，其对应A和B中的点即为控制点对。

$\left\{\begin{array}{c}h(A_i,B_j)=\underset{a_{\mathrm{ik}}\&Element;A_i}{\max }\underset{b_{\mathrm{jl}}\&Element;B_j}{\min }\left|\right|a_{\mathrm{ik}}-b_{\mathrm{jl}}\left|\right|\\ h(B_j,A_i)=\underset{b_{\mathrm{jl}}\&Element;B_j}{\max }\underset{a_{\mathrm{ik}}\&Element;A_i}{\min }\left|\right|b_{\mathrm{jl}}-a_{\mathrm{ik}}\left|\right|\end{array}\right.---\left(9\right)$

H(i，j)＝max{h(A_i，B_j)，h(B_j，A_i)}                            (10)

对多帧图像分别求取图中所有星点的实际坐标及与之相对应的预校正后坐标，将所得到的一组坐标不加区分地认为是从同一张图中检测得出的。然后按照每帧图像的视轴指向及视场角搜索星表列出可能落入每帧图像中的恒星，计算每颗星的理论坐标，记为理论坐标数组。将所得到的预校正后坐标数组和理论坐标数组按照HD方法进行匹配，去除匹配失败点，便得到一组星点的预校正后坐标和对应的理论坐标，进而可以得到各实际坐标(X_d，Y_d)及与之一一对应的理论坐标(X_u，Y_u)。

5.畸变参数计算

将所得到的多帧图像中星点的一组实际坐标及与之对应的理论坐标代入式(11)，利用最小二乘法求得校正参数a₀～a₉、b₀～b₉的值。在式(11)中，以图像中心为坐标原点，(X_u，Y_u)为校正后图像某点的坐标，(X_d，Y_d)为校正前图像某点的坐标。由于星载相机造成图像畸变的成因复杂，故采用三阶多项式拟合的方法进行校正，可以避免畸变模型不准确造成的校正误差，但需要的星点较多，因而提高了对输入图像的要求。

$\left\{\begin{array}{c}{X}_d=a_0+a_1{X}_u+a_2{Y}_u+a_3{X}_u^2+a_4{X}_u{Y}_u+a_5{Y}_u^2+a_6{X}_u^3+a_7{X}_u^2{Y}_u+a_8{X}_u{Y}_u^2+a_9{Y}_u^3\\ Y_d=b_0+b_1{X}_u+b_2{Y}_u+b_3{X}_u^2+b_4{X}_u{Y}_u+b_5{Y}_u^2+b_6{X}_u^3+b_7{X}_u^2{Y}_u+b_8{X}_u{Y}_u^2+b_9{Y}_u^3\end{array}\right.---\left(11\right)$

设单帧或多帧图像中恒星点的理论坐标数组和实际坐标数组分别为： ${\left(\begin{array}{cccc}{X}_{u1}& X_{u2}& ...& X_{\mathrm{uN}}\\ Y_{u1}& Y_{u2}& ...& Y_{\mathrm{uN}}\end{array}\right)}^T,{\left(\begin{array}{cccc}{X}_{d1}& X_{d2}& ...& X_{\mathrm{dN}}\\ Y_{d1}& Y_{d2}& ...& Y_{\mathrm{dN}}\end{array}\right)}^T,$ 则令A＝(a₀ a₁…a₉)^T，B＝(b₀ b₁…b₉)^T，

$\mathrm{\&Gamma;}=\left(\begin{array}{cccccccccc}1& X_{u1}& Y_{u1}& X_{u1}^2& X_{u1}{Y}_{u1}& Y_{u1}^2& X_{u1}^3& X_{u1}^2{Y}_{u1}& X_{u1}{Y}_{u1}^2& Y_{u1}^3\\ 1& X_{u2}& Y_{u2}& X_{u2}^2& X_{u2}{Y}_{u2}& Y_{u2}^2& X_{u2}^3& X_{u2}^2{Y}_{u2}& X_{u2}{Y}_{u2}^2& Y_{u2}^3\\ ...& ...& & & & & & & & \\ 1& X_{\mathrm{uN}}& Y_{\mathrm{uN}}& X_{\mathrm{uN}}^2& X_{\mathrm{uN}}{Y}_{\mathrm{uN}}& Y_{\mathrm{uN}}^2& X_{\mathrm{uN}}^3& X_{\mathrm{uN}}^2{Y}_{\mathrm{uN}}& X_{\mathrm{uN}}{Y}_{\mathrm{uN}}^2& Y_{\mathrm{uN}}^3\end{array}\right),$ ${\&aleph;}_d={\left(\begin{array}{cccc}{X}_{d1}& X_{d2}& ...& X_{\mathrm{dN}}\end{array}\right)}^T,$

Ψ_d＝(Y_d1 Y_d2…Y_dN)^T，Γ^＝[(Γ^TΓ)^-1Γ^T]。

则

B＝Γ^Ψ_d。

6.利用双线性插值生成校正后图像

将校正后图像某点坐标(x_u，y_u)代入式(11)，得到校正前图像上的对应点(x_d，y_d)。但(x_d，y_d)一般不为整数。若

代表下取整，则令：

则：

由于在太空中相机的内外参数不再会有剧烈的变化，所以通过流程1所保存的参数可以用于相机所拍摄的其它图像的校正，即流程3，除非认为需要再次通过流程1进行参数计算。而当星点过少等原因无法通过流程1得到满意的校正图像时，可以利用相机在地面的标定结果对图像进行校正，以降低图像的畸变程度。

图2为6幅视场角为2°×2°的不同视轴指向的标准模拟星空图像其视轴指向分别为：(5.1h，7°)、(5.15h，8°)、(5.2h，9°)、(5.3h，11°)、(5.45h，14°)、(5.5h，15°)，模拟卫星在运行中拍摄的序列图像。实验结果如图3和图4所示。图3为将图2加入径向、切向和薄棱镜畸变后的图像，其中恒星光斑部分用方框标出。将这6幅畸变图像读入文件，并利用图1中流程一的方法进行畸变校正，得到的6幅校正后的图像如图4所示，对比图2和图3可以看出，图3中存在明显的畸变，而图4与图2相比无明显区别，这表明了校正后图像与标准图像非常接近，校正效果显著。

令(X_si，Y_si)、(X_di，Y_di)分别为无畸变图像、有畸变图像上第i个星点的重心坐标，则可利用星点位置由式(13)计算图像的相对误差。

$E_r=\frac{\underset{i}{\mathrm{\&Sigma;}}\sqrt{{(X_{\mathrm{si}}-X_{\mathrm{di}})}^2+{(Y_{\mathrm{si}}-Y_{\mathrm{di}})}^2}}{\underset{i}{\mathrm{\&Sigma;}}\sqrt{X_{\mathrm{si}}^2+Y_{\mathrm{si}}^2}}\&times;100\%---\left(13\right)$

由于在实际处理中无畸变图像是无法得到的，所以，可以将星表匹配得到的理论坐标作为(X_si，Y_si)，进行校正效果的评价。

当相机视场角较大，其所拍到的一帧图像中就包含了不少于十颗恒星，此时可以作为本发明方法的特例，读取和处理图像时仅需对单帧图像进行操作，同样按照本发明所述的方法进行畸变参数计算和误差补偿。

7.副流程说明

图1中给出了主流程(流程一)的同时，还给出了备选的两个副流程，即流程二和流程三。由于相机在太空中运行过程中，其内部参数不会发生剧烈变化，所以不需要每次进行拍摄时都进行畸变参数的测算。而是可以采用将最近一次或几次测算的畸变参数保存，需要对相机拍摄的图像进行校正时仅需打开保存的畸变参数，即可通过具体实施方式6中所述的插值的方法得到校正后图像，即流程二所示，其优点在于基本不降低精度的同时简化了操作复杂度和计算量，可以实现图像的批处理。

由于实现流程一需要输入的图像或序列图像中共拍摄到十颗以上的恒星，可能有时无法达到，或由于其它原因造成不能实现流程一的处理，同时又没有保存的畸变参数可以使用。此时给出用于应急的流程三处理方法，即利用相机在随卫星升空前在地面上测算的畸变参数实现其所拍摄图像的畸变校正。该流程方法虽然在精度上相比于流程一和流程二更低，但流程三处理后的图像比相机拍摄的原始图像精度仍会有一定的提升，其具体实施步骤为：

(1)相机随卫星升空前，首先在地面固定在高精度转台上，利用标准模板按式(14)测算其参数k。

$\left\{\begin{array}{c}{X}_u=X_d+k{X}_d(X_d^2+Y_d^2)\\ Y_u=Y_d+k{Y}_d(X_d^2+Y_d^2)\end{array}\right.---\left(14\right)$

其中图像某参考点准确坐标为(X_u，Y_u)，存在畸变坐标为(X_d，Y_d)。

(2)设(x_u，y_u)为校正后图像上某点像素坐标，通过式(15)计算该点所对应的校正前图像点坐标(x_d，y_d)，

$x_d=x_u\frac{1}{6k(x_u^2+y_u^2)}\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}$

$-\frac{2}{\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}}$

$y_d=y_u\frac{1}{6k(x_u^2+y_u^2)}\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}---\left(15\right)$

$-\frac{2}{\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}}$

若

代表下取整，令：

则通过式(12)便计算出校正后图像每像素点的灰度值，得到校正后图像。

8.误差分析与补偿

考虑到相机光轴指向、星表信息以及图像中星点实际坐标的提取等方面均存在误差，会给校正带来一定影响。对于图像星点的提取，目前的研究较为成熟，能够较高精度地提取图像中星点的准确位置。星表误差主要源于岁差、自行、章动、光行差等造成星表随时间变化，而这些变化均可以通过将基本星表由标准历年到拍摄当天的换算进行补偿，因而也不会产生明显的误差。

光轴指向的误差主要源于卫星自身姿态误差和相机转台转角误差的累加，而前者为主要因素。当卫星姿态误差主要为漂移造成的系统误差时，可以认为在短时间内所拍摄的序列图像光轴指向误差基本相同。针对在多帧图像中光轴指向误差较小且基本相同的情况，提出补偿方法如下。由式(6)可知，当α₀、δ₀变化较小时，可通过式(16)估计偏差。当图像大小为512像素×512像素，视场角为2°×2°，计算可知光轴的赤经、赤纬每各有1′误差时，图像最大偏差约为6.0357像素，显然偏差不能忽略，必须进行补偿。并可得到，赤经每变化1′图像Y轴最大偏差量为0.0026像素；赤纬每变化1′图像X轴最大偏差量为0.0745像素。所以X、Y方向的偏差主要分别来源于赤经、赤纬的误差。

$E_{{\mathrm{\&alpha;}}_0{\mathrm{\&delta;}}_0}=\sqrt{{\left(\frac{\&PartialD;X}{{\&PartialD;\mathrm{\&alpha;}}_0}d{\mathrm{\&alpha;}}_0\right)}^2+{\left(\frac{\&PartialD;X}{{\&PartialD;\mathrm{\&delta;}}_0}d{\mathrm{\&delta;}}_0\right)}^2+{\left(\frac{\&PartialD;Y}{{\&PartialD;\mathrm{\&alpha;}}_0}d{\mathrm{\&alpha;}}_0\right)}^2+{\left(\frac{\&PartialD;Y}{{\&PartialD;\mathrm{\&delta;}}_0}d{\mathrm{\&delta;}}_0\right)}^2}---\left(16\right)$

下面进一步探讨光轴指向存在小的系统误差时图像各像素点偏差量的变化情况，若各像素点在X、Y方向的偏差量基本无变化，则可以认为α₀、δ₀的误差仅会造成图像整体的平移。通过式(6)计算可知，在上述图像大小和视场角条件下，相对于光轴指向准确的情况，α₀、δ₀各有1′误差时图像各像素点偏移量两两相差的最大值为：X方向0.1490像素，Y方向0.0758像素，差别很小，基本可以认为，图像上各像素点沿同一矢量偏移。图5为α₀、δ₀分别误差为1′，在X、Y方向像素点偏移量相差分别最大时光轴指向的情况下，图像各像素点的偏移情况，其中(a)为光轴指向误差导致X轴的偏差，(b)为光轴指向误差导致Y轴的偏差。从图5中可以看出，虽然图像上各点在X、Y方向均有偏移，但偏移量的差别很小，可以忽略。为光轴指向各有1′误差时图像X、Y轴的偏差，

图像准确的畸变模型可表示为式(17)，考虑到输入信息等方面的误差，能够计算得到带有误差的畸变模型估计式，表示为式(18)。

$\left\{\begin{array}{c}{X}_d=X_d(X_u,Y_u)\\ Y_d=Y_d(X_u,Y_u)\end{array}\right.---\left(17\right)$

$\left\{\begin{array}{c}{\hat{X}}_d={\hat{X}}_d(X_u,Y_u)\\ {\hat{Y}}_d={\hat{Y}}_d(X_u,Y_u)\end{array}\right.---\left(18\right)$

当校正误差主要由光轴指向系统误差(dα₀，dδ₀)造成，即式(18)可写为：

$\left\{\begin{array}{c}{\hat{X}}_d=X_d[X_u({\mathrm{\&alpha;}}_0+d{\mathrm{\&alpha;}}_0,{\mathrm{\&delta;}}_0+d{\mathrm{\&delta;}}_0),Y_u({\mathrm{\&alpha;}}_0+d{\mathrm{\&alpha;}}_0,{\mathrm{\&delta;}}_0+d{\mathrm{\&delta;}}_0)]\\ {\hat{Y}}_d=Y_d[X_u({\mathrm{\&alpha;}}_0+d{\mathrm{\&alpha;}}_0,{\mathrm{\&delta;}}_0+d{\mathrm{\&delta;}}_0),Y_u({\mathrm{\&alpha;}}_0+d{\mathrm{\&alpha;}}_0,{\mathrm{\&delta;}}_0+d{\mathrm{\&delta;}}_0)]\end{array}\right.---\left(19\right)$

由前面分析可知，光轴指向的系统误差对校正图像的影响为整体的移动，所以存在不随(X_u，Y_u)变化的常数误差量(ε_x，ε_y)使式(19)可写为式(20)的形式。

$\left\{\begin{array}{c}{\hat{X}}_d=X_d(X_u,Y_u)+{\mathrm{\&epsiv;}}_x\\ {\hat{Y}}_d=Y_d(X_u,Y_u)+{\mathrm{\&epsiv;}}_y\end{array}\right.---\left(20\right)$

由于(dα₀、dδ₀)未知，所以无法由式(19)直接计算误差大小。下面通过对误差量(ε_x，ε_y)的估计，补偿式(18)由于光轴指向系统误差所带来的误差影响，进而得到更精确的畸变校正公式。

设(X₀，Y₀)为图像的畸变中心，当X_d＝X₀，Y_d＝Y₀时，有X_u＝X₀，Y_u＝Y₀，即式(21)成立。

$\left\{\begin{array}{c}{X}_0=X_d(X_0,Y_0)\\ Y_0=Y_d(X_0,Y_0)\end{array}\right.---\left(21\right)$

将式(21)代入式(20)可得：

$\left\{\begin{array}{c}{\hat{X}}_d{|}_{X_u=X_0,Y_u=Y_0}=X_0+{\mathrm{\&epsiv;}}_x\\ {\hat{Y}}_d{|}_{X_u=X_0,Y_u=Y_0}=Y_0+{\mathrm{\&epsiv;}}_y\end{array}\right.---\left(22\right)$

将式(22)代入式(18)即可求得误差量的估计值，如式(23)所示，其中

即为(ε_x，ε_y)的估计值。

$\left\{\begin{array}{c}{\widehat{\mathrm{\&epsiv;}}}_x={\hat{X}}_d(X_0,Y_0)-X_0\\ {\widehat{\mathrm{\&epsiv;}}}_y={\hat{Y}}_d(X_0,Y_0)-Y_0\end{array}\right.---\left(23\right)$

基于误差量的估计，可以对式(18)进行误

差补偿，联立式(17)、式(18)、式(19)、式(23)可得：

$\left\{\begin{array}{c}{X}_d={\hat{X}}_d(X_u,Y_u)-{\hat{X}}_d(X_0,Y_0)+X_0\\ Y_d={\hat{Y}}_d(X_u,Y_u)-{\hat{Y}}_d(X_0,Y_0)+Y_0\end{array}\right.---\left(24\right)$

式(24)即为经过误差补偿的畸变校正模型。在本文方法中，将形式化表达式(18)用三阶多项式的拟合表达式(11)代换，并代入式(24)中可得式(25)。

$\left\{\begin{array}{c}{X}_d=a_1(X_u-X_0)+a_2(Y_u-Y_0)+a_3(X_u^2-X_0^2)+a_4(X_u{Y}_u-X_0{Y}_0)+a_5(Y_u^2-Y_0^2)\\ +a_6(X_u^3-X_0^3)+a_7(X_u^2{Y}_u-X_0^2{Y}_0)+a_8(X_u{Y}_u^2-X_0{Y}_0^2)+a_9(Y_u^3-Y_0^3)+X_0\\ Y_d=b_1(X_u-X_0)+b_2(Y_u-Y_0)+b_3(X_u^2-X_0^2)+b_4(X_u{Y}_u-X_0{Y}_0)+b_5(Y_u^2-Y_0^2)\\ +b_6(X_u^3-X_0^3)+b_7(X_u^2{Y}_u-X_0^2{Y}_0)+b_8(X_u{Y}_u^2-X_0{Y}_0^2)+b_9(Y_u^3-Y_0^3)+Y_0\end{array}\right.---\left(25\right)$

由于(X₀，Y₀)未知，在计算时可用相机在地面畸变校准时测定的畸变中心代替。当图像畸变中心与图像中心重合时，有X₀＝0，Y₀＝0，式(24)还可以进一步简化为

$\left\{\begin{array}{c}{X}_d=a_1{X}_u+a_2{Y}_u+a_3{X}_u^2+a_4{X}_u{Y}_u+a_5{Y}_u^2+a_6{X}_u^3+a_7{X}_u^2{Y}_u+a_8{X}_u{Y}_u^2+a_9{Y}_u^3\\ Y_d=b_1{X}_u+b_2{Y}_u+b_3{X}_u^2+b_4{X}_u{Y}_u+b_5{Y}_u^2+b_6{X}_u^3+b_7{X}_u^2{Y}_u+b_8{X}_u{Y}_u^2+b_9{Y}_u^3\end{array}\right.---\left(26\right)$

在利用插值法生成校正后图像时，用式(25)代替式(11)，从而可自适应地修正由于视轴指向输入不准确带来的误差。而在进行参数计算时，仍使用三阶多项式拟合模型(11)。

设dX、dY分别为光轴指向(α₀，δ₀)每各有1′误差时图像中心沿X、Y轴的偏移量，可通过式(6)求解偏微分方程并代入式(27)计算得到，则光轴指向的误差可由式(28)进行估算，从而进一步对卫星姿态漂移量给出参考，从而有助于卫星惯导系统的误差校正。

$\left\{\begin{array}{c}d{\mathrm{\&alpha;}}_0=\frac{{\mathrm{\&epsiv;}}_x}{\mathrm{dX}}\\ d{\mathrm{\&delta;}}_0=\frac{{\mathrm{\&epsiv;}}_y}{\mathrm{dY}}\end{array}\right.---\left(28\right)$

9.仿真实验

仿真实验结果如图6和图7所示。在此实验中，取输入视轴指向的赤经系统误差和赤纬系统误差均为-2′，图6为不考虑误差补偿情况下对图3进行校正后图像，图7为进行了误差补偿后对图3进行校正图像，图6、图7中6幅图像与图2中6幅图像对应。表1为图2、图6和图7中各星点的重心坐标。图8为校正结果的仿真图。由图8可知，当存在输入误差时，若不考虑误差补偿，则各星点与无畸变情况相差较大，不能满足精度要求；加入误差补偿后，校正后图像与标准图像间的偏差明显减小，校正效果良好。计算可知，校正后图像的相对误差率为3.91％，而利用本方法进行误差补偿后，校正图像的相对误差率为1.05％，相对误差下降显著，校正效果很明显。图9为无畸变a、校正前b、校正后c多帧图像中第一帧图像的对比图，为了便于观察，图中分别加入了与图像畸变程度相同的网格。

表1误差补偿星点坐标对比

无畸变图像星点重心坐标	未考虑误差补偿的坐标	进行误差补偿后的坐标
			(387.999，26.0015)	(396.456，7.4899)	(384.744，23.4044)

(119.093，31.9056)	(125.995，23.1015)	(116.516，29.9324)
			(88.0904，42.9082)	(95.5716，35.0239)	(86.0930，41.5065)
(186.00，128.00)	(194.072，120.479)	(185.321，126.974)
			(36.00，147.00)	(43.1725，137.993)	(31.8097，142.458)
(282.999，164.001)	(289.551，153.449)	(281.032，160.445)
			(135.907，169.093)	(144.961，163.017)	(135.471，168.607)
(86.00，190.00)	(92.9339，181.006)	(82.0155，185.511)
			(354.908，273.093)	(363.028，265.961)	(354.483，272.511)
(199.00，281.00)	(205.230，271.624)	(195.026，276.008)
			(354.094，336.907)	(363.173，330.848)	(353.910，336.459)
(498.00，432.00)	(506.105，424.333)	(495.866，429.794)

由式(27)、(28)可算得，光轴指向误差(dα₀，dδ₀)的估计值为(2.0836′，2.1002′)，与准确值(2′，2′)基本相同。在实现空间图像校正的同时，可以对卫星本身姿态漂移的校正提供有利的帮助。

Claims (9)

1、一种基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于步骤如下：

(1)读取的图像进行分割处理，提取图像中恒星光斑部分；

(2)通过计算各光斑的重心坐标，作为实际坐标；

(3)利用相机在地面标定时得到的畸变参数，对步骤(2)得到的实际坐标进行预校正，使之更接近无畸变情况下的坐标，从而提高星图匹配的成功率；

(4)利用单帧或多帧Hausdorff距离星图匹配方法将星载相机所拍摄的空间图像进行星图匹配，计算各实际坐标对应的理论坐标；

(5)利用实际坐标与理论坐标对应关系，计算图像的畸变参数；

(6)利用插值法得到校正后图像。

2、根据权利要求1所述的基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于：所述步骤(1)中的图像分割方法采用改进型大津阈值的方法，其具体流程为：

(1)令T_min为图像灰度均值，即令图像大小为M×N，图像上(i，j)点的灰度值为G(i，j)，则 $T_{\min }=\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}G(i,j)}{M\&times;N};$

(2)在灰度级从T_min～255区间中利用大津阈值选取最优阈值T，步骤为：

a.令p(k)为灰度级k的频率，即 $p\left(k\right)=\frac{\underset{G(i,j)=k}{\mathrm{\&Sigma;}}1}{M\&times;N},$ k∈[0，255]。

b.若t为分割阈值，将图像分割为两部分：{G(i，j)|G(i，j)≤t}和

{G(i，j)|G(i，j)＞t}，ω₀和ω₁分别为目标和背景部分比例；

${\mathrm{\&omega;}}_0=\underset{t<i\&le;255}{\mathrm{\&Sigma;}}p\left(i\right)$

${\mathrm{\&omega;}}_1=\underset{0\&le;i\&le;t}{\mathrm{\&Sigma;}}p\left(i\right)$

μ₀和μ₁为目标和背景部分的灰度均值：

${\mathrm{\&mu;}}_0=\frac{\underset{t<i\&le;255}{\mathrm{\&Sigma;}}\mathrm{ip}\left(i\right)}{{\mathrm{\&omega;}}_1\left(t\right)}$

${\mathrm{\&mu;}}_1=\frac{\underset{0\&le;i\&le;t}{\mathrm{\&Sigma;}}\mathrm{ip}\left(i\right)}{{\mathrm{\&omega;}}_0\left(t\right)}$

总均值μ＝ω₀(t)μ₀(t)+ω₁(t)μ₁(t)；

c.最佳阈值T可由式(1)计算得到：

$T=\mathrm{Arg}\underset{T_{\min }\&le;t\&le;255}{\mathrm{Max}}[{\mathrm{\&omega;}}_0\left(t\right){({\mathrm{\&mu;}}_0\left(t\right)-\mathrm{\&mu;})}^2+{\mathrm{\&omega;}}_1\left(t\right){({\mathrm{\&mu;}}_1\left(t\right)-\mathrm{\&mu;})}^2]---\left(1\right)$

令μ₀′和μ₁′为利用最佳阈值分割后图像目标和背景的灰度均值，σ₀和σ₁为目标和背景灰度标准差，考察公式(2)：μ₀′-μ₁′＞a(σ₀+σ₁)                (2)

是否成立，其中取a＝3～5，若成立则T即为所求的最佳阈值，若不成立则令T_min＝T_min+1并返回步骤(2)；

(3)遍历图像，若G(i，j)≤T则令G(i，j)＝0。

3、根据权利要求1所述的基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于：所述步骤(2)中的计算各光斑的重心坐标的计算方法采用重心法，具体如下：

(1)初始化与图像矩阵大小相同的标号矩阵S＝[0]、初始化记录搜索到的星点数量的变量int StarNum＝0；以及记录各星点重心坐标的数组float * StarPosition[2]；

(2)利用八连同域方法确定图像G中不为零像素的连通成分，将G中各自连通像素点的连同域标号记入矩阵S中的对应位置，同一连同域内像素点标号相同，不同连同域标号从1开始依次累加，设共搜索到N个连同域，

float*StarPosition[0]＝new float[N]；

float*StarPosition[1]＝new float[N]；

(3)初始化数组long Total[N+1][2]＝[0]；int Area[N+1]＝[0]；遍历S：

若S(i，j)＞0则

Total[S(i，j)][0]+＝G(i，j)×i；Total[S(i，j)][1]+＝G(i，j)×j；Area[S(i，j)]++；

(4)设定连同域面积双阈值T_large、T_small：

若Area[i]＜T_large且Area[i]＞T_small，则：

$\mathrm{StarPosition}\left[0\right]\left[\mathrm{StarNum}\right]=\frac{\mathrm{Total}\left[i\right]\left[0\right]}{\mathrm{Area}\left[i\right]};$

$\mathrm{StarPosition}\left[1\right]\left[\mathrm{StarNum}\right]=\frac{\mathrm{Total}\left[i\right]\left[1\right]}{\mathrm{Area}\left[i\right]};$

           StarNum++；

其中i＝1~N

至此将从图中检测到的星点数量和各星点的重心坐标分别存入StarNum和StarPosition中。

4、根据权利要求1所述的基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于：所述步骤(2)中的实际坐标进行预校正的方法如下：

设通过相机在地面标定的结果，图像某点存在畸变坐标(x_d，y_d)与准确坐标(x_u，y_u)之间的关系为：

$\left\{\begin{array}{c}{x}_u=x_d+k_1{x}_d(x_d^2+y_d^2)+p_1{x}_d({3x}_d^2+y_d^2)+2p_2{x}_d{y}_d\\ y_u=y_d+k_2{y}_d(x_d^2+y_d^2)+2p_1{x}_d{y}_d+p_2\mathrm{yd}(x_d^2+3y_d^2)\end{array}\right.---\left(4\right)$

公式(4)中各参数已测知，其中k₁、k₂为径向畸变系数，其作用使图像某点产生沿与畸变中心连线的径向方向产生畸变，p₁、p₂为切向畸变系数，其作用使图像沿畸变中心发生旋转，

则通过步骤(2)所得到的各星点的实际坐标(X_d，Y_d)，利用式(5)进行预校正，得到预校正后坐标

$\left\{\begin{array}{c}{\tilde{X}}_u=X_d+k_1{X}_d(X_d^2+Y_d^2)+p_1{X}_d({3X}_d^2+Y_d^2)+2p_2{X}_d{Y}_d\\ {\tilde{Y}}_u=Y_d+k_2{Y}_d(X_d^2+Y_d^2)+2p_1{X}_d{Y}_d+p_2\mathrm{Yd}(X_d^2+3Y_d^2)\end{array}\right.---\left(5\right).$

5、根据权利要求1所述的基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于：所述步骤(4)利用单帧或多帧Hausdorff距离星图匹配方法将星载相机所拍摄的空间图像进行星图匹配，计算各实际坐标对应的理论坐标的方法步骤如下：

(1)根据提供的相机拍摄时刻视轴指向及市场角搜索星表得到可能落入视场范围的各恒星的赤经、赤纬；

(2)根据赤经、赤纬(α₀，δ₀)及旋转角计算所有可能落入视场的恒星的理论坐标；

$\left\{\begin{array}{c}X=\frac{N_x}{2}\frac{1}{\tan ({\mathrm{FOV}}_x/2)}\frac{\cos \mathrm{\&delta;}\sin \left({\mathrm{\&alpha;}-\mathrm{\&alpha;}}_0\right)}{\sin \mathrm{\&delta;}\sin {\mathrm{\&delta;}}_0+\cos \mathrm{\&delta;}\cos {\mathrm{\&delta;}}_0\cos (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)}\\ Y=\frac{N_y}{2}\frac{1}{\tan ({\mathrm{FOV}}_y/2)}\frac{\sin \mathrm{\&delta;}\cos {\mathrm{\&delta;}}_0-\cos {\mathrm{\&delta;}\sin \mathrm{\&delta;}}_0\cos (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)}{\sin \mathrm{\&delta;}\sin {\mathrm{\&delta;}}_0+\cos \mathrm{\&delta;}\cos {\mathrm{\&delta;}}_0\cos (\mathrm{\&alpha;}-{\mathrm{\&alpha;}}_0)}\end{array}\right.$

$\left\{\begin{array}{c}{X}^{\&prime;}=X\cos \mathrm{\&psi;}+Y\sin \mathrm{\&psi;}\\ Y^{\&prime;}=-X\sin \mathrm{\&psi;}+Y\cos \mathrm{\&psi;}\end{array}\right.$

(3)将理论坐标与预校正后坐标通过Hausdorff距离方法建立对应关系，并去除匹配失败点；

(4)通过实际坐标与预校正后坐标的对应关系，得到实际坐标与理论坐标的一一对应关系，即各匹配成功点对。

6、根据权利要求1所述的基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于：所述步骤(5)利用实际坐标与理论坐标对应关系，计算图像的畸变参数的方法如下：

将实际坐标数组与理论坐标数组代入三阶多项拟和表达式中，并通过最小二乘法求得式中各系数a₀～a₉、b₀～b₉，即为所需的畸变参数，将畸变参数保存为文件，其中(X_d，Y_d)为实际坐标，(X_u，Y_u)为理论坐标，

$\left\{\begin{array}{c}{X}_d=a_0+a_1{X}_u+a_2{Y}_u+a_3{X}_u^2+a_4{X}_u{Y}_u+a_5{Y}_u^2+a_6{X}_u^3+a_7{X}_u^2{Y}_u+a_8{X}_u{Y}_u^2+a_9{Y}_u^3\\ Y_d=b_0+b_1{X}_u+b_2{Y}_u+b_3{X}_u^2+b_4{X}_u{Y}_u+b_5{Y}_u^2+b_6{X}_u^3+b_7{X}_u^2{Y}_u+b_8{X}_u{Y}_u^2+b_9{Y}_u^3\end{array}\right.$

7、根据权利要求1所述的基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于：所述步骤(6)中利用插值法得到校正后图像的方法如下：

设(x_u，y_u)为校正后图像上某点像素坐标，通过下式计算该点所对应的校正前图像点坐标(x_d，y_d)，

$\left\{\begin{array}{c}{x}_d=a_0+a_1{x}_u+a_2{y}_u+a_3{x}_u^2+a_4{x}_u{y}_u+a_5{y}_u^2+a_6{x}_u^3+a_7{x}_u^2{y}_u+a_8{x}_u{y}_u^2+a_9{y}_u^3\\ y_d=b_0+b_1{x}_u+b_2{y}_u+b_3{x}_u^2+b_4{x}_u{y}_u+b_5{y}_u^2+b_6{x}_u^3+b_7{x}_u^2{y}_u+b_8{x}_u{y}_u^2+b_9{y}_u^3\end{array}\right.$

若

代表下取整，则令：

则：

由此便计算得到校正后图像各像素点的灰度值。

8、一种基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于步骤如下：

将权利要求6中所述的畸变参数文件调出，并按照权利要求7所述方法将待校正图像进行畸变校正，得到校正后图像。

9、一种基于空间图像的Hausdorff距离星图匹配畸变校正方法，其特征在于步骤如下：

(1)相机随卫星升空前，首先在地面固定在高精度转台上，利用标准模板按下式测算其参数k。

$\left\{\begin{array}{c}{X}_u=X_d+k{X}_d(X_d^2+Y_d^2)\\ Y_u=Y_d+k{Y}_d(X_d^2+Y_d^2)\end{array}\right.$

其中图像某参考点准确坐标为(X_u，Y_u)，存在畸变坐标为(X_d，Y_d)。

(2)设(x_u，y_u)为校正后图像上某点像素坐标，通过下式计算该点所对应的校正前图像点坐标(x_d，y_d)，

$x_d=x_u\frac{1}{6k(x_u^2+y_u^2)}\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}$

$-\frac{2}{\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}}$

$y_d=y_u\frac{1}{6k(x_u^2+y_u^2)}\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}$

$-\frac{2}{\sqrt[3]{(108+12\sqrt{3}\sqrt{\frac{4+27k(x_u^2+y_u^2)}{k(x_u^2+y_u^2)}})k^2{(x_u^2+y_u^2)}^2}}$

若

代表下取整，令：

则：

便计算出校正后图像每像素点的灰度值，得到校正后图像。

Images