CN105205783A - 一种联合先验估计的sar图像盲超分辨率重建方法 - Google Patents

Abstract

本发明公开了一种联合先验估计的SAR图像盲超分辨率重建方法，首先获取同场景时相相近的多帧SAR图像，对其配准并融合，形成单帧SAR图像作为重建的初始图像；然后设定噪声模型并估计其参数值，利用特征驱动法确定原始SAR图像的先验模型；接着将估计量Θ的后验分布p(Θ|g)近似为q(Θ)，通过最小化KL散度求得q⁰(f)、q¹(h)，估计h值并利用特征驱动法确定h的先验模型；最后利用KL散度和变分逼近法迭代估计出f，h，α_im和α_h。本发明在SAR图像盲重建的同时，联合先验估计，先验模型参数估计与图像重建及PSF估计相耦合迭代，所得先验符合SAR的物理统计特征，更准确地刻画图像和PSF，进而提高SAR图像盲超分辨率重建的准确性。

Description

一种联合先验估计的SAR图像盲超分辨率重建方法

技术领域

本发明涉及一种联合先验估计的SAR图像盲超分辨率重建方法，属于遥感图像处理技术领域。

背景技术

合成孔径雷达(SyntheticApertureRadar，SAR)是一种全天时、全天候的微波遥感成像雷达。在SAR成像技术中，空间分辨率作为衡量图像质量的重要指标一直都备受关注。然而，SAR所处工作环境复杂、恶劣，受到诸如硬件系统限制、平台非理想运动、成像条件不理想及系统相干斑噪声等因素影响，所获雷达图像的质量和空间分辨率难以满足科研和应用需求。

提高SAR图像的空间分辨率方法有三：一、直接改善雷达系统的天线长度、合成孔径长度和扩展信号带宽，但此法受限于硬件技术瓶颈和开发周期，且代价高；二、利用SAR成像系统的原始回波数据(非图像数据)，通过超分辨率成像方法获取高分辨率的SAR图像，但由于信息保密等因素，原始回波数据有时无法获得，只能得到成像系统处理后甚至是退化的图像数据；三、利用得到的SAR图像数据，通过超分辨率重建方法提高SAR图像的分辨率，此法实现的受限条件较少，是突破技术壁垒和信息封锁的可行方法。

对于典型的SAR图像超分辨率重建，首先，独立于图像重建之前，要确定SAR成像中的模糊函数也即点扩展函数(PSF)。获得途径有三：一、实验获取：收集一个或更多的点源图像，综合平均后，得到PSF的模型表示，但实验的开展条件难以保证；二、理论建模：根据SAR的特定物理成像、退化机理，经数学分析推导，将PSF设定为散焦、高斯或巴特沃思等模型，但特定SAR的成像、退化机理有时难以获得；三、直接凭经验假设：根据已有成像经验，直接假定PSF，但这种假定随机性大，PSF准确性低，且不随重建过程自适应调整以逼近真实值。

与典型的SAR图像超分辨率重建不同，一种将PSF估计与图像重建联合实现的盲超分辨率重建，可以不必囿于PSF实验和理论获取的局限，更省去了实验开展和理论建模的时间，适用于紧急情况下的SAR图像重建。另外，盲超分辨率重建的PSF可随重建迭代过程不断优化以逼近真实值，进而提高SAR图像超分辨率重建精度。

但是，在盲超分辨率重建过程中，图像和PSF的先验模型及其参数一般都事先假定，没有针对SAR的特点严格估计或导出，不能准确刻画SAR图像和系统PSF的统计特征，致使一般盲超分辨率难以精确重建SAR图像，给SAR图像的检测、跟踪、识别等应用造成困难，容易造成错误判断。

发明内容

现有的SAR图像盲超分辨率重建，独立设定的先验模型及其参数不能根据SAR系统和图像的自身特点刻画其统计特征，进而难以准确重建SAR图像。为了克服上述不足，本发明提供一种SAR图像盲超分辨率重建方法，该方法不仅能联合实现PSF估计与图像重建，而且在PSF估计与图像重建过程中耦合先验估计，所得先验模型及其参数更符合SAR的物理和统计特征，同时也使SAR图像重建的结果更加精确。

本发明解决其技术问题所采用的技术方案如下：

一种联合先验估计的SAR图像盲超分辨率重建方法，包括以下步骤：

1)输入同一场景时相相近的多帧SAR图像；

2)将输入的多时相SAR图像进行空间配准；

3)根据配准信息，将几何校正后的SAR图像融合初始化；

4)设定噪声模型，并确定模型参数；

5)确定原始图像的先验模型；

6)为了利用KL散度导出估计量Θ的后验分布p(Θ|g)，将p(Θ|g)近似为q(Θ)，所述估计量Θ包括原始图像f、点扩展函数h、原始图像f先验模型参数α_im，即以q(f)作为p(f|g)的近似，q(h)作为p(h|g)的近似，q(α_im)作为p(α_im|g)的近似；

KL散度的解析式为：

$C_{KL}\left(q\right(\mathrm{\Θ}\left)\right|\left|p\right(\mathrm{\Θ}|g\left)\right)=\underset{\mathrm{\Θ}}{\∫}q\left(\mathrm{\Θ}\right)\log \left(q\right(\mathrm{\Θ})/p(\mathrm{\Θ}|g\left)\right)d\mathrm{\Θ};$

7)估计点扩展函数h值，并确定点扩展函数h的先验模型：

8)以q(α_h)作为p(α_h|g)的近似，给定点扩展函数先验模型参数α_h的初始分布q⁰(α_h)为，其中，点扩展函数先验模型为高斯模型时，点扩展函数先验模型为拉普拉斯模型时，

9)利用KL散度和变分逼近法，循环迭代估计量，原始图像f、点扩展函数h、原始图像f先验模型参数α_im、点扩展函数h先验模型参数α_h，得到SAR图像超分辨率重建的最终结果输出。

前述的步骤2)进行空间配准包括，

a)对各个SAR图像进行几何校正，步骤如下：

2-1)确定输入图像上的像点坐标(x,y)和几何校正后的图像上的像点坐标(X,Y)之间的变换公式(X,Y)＝F(x,y)；

2-2)根据变换公式，逐行逐列计算输入图像上每个像素在校正后的图像上的像点坐标；

2-3)如果校正后的图像上的像点坐标不是整数值，则通过内插的方法获得整数像点位置的灰度值，最后得到几何校正图像；

b)计算经几何校正后的图像间的相对位置信息。

前述的步骤3)图像融合初始化的方法为：根据几何校正后的图像间的相对位置信息，将几何校正后的所有SAR图像置于一个控制网格内，并用插值法使网格内不均匀分布的像素达到均匀分布，形成单帧SAR图像g，作为超分辨率重建的初始图像。

前述的步骤4)的噪声模型为：

将单帧SAR图像g中噪声设定为零均值、方差为β^-1的高斯噪声，则单帧SAR图像g的条件概率分布p(g|f,h)为： $p\left(g\right|f,h)\∝{\mathrm{\β}}^{N/2}\exp \left\{-\frac{1}{2}\mathrm{\β}|\lg -h*f|{|}_2^2\right\},$

其中，f表示原始图像，h表示点扩展函数，*表示卷积运算符，N表示原始图像中的像素数，||·||₂表示2-范数运算符，

噪声模型的参数β^-1为：截取单帧SAR图像g中亮度均匀的区域，以区域中像素值的方差作为β^-1。

前述的步骤5)确定原始图像的先验模型采用特征驱动法，从高斯模型 $p_G\left(f\right|{\mathrm{\α}}_{im})={({\mathrm{\α}}_{im}/2\mathrm{\π})}^{N/2}\exp \left\{-\frac{1}{2}{\mathrm{\α}}_{im}\left|\right|Cf|{|}_2^2\right\}$ 和拉普拉斯模型 $p_L\left(f\right|{\mathrm{\α}}_{im})={({\mathrm{\α}}_{im}/2)}^N\exp \left\{-{\mathrm{\α}}_{im}\left|\right|Cf|{|}_1^1\right\}$ 中确定一种作为原始图像f的先验模型，

其中，α_im表示原始图像先验模型参数，服从均匀分布p(α_im)＝constant，C表示已给定的高通算子，||·||₂表示2-范数运算符，||·||₁表示1-范数运算符，具体步骤如下：

5-1)设原始图像f的先验模型为高斯模型且其估计为则α_im的最大似然估计为，将和代入高斯模型p_G(f|α_im)，可计算出

5-2)设原始图像f的先验模型为拉普拉斯模型且其估计为则α_im的最大似然估计为，将和代入拉普拉斯模型p_L(f|α_im)，可计算出

5-3)比较和若前者的值较大，则原始图像f的先验模型确定为高斯模型；反之，原始图像f的先验模型确定为拉普拉斯模型。

前述的步骤7)中，

通过最小化KL散度估计点扩展函数h值，并采用特征驱动法确定点扩展函数的先验模型，包括：

7-1)将h的初始值h⁰设定为一椭圆抛物面，则初始分布q⁰(h)为，q⁰(h＝h⁰)＝1，表示h＝h⁰时，q⁰(h)＝1；给定原始图像先验模型参数α_im的初始分布q⁰(α_im)为，其中，原始图像f的先验模型为高斯模型时原始图像f的先验模型为拉普拉斯模型时

7-2)通过最小化KL散度，求得原始图像f的初始分布q⁰(f)， $q^0\left(f\right)=\underset{q\left(f\right)}{\mathrm{argmin}}{C}_{KL}\left(q^0\right({\mathrm{\α}}_{im}\left)q\right(f\left)q^0\right(h\left)\right|\left|p\right({\mathrm{\α}}_{im},f,h|g\left)\right),$ 求得点扩展函数h的第一次迭代分布q¹(h)， $q^1\left(h\right)=\underset{q\left(h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^0\right({\mathrm{\α}}_{im}\left)q^0\right(f\left)q\right(h\left)\right|\left|p\right({\mathrm{\α}}_{im},f,h|g\left)\right),$ 进而以h的均值作为初始估计h¹，其中，E(·)表示均值运算符；

7-3)采用特征驱动法确定点扩展函数h的先验模型，从高斯模型 $p_G\left(h\right|{\mathrm{\α}}_h)={({\mathrm{\α}}_h/2\mathrm{\π})}^{M/2}\exp \left\{-\frac{1}{2}{\mathrm{\α}}_h\left|\right|Ch|{|}_2^2\right\}$ 和拉普拉斯模型 $p_L\left(h\right|{\mathrm{\α}}_h)={({\mathrm{\α}}_h/2)}^M\exp \left\{-{\mathrm{\α}}_h\left|\right|Ch|{|}_1^1\right\}$ 中确定一种作为点扩展函数h的先验模型，其中，α_h表示点扩展函数先验模型参数，服从均匀分布p(α_h)＝constant，M表示点扩展函数h的元素数，步骤如下：

7-3-1)设点扩展函数h的先验模型为高斯模型且其估计为则α_h的最大似然估计为， ${\widehat{\mathrm{\α}}}_h=M/\left|\right|{\mathrm{Ch}}^1|{|}_2^2,$ 由此计算

7-3-2)设点扩展函数h的先验模型为拉普拉斯模型且其估计为则α_h的最大似然估计为， ${\widehat{\mathrm{\α}}}_h=M/\left|\right|{\mathrm{Ch}}^1|{|}_1^1;$ 由此计算

7-3-3)比较和若前者的值较大，则点扩展函数h的先验模型确定为高斯模型；反之，点扩展函数h的先验模型确定为拉普拉斯模型。

前述的步骤9)中，循环迭代步骤如下：

9-1)首先确保原始图像先验模型参数α_im的初始分布q⁰(α_im)，点扩展函数先验模型参数α_h的初始分布q⁰(α_h)，原始图像f的初始分布q⁰(f)，点扩展函数h的第一次迭代分布q¹(h)已设定或已求得，设迭代序数k＝1；

9-2)设定终止迭代下界ε；

9-3)根据下式求得原始图像先验模型参数α_im和点扩展函数h先验模型参数α_h的后验分布：

$q^k\left({\mathrm{\α}}_{im}\right)=\underset{q\left({\mathrm{\α}}_{im}\right)}{\mathrm{argmin}}{C}_{KL}\left(q\right({\mathrm{\α}}_{im}\left)q^{k-1}\right({\mathrm{\α}}_h\left)q^{k-1}\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\α}}_{im},{\mathrm{\α}}_h,f,h|g\left)\right),$

$q^k\left({\mathrm{\α}}_h\right)=\underset{q\left({\mathrm{\α}}_h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^{k-1}\right({\mathrm{\α}}_{im}\left)q\right({\mathrm{\α}}_h\left)q^{k-1}\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\α}}_{im},{\mathrm{\α}}_h,f,h|g\left)\right);$

9-4)将α_im估计为 ${\mathrm{\α}}_{im}^k=E{\left({\mathrm{\α}}_{im}\right)}_{q^k\left({\mathrm{\α}}_{im}\right)},$ 将α_h估计为 ${\mathrm{\α}}_h^k=E{\left({\mathrm{\α}}_h\right)}_{q^k\left({\mathrm{\α}}_h\right)};$

9-5)根据下式求得原始图像的后验分布：

$q^k\left(f\right)=\underset{q\left(f\right)}{\mathrm{argmin}}{C}_{KL}\left(q^k\right({\mathrm{\α}}_{im}\left)q^k\right({\mathrm{\α}}_h\left)q\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\α}}_{im},{\mathrm{\α}}_h,f,h|g\left)\right),$

9-6)将原始图像f估计为：

9-7)根据下式求得点扩展函数的后验分布：

$q^{k+1}\left(h\right)=\underset{q\left(h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^k\right({\mathrm{\α}}_{im}\left)q^k\right({\mathrm{\α}}_h\left)q^k\right(f\left)q\right(h\left)\right|\left|p\right({\mathrm{\α}}_{im},{\mathrm{\α}}_h,f,h|g\left)\right)$

9-8)将点扩展函数h估计为：

9-9)如果成立，则将f^k，h^k+1作为SAR图像超分辨率重建的最终结果输出；否则，令k＝k+1，返回步骤9-3)，直至 $\left|\right|f^k-f^{k-1}|{|}_2^2/|\left|f^{k-1}\right|{|}_2^2<\mathrm{\&epsiv;}$ 成立。

前述的迭代下界ε取10^-4。

本发明的有益效果是，可以在重建SAR图像和估计PSF的同时，联合进行先验估计，先验模型的确定采用特征驱动方式，且模型参数估计与盲超分辨率重建的图像重建及PSF估计相耦合迭代，所得先验针对SAR的物理统计特征，更准确地刻画图像和PSF，进而提高SAR图像盲超分辨率重建的准确性，可为SAR图像的检测、跟踪、识别等应用提供有力保障。

附图说明

图1是本发明方法的流程图；

图2是图像几何校正示意图；

图3是图像融合初始化示意图；

图4是特征驱动法确定先验模型的流程图；

图5是利用KL散度和变分逼近法的迭代流程图。

具体实施方式

现结合附图和具体实施方式对本发明作进一步详细说明。

如图1所示，本发明提供一种联合先验估计的SAR图像盲超分辨率重建方法如下：

首先获取同一场景时相相近的多帧SAR图像，对其进行空间配准，并根据配准信息将所有SAR图像融合，形成单帧SAR图像g，作为超分辨率重建的初始图像；然后设定噪声模型，确定其参数β值，并利用特征驱动法确定原始图像的先验模型；接着将估计量Θ(包括原始SAR图像f、PSFh、原始图像f先验模型参数α_im和PSFh先验模型参数α_h)的后验分布p(Θ|g)近似为q(Θ)，通过最小化KL散度求得q⁰(f)、q¹(h)，估计h值，并利用特征驱动法确定PSF的先验模型；最后，利用KL散度，通过变分逼近法，循环迭代估计出f，h，α_im，和α_h。

具体如下：

步骤1，输入同一场景时相相近的多帧SAR图像；

步骤2，配准：将输入的多时相SAR图像进行空间配准，包括，

a)对各个图像进行如图2所示的几何校正，消除图像上的几何误差、得到正射图像，操作过程如下：

2-1)确定输入图像上的像点坐标(x,y)和几何校正后的图像上的像点坐标(X,Y)之间的变换公式(X,Y)＝F(x,y)。

2-2)根据变换公式，逐行逐列计算输入图像上每个像素在校正后的图像上的像点坐标。

2-3)如果校正后的图像上的像点坐标不是整数值，则通过内插的方法获得整数像点位置的灰度值，最后得到几何校正图像。

b)计算经几何校正后的图像间的相对位置信息；

步骤3，进行如图3所示的图像融合初始化：根据位置信息，将几何校正后的所有SAR图像置于一个控制网格内，并用插值法使网格内不均匀分布的像素达到均匀分布，形成单帧SAR图像g，作为超分辨率重建的初始图像，其像素密度比输入图像的像素密度高；

步骤4，确定噪声及其参数：将单帧SAR图像g中噪声设定为零均值、方差为β^-1的高斯噪声，则单帧SAR图像g的条件概率分布p(g|f,h)为：

$p\left(g\right|f,h)\&Proportional;{\mathrm{\&beta;}}^{N/2}\exp \left\{-\frac{1}{2}\mathrm{\&beta;}\left|\right|g-h*f|{|}_2^2\right\},$

其中，f表示原始图像，h表示点扩展函数(PSF)，*表示卷积运算符，N表示原始图像中的像素数，||·||₂表示2-范数运算符，截取单帧SAR图像g中亮度较为均匀区域，以区域中像素值的方差作为β^-1；

步骤5，确定原始图像的先验模型：从高斯模型 $p_G\left(f\right|{\mathrm{\&alpha;}}_{im})={({\mathrm{\&alpha;}}_{im}/2\mathrm{\&pi;})}^{N/2}\exp \left\{-\frac{1}{2}{\mathrm{\&alpha;}}_{im}\left|\right|Cf|{|}_2^2\right\}$ 和拉普拉斯模型 $p_L\left(f\right|{\mathrm{\&alpha;}}_{im})={({\mathrm{\&alpha;}}_{im}/2)}^N\exp \left\{-{\mathrm{\&alpha;}}_{im}\left|\right|Cf|{|}_1^1\right\}$ 中确定一种作为原始图像f的先验模型，

其中，α_im表示原始图像先验模型参数，服从均匀分布p(α_im)＝constant，C表示已给定的高通算子，||·||₂表示2-范数运算符，||·||₁表示1-范数运算符，原始图像的先验模型具体确定方法为如图4所示的特征驱动法，步骤如下：

步骤5-1，设原始图像f的先验模型为高斯模型且其估计为则α_im的最大似然估计为，将和代入高斯模型p_G(f|α_im)，可计算出

步骤5-2，设原始图像f的先验模型为拉普拉斯模型且其估计为则α_im的最大似然估计为将和代入拉普拉斯模型p_L(f|α_im)，可计算出

步骤5-3，比较和若前者的值较大，则原始图像f的先验模型确定为高斯模型；反之，原始图像f的先验模型确定为拉普拉斯模型；

步骤6，为了利用KL(Kullback–Leibler)散度导出估计量Θ的后验分布p(Θ|g)，将p(Θ|g)近似为q(Θ)，即以q(f)作为p(f|g)的近似，q(h)作为p(h|g)的近似，q(α_im)作为p(α_im|g)的近似；KL散度的解析式为：

$C_{KL}\left(q\right(\mathrm{\&Theta;}\left)\right|\left|p\right(\mathrm{\&Theta;}|g\left)\right)=\underset{\mathrm{\&Theta;}}{\&Integral;}q\left(\mathrm{\&Theta;}\right)\log \left(q\right(\mathrm{\&Theta;})/p(\mathrm{\&Theta;}|g\left)\right)d\mathrm{\&Theta;};$

步骤7，估计点扩展函数h值，并据此确定PSF的先验模型：

步骤7-1，将h的初始值h⁰设定为一椭圆抛物面，以更加符合SAR的PSF特性，则初始分布q⁰(h)为，q⁰(h＝h⁰)＝1，表示h＝h⁰时，q⁰(h)＝1；给定原始图像先验模型参数α_im的初始分布q⁰(α_im)，其中，原始图像f的先验模型为高斯模型时原始图像f的先验模型为拉普拉斯模型时

步骤7-2，通过最小化KL散度，求得原始图像f的初始分布q⁰(f)， $q^0\left(f\right)=\underset{q\left(f\right)}{\mathrm{argmin}}{C}_{KL}\left(q^0\right({\mathrm{\&alpha;}}_{im}\left)q\right(f\left)q^0\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},f,h|g\left)\right),$ 求得点扩展函数h第一次迭代的分布q¹(h)， $q^1\left(h\right)=\underset{q\left(h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^0\right({\mathrm{\&alpha;}}_{im}\left)q^0\right(f\left)q\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},f,h|g\left)\right),$ 进而以h的均值作为初始估计h¹，其中，E(·)表示均值运算符；

步骤7-3，根据初始估计h¹确定PSF的先验模型，方法同步骤5：从高斯模型 $p_G\left(h\right|{\mathrm{\&alpha;}}_h)={({\mathrm{\&alpha;}}_h/2\mathrm{\&pi;})}^{M/2}\exp \left\{-\frac{1}{2}{\mathrm{\&alpha;}}_h\left|\right|Ch|{|}_2^2\right\}$ 和拉普拉斯模型中确定一种作为PSF的先验模型，其中α_h表示PSF先验模型参数，服从均匀分布p(α_h)＝constant，M表示h的元素数，具体确定方法为特征驱动法，步骤如下：

步骤7-3-1，设h的先验模型为高斯模型且其估计为则α_h的最大似然估计为， ${\widehat{\mathrm{\&alpha;}}}_h=M/\left|\right|{\mathrm{Ch}}^1|{|}_2^2,$ 由此计算

步骤7-3-2，设h的先验模型为拉普拉斯模型且其估计为则α_h的最大似然估计为， ${\widehat{\mathrm{\&alpha;}}}_h=M/\left|\right|{\mathrm{Ch}}^1|{|}_1^1,$ 由此计算

步骤7-3-3，比较和若前者的值较大，则h的先验模型确定为高斯模型；反之，h的先验模型确定为拉普拉斯模型；

步骤8，以q(α_h)作为p(α_h|g)的近似，给定PSF先验模型参数α_h的初始分布q⁰(α_h)为，其中，PSF先验模型为高斯模型时PSF先验模型为拉普拉斯模型时 ${\widehat{\mathrm{\&alpha;}}}_h=M/\left|\right|{\mathrm{Ch}}^1|{|}_1^1;$

步骤9，利用KL散度和变分逼近法，如图5所示，循环迭代估计f，h，α_im和α_h，迭代序数k＝1,2,...，具体按如下顺序循环迭代：

(1)首先确保q⁰(α_im)，q⁰(α_h)，q⁰(f)，q¹(h)已设定或已求得，设迭代序数k＝1；

(2)设定终止迭代下界ε，本发明中ε取10^-4；

(3)根据下式求得原始图像先验模型参数α_im和h先验模型参数α_h的后验分布：

$q^k\left({\mathrm{\&alpha;}}_{im}\right)=\underset{q\left({\mathrm{\&alpha;}}_{im}\right)}{\mathrm{argmin}}{C}_{KL}\left(q\right({\mathrm{\&alpha;}}_{im}\left)q^{k-1}\right({\mathrm{\&alpha;}}_h\left)q^{k-1}\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right),$

$q^k\left({\mathrm{\&alpha;}}_h\right)=\underset{q\left({\mathrm{\&alpha;}}_h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^{k-1}\right({\mathrm{\&alpha;}}_{im}\left)q\right({\mathrm{\&alpha;}}_h\left)q^{k-1}\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right);$

(4)将α_im估计为 ${\mathrm{\&alpha;}}_{im}^k=E{\left({\mathrm{\&alpha;}}_{im}\right)}_{q^k\left({\mathrm{\&alpha;}}_{im}\right)},$ 将α_h估计为 ${\mathrm{\&alpha;}}_h^k=E{\left({\mathrm{\&alpha;}}_h\right)}_{q^k\left({\mathrm{\&alpha;}}_h\right)};$

(5)根据下式求得原始图像的后验分布：

$q^k\left(f\right)=\underset{q\left(f\right)}{\mathrm{argmin}}{C}_{KL}\left(q^k\right({\mathrm{\&alpha;}}_{im}\left)q^k\right({\mathrm{\&alpha;}}_h\left)q\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right),$

(6)将原始图像f估计为：

(7)根据下式求得PSF的后验分布：

$q^{k+1}\left(h\right)=\underset{q\left(h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^k\right({\mathrm{\&alpha;}}_{im}\left)q^k\right({\mathrm{\&alpha;}}_h\left)q^k\right(f\left)q\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right)$

(8)将PSF估计为：

(9)如果成立，则将f^k，h^k+1作为SAR图像超分辨率重建的最终结果输出；否则，令k＝k+1，返回步骤(3)，直至成立。

Claims (8)

1.一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，包括以下步骤：

1)输入同一场景时相相近的多帧SAR图像；

2)将输入的多时相SAR图像进行空间配准；

3)根据配准信息，将几何校正后的SAR图像融合初始化；

4)设定噪声模型，并确定模型参数；

5)确定原始图像的先验模型；

6)为了利用KL散度导出估计量Θ的后验分布p(Θ|g)，将p(Θ|g)近似为q(Θ)，所述估计量Θ包括原始图像f、点扩展函数h、原始图像f先验模型参数α_im，即以q(f)作为p(f|g)的近似，q(h)作为p(h|g)的近似，q(α_im)作为p(α_im|g)的近似；

KL散度的解析式为：

$C_{KL}\left(q\right(\mathrm{\&Theta;}\left)\right|\left|p\right(\mathrm{\&Theta;}|g\left)\right)=\underset{\mathrm{\&Theta;}}{\&Integral;}q\left(\mathrm{\&Theta;}\right)log\left(q\right(\mathrm{\&Theta;})/p(\mathrm{\&Theta;}|g\left)\right)d\mathrm{\&Theta;};$

7)估计点扩展函数h值，并确定点扩展函数h的先验模型：

8)以q(α_h)作为p(α_h|g)的近似，给定点扩展函数先验模型参数α_h的初始分布q⁰(α_h)为，其中，点扩展函数先验模型为高斯模型时，点扩展函数先验模型为拉普拉斯模型时，

9)利用KL散度和变分逼近法，循环迭代估计量，原始图像f、点扩展函数h、原始图像f先验模型参数α_im、点扩展函数h先验模型参数α_h，得到SAR图像超分辨率重建的最终结果输出。

2.根据权利要求1所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述步骤2)进行空间配准包括，

a)对各个SAR图像进行几何校正，步骤如下：

2-1)确定输入图像上的像点坐标(x,y)和几何校正后的图像上的像点坐标(X,Y)之间的变换公式(X,Y)＝F(x,y)；

2-2)根据变换公式，逐行逐列计算输入图像上每个像素在校正后的图像上的像点坐标；

2-3)如果校正后的图像上的像点坐标不是整数值，则通过内插的方法获得整数像点位置的灰度值，最后得到几何校正图像；

b)计算经几何校正后的图像间的相对位置信息。

3.根据权利要求1所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述步骤3)图像融合初始化的方法为：根据几何校正后的图像间的相对位置信息，将几何校正后的所有SAR图像置于一个控制网格内，并用插值法使网格内不均匀分布的像素达到均匀分布，形成单帧SAR图像g，作为超分辨率重建的初始图像。

4.根据权利要求3所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述步骤4)的噪声模型为：

将单帧SAR图像g中噪声设定为零均值、方差为β^-1的高斯噪声，则单帧SAR图像g的条件概率分布p(g|f,h)为： $p\left(g\right|f,h)\&Proportional;{\mathrm{\&beta;}}^{N/2}\exp \{-\frac{1}{2}\mathrm{\&beta;}||g-h*f|{|}_2^2\},$

其中，f表示原始图像，h表示点扩展函数，*表示卷积运算符，N表示原始图像中的像素数，||·||₂表示2-范数运算符，

噪声模型的参数β^-1为：截取单帧SAR图像g中亮度均匀的区域，以区域中像素值的方差作为β^-1。

5.根据权利要求1所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述步骤5)确定原始图像的先验模型采用特征驱动法，从高斯模型 $p_G\left(f\right|{\mathrm{\&alpha;}}_{im})={({\mathrm{\&alpha;}}_{im}/2\mathrm{\&pi;})}^{N/2}\exp \{-\frac{1}{2}{\mathrm{\&alpha;}}_{im}|\left|Cf\right|{|}_2^2\}$ 和拉普拉斯模型 $p_L\left(f\right|{\mathrm{\&alpha;}}_{im})={({\mathrm{\&alpha;}}_{im}/2)}^N\exp \{-{\mathrm{\&alpha;}}_{im}|\left|Cf\right|{|}_1^1\}$ 中确定一种作为原始图像f的先验模型，

其中，α_im表示原始图像先验模型参数，服从均匀分布p(α_im)＝constant，C表示已给定的高通算子，||·||₂表示2-范数运算符，||·||₁表示1-范数运算符，具体步骤如下：

5-1)设原始图像f的先验模型为高斯模型且其估计为则α_im的最大似然估计为， ${\widehat{\mathrm{\&alpha;}}}_{im}=N/\left|\right|Cg|{|}_2^2,$ 将和代入高斯模型p_G(f|α_im)，可计算出

5-2)设原始图像f的先验模型为拉普拉斯模型且其估计为则α_im的最大似然估计为，将和代入拉普拉斯模型p_L(f|α_im)，可计算出

5-3)比较和若前者的值较大，则原始图像f的先验模型确定为高斯模型；反之，原始图像f的先验模型确定为拉普拉斯模型。

6.根据权利要求1所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述步骤7)中，

通过最小化KL散度估计点扩展函数h值，并采用特征驱动法确定点扩展函数的先验模型，包括：

7-1)将h的初始值h⁰设定为一椭圆抛物面，则初始分布q⁰(h)为，q⁰(h＝h⁰)＝1，表示h＝h⁰时，q⁰(h)＝1；给定原始图像先验模型参数α_im的初始分布q⁰(α_im)为，其中，原始图像f的先验模型为高斯模型时原始图像f的先验模型为拉普拉斯模型时

7-2)通过最小化KL散度，求得原始图像f的初始分布q⁰(f)为， $q^0\left(f\right)=\underset{q\left(f\right)}{\arg \min }{C}_{KL}\left(q^0\right({\mathrm{\&alpha;}}_{im}\left)q\right(f\left)q^0\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},f,h|g\left)\right),$ 求得点扩展函数h的第一次迭代分布q¹(h)为， $q^1\left(h\right)=\underset{q\left(h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^0\right({\mathrm{\&alpha;}}_{im}\left)q^0\right(f\left)q\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},f,h|g\left)\right),$ 进而以h的均值作为初始估计h¹，其中，E(·)表示均值运算符；

7-3)采用特征驱动法确定点扩展函数h的先验模型，从高斯模型 $p_G\left(h\right|{\mathrm{\&alpha;}}_h)={({\mathrm{\&alpha;}}_h/2\mathrm{\&pi;})}^{M/2}\exp \{-\frac{1}{2}{\mathrm{\&alpha;}}_h|\left|Ch\right|{|}_2^2\}$ 和拉普拉斯模型中确定一种作为点扩展函数h的先验模型，其中，α_h表示点扩展函数先验模型参数，服从均匀分布p(α_h)＝constant，M表示点扩展函数h的元素数，步骤如下：

7-3-1)设点扩展函数h的先验模型为高斯模型且其估计为则α_h的最大似然估计为， ${\widehat{\mathrm{\&alpha;}}}_h=M/\left|\right|{\mathrm{Ch}}^1|{|}_2^2,$ 由此计算

7-3-2)设点扩展函数h的先验模型为拉普拉斯模型且其估计为则α_h的最大似然估计为， ${\widehat{\mathrm{\&alpha;}}}_h=M/\mathrm{||}{\mathrm{Ch}}^1{\mathrm{||}}_1^1,$ 由此计算

7-3-3)比较和若前者的值较大，则点扩展函数h的先验模型确定为高斯模型；反之，点扩展函数h的先验模型确定为拉普拉斯模型。

7.根据权利要求1所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述步骤9)中，循环迭代步骤如下：

9-1)首先确保原始图像先验模型参数α_im的初始分布q⁰(α_im)，点扩展函数先验模型参数α_h的初始分布q⁰(α_h)，原始图像f的初始分布q⁰(f)，点扩展函数h的第一次迭代分布q¹(h)已设定或已求得，设迭代序数k＝1；

9-2)设定终止迭代下界ε；

9-3)根据下式求得原始图像先验模型参数α_im和点扩展函数h先验模型参数α_h的后验分布：

$q^k\left({\mathrm{\&alpha;}}_{im}\right)=\underset{q\left({\mathrm{\&alpha;}}_{im}\right)}{\mathrm{argmin}}{C}_{KL}\left(q\right({\mathrm{\&alpha;}}_{im}\left)q^{k-1}\right({\mathrm{\&alpha;}}_h\left)q^{k-1}\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right),$ $q^k\left({\mathrm{\&alpha;}}_h\right)=\underset{q\left({\mathrm{\&alpha;}}_h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^{k-1}\right({\mathrm{\&alpha;}}_{im}\left)q\right({\mathrm{\&alpha;}}_h\left)q^{k-1}\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right);$

9-4)将α_im估计为 ${\mathrm{\&alpha;}}_{im}^k=E{\left({\mathrm{\&alpha;}}_{im}\right)}_{q^k\left({\mathrm{\&alpha;}}_{im}\right)},$ 将α_h估计为 ${\mathrm{\&alpha;}}_h^k=E{\left({\mathrm{\&alpha;}}_h\right)}_{q^k\left({\mathrm{\&alpha;}}_h\right)};$

9-5)根据下式求得原始图像的后验分布：

$q^k\left(f\right)=\underset{q\left(f\right)}{\mathrm{argmin}}{C}_{KL}\left(q^k\right({\mathrm{\&alpha;}}_{im}\left)q^k\right({\mathrm{\&alpha;}}_h\left)q\right(f\left)q^k\right(h\left)\right|\left|p\right({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\left)\right),$

9-6)将原始图像f估计为：

9-7)根据下式求得点扩展函数的后验分布：

$q^{k+1}\left(h\right)=\underset{q\left(h\right)}{\mathrm{argmin}}{C}_{KL}\left(q^k\right({\mathrm{\&alpha;}}_{im}\left)q^k\right({\mathrm{\&alpha;}}_h\left)q^k\right(f\left)q\right(h)\left|\right|p\left({\mathrm{\&alpha;}}_{im},{\mathrm{\&alpha;}}_h,f,h|g\right))$

9-8)将点扩展函数h估计为：

9-9)如果成立，则将f^k，h^k+1作为SAR图像超分辨率重建的最终结果输出；否则，令k＝k+1，返回步骤9-3)，直至 $\left|\right|f^k-f^{k-1}|{|}_2^2/|\left|f^{k-1}\right|{|}_2^2<\mathrm{\&epsiv;}$ 成立。

8.根据权利要求7所述的一种联合先验估计的SAR图像盲超分辨率重建方法，其特征在于，所述迭代下界ε取10^-4。