1. 暨南大学 计算机科学系, 广东 广州 510632;
2. 暨南大学 中法天体测量、 动力学与空间科学联合实验室, 广东 广州 510632
收稿日期: 2015-09-39; 修订日期: 2015-11-23
基金项目: 国家自然科学基金(11273014,U1431227)资助.
作者简介: 傅夏青,女,硕士.研究方向:图像处理.Email:cherishfu2015@sina.com
Residuals and Analysis of Geometric Distortion Correction of CCD Images Based on Geometric Distortion Pattern
1. Department of Computer Science, Jinan University, Guangzhou 510632, China;
2. Sino-French Joint Laboratory for Astrometry, Dynamics and Space Science, Jinan University, Guangzhou 510632, China
在天文观测中，望远镜光学成像系统经常存在几何扭曲，导致观测得到的天体位置存在几何扭曲误差。如哈勃空间望远镜(Habble Space Telescope,HST)的第二代广域和行星照相机(Wide-Field Planetary Camera,WFPC2)拍摄的CCD图像存在严重的几何扭曲。其非常小视场(WFPC2视场中的宽场仅为80×80 arcsec2)的几何扭曲最大可达5 pixel。这种由光学成像系统导致的星像位置移动，就是几何扭曲(Geometric Distortion,GD)。文[1]通过提高线性项的精度提出了一种改正几何扭曲的方案，求解了哈勃空间望远镜的几何扭曲，使得哈勃空间望远镜拍摄的图像在行星相机视场精确到0.02 pixel，在宽场视场精确到0.01 pixel。解决了几何扭曲问题后，哈勃空间望远镜才真正发挥了其作为空间设备在天体测量上的潜能。在地面望远镜CCD成像观测中,也同样存在视场几何扭曲的影响。文[2, 3, 4, 5, 6]对地基天文望远镜几何扭曲的求解做了系统的研究，求解了云南天文台2.4 m和1 m望远镜不同滤光片的几何扭曲模型并做了实际应用。本文基于已知的几何扭曲模型改正观测的CCD图像，并对其中产生的星像位置几何扭曲残差进行了仿真研究。
1 观测资料和几何扭曲模型
1.1 观测资料
本文所采用的图像是2013年2月5日晚使用云南天文台2.4 m望远镜配备的云南暗弱天体光谱及成像仪(Yunnan Faint-Object Spectrog-raph and Camera,YFOSC)B滤光片观测近地小行Apophis所得图像中的58帧。因为云南暗弱天体光谱及成像仪图像的四周有黑框，所以先对原始图像进行适当的裁剪，裁剪后的图像大小为1 900×1 900像素。
1.2 几何扭曲模型
文[2, 3]提出了一种简便有效的测定望远镜CCD视场几何扭曲的新方法，并应用到云南天文台1 m和2.4 m望远镜抖动观测疏散星团所得的图像，求解了图像的几何扭曲模型。经扭曲改正后，视场中星像的定位精度有了明显提高[5, 6]。本文使用的几何扭曲模型正是用这种方法求解得到的。
如图 1，将大小为1 900×1 900像素的图像划分为19×19个大小为100×100像素的子区域，黑点表示子区域的中心像素点(格点)。图 2为求解得到的云南天文台2.4 m望远镜B滤光片的几何扭曲矩阵的部分数值。矩阵大小为19×19，其中每一个数值代表对应格点位置的几何扭曲量。对于该几何扭曲模型数值对应的精度，文[2]做了相关研究，结果显示绝大部分几何扭曲量数值对应的标准差小于0.05 pixel，并且各个标准差的分布是相对随机的，与图像坐标没有显著的相关性。图 3是将几何扭曲矩阵放大200倍后用矢量方式绘出。从图 3可以看到，几何扭曲在图像的上面角落位置比较大。
2 几何扭曲校正
已知的几何扭曲以无扭曲位置为参考。如图 4，点P(xP,yP)受到几何扭曲后移动到P′(xP′,yP′)，记为P→P′，则P点的几何扭曲为$\overrightarrow {{\text{d}}\left( P \right)} = \left( {{\text{d}}x\left( p \right),{\text{d}}y\left( p \right)} \right)$，如(1)式。
|
$\left\{ \begin{gathered}
{\text{d}}x\left( p \right) = {x_{p'}} - {x_p} \hfill \\
{\text{d}}y\left( p \right) = {y_{p'}} - {x_p} \hfill \\
\end{gathered} \right.$
|
(1)
|
由已知的几何扭曲矩阵，图像中任意位置处的几何扭曲可由离该位置最近的4个格点的几何扭曲量用双线性插值的方法求出。
对图像进行扭曲校正的过程如图 5。黑点表示一幅图像中不受扭曲的整像素点，受到几何扭曲后偏移到了空心点的位置。如整像素点P(i,j)受到扭曲后移动到P′(x,y)，那么用(2)式求出P′的位置。在扭曲图像中插值求出P′位置的灰度值，并赋值给像素P。这样就恢复了P像素点无扭曲的灰度值。对图像中的每一个像素应用此过程，就可以得到一幅扭曲校正后的图像。
|
$\left\{ \begin{gathered}
x = i + {\text{d}}x\left( p \right) \hfill \\
y = j + {\text{d}}y\left( p \right) \hfill \\
\end{gathered} \right.$
|
(2)
|
3 实验结果
对观测的CCD图像搜星得到存在几何扭曲的星像坐标，记为B(xgd,ygd)。用(3)式计算其扣除几何扭曲后的星像坐标C(x,y)。式中C的位置未知，其扭曲量也未知，所以无法直接计算C的位置。实际上相邻位置的几何扭曲几乎相等(实际求解的几何扭曲模型最大扭曲量约为±2 pixel)。所以先用B位置的几何扭曲近似C位置的几何扭曲，求出C的近似位置，然后再用迭代的方法计算C更精确的位置。用(4)式将两组坐标相减，得到星像的几何扭曲。图 6中灰色点为星像位置的几何扭曲分别随星像位置的x、y坐标的分布情况。
|
$\left\{ \begin{gathered}
x = {x_{{\text{gd}}}} + {\text{d}}x\left( p \right) \hfill \\
y = {y_{{\text{gd}}}} + {\text{d}}y\left( p \right) \hfill \\
\end{gathered} \right.$
|
(3)
|
|
$\left\{ \begin{gathered}
\Delta x = {x_{{\text{gd}}}} - x \hfill \\
\Delta y = {y_{{\text{gd}}}} - y \hfill \\
\end{gathered} \right.$
|
(4)
|
对扭曲校正后的CCD图像搜星得到的星像中心，记为A(xgdc,ygdc)，与C(x,y)相减，得到扭曲校正后星像位置的几何扭曲残差，即(5)式。图 6中黑色点为星像位置的几何扭曲残差分别随x、y坐标的分布情况。将图 6中黑色点表示的几何扭曲残差量程缩小到± 0.05 pixel ，得到图 7。图 8为星像的几何扭曲残差放大200倍后的矢量分布图。
|
$\left\{ \begin{gathered}
\delta x = {x_{{\text{gdc}}}} - x \hfill \\
\delta y = {y_{{\text{gdc}}}} - y \hfill \\
\end{gathered} \right.$
|
(5)
|
从图 6可以看到，图像做扭曲改正前，星像的几何扭曲最大可达1.7 pixel。扭曲改正后，星像的几何扭曲残差基本分布在0附近。从图 8可以看到，星像的几何扭曲残差分布的均值为0.008 pixel，标准差为0.009 pixel，并且与图 3相比明显消除了几何扭曲的区域系统性。在图 7中将残差量程缩小后，可以看到几何扭曲残差大部分在± 0.02 pixel内，但发现几何扭曲残差的分布呈现一定的系统趋势。出现这个现象，可能的一个原因是几何扭曲模型在图像边缘处的变化比较大，造成图像改正后的几何扭曲残差也比较大。为了进一步研究影响几何扭曲残差的因素，本文做了星像几何扭曲的仿真实验。
4 仿真实验
4.1 仿真实验设计
如图 9，(a)为用一维连续高斯函数模拟的星像f(x)，表达式为(6)式，xc是模拟星像的中心。(b)为用正弦函数模拟的扭曲函数模型g(x)，表达式为(7)式。图中为了突出星像扭曲的效果，在对扭曲函数绘图时缩小了扭曲函数的周期。实际上，一个星像范围内的扭曲变化非常缓慢。(c)为对模拟星像函数添加泊松噪声，再根据扭曲函数添加扭曲后得到的连续星像函数。根据扭曲前后模拟星像函数的积分面积不变，由(8)式可以推导出受扭曲后的连续星像函数h(t)，如(9)式。
|
$f\left( x \right) = B + H{e^{\frac{{ - {{\left( {x - {x_{\text{c}}}} \right)}^2}}}{{2{\sigma ^2}}}}}$
|
(6)
|
|
$g\left( x \right) = A\sin \left( {\frac{{2\pi }}{T}x} \right)$
|
(7)
|
|
$\left\{ \begin{gathered}
x + g\left( x \right) = t \hfill \\
f\left( x \right){\text{d}}x = h\left( t \right){\text{d}}t \hfill \\
\end{gathered} \right.$
|
(8)
|
|
$h\left( t \right) = \frac{{f\left( x \right)}}{{1 + g'\left( x \right)}}$
|
(9)
|
根据文[7]对不同星像定心算法的比较和精度分析表明，高斯拟合法是一种精度相对较高的定心算法，且星像分布与高斯函数相似程度越高，则测量精度越高。实际求解的几何扭曲模型最大扭曲量约为2 pixel，扭曲后的星像分布函数非常接近高斯函数。因此采用高斯拟合法对受扭曲后的星像定心。具体地，对受扭曲后的连续星像函数在其± 2.5σ像素范围内采样(σ是高斯函数的标准差)。设此范围内覆盖的像素区间为(i,i+n)。即对(i,i+1)、(i+1,i+2)…(i+n-1,i+n)区间分别积分，结果得到一组离散数据点，用图 9(c)中的点表示。对这些离散数据点用(10)式做最小二乘拟合，求出受扭曲后的星像中心xc′。实验中对每个像素区间积分时，采用两点高斯数值积分[8]。
|
$H\left( t \right) = B + H'{e^{\frac{{ - {{\left( {t - {x_{{\text{c'}}}}} \right)}^2}}}{{2\sigma {'^2}}}}}$
|
(10)
|
(6)~(10)式中，B为背景值，实验中设定星像扭曲前后背景值不变；H为星像亮度值；A是扭曲函数幅度，实验中参考真实资料求解得到的扭曲矩阵，将A设为2 pixel；T是扭曲函数的周期，T越小，扭曲变化越快；H′、xc′、σ′ 为对离散点做最小二乘拟合后求得的高斯函数的参数。
4.2 仿真实验结果
设图像的宽度是1 024 pixel，在[0, 1023]区间内随机生成500个一维模拟星像函数。由(11)式计算几何扭曲残差，实验结果如图 9。
|
$\delta {\text{d = }}{x_{\text{c}}}' - \left[{{x_{\text{c}}} + g\left( {{x_{\text{c}}}} \right)} \right]$
|
(11)
|
图 10中(a)表示设定扭曲函数的周期T=512 pixel，模拟星像的全峰半宽分别为13、9、4 pixel时，几何扭曲残差随星像x坐标的分布情况。表 1是它的数值统计。(b)表示设定模拟星像的FWHM=9 pixel，扭曲函数的周期分别为256、512、1 024 pixel时，几何扭曲残差随x坐标的分布情况。表 2是它的数值统计。
表 1 相同周期不同全峰半宽下扭曲残差数值统计
Table 1 The GD residuals statistics in the same period but with different FWHM
|
FWHM(T=512)/pixel |
Max/pixel |
SD/pixel |
| 4 |
0.005 9 |
0.001 86 |
| 9 |
0.010 1 |
0.002 76 |
| 13 |
0.010 3 |
0.004 00 |
表 2 相同全峰半宽不同周期下扭曲残差数值统计
Table 2 The GD residuals statistics in the same FWHM but in different periods
|
T(FWHM=9)/pixel |
Max/pixel |
SD/pixel |
| 256 |
0.015 8 |
0.005 64 |
| 512 |
0.008 7 |
0.002 72 |
| 1 024 |
0.007 5 |
0.002 62 |
从图 10(a)和表 1可以看到，当扭曲函数的周期相同时，星像的全峰半宽越大，几何扭曲残差越大。例如在固定周期T=512 pixel情况下，当FWHM=4 pixel时，几何扭曲残差最大为0.006 pixel，标准差为0.001 86 pixel；当FWHM=9 pixel时，几何扭曲残差最大增至0.01 pixel，标准差增大到0.002 76 pixel；当FWHM=13 pixel时，几何扭曲残差最大值虽然几乎没有增大，但标准差达到了0.004 pixel。从图 10(b)和表 2可以看到，当星像的全峰半宽相同时，扭曲函数的周期越小，几何扭曲残差越大。例如在固定FWHM=9 pixel情况下，当T=1 024 pixel时，几何扭曲残差最大为0.007 5 pixel；当T=256 pixel时，几何扭曲残差最大达到0.015 8 pixel。也就是说，星像的全峰半宽越大，扭曲函数的周期越小，几何扭曲残差越大。因为星像的全峰半宽越大，扭曲函数的周期越小，都将导致一个星像范围内所包含的扭曲变化越大。从图 10还可以看到，几何扭曲残差的变化趋势跟扭曲函数的变化一致。图 7中几何扭曲残差的分布在图像的四周比较大，中间部分则较小。对应图 3的几何扭曲模型，几何扭曲在图像的角落处变化较快，相当于仿真实验中的扭曲函数的周期T较小的情况；几何扭曲在图像的中间部分变化较小，相当于仿真实验中扭曲函数的周期T较大的情况。从而导致了图 7中分布在图像四周的星像几何扭曲残差较大，分布在图像中间的星像几何扭曲残差较小。
5 结论
在已经求解出几何扭曲模型的情况下，对CCD图像进行几何扭曲校正。校正后，星像位置在X、Y两个方向的几何扭曲残差约在±0.02 pixel内。对星像位置的几何扭曲残差进行仿真研究，分别测试了星像在相同的扭曲函数周期下不同全峰半宽，以及星像在相同的全峰半宽下不同扭曲函数周期的几何扭曲残差的大小和变化规律。结果发现星像的全峰半宽越大，扭曲函数周期越小即扭曲模型变化越快，星像位置的几何扭曲残差越大。这个结论解释了CCD图像经几何扭曲校正后分布在图像四周的星像其几何扭曲残差较大，而分布在图像中间的星像其几何扭曲残差较小的现象。
