多次曝光图像的PSF反卷积以及欠采样图像的反混叠技术II：LSPDAA的制作方法

1技术领域

在图像的点扩散函数(psf)反卷积领域，人们发展了很多方法：比如基于最大似然估计的lucy-richardson图像psf反卷积技术和在该技术基础上加入的正则化手段，比如基于最小均方差或最小二乘法原理的wienner滤波psf反卷积技术，再比如基于熵最大化的psf反卷积技术等等。以上技术所涉及到的图像反卷积都是单幅图像的psf反卷积。但在某些领域往往需要多次曝光，比如视场中有极亮的源，单次长曝光会导致亮源像的过饱和，短曝光又看不到背景信息，因此需要多次曝光再后期叠加，再比如天文观测领域，遥远星系的光子流量有限，需要深度曝光，但往往望远镜的姿态控制问题会导致深度曝光具有很大难度，因此多次曝光也是必要的选项。有了多次曝光图像，那对这些图像的psf反卷积和对欠采样图像的反混叠技术就是处理这些图像的必不可少的手段。本专利技术是一系列关于“多次曝光图像的psf反卷积以及欠采样图像的反混叠技术”的第二项专利，是基于最小二乘法的psf反卷积反混叠技术(leastsquarepsfdeconvolutionandanti-aliasing，lspdaa)；该系列第一项专利是《多次曝光图像的psf反卷积以及欠采样图像的反混叠技术i：psfdelrsa》。

2

背景技术：

在图像成像领域，像素点是采集光子(或其他粒子)和组成图像的基本单元。实际中像素点数目是有限的，如此一来图像分辨率不仅受限于光学分辨率还受限于图像采样率，也就是所谓的欠采样，欠采样会导致信号混叠。另一方面由于仪器的光路中支架，以及仪器本身的口径会导致光学衍射，形成形如爱里斑机制那样的衍射图案，也就是点扩散函数(psf)，psf会导致图像细节模糊。信号混叠和psf效应都是对图像高频端(也就是细节)起作用的，二者联合起来再加上噪音就是我们观测到的图像，因此我们最终的观测图像会被抹掉许多高频的细节。实际中这些细节往往才是我们关注的，比如监控下嫌疑人的面目，比如遥远星系的形态，旋涡星系的旋臂，再比如相控阵雷达下敌机的特征等等，关注细节也是增加像素或者信号采集单元的初衷。所以我们需要一种高效的技术手段从信号混叠兼psf退化的并被噪音污染的多次曝光图像中来提取细节信息，本专利技术使用最小二乘法来作为限制目标图像与真实图像差异的手段，使观测图像和模拟观测(就是用与观测网格相同的网格对目标网格进行降采样，shift，旋转，像场畸变等等)图像的(噪音加权的)均方差最小化。

3

技术实现要素：

3.1基于最小二乘法的psf反卷积反混叠技术-lspdaa

有nimg幅经历了psf，欠采样和加噪音的观测图像，每一幅尺寸是lx，×ly，那么我们可以用ik，i，j来描述所有观测像素(其中ik，i，j∈ds，ds是所有nimg×lx×ly个观测像素的值)。这nimg幅观测图像要叠加到达到临界采样率的高分辨率目标网格m，n上，得到的在目标网格上的值是om，n，而观测像素ik，i，j与目标网格(m，n)像素重合的面积记为wk，i，j，m，n.描述观测像素与模拟观测像素之间拟合程度的χ²可以写为：

其中σk，i，j是第(k，i，j)观测像素处的噪音大小，是模拟观测在第(k，i，j)观测像素处的值，这里mp空间是第(k，i，j)观测像素通过psfk，pq，s，t以及wk，i，j，p，q影响到的所有目标网格上的像素。要求最小化的χ²，那么就要让公式1对目标网格像素om，n求偏导数并令其为零：

我们以上标r表示第r次迭代，那么对于目标网格上第(m，n)像素而言公式2可以展开为：

若令和则化简后可以得到第(m，n)像素在第r+1次迭代后的值是：

很明显等式右边第一项是一个差异改正项，是加性的。可以看到迭代中带有上标r的是需要重新算的，而ψk，i，j，m，n是不需要重新计算的。并且可以看出其实是上一次的迭代结果重新卷积psf并模拟观测后的模拟观测像素值。

值得注意的是我们之前开发的psfdelrsa技术是乘性改正项的，而本发明lspdaa则是加性改正项的。另一个区别是lspdaa使用了噪音加权，后续工作中我们会在psfdelrsa中也考虑噪音加权。

3.2结果与比较为了跟实际情况接近，我们每次曝光使用的psf都是不同的，这里使用著名的misslena图像进行测试。如图1中只有右下图是低于奈奎斯特临界采样率1/2×1/2倍的欠采样图像也就是观测图像，其余三幅图都是在目标网格上达到临界采样率的。

实例一：我们对临界分辨率511×511的misslena真实图像进行1/2×1/2倍分辨率的10次低采样率观测，包括不同的位移旋转。图1中，左上是misslena真实图像，而右上是左上的真实图像卷积了一个psf，左下图是右上图再加上柏松噪音的结果。右下是10副观测图像中的一副，4倍降采样。

图2中我们展示了本发明lspdaa和psfdelrsa的psf反卷积和反混叠能力。左上是misslena真实图像，右上是10幅形如1中右下那样的观测图像的直接shiftandadd(sa或drizzle)叠加(snr＝21.816db；psnr＝1.602db)，可见drizzle既没有反卷积psf也没有很好地反混叠。左下是psfdelrsa在30次迭代中信噪比snr最高的一次(第5次，snr＝25.137db；psnr＝8.243db)，右下图是lspdaa在30次迭代中信噪比snr最高的一次(第29次，snr＝26.550db；psnr＝11.069db)。可以看到lspdaa与psfdelrsa能恢复大部分细节，包括双层帽檐，睫毛，脑后的小花，边缘更锐利。但随着迭代次数增加图像更加锐利，噪音也逐渐增加，从而导致信噪比下降。这是基于lucy-richardson迭代所共有的特征，经过测试如果在迭代途中加入滤波，峰值信噪比会提高2-3db。lspdaa边界上的不连续性是我们留的缓冲区，信噪比统计是所有图像都不计这个缓冲区的。

图3中我们展示了drizzle，psfdelrsa和lspdaa的结果与真实图像的残差(resid-ual)。左上是drizzle的残差图，右上是lspdaa第5次迭代的残差，左下是psfdelrsa的最好拟合的残差，右下是lspdaa的最好拟合的残差。lena图像越清晰，说明残留信号越多，可以看到drizzle是最差的，lspdaa第5次迭代的信噪比与psfdelrsa的最好拟合的信噪比是相当的，不同的是前者在残差中留下了很多信号而后者留下了很多在平滑区域的噪音。当然右下的lspdaa的最好拟合残差，信号残留是最少的。

图4和图5分别对比了psfdelrsa(划线)和lspdaa(实线)在psnr和snr指标中的优劣。很明显可以看到psfdelrsa以较快的速度达到信噪比最高点，而后快速下降。而lspdaa几乎一直保持了上升的势头过了20次之后上升缓慢，第29次达到峰值。总之，我们的lspdaa在峰值信噪比psnr上，比没有使用lspdaa技术的sna的结果提高了将近10个db，比此前开发的psfdelrsa也高了将近3个db。

3.3计算复杂度

psfdelrsa的计算复杂度其实是lucy-richardson迭代与fidrizzle之和，因为涉及psf反卷积所以每一步都要有一次psf“卷积”和“相关”的操作，这恰恰是lucy-richardson迭代带来的计算量。lspdaa的计算复杂度取决于psf的大小，因为大尺寸的psf会影响更多的目标像素，从而加大mp空间的尺寸，导致求和项数增加。实际计算中如实例一使用的10个psf都是15×15的，对于相同的迭代次数，运算速度比psfdelrsa快50％左右。如果psf尺寸增加到25×25，二者速度是相当的。

3.4讨论和结论

lspdaa相对于psfdelrsa的优势：

1.lspdaa能够比psfdelrsa获得更好的图像重建效果。

2.比psfdelrsa更好地抑制了噪音。

3.lspdaa并不要求观测的像素计数值非负，这对于扣过背景的图像而言是一大利好，因为减过背景的图像在某些像素往往会出现负值的计数。

lspdaa尚存的缺陷：比较消耗计算资源，收敛较慢；编程比较复杂，要求较高的编程能力。

在不久的将来很多新的望远镜设备要投入观测中来：比如nasa的widefieldinfraredsurveytelescope(wfirst)，欧空局的euclid，美国国家自然科学基机会资助的largesynopticsurveytelescope(lsst)和中国的天宫望远镜(thechinesespacestationopticaltelescope，cssot)。届时会产生海量的数据，如何快速地高保真地处理这些数据将成为当务之急。实际上这恰好给了多次曝光图像重建技术很大的用武之地。

天文领域：

1.对于跟踪姿态保持不够好的望远镜，可以使用多次曝光技术，缩短曝光时间，这样就会大大降低由于姿态不稳引起的图像轮廓模糊，然后再用多次曝光图像重建技术重建高分辨率观测图像。

2.无论是空间望远镜还是地面望远镜，都可以选择好的时段对准被观测源进行多次曝光，然后用多次曝光图像重建技术来重建源的高分辨率图像。这就避免了不良天气条件或者空间望远镜被地球月球等遮挡的因素影响。

3.由于天文观测中，暗弱源往往占大多数，有了高效可靠的多次曝光图像重建技术就可以使用这些数目可观的以前只能舍弃的暗弱源来做更深空间的天文研究。

4.只要知道拍摄望远镜的硬件情况，位置等参数我们就可以在多次曝光图像重建技术的协助下对同一个源的来自不同望远镜不同历史时期的观测数据进行处理并生成高分辨率图像，从而达到充分利用历史数据的目的。

5.天文观测在有些情况下进行多次曝光是必须的，因为在有多个信号源而且亮度相差很大的情况下，为了增强弱源的信号强度同时又要保证强源不饱和，就必须采用多次曝光，这时候再用我们的多次曝光图像重建技术快有有效地重建具有更高保真度的源的图像就很必要了。

数字图像监控领域，lspdaa多次曝光图像重建技术能够利用连续的多帧录像画面重建得到被监控事物的高清图像。

微观物理学，微生物，医学成像等领域可以在像素分辨率有限的情况下用多次曝光结合多次曝光图像重建技术得到病毒，细菌或者有机分子的高分辨率图像。

4具体实施方式

我们用c语言开发了lspdaa软件，软件会自动读取图像fits头文件的一些参数(如位置，旋转，ccd形变等)并重建高分辨率的源图像。

致谢：本发明受国家973项目(no.2015cb857003，no.2015cb857000，no.2013cb834900)，国家基金委项目(no.11333008，no.11233005，no.11273061)，江苏省杰出青年项目(no.bk20140050)和中科院宇宙学结构先导项目(no.xdb09010000)的支持，发明者在此对这些支持项目表示感谢。

5附图说明

本专利共有5幅附图，用来在视觉表观和定量方面展示psfdelrsa的效果。附图说明如下：

1.实例一图1：左上是misslena真实图像，而右上是左上的真实图像卷积了一个psf，左下图是右上图再加上柏松噪音的结果。右下是10副观测图像中的一副，被降采样4倍。

2.实例一图2：左上是misslena真实图像，右上是10幅形如1中右下那样的观测图像的直接shiftandadd(sa或drizzle)叠加(snr＝21.816db；psnr＝1.602db)，左下是psfdelrsa在30次迭代中信噪比snr最高的一次(第5次，snr＝25.137db；psnr＝8.243db)，右下图是lspdaa在30次迭代中信噪比snr最高的一次(第29次，snr＝26.550db；psnr＝11.069db)。

3.实例一图3：展示了drizzle，psfdelrsa和lspdaa的结果与真实图像的残差(resid-ual)。左上是drizzle的残差图，右上是lspdaa第5次迭代的残差，左下是psfdelrsa的最好拟合的残差，右下是lspdaa的最好拟合的残差。

4.实例一图4：psfdelrsa(划线)和lspdaa(实线)重建图像的峰值信噪比psnr随迭代次数增加的变化情况。

5.实例一图5：psfdelrsa(划线)和lspdaa(实线)重建图像的信噪比snr随迭代次数增加的变化情况。