基于混合高斯模型与稀疏约束的PSF估计方法

技术领域

本发明涉及光学图像复原技术领域，尤其涉及一种基于混合高斯模型与稀疏约束的PSF估计方法。

背景技术

根据地物成像特性估计图像的点扩散函数Point Spread Function(PSF)的方法可最佳的反映图像的降质特性，并能实现对整幅图像降质模型的最优估计。结合图像成像的降质特性，基于实际地物成像特点估计PSF的手段非常适合应用于图像处理中降质模型的估计问题，例如：基于PSF的图像滤波去噪和图像超分辨率重建、图像复原等。

相应的，基于光学图像中地物成像特点估计PSF的方法被陆续提出。Ruiz C P提出了一种利用SPOT图像中线状目标信息及逆滤波器估计PSF并进行图像复原的方法(参考：Ruiz CP，López F J A.“Restoring SPOT images using PSF-derived deconvolutionfilters，”International Journal of Remote Sensing，vol.12，no.23，pp.2379-2391.2002)，由于该方法通过图像中线状地物实际成像估计PSF，因此有好的效果，缺点是要求图像要具有比较高的信噪比。刘正军等人提出了利用图像中的特定线状目标结合拟合的方式估计PSF的方法(参考：刘正军，王长耀，骆成凤.“CBERS-1 PSF估计与图像复原”.遥感学报，第8卷第3期，pp.234-238.2004)，首先借助图像中的特定线状目标信息(如桥梁、堤坝等)，采用正交的线扩展函数Linear Spread Function(LSF)叉积经验拟合的方式提取图像中的点扩散函数，然后利用该点扩散函数，结合频域维纳滤波器求解去图像模糊的空域反卷积算子。该方法效果较好，但该方法是在假设点扩展函数是线性可分解，且相移引起的距离向和方位向的不对称性很小的条件下估计的，并没有考虑图像实际在距离向和方位向上点扩散函数是否对称。

上述的成像方法虽然有较好的效果，但都存在局限性，且都以线状地物为初始模型的估计条件。

发明内容

针对现有技术的缺陷，本发明提供一种基于混合高斯模型与稀疏约束的PSF估计方法，利用图像中的面状地物目标进行PSF估计，实现了从成像的列和行两个方向估计2D的PSF模板，并且方法求解效率高。

一种基于混合高斯模型与稀疏约束的PSF估计方法，包括如下步骤：

步骤S1：读取整幅图像，在图像中选取成像较好、周围没有过多干扰的面状目标区域，根据成像区域的大小，选择合适的窗口截取目标区域图像；

步骤S2：读取面状目标区域图像，提取其像素值；

步骤S3：建立初始二维高斯混合模型，拟合目标区域的像素值；

步骤S4：根据行方向和列方向的初步拟合结果，更新先验分布参数，包括高斯分量均值、标准差、权重和个数，再合并或分裂高斯分量，最终生成或删除高斯分量；

步骤S5：利用得到的二维高斯混合模型，建立图像的初始估计的PSF混合高斯模板；

步骤S6：将初始估计的PSF混合高斯模板代入已知的稀疏约束图像复原模型，进一步估计PSF混合高斯模板；

步骤S7：采用迭代加权最小二乘法(IRLS，Iteratively reweighted leastsquares)求解更新初始估计的PSF混合高斯模板，并结合PSF的非负和归一约束进行调整；用迭代次数阈值N判断迭代终止条件，得到新估计的PSF混合高斯模板，在迭代中，使用共轭梯度迭代解决内部IRLS系统的低尺度精度问题；

步骤S8：稀疏约束图像复原模型求解后生成最终估计的PSF混合高斯模板。

进一步地，所述步骤S1中的输入光学图像为g(x，y)，其表达式如式(1)所示；

其中，g(x，y)为观测到的图像，h(x，y)为成像系统的综合点扩散函数即PSF，符号表示图像空间域的卷积操作，f(x，y)表示没有降质的目标图像，即待求解的图像，n(x，y)为噪声因子，截取的小面状目标区域图像为g′(x，y)，x和y为图像的行列值。

进一步地，所述步骤S2中的面状目标区域图像中的任一行或一列的像素值矩阵为X＝{F，R}^T，沿行和列方向的表达式分别为F＝(F_i；i＝1，2，…，n)和R＝(R_j；j＝1，2，…，w)，其中，F和R分别为面状目标区域的沿行和列方向的像素值向量，F_i为行方向第i个向量，R_j为列方向第j个向量，n表示行方向的向量个数，w表示列方向的向量个数。

进一步地，所述步骤S3中的二维高斯混合模型，其概率密度函数如式(2)所示；

其中，P(X)表示二维高斯混合模型的概率密度函数，N(x；ε_q，∑_q)表示子分布的概率密度函数，q＝1，2，…，Q，q表示像素值矩阵X中第几个向量，Q表示像素值矩阵X中总的向量个数，每个子分布均为二维的混合高斯概率分布，其概率密度函数如式(3)所示；

其中，ε_q表示像素值矩阵X中行向量和列向量混合高斯分布的第q个均值矩阵，ε_q＝(μ，μ′)^T，μ和μ′分别表示像素值矩阵X中行向量和列向量的混合高斯函数的均值向量，其中μ＝μ₁，μ₂，…，μ_m，μ′＝μ′₁，μ′₂，…，μ′_m′，μ_k表示μ中某一行方向的混合高斯分布的第k个均值向量，μ′_k′表示μ′中某一列方向的混合高斯分布的第k′个均值向量，k＝1，2，…，m，k′＝1，2，…，m′，∑_q表示X矩阵中行向量和列向量混合高斯分布的第q个协方差矩阵；

ρ_q表示第q个混合权重矩阵，ρ_q＝(ρ，ρ′)^T，ρ和ρ′分别表示像素值矩阵X中行向量和列向量的混合高斯函数的权重向量，其中ρ＝ρ₁，ρ₂，…，ρ_m，ρ′＝ρ′_1′，ρ′_2′，…，ρ′_m′，ρ_k表示ρ中某一行方向的混合高斯分布的第k个权重向量，ρ′_k′表示ρ′中某一列方向的混合高斯分布的第k′个权重向量，k＝1，2，…，m，k′＝1，2，…，m′，ρ_k和ρ′_k′分别满足如式(4)和式(5)所示的条件：

在行方向的混合高斯密度函数如式(6)所示；

其中，p(F)表示行方向的混合高斯密度函数，μ_k和σ_k分别表示行方向高斯分量的第k个均值和标准差，m为行方向高斯分量的个数；

在列方向的混合高斯密度函数如式(7)所示；

其中，p(R)表示列方向的混合高斯密度函数，μ′_k′和σ′_k′分别表示列方向高斯分量的第k′个均值和标准差，m′为列方向高斯分量的个数。

进一步地，所述步骤S4中的先验分布参数，其中行方向的高斯分量均值μ的联合概率密度函数如式(8)所示；

其中，p(μ)为先验分布参数中行方向的高斯分量均值μ的联合概率密度函数，L为高斯分量分布区间，μ_k(k＝1，2，…，m)为[0，L]上的独立均匀分布，将高斯分量的均值按照增序排列为0＜μ₁＜μ₂＜...＜μ_m＜L；

先验分布参数中行方向的高斯分量标准差σ的联合概率密度函数如式(9)所示；

其中，p(σ)为先验分布参数中行方向的高斯分量标准差σ的联合概率密度函数，σ＝σ₁，σ₂，…，σ_k，…，σ_m，σ_k(k＝1，2，…，m)为服从μ_σ和σ_σ的正态分布且相互独立；

先验分布参数中行方向的高斯分量权重ρ的概率密度函数如式(10)所示；

其中，p(ρ)为先验分布参数中行方向的高斯分量权重ρ的概率密度函数，ρ_k(k＝1，2，…，m)服从[0，1]上的独立均匀分布；

先验分布参数中行方向的高斯分量个数m的概率密度函数如式(11)所示；

其中，λ为行方向泊松分布的均值，2≤m≤m_max，m_max为行方向最大高斯分量数；

根据贝叶斯定理，得到在数据F条件下高斯分量的均值μ、标准差σ、权重ρ、个数m的联合后验概率密度函数如式(12)所示；

其中，y_i和分别为行方向像素值的真值和拟合值，C表示常数；

同理，得到在数据R条件下高斯分量的均值μ′、标准差σ′、权重ρ′、个数m′的联合后验概率密度函数如式(13)所示；

其中，y′_i和分别为列方向像素值的真值和拟合值，λ′为列方向泊松分布的均值，C′表示常数。

进一步地，所述步骤S5中的初始PSF混合高斯模板生成过程为，根据贝叶斯定理，计算初始PSF混合高斯模板的各个高斯分量值，具体方法为：

步骤S5.1：按步骤4中的各个参数分布公式，计算更新行方向和列方向高斯分量均值μ和μ′、高斯分量标准差σ和σ′、高斯分量权重ρ和ρ′及个数m和m′；

步骤S5.2：利用行和列方向的高斯分量参数，建立空间格网，生成新的二维高斯混合模型；

步骤S5.3：利用步骤S5.2的新的二维高斯混合模型生成图像域的初始估计的PSF混合高斯模板H。

进一步地，所述步骤S6中的稀疏约束图像复原模型，其模型函数表达式如式(14)所示；

其中，f为没有降质的目标图像，g为观测到的图像，即降质后的图像，H为初始估计的PSF混合高斯模板，符号表示卷积，γ与β是正则化参数，用于平衡各正则项之间的权重，||·||₁和||·||₂分别表示1范数和2范数。

由上述技术方案可知，本发明的有益效果在于：本发明提供的基于混合高斯模型与稀疏约束的PSF估计方法，能够快速、可靠的估计图像成像时的PSF，在处理过程中没有近似，并且主要操作在图像空间域内进行，利用图像自身成像特性估计PSF模板的初值，无需反复迭代计算，计算速度快，估计出的PSF模板能够用于图像复原操作，能够有效提升图像的质量，能实现可靠、高效的PSF估计。

附图说明

图1为本发明实施例提供的基于混合高斯模型与稀疏约束的PSF估计方法中选择的25×25方形区域样本示意图；

图2为图1所示的方形区域样本的3D示意图；

图3为本发明实施例提供的基于混合高斯模型与稀疏约束的PSF估计方法流程图；

图4为本发明实施例提供的仿真数据区域第15行方向的混合高斯模型拟合示意图；

图5为本发明实施例提供的仿真数据区域第15列方向的混合高斯模型拟合示意图；

图6为本发明实施例提供的估计出的PSF混合高斯模板2D示意图；

图7为本发明实施例提供的估计出的PSF混合高斯模板3D示意图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

从观测到的整幅图像中选择成像效果好的方形区域作为本实施例中PSF初始估计的样本，如图1和图2所示，为小窗口截取的25×25方形区域样本，从图2所示的三维示意图中能明确得知行方向和列方向的像素值。

为了对图像成像时的PSF进行估计，本实施例提供一种基于混合高斯模型与稀疏约束的PSF估计方法，如图3所示，包括步骤如下：

步骤S1：读取观测到的整幅光学图像g(x，y)，其表达式为：

其中，g(x，y)为观测到的图像，h(x，y)为成像系统的综合点扩散函数即PSF，符号表示图像空间域的卷积操作，f(x，y)表示没有降质的目标图像，即待求解的图像，n(x，y)为噪声因子，x和y为图像的行列值。

在图像g(x，y)中选取成像较好、周围没有过多干扰的面状目标区域，根据成像区域的大小，选择合适的窗口截取目标区域图像g′(x，y)，如图1所示；

步骤S2：读取目标区域图像，提取其像素值，如图2所示。面状目标区域图像中的任一行或一列的像素值矩阵为X＝{F，R}^T，沿行和列方向的表达式分别为F＝(F_i；i＝1，2，…，n)和R＝(R_j；j＝1，2，…，w)，其中，F和R分别为面状目标区域的沿行和列方向的像素值向量，F_i为行方向第i个向量，R_j为列方向第j个向量，n表示行方向的向量个数，w表示列方向的向量个数。

步骤S3：建立二维高斯混合模型，拟合目标区域的像素值。

二维高斯混合模型的概率密度函数如式(2)所示；

其中，P(X)表示二维高斯混合模型的概率密度函数，N(X；ε_q，∑_q)表示子分布的概率密度函数，q＝1，2，…，Q，q表示像素值矩阵X中第几个向量，Q表示像素值矩阵X中总的向量个数，每个子分布均为二维的混合高斯概率分布，其概率密度函数如式(3)所示；

其中，ε_q表示像素值矩阵X中行向量和列向量混合高斯分布的第q个均值矩阵，ε_q＝(μ，μ′)^T，μ和μ′分别表示像素值矩阵X中行向量和列向量的混合高斯函数的均值向量，其中μ＝μ₁，μ₂，…，μ_m，μ′＝μ′₁，μ′₂，…，μ′_m′，μ_k表示μ中某一行方向的混合高斯分布的第k个均值向量，μ′_k′表示μ′中某一列方向的混合高斯分布的第k′个均值向量，k＝1，2，…，m，k′＝1，2，…，m′，∑_q表示X矩阵中行向量和列向量混合高斯分布的第q个协方差矩阵；

ρ_q表示第q个混合权重矩阵，ρ_q＝(ρ，ρ′)^T，ρ和ρ′分别表示像素值矩阵X中行向量和列向量的混合高斯函数的权重向量，其中ρ＝ρ₁，ρ₂，…，ρ_m，ρ′＝ρ′_1′，ρ′_2′，…，ρ′_m′，ρ_k表示ρ中某一行方向的混合高斯分布的第k个权重向量，ρ′_k′表示ρ′中某一列方向的混合高斯分布的第k′个权重向量，k＝1，2，…，m，k′＝1，2，…，m′，ρ_k和ρ′_k′分别满足如式(4)和式(5)所示的条件。

在行方向的混合高斯密度函数如式(6)所示；

其中，p(F)表示行方向的混合高斯密度函数，μ_k和σ_k分别表示行方向高斯分量的第k个均值和标准差，m为行方向高斯分量的个数。

在列方向的混合高斯密度函数如式(7)所示；

其中，p(R)表示列方向的混合高斯密度函数，μ′_k′和σ′_k′分别表示列方向高斯分量的第k′个均值和标准差，m′为列方向高斯分量的个数。

本实施例中，仿真数据区域第15行方向和第15列方向的混合高斯模型拟合如图4和图5所示，由图中可知，行方向和列方向均为双峰高斯模型拟合，其中行方向拟合得到的的两个高斯分量均值为μ₁＝5.316，μ₂＝14.39，行方向混合高斯分布的协方差矩阵为∑₁＝[4.684>1′＝3.995，μ′_2′＝13.7，列方向混合高斯分布的协方差矩阵为∑′_1′＝[3.008>

步骤S4：更新先验分布参数，包括行方向和列方向高斯分量均值μ和μ′、高斯分量标准差σ和σ′、高斯分量权重ρ和ρ′及个数m和m′，再合并或分裂高斯分量，最终生成或删除高斯分量。

先验分布参数中行方向的高斯分量均值μ的联合概率密度函数如式(8)所示；

其中，p(μ)为先验分布参数中行方向的高斯分量均值μ的联合概率密度函数，L为高斯分量分布区间，μ_k(k＝1，2，…，m)为[0，L]上的独立均匀分布，将高斯分量的均值按照增序排列为0＜μ₁＜μ₂＜...＜μ_m＜L；

先验分布参数中行方向的高斯分量标准差σ的联合概率密度函数如式(9)所示；

其中，p(σ)为先验分布参数中行方向的高斯分量标准差σ的联合概率密度函数，σ＝σ₁，σ₂，…，σ_k，…，σ_m，σ_k(k＝1，2，…，m)为服从μ_σ和σ_σ的正态分布且相互独立；

先验分布参数中行方向的高斯分量权重ρ的概率密度函数如式(10)所示；

其中，p(ρ)为先验分布参数中行方向的高斯分量权重ρ的概率密度函数，ρ_k(k＝1，2，…，m)服从[0，1]上的独立均匀分布；

先验分布参数中行方向的高斯分量个数m的概率密度函数如式(11)所示；

其中，λ为行方向泊松分布的均值，2≤m≤m_max，m_max为行方向最大高斯分量数；

根据贝叶斯定理，得到在数据F条件下(即行方向上)高斯分量的均值μ、标准差σ、权重ρ、个数m的联合后验概率密度函数如式(12)所示；

其中，y_i和分别为行方向像素值的真值和拟合值，C表示常数。

同理，得到在数据R条件下高斯分量的均值μ′、标准差σ′、权重ρ′、个数m′的联合后验概率密度函数如式(13)所示；

其中，y′_i和分别为列方向像素值的真值和拟合值，λ′为列方向泊松分布的均值，C′表示常数。

步骤S5：利用得到的二维高斯混合模型，建立图像的初始PSF混合高斯模板。

根据贝叶斯定理，计算初始PSF混合高斯模板的各个高斯分量值，具体方法为：

步骤S5.1：按步骤4中的各个参数分布公式，计算更新行方向和列方向高斯分量均值μ和μ′、高斯分量标准差σ和σ′、高斯分量权重ρ和ρ′及个数m和m′；

步骤S5.2：利用行和列方向的高斯分量参数，建立空间格网，生成新的二维高斯混合模型；

步骤S5.3：利用步骤S5.2生成的新的二维高斯混合模型生成图像域的初始估计的PSF混合高斯模板H。

步骤S6：将初始PSF混合高斯模板H代入到稀疏约束图像复原模型，求解最终估计的PSF混合高斯模板。稀疏约束图像复原模型函数表达式如式(14)所示；

其中，f为没有降质的目标图像，g为观测到的图像，即降质后的图像，H为初始估计的PSF混合高斯模板，符号表示卷积，γ与β是正则化参数，用于平衡各正则项之间的权重，||·||₁和||·||₂分别表示1范数和2范数。

步骤S7：采用迭代加权最小二乘法(IRLS，Iteratively reweighted leastsquares)求解更新PSF混合高斯模板，并结合PSF的非负和归一约束进行调整；用迭代次数阈值N判断迭代终止条件，当迭代求解次数大于设定次数N时，迭代计算结束，得到新估计的PSF混合高斯模板H′矩阵。在迭代中，使用共轭梯度迭代解决内部IRLS系统的低尺度精度问题。

步骤S8：稀疏约束图像复原模型求解后，利用新估计的H′矩阵像素值生成图像域的新的PSF混合高斯模板，即最终估计的PSF混合高斯模板，如图6和图7所示。

本发明提供的基于混合高斯模型与稀疏约束的PSF估计方法，能够快速、可靠的估计图像成像时的PSF，在处理过程中没有近似，并且主要操作在图像空间域内进行，无需反复迭代计算，因此可实现可靠、高效的PSF估计。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。