CN103761712B - 基于自适应光学系统点扩散函数重建的图像盲卷积方法 - Google Patents

Abstract

本发明属于自适应光学图像处理领域，具体是一种基于自适应光学系统点扩散函数重建的图像盲卷积方法。该方法分析了成像过程并将光学传递函数分为两个主要部分：由望远镜系统的象差，重力变形和热变形等所造成的准静态的光学传递函数与自适应光学系统校正之后的残差所建立的动态光学传递函数。通过测量数据计算出这两个光学传递函数，并将其作为基础计算得到盲卷积的先验点扩散函数。利用估计得到的光学系统的点扩散函数作为盲卷积优化点扩散函数的先验分布，以图像总的梯度和作为评价函数重要指标，同时利用相关的大气湍流参数确定可重建部分图像的大小，通过迭代可以快速获得高质量的图像。

Description

基于自适应光学系统点扩散函数重建的图像盲卷积方法

技术领域

本发明属于自适应光学图像处理领域，具体是一种基于自适应光学系统点扩散函数重建的图像盲卷积方法。

背景技术

地面通过望远镜在对天体进行观察，成像过程中会受到多种外部作用而影响成像质量，包括：大气湍流作用，望远镜的重力变形和热变形等。主动光学技术可以部分的校正望远镜主镜的重力变形和热变形，但是对于没有配备主动光学的望远镜，镜面在进行观测时的变形往往难以校正和测量。自适应光学技术主要用来校正大气湍流对于成像的影响，但是受限于波前测量精度，变形镜校正自由度，校正速度，系统的光学设计等，往往无法完全校正大气湍流影响。此外，由于自适应光学系统的测量光路和科学相机成像光路的结构不是完全相同的，非共光路像差对于成像也有很严重影响。

为了进一步提高图像质量，文献《Deconvolutionfromwave-frontsensing:anewtechniqueforcompensatingturbulence—degradedimages》提出了使用波前测量信息对图像进行解卷积处理，来复原图象的方法，但是受限于波前测量精度，这种方法实用性不佳。

专利申请号为：201210246852.X的专利曾提出过使用基于非迭代自解卷积图像盲卷积方法的图像恢复技术，该方法通过选取质量较好的帧图像作为基础，并结合波前测量结果确定解卷积中的可调参数。通过累加得到系统的一个接近长曝光的光学传递函数，并将这一曝光函数作为解卷积时可调参数的参考，通过一次的解卷积计算去除自适应光学系统校正残差，以达到提高图像质量的目的。可见，对于该方法来说，波前测量得到的残差与成像光路的残差差距越小，则使用这一方法经过处理后的图像质量越好。但是，由于波前探测器和变形镜的物理结构等原因，直接使用波前探测器的测量结果而不考虑变形镜的结构估计得到的测量残差，对于高频误差的估计往往是不足的。更为重要的是，由于波前测量系统和科学相机所在的为不同的光路，由光路不同造成的非共光路象差对于最终图像质量影响也比较大。此外，由于望远镜的准静态象差(望远镜的重力变形和热变形)的影响，主镜拼接结构，以及望远镜的副镜和支撑副镜的部分的衍射的影响，也会使最终的图像质量下降，在图像重建的过程，这一部分也要考虑。因此想要能够较好的恢复图像，要求能够以高的精度估计出整个光学系统的影响，即系统的点扩散函数(PSF)。之后以这一PSF为基础，可以通过反卷积，提高图像质量。

文献《FirstSuccessfulAdaptiveOpticsPSFReconstructionatW.M.KeckObservatory》报告了在Keck望远镜上估计自适应光学系统PSF的工作。他们通过使用变形镜的控制信息以及差分传感器，估计出了大气湍流参数和系统的静态象差等。在对于自适应光学系统的PSF估计的过程中，将光学系统的点扩散函数的估计分为两个部分：系统的静态象差部分和系统的动态象差部分。对于系统的动态象差，通过使用变形镜信号进行了估计；在估计系统的静态象差时通过使用导星和差分传感进行了测量。但是，结果表明估计出的PSF和真实的PSF还有一些差距。因此，如果直接使用这样的方法得到的PSF进行反卷积，其结果还是不理想的。

在PSF未知或者了解不多的情况下，盲卷积是一种使用比较广泛的后期图像处理方法，这种图像盲卷积方法将成像过程用以下的表达式进行描述：g(x，y)为观测到的图像，f(x，y)为真实目标；h(x，y)表示光学系统的点扩散函数；用来描述包括大气湍流；望远镜以及自适应光学系统等对于成像的影响；n(x，y)表示系统的噪声；表示卷积操作。参考先验条件，通过优化图像盲卷积方法迭代逐渐改变PSF和图像，可以达到恢复图像的结果。文献《Totalvariationblinddeconvolution》提出了基于全局变分的盲卷积图像盲卷积方法，这种图像盲卷积方法被广泛的应用在众多的图像处理任务中。该图像盲卷积方法最小化如下评价函数，通过迭代得到优化后的图像：上式中，f为目标函数，k为点扩散函数，u为优化之后的图像，z为原始图像，α₁和α₂为可以调整的正常数。这一优化图像盲卷积方法的收敛速度和最终的结果，依赖于点扩散函数的初始值和评价函数。目前所出现的图像盲卷积方法，一般采用的是人工经验确定的点扩散函数，很难实现对于点扩散函数的准确估计，导致收敛速度慢，且处理后的图像效果不佳。

发明内容

本发明所要解决的技术问题是：如何能更好地提高经过自适应光学系统校正之后的天文望远镜的图像质量。

本发明所采用的技术方案是：基于自适应光学系统点扩散函数重建的图像盲卷积方法，包括

时钟电路：为望远镜系统提供时钟服务；

第一存储器：连接望远镜系统主镜面形测量信号传入接口，实时记录望远镜主镜面形状随时间变化的信号φ_M(t)；

第二存储器：连接自适应光学系统波前探测器探测信号传入接口，实时记录波前探测器随时间变化的测量信号φ_wfs(t)；

第三存储器：连接望远镜科学相机曝光信号传入接口，实时记录天文图像曝光时间T,(T为一个离散数据)；

第四存储器：连接望远镜系统的自适应光学系统控制信号连接接口，实时记录变形镜(变形镜是自适应光学系统的一部分)随时间变化的控制信号φ_dm(t)；

第五存储器：连接望远镜科学相机，用来记录曝光图像img(t)；

核心处理器：连接第一存储器、第二存储器、第三存储器、第四存储器、第五存储器和时钟电路；

核心处理器控制存储器进行信号的存取工作。

核心处理器控制存储器的工作步骤如下：

步骤一，根据科学需求，选择需要处理的图象的曝光时刻t₀(t₀为T的一个值)和曝光时长Δt(t₀为T的一段时间)，从第一存储器中读取φ_M(Δt)，从第二存储器读取φ_wfs(Δt)，从第四存储器读取φ_dm(Δt)，从第五存储器中读取img(t)，根据φ_wfs(Δt)和φ_dm(Δt)计算出大气相干长度r0和外尺度大小L0随时间变化的序列，通过统计选择概率最大的相干长度和外尺度，计算出大气湍流的相干角度，从img(t)里选择包含感兴趣的目标大小为相干角度之内的原始图像subimg(Δt)作为初始的需要处理的图像；

步骤二，主镜影响的光学传递函数OTF_static(Δt)＝A·exp(-i·φ_M(Δt))·P·OTF_sec(Δt)，A为描述主镜的口径和拼接结构的光学传递函数，P根据主镜和副镜之间结构确定的光线自由传递函数，OTF_sec(Δt)为副镜口径和副镜支架确定的副镜光学传递函数；

步骤三，自适应光学系统的矫正残差的光学传递函数OTF_dynamical(Δt)＝OTF_wfs(Δt)·OTF_dm(Δt)(参考论文Imagereconstructionbymeansofwave—frontsensormeasurementsinclosed—loopadaptiveoptics)，波前探测器测量残差光学传递函数OTF_wfs(Δt)＝exp(-i·φ_wfs(Δt))，变形镜响应函数光学传递函数OTF_dm(Δt)＝exp(-i·φ_dm(Δt)·h_dm)，变形镜响应函数h_dm是由变形镜物理结构决定的，一般由厂家直接给出，总光学传递函数OTF(Δt)＝OTF_static(Δt)·OTF_dynamical(Δt)；

步骤四，利用MATLAB中的函数OTF2PSF，把光学传递函数计算为系统点扩散函数PSF(Δt)＝OTF2PSF{OTF(Δt)}，Nor()为归一化函数，将点扩散函数总能量归一化；

步骤五，计算初始图像的整体梯度范数和确定优化的参数λ的初始值λ₀＝0.2；

步骤六，利用点扩散函数的初始值PSF_all，可以通过维纳逆滤波得到经过该点扩散函数重建之后的图像img_all(Δt)；

步骤七，计算此时的评价函数值：

$F\left({\mathrm{img}}_{all}\right(\mathrm{\Δ}t),{\mathrm{PSF}}_{all})={\mathrm{\λ}}_0||{\mathrm{PSF}}_{all}\⊗{\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right)-subimg\left(\mathrm{\Δ}t\right)|{|}^2-log\left(\mathrm{\Σ}\right|\▿\left({\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right)\right)\left|\right);$

步骤八，根据评价函数值，随机地改变点扩散函数的初始值PSF_all局部的大小(改变点扩散函数中高频部分的灰度值，以达到改变点扩散函数形状的效果)，得到点扩散函数PSF₁和利用此点扩散函数通过维纳逆滤波得到图像img₁(Δt)，计算此时的评价函数值：

$F\left({\mathrm{img}}_1\right(\mathrm{\Δ}t),{\mathrm{PSF}}_1)={\mathrm{\λ}}_0\left|\right|{\mathrm{PSF}}_1\⊗{\mathrm{img}}_1\left(\mathrm{\Δ}t\right)-subimg\left(\mathrm{\Δ}t\right)|{|}^2-log\left(\mathrm{\Σ}\right|\▿\left({\mathrm{img}}_1\left(\mathrm{\Δ}t\right)\right)\left|\right)$

步骤九，比较img_all(Δt)和img₁(Δt)的整体梯度范数和之差，若大于0.001，则改变优化参数λ的大小，若整体梯度范数和之差增加，则减小0.01，否则增加0.01，返回步骤八反复迭代，直到整体梯度范数和之差小于等于0.001。

本发明的有益效果是：通过记录自适应光学系统的相关参数来确定大气的参数，用以确定点扩散函数一致的区域，以避免由于图像的点扩散函数不一致对于真实图像的恢复的影响。在对于点扩散函数初值的估计中，参考了望远镜的结构，自适应光学系统的变形镜和波前探测器的物理结构，以及望远镜主镜的测量信号和自适应光学系统的控制，测量信号。这样重建的点扩散函数与物理实际模型联系了起来，从而使得初始的点扩散函数与真实情况下的点扩散函数更加接近，大大减小了优化的时间，避免了局部最优对于解的干扰。通过迭代，以初始的点扩散函数为基础，根据评价函数，初步变化点扩散函数，从而得到最终的点扩散函数。这样的迭代，可以部分地消除非共光路象差对于最终图像的影响。

附图说明

图1是本发明的电路连接示意图；

图2是本发明的流程示意图。

具体实施方式

本发明通过对于望远镜，自适应光学的控制和测量信息进行记录，并利用记录信息进行点扩散函数的重建和恢复图象尺寸的选择。之后利用重建的点扩散函数作为初值，利用图像的全变分范数和与恢复后图像的最大相似度作为评价函数，经过优化获得质量提高的图像。该图像盲卷积方法使用了更加接近实际的点扩散函数作为初值，通过迭代部分的消除了非共光路象差和自适应光学系统残差对于图象的影响。整个图像盲卷积方法收敛速度快，结果可靠。

基于自适应光学系统点扩散函数重建的图像盲卷积方法，包括

时钟电路：为望远镜系统提供时钟服务；

第一存储器：连接望远镜系统主镜面形测量信号传入接口，实时记录望远镜主镜面形状随时间变化的信号φ_M(t)；

第二存储器：连接自适应光学系统波前探测器探测信号传入接口，实时记录波前探测器随时间变化的测量信号φ_wfs(t)；

第三存储器：连接望远镜科学相机曝光信号传入接口，实时记录天文图像曝光时间T,(T为一个离散数据)；

第四存储器：连接望远镜系统的自适应光学系统控制信号连接接口，实时记录变形镜(变形镜是自适应光学系统的一部分)随时间变化的控制信号φ_dm(t)；

第五存储器：连接望远镜科学相机，用来记录曝光图像img(t)；

核心处理器：连接第一存储器、第二存储器、第三存储器、第四存储器、第五存储器和时钟电路；

核心处理器控制存储器进行信号的存取工作。

核心处理器控制存储器的工作步骤如下：

步骤一，根据科学需求，选择需要处理的图象的曝光时刻t₀(t₀为T的一个值)和曝光时长Δt(t₀为T的一段时间)，从第一存储器中读取φ_M(Δt)，从第二存储器读取φ_wfs(Δt)，从第四存储器读取φ_dm(Δt)，从第五存储器中读取img(t)，根据φ_wfs(Δt)和φ_dm(Δt)计算出大气相干长度r₀和外尺度大小L₀随时间变化的序列，通过统计选择概率最大的相干长度和外尺度，计算出大气湍流的相干角度，从img(t)里选择包含感兴趣的目标大小为相干角度之内的原始图像subimg(Δt)作为初始的需要处理的图像，感兴趣的目标是指预备要处理的目标，是指人对一块图像感觉有研究意义，需要仔细观察；

步骤二，主镜影响的光学传递函数OTF_static(Δt)＝A·exp(-i·φ_M(Δt))·P·OTF_sec(Δt)，A为描述主镜的口径和拼接结构的光学传递函数，P根据主镜和副镜之间结构确定的光线自由传递函数，OTF_sec(Δt)为副镜口径和副镜支架确定的副镜光学传递函数；

步骤三，自适应光学系统的矫正残差的光学传递函数OTF_dynamical(Δt)＝OTF_wfs(Δt)·OTF_dm(Δt)，波前探测器测量残差光学传递函数OTF_wfs(Δt)＝exp(-i·φ_wfs(Δt))，变形镜响应函数光学传递函数OTF_dm(Δt)＝exp(-i·φ_dm(Δt)·h_dm)，变形镜响应函数h_dm是由变形镜物理结构决定的，一般由厂家直接给出，总光学传递函数OTF(Δt)＝OTF_static(Δt)·OTF_dynamical(Δt)；

步骤四，利用MATLAB中的函数OTF2PSF，把光学传递函数计算为系统点扩散函数PSF(Δt)＝OTF2PSF{OTF(Δt)}，Nor()为归一化函数，将点扩散函数总能量归一化；

步骤五，计算初始图像的整体梯度范数和确定优化的参数λ的初始值λ0＝0.2；

步骤六，利用点扩散函数的初始值PSF_all，可以通过维纳逆滤波得到经过该点扩散函数重建之后的图像img_all(Δt)；

步骤七，计算此时的评价函数值：

$F\left({\mathrm{img}}_{all}\right(\mathrm{\Δ}t),{\mathrm{PSF}}_{all})={\mathrm{\λ}}_0||{\mathrm{PSF}}_{all}\⊗{\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right)-subimg\left(\mathrm{\Δ}t\right)|{|}^2-log\left(\mathrm{\Σ}\right|\▿\left({\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right)\right)\left|\right);$

步骤八，根据评价函数值，随机地改变点扩散函数的初始值PSF_all局部的大小(改变点扩散函数中高频部分的灰度值，以达到改变点扩散函数形状的效果)，得到点扩散函数PSF₁和利用此点扩散函数通过维纳逆滤波得到图像img₁(Δt)，计算此时的评价函数值：

$F\left({\mathrm{img}}_1\right(\mathrm{\Δ}t),{\mathrm{PSF}}_1)={\mathrm{\λ}}_0\left|\right|{\mathrm{PSF}}_1\⊗{\mathrm{img}}_1\left(\mathrm{\Δ}t\right)-subimg\left(\mathrm{\Δ}t\right)|{|}^2-log\left(\mathrm{\Σ}\right|\▿\left({\mathrm{img}}_1\left(\mathrm{\Δ}t\right)\right)\left|\right)$

步骤九，比较img_all(Δt)和img₁(Δt)的整体梯度范数和之差，若大于0.001，则改变优化参数λ的大小，若整体梯度范数和之差增加，则减小0.01，否则增加0.01，返回步骤八反复迭代，直到整体梯度范数和之差小于等于0.001。

我们定义对比度B：式子中，I_max和I_min为图像的(灰度值)最大值和最小值。计算源图像和复原图像可以得到，处理前图像对比度为0.3927，处理之后图像对比度为0.9824。综上所述，本发明能够有效的提高自适应光学系统的成像质量。

Claims (1)

1.基于自适应光学系统点扩散函数重建的图像盲卷积方法，其特征在于：包括时钟电路：为望远镜系统提供时钟服务；

第一存储器：连接望远镜系统主镜面形测量信号传入接口，实时记录望远镜主镜面形状随时间变化的信号φ_M(t)；

第二存储器：连接自适应光学系统波前探测器探测信号传入接口，实时记录波前探测器随时间变化的测量信号φ_wfs(t)；

第三存储器：连接望远镜科学相机曝光信号传入接口，实时记录天文图像曝光时间T；

第四存储器：连接望远镜系统的自适应光学系统控制信号连接接口，实时记录变形镜(变形镜是自适应光学系统的一部分)随时间变化的控制信号φ_dm(t)；

第五存储器：连接望远镜科学相机，用来记录曝光图像img(t)；

核心处理器：连接第一存储器、第二存储器、第三存储器、第四存储器、第五存储器和时钟电路；

核心处理器控制存储器进行信号的存取工作，具体工作步骤如下：

步骤一，根据科学需求，选择需要处理的图象的曝光时刻t₀和曝光时长Δt，从第一存储器中读取φ_M(Δt)，从第二存储器读取φ_wfs(Δt)，从第四存储器读取φ_dm(Δt)，从第五存储器中读取img(t)，根据φ_wfs(Δt)和φ_dm(Δt)计算出大气相干长度r₀和外尺度大小L₀随时间变化的序列，通过统计选择概率最大的相干长度和外尺度，计算出大气湍流的相干角度，从img(t)里选择包含感兴趣的目标大小为相干角度之内的原始图像subimg(Δt)作为初始的需要处理的图像；

步骤二，计算主镜影响的光学传递函数

OTF_static(Δt)＝A·exp(-i·φ_M(Δt))·P·OTF_sec(Δt)，A为描述主镜的口径和拼接结构的光学传递函数，P根据主镜和副镜之间距离确定的光线自由传递函数，OTF_sec(Δt)为副镜口径和副镜支架确定的副镜光学传递函数；

步骤三，计算自适应光学系统的矫正残差的光学传递函数OTF_dynamical(Δt)＝OTF_wfs(Δt)·OTF_dm(Δt)，波前探测器测量残差光学传递函数OTF_wfs(Δt)＝exp(-i·φ_wfs(Δt))，变形镜响应函数光学传递函数OTF_dm(Δt)＝exp(-i·φ_dm(Δt)·h_dm)，变形镜响应函数h_dm是由变形镜物理结构决定，OTF(Δt)＝OTF_static(Δt)·OTF_dynamical(Δt)；

步骤四，利用MATLAB中的函数OTF2PSF，把光学传递函数计算为系统点扩散函数PSF(Δt)＝OTF2PSF{OTF(Δt)}，Nor()为归一化函数，将点扩散函数总能量归一化；

步骤五，计算初始图像的整体梯度范数和确定优化的参数λ的初始值λ₀＝0.2；步骤六，利用点扩散函数的初始值PSF_all，可以通过维纳逆滤波得到经过该点扩散函数重建之后的图像img_all(Δt)；

步骤七，计算此时的评价函数值：

$F\left({\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right),{\mathrm{PSF}}_{all}\right)={\mathrm{\λ}}_0\left|\right|{\mathrm{PSF}}_{all}\⊗{\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right)-subimg\left(\mathrm{\Δ}t\right)|{|}^2-\log \left(\mathrm{\Σ}|\▿\left({\mathrm{img}}_{all}\left(\mathrm{\Δ}t\right)\right)|\right);$

步骤八，根据评价函数值，随机地改变点扩散函数的初始值PSF_all局部的大小，即改变点扩散函数中高频部分的灰度值，以达到改变点扩散函数形状的效果，得到点扩散函数PSF₁和利用此点扩散函数通过维纳逆滤波得到图像img₁(Δt)，计算此时的评价函数值：

$F\left({\mathrm{img}}_1\left(\mathrm{\Δ}t\right),{\mathrm{PSF}}_1\right)={\mathrm{\λ}}_0\left|\right|{\mathrm{PSF}}_1\⊗{\mathrm{img}}_1\left(\mathrm{\Δ}t\right)-subimg\left(\mathrm{\Δ}t\right)|{|}^2-\log \left(\mathrm{\Σ}|\▿\left({\mathrm{img}}_1\left(\mathrm{\Δ}t\right)\right)|\right);$ 步骤九，比较img_all(Δt)和img₁(Δt)的整体梯度范数和之差，若大于0.001，则改变优化参数λ的大小，若整体梯度范数和之差增加，则减小0.01，否则增加0.01，返回步骤八反复迭代，直到整体梯度范数和之差小于等于0.001。