综合孔径干涉望远镜的高分辨率图像重建

综合孔径干涉望远镜的高分辨率图像重建
高莉莉;刘忠;金振宇
【摘要】利用干涉望远镜成像,可以获取最长基线对应的高频率信息,但往往只能获得部分频域覆盖.为了获得尽可能多的频率信息,可以先通过孔径排布变换,进行更充分的频率采样,再经过干涉图像合成,得到含有完备频率信息的目标高分辨率重建像.介绍了综合孔径干涉望远镜的高分辨率图像重建工作,重点讨论了孔径旋转条件下的干涉图像对齐和合成问题,并成功实施了天文目标的干涉成像观测实验,获得了有完备频率信息的目标高分辨率重建像.
【期刊名称】《天文研究与技术－国家天文台台刊》
【年(卷),期】2009(006)004
【总页数】7页(P327-333)
【关键词】综合孔径干涉望远镜;图像重建;频域覆盖;选代位移叠加法
【作者】高莉莉;刘忠;金振宇
【作者单位】中国科学院国家天文台云南天文台,云南,昆明,650011;中国科学院国家天文台云南天文台,云南,昆明,650011;中国科学院国家天文台云南天文台,云南,昆明,650011
【正文语种】中文
【中图分类】TN911.73
望远镜的口径越大，可以获取的空间频率越高。

然而大口径望远镜的制造面临着技
术、资金等诸多难题。

天文光干涉成像技术的出现为解决这个问题提供了一种方案，利用干涉成像，可以获取最长基线对应的高频率信息。

目前主要的干涉望远镜和光学干涉阵有Keck Interferometer、VLTI、LBT、COAST、CHARA等[1]，其中CHARA的基线最长为350m。

这类望远镜或望远镜阵列的成像特点是一次曝光只能获得部分频域覆盖。

为了获得尽可能多的频率信息，可以通过掩模旋转、孔径排布变换等方法，来进行更充分的频率采样，从而得到含有完备频率信息的目标高分辨率重建像。

由于受大气湍流的影响，地基望远镜的长曝光分辨率远远小于衍射极限分辨率，衍射极限分辨率由λ/r0决定，其中r0为大气视宁度参数[2-3]。

为了
获取清晰的含有部分频率信息的目标干涉图，可以采用自适应光学或高分辨统计重建技术等方法[4-11]来消除大气影响。

中国从上世纪90年代陆续开展了干涉望远镜的相关技术研究。

包括光束合成器，图像重建，光学延迟线等。

其中图像重建是干涉望远镜成功获取天体高分辨率图像的关键技术，也是“斐索型干涉望远镜原理研究” 这个国家自然科学基金重点项
目的核心研究内容。

本文针对该问题展开了深入研究，并进行了相关天文实验。

1 综合孔径干涉望远镜高分辨率图像重建原理
1.1 综合孔径干涉望远镜的频域覆盖特性
与全孔径望远镜相比，干涉望远镜的频域覆盖是不完备的[12]。

可以通过观察望远镜的传递函数(OTF)来研究其频域覆盖特性。

图1显示了由3个子孔径组成的最长基线为1m的干涉望远镜及其传递函数，和直径为1m的全孔径望远镜及其传递
函数。

从图中可以看出，干涉望远镜的频率覆盖是不完备的，是对目标频域信息的一种不完备采样。

要想得到完备频域覆盖，可以通过掩模旋转或孔径排布变换等方法来实现。

如果把图1(a)所示的干涉望远镜旋转3次，再把3次得到的传递函数
相加后就可以得到具有更完整频域覆盖的合成传递函数:
(1)
其中OTF(u,v)表示第i次旋转后具有部分频率覆盖的干涉望远镜的传递函数，OTF(u,v)表示具有更完整频域覆盖的合成传递函数，如图2所示。

(a) (b) (c) (d)
图1 图1(a)干涉望远镜的孔径排布图 (b)干涉望远镜的调制传递函数
(c)全孔径望远镜孔径示意图 (d)全孔径望远镜的调制传递函数
Fig.1 (a)Aperture configuration of the interferometry telescope (b)OTF of the interferometry telescope
(c)The full aperture of the telescope (d)OTF of the telescope with the full aperture transparent
图2 干涉望远镜的合成传递函数Fig.2 Combined OTF of the interferometry telescope
1.2 综合孔径干涉望远镜高分辨率图像重建
综合孔径干涉望远镜的高分辨率图像重建主要包括两部分。

第一部分是消除大气影响，获取清晰的含有部分频率信息的目标干涉图；第二部分是把干涉图合成，得到含有完备频率信息的目标高分辨率重建像。

本文采用迭代位移叠加法(ISA)[11]对单组斑点图实施高分辨率重建。

第i次观测时，ISA的统计传递函数记为Hi(u,v)，对应的点扩散函数hi(x,y)可由单星斑点图统计
平均得：
(2)
其中，hit(x-xmt,y-ymt)表示对t时刻的单星斑点图进行位移；xmt和ymt分别表示x方向和y方向的位移量。

目标经过大气望远镜系统后的衰减像记为ii(x,y)，可
由目标斑点图统计平均得：
(3)
其中，iit(x-xmt,y-ymt)表示对t时刻的目标斑点图进行位移；xmt和ymt分别表示x方向和y方向的位移量。

实际成像过程是，目标星卷积大气望远镜综合系统点扩散函数得到衰减的目标像。

所以，用ISA统计点扩散函数对目标衰减像退卷积，就可得到单组目标干涉图，
在频域表示为：
(4)
其中，i(u,v)为第i组含有部分频率信息的目标干涉图的频谱；Hi(u,v)为ISA的统
计传递函数；Ii(u,v)为目标经过大气望远镜系统后的衰减像的频谱；(u,v)是频域的两个变量。

干涉望远镜一次曝光的频域覆盖不完备，如果把多次旋转得到的频域覆盖相加，就可以得到较为完备的频域覆盖。

同样，用ISA算法重建出来的单组干涉图也只含
有部分频率信息，假设各组干涉图之间没有位移偏差，那么直接相加就可以得到含有完备频率信息的合成干涉图，在频域表示为：
i(u,v)
(5)
其中，(u,v)是含有完备频率信息的合成干涉图的频谱，是目标频谱的估计。

经过上述过程，得到的频谱(u,v)是对目标频谱在频域覆盖范围内截断表示的一个
估计，傅里叶逆变换后得到的目标估计(x,y)存在趾效应，同时由于高频信噪比不足，会导致高频噪声放大现象。

通常会卷积一个最长基线B对应的全孔径望远镜的点
扩散函数psf0(x,y)(爱里斑)来对图像进行切趾，并减小高频噪声放大的影响，这样
做也是为了方便与全孔径望远镜的成像结果进行对比。

i0(x,y)=(x,y)⊗psf0(x,y)
(6)
其中，(x,y)为目标估计；i0(x,y)为等效全孔径望远镜的综合孔径干涉望远镜重建像。

2 干涉图像合成
2.1 观测实验介绍
文[13]作者曾于2002年在实验室里用模板干涉法模拟干涉成像。

类似地，可以在云南天文台1m镜上放置掩模，研究综合孔径干涉望远镜的图像重建工作。

假设
掩模旋转3次可以达到完备频率覆盖，同时考虑1m镜的实际情况，可以提出以
下设计方案：3个直径为的330mm的圆形子孔径，子孔径圆心呈直角三角形排布在直径为670mm的圆上(如图1(a)),最长基线可达1000mm。

再分别把掩模旋转60°和120°，共得到3种排布方式:m1(x,y)、m2(x,y)、m3(x,y)。

2008年11月
14日，我们对几组双星进行了观测。

采用的设备是云南天文台1m RCC反射式望远镜，并在入瞳放置制作好的掩模；CCD为ANDOR ixonEM+ DU897-EI，曝光时间8ms；观测波长为580nm，带宽50nm；大气视宁度参数r0约为100mm；斑点图大小为256×256；每个目标共拍摄3组，每组对应一种掩模方式，而且每组包含目标星斑点图和点源单星斑点图各3000幅，光路如图3。

图3 云南天文台1m镜模板干涉实验光路图Fig.3 Layout of the experiment of pupil-plane interference with the 1 meter optical telescope in the Yunnan Observatory
首先要消除大气湍流对望远镜成像的影响，获取清晰的含有部分频率信息的目标干涉图。

上述掩膜经3次旋转，获得3组不同的孔径排布，目标双星BU989通过不同的孔径排布成像，得到各组目标斑点图，经ISA算法处理后，可重建出各组对
应的目标干涉图，见图4。

2.2 干涉图像对齐
在给出公式(5)时，假设各组干涉图像之间没有位移偏差，然而实际上，由于各组
观测是独立实施的，彼此存在指向精度，迭代计算等误差，导致各组干涉图之间有位移偏差。

所以，需要先来解决干涉图像对齐问题。

(a) (b) (c)
图4 BU989含有部分频率信息的目标干涉图，(a)第1组，(b)第2组，(c)第3组Fig.4 Interferograms of the double-star system BU989 with partial U-V coverage
(a)The first group (b)The second group (c)The third group
观察3组干涉图对应的3个掩模的传递函数如图5(a)、(b)、(c)，发现它们存在公共频率部分(图(d))，即330mm的子孔径对应的频率范围。

理论上，3组干涉图的频谱在这个范围内应该一样，所以可以在频域寻找它们的相关极大值，从而找到其在空域的偏移量。

具体做法如下：
(a) (b) (c) (d)
图5 (a)掩膜1的频域覆盖 (b)掩膜2的频域覆盖 (c)掩膜3的频域覆盖 (d)公共频
率覆盖范围
Fig.5 (a)U-V coverage of mask 1 (b)U-V coverage of mask 2
(c)U-V coverage of mask 3 (d)Common U-V coverage
首先找出3个传递函数的公共频率部分(图d)做一个滤波器，然后对第2幅图
2(x,y)进行滤波得到f 2(x,y)，接着分别让3幅图跟f 2(x,y)互相关，并按极大值位移，就可以得到3幅对齐的目标干涉图，记为i(x+δx,y+δy),其中(δx,δy)为位移量。

类似地，用该方法也可求出对齐的目标衰减像ii(x+δx,y+δy)。

另外，在这里选择对第2幅图进行滤波，是因为无论从拍摄时间，还是从掩模旋转角度来看，第2幅都是最接近其他两幅的，实验也证明，其重建效果最好。

如果参数n大于3，可能要涉及各组干涉条纹移动的相位递推问题，这里就不详细讨论了。

2.3 干涉图像合成
在干涉图像对齐的基础上讨论下面两种干涉图像合成法。

(1)先退卷积后相加
先对每组目标衰减像退卷积，重建出含有部分频率信息单组干涉图i(x,y)(图4)，再对齐相加得到含有完备频率信息的干涉图(x,y)。

由公式(4)和(5)可得:
(7)
其中F-1[]表示傅里叶逆变换。

此时再按照公式(6)，卷积1m全孔径望远镜的点扩散函数后，可以最终得到含有完备频率信息的目标高分辨率重建像。

图6给出了用该方法得到的目标BU989的高分辨率重建像。

观察频谱图发现其存在严重的中频噪声放大现象。

原因是单组干涉图的频率覆盖不完备，尤其是中频率信息不足，直接退卷积就会产生中频噪声放大现象。

针对这个问题，一方面可以采用对斑点图进行滤波[14]的方法加以解决，另一方面也可以通过研究孔径排布变换，优化频率覆盖[15-16]，从而减小噪声放大的影响。

(a) (b) (c)
图6 目标BU989的高分辨率重建像 (a)二维图 (b)三维图 (c)频谱的模
Fig.6 High-resolution reconstruction of the image of the double-star system BU989
(a)Two-dimensional image (b)Three-dimensional plot of the image fluxes
(c)Power of the U-V Spectrum
(2)先相加后退卷积
另外一种方法是先把3组目标衰减像ii(x,y)对齐相加，再进行退卷积，从而得到含有完备频率信息的目标估计(x,y)：
(8)
其中F[]表示傅里叶变换。

此时再卷积1m全孔径望远镜的点扩散函数，得到最终的目标高分辨率重建像。

用这种方法对另一目标BU720实施重建，最终的结果如图7。

(a) (b) (c)
图7 目标BU720的高分辨率重建象(a)二维图(b)三维图(c)频谱的模
Fig.7 High-resolution reconstruction of the image of double-stars system BU720,
(a)Two-dimensional image (b)Three-dimensional plot of the image fluxes (c)Power of the U-V spectrum
Lawson等人在2002年的模板干涉实验就采用了这种合成算法。

该算法通过斑点图位移叠加像相加，可以使中频信息丰厚一些，再退卷积就避免了噪声放大问题。

但同时可能会引入视宁度不匹配的问题，为了减小这种影响，在观测时应该尽可能缩短拍摄周期，减小天顶角的变化范围。

3 实验结果讨论
本实验观测了多颗双星，重建出来的角间距和与WDS星表[17]里的数据相对比，有很好的一致性，如表1所示。

图8给出了目标高分辨率重建像的二维图，以及在波长580nm，带宽50nm时观测到的强度比，均采用第二种合成方法，开窗
显示值从零开始。

由重建结果可以看出，本实验成功实施了天文目标的干涉成像观测实验，获得了一批有完备频率信息的目标高分辨率重建像。

实际观测过程中还要注意以下一些问题：噪声、平场和暗流的处理[18]；大气视宁度的变化范围；掩模方向的一致性；增益固定；防止其他光干扰等。

表1 部分双星的重建结果与WDS星表数据相对比Table 1 Comparison between the r esults reconstructed by“ISA”and the results from WDS catalogue双星BU720HEI42BU989重建结果0 49″0 29″0 22″查星表结果0 4″～0 5″0 2″～0 3″0 2″～0 3″
(a)强度比为2.5∶1 (b)强度比为3.25∶1 (c)强度比为1.31∶1
图8 目标高分辨率重建象(a)BU720，(b)HEI42,(c)BU989
Fig.8 Reconstructed images for(a)BU720，(b)HEI42, and(c)BU989
4 结语
本文介绍了综合孔径干涉望远镜的图像重建工作和相关的天文实验。

重点讨论了孔径旋转条件下的干涉图像合成算法，包括对齐和合成两部分。

其中干涉图像对齐采用的滤波相关法很好地实现了空域对齐。

干涉图像合成部分介绍了两种算法，在合成过程中存在频率噪声放大和视宁度不匹配等问题，经分析可以采用孔径排布变换研究，维纳滤波，减小天顶角的变化范围，缩短拍摄周期等方法来减小其带来的不良影响，从而最终实现综合孔径干涉望远镜的高分辨率图象重建。

参考文献:
[1] A Labeyrie，S G Lipson，P Nisenson. An introduction to optical stellar interferometry[M].Cambridge University Press，2006.
[2] D L Fried.Optical resolution through a randomly inhomogeneous medium for very long and very short exposures[J].J.Opt.Soc.Am.1966, 56
(10):1372-1379.
[3] 刘忠，邱耀辉. 大气象质衰减参数[J]. 云南天文台台刊，1997,3:75-81.
[4] F Roddier. Adaptive Optics in Astronomy[M]. Cambridge University Press，1999.
[5] D Y Gezari， A Labeyrie， P V Stachnik. Speckle interferometry：Diffractionlimited measurements of nine stars with the 200-inch telescope[J].Astrophys，1972，173：L1-L5.
[6] J E Baldwin，C A Haniff，C D Mackay ,et al.Closure phase in highresolution optical imaging[J]. Nature,1986,320:595-597.
[7] K T Knox，B J Thompson.Recovery of images from atmospherically degraded short-exposure photographs[J].Astron，1974，193 (1): L45-L48.
[8] G Weigelt. Modified speckle interferometry:speckle
masking[J]mun,1977，21 (1):55-59.
[9] R H T Bates, W R Fright,F M Cady,et al.Speckle processing,shift-and-add, and compensating for instrument aberration[J].SPIE,1981,373:197-202. [10] J C Christou,E K Hege,J D Freeman et al.Images from astronomical speckle data: weighted shift-and-add analysis[J].SPIE,1985,556: 255-262. [11] Z Liu，Y Qiu，R W Lu.Reconstruction of Video Images Through Turbulent Atmosphere[J].SPIE, 1998，3561：326-331.
[12] 王海涛，周必方. 光学综合孔径干涉成像技术[J]. 光学精密工程，2002，
10(5):434-442.
[13] Lawson P R ，J E Baldwin，D Wilson.Desktop interferometer for optical synthesis imaging[J]. Proc.SPIE，2002，4838:404.
[14] 王大勇，伏西洋，郭红锋，等. 光学稀疏孔径系统的成像及其图像复原[J]. 光
子学报, 2005,34(10):1557-1560.
[15] 陈海亭，江月松，钟宇. 二维圆周光综合孔径阵的优化排列及其成像特性研究[J]. 光学学报，2005，25(12)：1616-1621.
[16] 陈欣扬. Fizeau型天文光学综合孔径望远镜的若干关键技术研究[D]. 中国科学院研究生院博士学位论文，2006.
[17] http://catserver.ing.iac.es/binary/html/catalogue/
binary_catalogue.html.
[18] 刘忠，邱耀辉，楼柯，等. 天文斑点成像中的数据预处理[J]. 云南天文台台刊，1997,4：42-48.。