CN105184762A - 同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法 - Google Patents

Abstract

本申请公开了一种同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，包括步骤：读取同步轨道卫星时间欠采样的低分辨率序列图像IL1；利用已有的低分辨率图像IL1序列进行建模，得到符合其基本变化规律的模型方程，通过得到的模型方程推导出卡尔曼滤波状态方程和量测方程，根据离散卡尔曼滤波基本方程进行预测，得到时间信息充分的同步轨道卫星姿态变化的低分辨率图像序列IL2；采用泰勒级数法对所述低分辨率图像序列IL2运动估计结果进行图像配准；采用时间信息标志的凸集投影POCS算法对低分辨率图像序列IL2进行空域超分辨率重建。本发明利用欠采样量测信息约束下的低分辨率图像序列重建出超分辨率的高清图像，解决了姿态量测信息欠采样的影响。

Description

同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法

技术领域

本发明涉及一种基于星载画幅式相机姿态量测信息欠采样条件下的多画面超分辨率图像重建技术的方法，适用于提高空间遥感成像分辨率，提升影像质量的遥感应用领域，尤其适用于地球同步轨道卫星的凝视成像模式，在获取的姿态量测信息更新速率低于相机输出画面速率的情况下，利用欠采样量测信息约束下的低分辨率图像序列重建出超分辨率的高清图像，具体地说，是涉及一种同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法。

背景技术

地球同步轨道卫星搭载的画幅式相机采用对地凝视成像的观测机制，具有对同一地区持续获取图像的能力，时间分辨率大大提高。然而，在图像获取过程中，诸多因素如光学系统的像差、大气的扰动、目标与成像器件之间的相对运动、离散采样和系统噪声等等会导致影像退化而降低了分辨率。而且受CCD成像器件像元尺寸的限制，离散采样间距过大造成的欠采样效应也会导致图像空间分辨力的降低。如何在不改变现有探测器的前提下提高光电成像系统整体指标，获得具有更高空间分辨率的图像，是超分辨重建技术旨在解决的问题。利用多幅低分辨率图像，可以有效降低空间离散采样间距，实现图像在空间上的过采样，从而重建出更高分辨率的图像。该技术可以在不改变现有成像设备的前提下，克服成像系统的固有分辨率局限，提高现有数字图像处理系统的性能，具有经济实用的特点，应用前景十分广阔。

传统的图像过采样技术主要围绕着单帧图像的最邻近插值、双线性插值、立方样条插值等方法。其仅利用了图像的空间信息，改变了图像的视觉效果，但没有提高图像的信息量，得到的图像分辨率有很大的限制，并不能满足实际应用要求。所以目前的研究主要集中在多帧图像的重建，但是如果仅仅只是对光电系统拍摄的序列图像进行处理，序列图像之间的相对位移具有不确定性，这样重建出来的图像，其分辨率并不能有很大的改观。所以应用新的信息获取及处理技术实现超分辨成像，在光学系统参数和探测器基本不变的情况下，提高光电成像系统分辨能力和识别概率的方法成为当前最具潜力的技术手段。其有效途径之一是通过特定的光电成像手段获取内容相关但存在相对位移的序列图像，通过有针对性软件重建算法，重建出更高分辨率的图像。这样重建的图像，包含了多帧图像序列的信息，大大增加了图像信息内容，因此可以真正实现在现有的探测器及光学系统不变的前提下，有效提高了成像分辨力。

图1为当今主要使用的超分辨率图像重建方法步骤，而空域超分辨率重建算法中的大多算法(如迭代反向投影、凸集投影、最大后验概率估计等)把插值和消除模糊与噪声这两个环节综合成了一个过程。

当今基于低分辨率图像序列的超分辨重建技术需要对每幅图像估计相对于基准图像的运动位移，普通的航摄相机需要在拍摄每一幅图片的时候都给出当时的相机姿态信息，用于帧间相对位移估计和后续的几何校正等处理。但对于同步轨道卫星的姿态量测信息的获取时间间隔远长于拍摄的时间间隔，相对于输出图像帧频处于欠采样状态，无法采用传统航摄相机的处理手段来进行超分辨率重建。

发明内容

为解决上述技术问题，本发明提供了一种同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，包括步骤：

读取同步轨道卫星时间欠采样的低分辨率序列图像I_L1；

视同步轨道卫星姿态变化为简谐运动，利用已有的低分辨率图像I_L1序列进行建模，得到符合其基本变化规律的模型方程，然后通过得到的模型方程推导出卡尔曼滤波状态方程和量测方程，接着根据离散卡尔曼滤波基本方程进行预测，得到时间信息充分的同步轨道卫星姿态变化的低分辨率图像序列I_L2；

采用泰勒级数法对所述低分辨率图像序列I_L2运动估计结果进行图像配准；

采用时间信息标志的凸集投影POCS算法对低分辨率图像序列I_L2进行空域超分辨率重建。

优选地，所述步骤视同步轨道卫星姿态变化为简谐运动，利用已有的低分辨率图像I_L1序列进行建模，得到符合其基本变化规律的模型方程，然后通过得到的模型方程推导出卡尔曼滤波状态方程和量测方程，接着根据离散卡尔曼滤波基本方程进行预测，得到时间信息充分的同步轨道卫星姿态变化的低分辨率图像序列I_L2，进一步为，

同步轨道卫星姿态变化为两维变化，将相机和星体设为一个完整刚体，仅考虑其振动的低频信息，根据下式建立一个反映同步轨道卫星运动变化的低阶模型，得到同步轨道卫星姿态变化引起的完整图像序列I_L2rough：

X(t)指图像x方向随时间变化发生的位移，A指图像在x方向的振幅，ω_x指图像在x方向的振动频率，指图像在x方向的相位，Y(t)指图像y方向随时间变化发生的位移，B指图像在y方向的振幅，ω_y指图像在y方向的振动频率，指图像在y方向的相位；

通过二维卡尔曼滤波预测法递推得到较高精度，

低分辨率图像离散系统的状态方程：

X(x_k,y_k)＝Φ(k|k-1)X(x_k-1,y_k-1)+G(k)w(k)

低分辨率图像离散系统的量测方程：

Z(x_k,y_k)＝H(k)X(x_k,y_k)+v(k)

其中，Z(x_k,y_k)为k时刻的量测矩阵，w(k)和v(k)分别为状态噪声和量测噪声，均为互不相关的高斯白噪声向量序列，其协方差矩阵分别为Q(k)和R(k)；Φ(k)∈R^n×n，G(k)∈R^n×p和H(k)∈R^m×n分别为状态转移矩阵，输入转移矩阵和量测转移矩阵，这里的n和p指矩阵的维度，所述量测转移矩阵指测量得到的已知值；

已知状态方程X(x_k,y_k)＝Φ(k|k-1)X(x_k-1,y_k-1)+G(k)w(k)和所述量测方程Z(x_k,y_k)＝H(k)X(x_k,y_k)+v(k)，假设现在的系统状态是k，则根据系统的模型，可以基于系统的上一状态k-1而预测出现在的状态：

X(x_k,y_k|x_k-1,y_k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)

其中X(x_k-1,y_k-1|x_k-1,y_k-1)代表上一状态最优结果。

则X(x_k,y_k|x_k-1,y_k-1)对应的协方差：

P(k|k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)Φ^T(k|k-1)+G(k-1)Q(k-1)G^T(k-1)

得到现在状态k的预测结果，再收集现在状态的量测值，结合这二者，得到现在状态k的最优化估算值X(x_k,y_k|x_k,y_k)：

X(x_k,y_k|x_k,y_k)＝X(x_k,y_k|x_k-1,y_k-1)+K(k)[Z(x_k,y_k)-H(k)X(x_k,y_k|x_k-1,y_k-1)]

其中K(k)为卡尔曼增益：

K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1

再更新k状态下X(x_k,y_k|x_k,y_k)的协方差：

P(k|k)＝(I-K(k)H)P(k|k-1)

其中I为1的矩阵，对于单模型单测量，I＝1。

优选地，所述步骤采用泰勒级数法对所述低分辨率图像序列I_L2运动估计结果进行图像配准，进一步为，

令参考图像为f，待配准的图像为g，假设相互之间的位移为a和b，将图像g通过泰勒级数展开，并忽略高次项得:

$g(x,y)=f(x+a,y+b)\≈f(x,y)+a\frac{\∂f}{\∂x}+b\frac{\∂f}{\∂b}$

定义误差函数E为:

$E=\mathrm{\Σ}\[f(x,y)+a\frac{\∂f}{\∂x}+b\frac{\∂f}{\∂b}-g(x,y)\]$

则配准的过程为：

$a,b=\arg \underset{(a,b)}{min}E.$

优选地，所述步骤采用时间信息标志的凸集投影POCS算法对低分辨率图像序列I_L2进行空域超分辨率重建，进一步为，

对于每一个低分辨率图像的像素y(m₁,m₂,t)定义如下凸约束集合：

$C_{(m_1,m_2,t)}=\{f(n_1,n_2,k):|r^{\left(f\right)}(m_1,m_2,t)|\≤{\mathrm{\δ}}_0(m_1,m_2,t)\}$

其中y(m₁,m₂,t)表示t时间的低分辨率图像的(m₁,m₂)像素值，f(n₁,n₂,k)表示高分辨率图像(n₁,n₂)的像素值，r^(f)(m₁,m₂,t)为与y(m₁,m₂,t)有关的残差，δ₀(m₁,m₂,t)是阈值，

$r^{\left(f\right)}(m_1,m_2,t)\≈y(m_1,m_2,t)-\underset{n_1,n_2}{\Σ}f(n_1,n_2,k)\·h(m_1,m_2;n_1,n_2,t)$

f(n₁,n₂,k)向凸集多次投影，最终可确定改点处的像素值，得到矩阵：

$\begin{array}{c}{P}_{(m_1,m_2,t)}f(n_1,n_2,k)=f(n_1,n_2,k)+\\ \left\{\begin{array}{c}\frac{r^{\left(f\right)}(m_1,m_2,k)-{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)>{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ 0,-{\mathrm{\&sigma;}}_0(m_1,m_2,t)<r^{\left(f\right)}(m_1,m_2,k)<{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ \frac{r^{\left(f\right)}(m_1,m_2,k)+{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)<-{\mathrm{\&sigma;}}_0(m_1,m_2,t)\end{array}\right.\end{array}.$

与现有技术相比，本发明所述的同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，达到了如下效果：

超分辨率图像重建成功的关键是低分辨率图像序列含有非冗余的互补信息。从统计学意义上讲，上述的互补信息称为样本信息，它们为超分辨率图像重建过程提供了数据样本，互补性的样本信息越充分，增加的信息量越大，超分辨率重建的效果越好。与普通航摄相机不同，地球同步轨道卫星搭载的画幅式相机无法获取与相机拍摄同步的姿态量测信息，因此无法按照传统相机共面方程的形式做变形校正；而且地面场景复杂，成像画幅大，没有姿态量测信息约束的盲重建方法存在着运动估计偏差大，处理时间长等缺点。本发明的优点在于利用欠采样量测信息约束下的低分辨率图像序列重建出超分辨率的高清图像，其中针对星体姿态的周期性变化规律开展简谐振动建模和利用二维卡尔曼滤波的方法完成对每帧图像运动位移的估计，解决了姿态量测信息欠采样的影响，完成了序列图像的时-空配准，保证了重建图像分辨率提升的质量。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为传统超分辨率重建的方法流程图；

图2为本发明基于同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法流程图；

图3为实施例3中图像处理前灰度直方图；

图4为实施例3中图像处理后灰度直方图。

具体实施方式

如在说明书及权利要求当中使用了某些词汇来指称特定组件。本领域技术人员应可理解，硬件制造商可能会用不同名词来称呼同一个组件。本说明书及权利要求并不以名称的差异来作为区分组件的方式，而是以组件在功能上的差异来作为区分的准则。如在通篇说明书及权利要求当中所提及的“包含”为一开放式用语，故应解释成“包含但不限定于”。“大致”是指在可接收的误差范围内，本领域技术人员能够在一定误差范围内解决所述技术问题，基本达到所述技术效果。此外，“耦接”一词在此包含任何直接及间接的电性耦接手段。因此，若文中描述一第一装置耦接于一第二装置，则代表所述第一装置可直接电性耦接于所述第二装置，或通过其他装置或耦接手段间接地电性耦接至所述第二装置。说明书后续描述为实施本发明的较佳实施方式，然所述描述乃以说明本发明的一般原则为目的，并非用以限定本发明的范围。本发明的保护范围当视所附权利要求所界定者为准。

以下结合附图对本发明作进一步详细说明，但不作为对本发明的限定。

实施例1：

本实施例提供一种同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，包括步骤：

步骤101：读取时间欠采样的低分辨率序列图像I_L1；由于同步轨道卫星获取的姿态量测信息的频率低于相机拍摄图像的频率，因此获取每幅低分辨率图像的同时，无法根据卫星姿态量测信息建立的成像公式来估算相邻两幅图像间的相对位移，从而也严重影响后续超分辨重建算法的效果。因此首先对低分辨序列图像建立相对运动模型。

步骤102：由于序列图像中相对于基准图像的偏移是由同步轨道卫星的姿态稳定度造成的，所以各个图像之间的偏移基本是亚像素级的，可以根据星体姿态变化简谐运动模型和二维卡尔曼滤波相结合的混合算法对获取低分辨率图像IL1序列中各个时间采样点的相对位移进行运动估计，以消除量测信息的欠采样效应，得到低分辨率图像IL2；

步骤102对序列图像间的运动估计所采用的根据星体姿态变化简谐运动模型和二维卡尔曼滤波相结合的混合算法，其基本思路为：先视卫星姿态变化为简谐运动，利用已有的低分辨率图像IL1序列进行建模，得到符合其基本变化规律的模型方程，然后通过得到的模型方程推导出卡尔曼滤波的状态方程和量测方程，接着根据离散卡尔曼滤波基本方程进行预测，得到时间信息充分的卫星姿态变化的图像序列IL2，具体方法如下：

卫星姿态变化为两维变化，因为在本申请中将相机和星体考虑为一个完整刚体，所以只考虑其振动的低频信息，先根据下式建立一个能反映其运动变化的低阶模型，得到基本的卫星姿态变化引起的完整图像序列I_L2rough：

X(t)指图像x方向随时间变化发生的位移，A指图像在x方向的振幅，ω_x指图像在x方向的振动频率，指图像在x方向的相位，Y(t)指图像y方向随时间变化发生的位移，B指图像在y方向的振幅，ω_y指图像在y方向的振动频率，指图像在y方向的相位。

根据步骤101可知已得到基本的卫星姿态变化引起的完整图像序列I_L2rough，但低阶模型精度明显存在一定的误差，需要通过二维卡尔曼滤波预测法可以不断修改预测权值的优点来递推得到较高精度。

低分辨率图像离散系统的状态方程：

X(x_k,y_k)＝Φ(k|k-1)X(x_k-1,y_k-1)+G(k)w(k)

低分辨率图像离散系统的量测方程：

Z(x_k,y_k)＝H(k)X(x_k,y_k)+v(k)

其中，w(k)和v(k)分别为状态噪声和量测噪声，均为互不相关的高斯白噪声向量序列，其协方差矩阵分别为Q(k)和R(k)；Φ(k)∈R^n×n，G(k)∈R^n×p和H(k)∈R^m×n分别为状态转移矩阵，输入矩阵和量测矩阵。

由于已知基本条件状态和量测方程，于是可继续以下的递推方程。

假设现在的系统状态是k，则根据系统的模型，可以基于系统的上一状态而预测出现在的状态：

X(x_k,y_k|x_k-1,y_k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)

其中X(x_k-1,y_k-1|x_k-1,y_k-1)代表上一状态最优结果。

则X(x_k,y_k|x_k-1,y_k-1)对应的协方差：

P(k|k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)Φ^T(k|k-1)+G(k-1)Q(k-1)G^T(k-1)

现在我们有现在状态的预测结果，再收集现在状态的量测值，结合这二者，即可得到现在状态k的最优化估算值X(x_k,y_k|x_k,y_k)：

X(x_k,y_k|x_k,y_k)＝X(x_k,y_k|x_k-1,y_k-1)+K(k)[Z(x_k,y_k)-H(k)X(x_k,y_k|x_k-1,y_k-1)]

其中K(k)为卡尔曼增益：

K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1

为了使卡尔曼滤波器不断运行下去直到系统结束，我们还要更新k状态下X(x_k,y_k|x_k,y_k)的协方差：

P(k|k)＝(I-K(k)H)P(k|k-1)

其中I为1的矩阵，对于单模型单测量，I＝1。

步骤103：多帧的低分辨率图像之间必然会存在一些相对平移、旋转、缩放的差异。这些差异可能是在拍摄中由于镜头的抖动和焦距的变化引起的，如果在这些运动不明确的情况下进行盲重建，将会使重建后的图像变得更加模糊。所以，需对步骤102获得的新的低分辨率图像序列运动估计结果进行图像配准，再将配准之后的图像进行重建才会取到更好的效果；

本步骤中所采用的配准方法为泰勒级数法，其相关函数公式为：

令参考图像为f，待配准的图像为g，假设相互之间的位移为a和b，将图像g通过泰勒级数展开,并忽略高次项得:

$g(x,y)=f(x+a,y+b)\&ap;f(x,y)+a\frac{\&part;f}{\&part;x}+b\frac{\&part;f}{\&part;b}$

定义误差函数E为:

$E=\mathrm{\&Sigma;}\&lsqb;f(x,y)+a\frac{\&part;f}{\&part;x}+b\frac{\&part;f}{\&part;b}-g(x,y)\&rsqb;$

则配准的过程就是：

$a,b=\arg \underset{(a,b)}{min}E$

步骤104：对步骤102获得的新的低分辨率图像序列进行空域超分辨率重建，采用时间信息标志的凸集投影(POCS)算法，其中的一个运算参数由步骤103获得。

本步骤中所采用时间信息标志的POCS算法。

对于每一个低分辨率图像的像素y(m₁,m₂,t)定义如下凸约束集合：

$C_{(m_1,m_2,t)}=\{f(n_1,n_2,k):|r^{\left(f\right)}(m_1,m_2,t)|\&le;{\mathrm{\&delta;}}_0(m_1,m_2,t)\}$

其中y(m₁,m₂,t)表示t时间的低分辨率图像的(m₁,m₂)像素值，f(n₁,n₂,k)表示高分辨率图像(n₁,n₂)的像素值，r^(f)(m₁,m₂,t)为与y(m₁,m₂,t)有关的残差，δ₀(m₁,m₂,t)是阈值。

$r^{\left(f\right)}(m_1,m_2,t)\&ap;y(m_1,m_2,t)-\underset{n_1,n_2}{\&Sigma;}f(n_1,n_2,k)\&CenterDot;h(m_1,m_2;n_1,n_2,t)$

f(n₁,n₂,k)向凸集多次投影，最终可确定改点处的像素值，写成矩阵形式：

$\begin{array}{c}{P}_{(m_1,m_2,t)}f(n_1,n_2,k)=f(n_1,n_2,k)+\\ \left\{\begin{array}{c}\frac{r^{\left(f\right)}(m_1,m_2,k)-{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)>{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ 0,-{\mathrm{\&sigma;}}_0(m_1,m_2,t)<r^{\left(f\right)}(m_1,m_2,k)<{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ \frac{r^{\left(f\right)}(m_1,m_2,k)+{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)<-{\mathrm{\&sigma;}}_0(m_1,m_2,t)\end{array}\right.\end{array}$

实施例2：

如图2所示，本发明的具体实施方法如下：

(1)首先读取时间欠采样的低分辨率序列图像I_L1；

(2)卫星姿态变化为两维变化，因为在本专利中将相机和星体考虑为一个完整刚体，所以只考虑其振动的低频信息，先根据下式建立一个能反映其运动变化的低阶模型，得到基本的卫星姿态变化引起的完整图像序列I_L2rough，

X(t)指图像x方向随时间变化发生的位移，A指图像在x方向的振幅，ω_x指图像在x方向的振动频率，指图像在x方向的相位，Y(t)指图像y方向随时间变化发生的位移，B指图像在y方向的振幅，ω_y指图像在y方向的振动频率，指图像在y方向的相位；

(3)因为低阶模型精度明显存在一定的误差，需要通过二维卡尔曼滤波预测法可以不断修改预测权值的优点来递推得到较高精度。

低分辨率图像离散系统的状态方程：

X(x_k,y_k)＝Φ(k|k-1)X(x_k-1,y_k-1)+G(k)w(k)低分辨率图像离散系统的量测方程：

Z(x_k,y_k)＝H(k)X(x_k,y_k)+v(k)

其中，Z(x_k,y_k)为k时刻的量测矩阵，w(k)和v(k)分别为状态噪声和量测噪声，均为互不相关的高斯白噪声向量序列，其协方差矩阵分别为Q(k)和R(k)；Φ(k)∈R^n×n，G(k)∈R^n×p和H(k)∈R^m×n分别为状态转移矩阵，输入矩阵和量测矩阵。

(4)由于已知基本条件状态(状态方程)和量测方程，于是可继续以下的递推方程。假设现在的系统状态是k，则根据系统的模型，可以基于系统的上一状态k-1而预测出现在的状态：

X(x_k,y_k|x_k-1,y_k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)

其中X(x_k-1,y_k-1|x_k-1,y_k-1)代表上一状态最优结果。

(5)则X(x_k,y_k|x_k-1,y_k-1)对应的协方差：

P(k|k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)Φ^T(k|k-1)+G(k-1)Q(k-1)G^T(k-1)(6)

现在我们有现在状态的预测结果，再收集现在状态的量测值，结合这二者，即可得到现在状态k的最优化估算值X(x_k,y_k|x_k,y_k)：

X(x_k,y_k|x_k,y_k)＝X(x_k,y_k|x_k-1,y_k-1)+K(k)[Z(x_k,y_k)-H(k)X(x_k,y_k|x_k-1,y_k-1)]

其中K(k)为卡尔曼增益：

K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1

(7)为了使卡尔曼滤波器不断运行下去直到系统结束，还要更新k状态下X(x_k,y_k|x_k,y_k)的协方差：

P(k|k)＝(I-K(k)H)P(k|k-1)

其中I为1的矩阵，对于单模型单测量，I＝1。

(8)对时间信息充分的低分辨率图像序列I_L2进行图像的空间配准，其中，采用的配准方法为泰勒级数法，其相关函数公式为：

令参考图像为f，待配准的图像为g，假设相互之间的位移为a和b，将图像g通过泰勒级数展开,并忽略高次项得:

$g(x,y)=f(x+a,y+b)\&ap;f(x,y)+a\frac{\&part;f}{\&part;x}+b\frac{\&part;f}{\&part;b}$

定义误差函数E为:

$E=\mathrm{\&Sigma;}\&lsqb;f(x,y)+a\frac{\&part;f}{\&part;x}+b\frac{\&part;f}{\&part;b}-g(x,y)\&rsqb;$

则配准的过程就是：

$a,b=\arg \underset{(a,b)}{min}E$

(9)对低分辨率图像序列IL2进行空域超分辨率重建，改进的时间信息标志的POCS算法。

对于每一个低分辨率图像的像素y(m₁,m₂,t)定义如下凸约束集合：

$C_{(m_1,m_2,t)}=\{f(n_1,n_2,k):|r^{\left(f\right)}(m_1,m_2,t)|\&le;{\mathrm{\&delta;}}_0(m_1,m_2,t)\}$

其中y(m₁,m₂,t)表示t时间的低分辨率图像的(m₁,m₂)像素值，f(n₁,n₂,k)表示高分辨率图像(n₁,n₂)的像素值，r^(f)(m₁,m₂,t)为与y(m₁,m₂,t)有关的残差，δ₀(m₁,m₂,t)是阈值。

$r^{\left(f\right)}(m_1,m_2,t)\&ap;y(m_1,m_2,t)-\underset{n_1,n_2}{\&Sigma;}f(n_1,n_2,k)\&CenterDot;h(m_1,m_2;n_1,n_2,t)$

(10)f(n₁,n₂,k)向凸集多次投影，最终可确定改点处的像素值，写成矩阵形式：

$\begin{array}{c}{P}_{(m_1,m_2,t)}f(n_1,n_2,k)=f(n_1,n_2,k)+\\ \left\{\begin{array}{c}\frac{r^{\left(f\right)}(m_1,m_2,k)-{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)>{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ 0,-{\mathrm{\&sigma;}}_0(m_1,m_2,t)<r^{\left(f\right)}(m_1,m_2,k)<{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ \frac{r^{\left(f\right)}(m_1,m_2,k)+{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)<-{\mathrm{\&sigma;}}_0(m_1,m_2,t)\end{array}\right.\end{array}$

实施例3：

在实施例1和实施例2的基础上，本实施例为应用实施例。

已知某相机在0s，5s，10s，25s，35s的拍摄图像，根据拍摄图像之间的位移关系建立一个反映其振动变化的低阶模型作为下步的状态转移矩阵，如下，

$\left\{\begin{array}{c}X\left(t\right)=5\sin (2t+3)\\ Y\left(t\right)=7sin(2t+3)\end{array}\right.$

根据卡尔曼二维滤波得到最终图像序列变化模型：

x_t+0.5281x_t-1-0.8861x_t-2-0.01886x_t-3+0.8568x_t-4+0.2776x_t-5＝ε_t+0.299ε_t-1-1.117ε_t-2+0.2802ε_t-3+0.9399ε_t-4

y_t+0.2337y_t-1-1.114y_t-2-0.9148y_t-3+0.01707y_t-4+0.7262y_t-5+0.1876y_t-6＝ε_t+0.2459ε_t-1-1.322ε_t-2-1.096ε_t-3-0.5922ε_t-4+0.9435ε_t-5+0.5864ε_t-6+0.4674ε_t-7-0.1208ε_t-8-0.6214ε_t-9

根据上面两个式子得到0s，5s，10s，15s，20s，25s，30s，35s的图像数据，然后将这几幅图进行配准，最后将配准后的图利用改进的pocs算法进行超分辨重建，得到最后高分辨率的图像。

表1中展示了对图像处理前后的量化指标对比，从表1中可以看出图像经本发明的方法处理后，方差和信噪比都有明显提高。

表1图像处理前后方差与信噪比的对比：

从图3和图4的对比中可以看出，利用本发明的方法图像处理前后的灰度值变化很大。

上述说明示出并描述了本发明的若干优选实施例，但如前所述，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (4)

1.一种同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，其特征在于，包括步骤：

读取同步轨道卫星时间欠采样的低分辨率序列图像I_L1；

视同步轨道卫星姿态变化为简谐运动，利用已有的低分辨率图像I_L1序列进行建模，得到符合其基本变化规律的模型方程，然后通过得到的模型方程推导出卡尔曼滤波状态方程和量测方程，接着根据离散卡尔曼滤波基本方程进行预测，得到时间信息充分的同步轨道卫星姿态变化的低分辨率图像序列I_L2；

采用泰勒级数法对所述低分辨率图像序列I_L2运动估计结果进行图像配准；

采用时间信息标志的凸集投影POCS算法对低分辨率图像序列I_L2进行空域超分辨率重建。

2.根据权利要求1所述的同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，其特征在于，所述步骤视同步轨道卫星姿态变化为简谐运动，利用已有的低分辨率图像I_L1序列进行建模，得到符合其基本变化规律的模型方程，然后通过得到的模型方程推导出卡尔曼滤波状态方程和量测方程，接着根据离散卡尔曼滤波基本方程进行预测，得到时间信息充分的同步轨道卫星姿态变化的低分辨率图像序列I_L2，进一步为，

同步轨道卫星姿态变化为两维变化，将相机和星体设为一个完整刚体，仅考虑其振动的低频信息，根据下式建立一个反映同步轨道卫星运动变化的低阶模型，得到同步轨道卫星姿态变化引起的完整图像序列I_L2rough：

X(t)指图像x方向随时间变化发生的位移，A指图像在x方向的振幅，ω_x指图像在x方向的振动频率，指图像在x方向的相位，Y(t)指图像y方向随时间变化发生的位移，B指图像在y方向的振幅，ω_y指图像在y方向的振动频率，指图像在y方向的相位；

通过二维卡尔曼滤波预测法递推得到较高精度，

低分辨率图像离散系统的状态方程：

X(x_k,y_k)＝Φ(k|k-1)X(x_k-1,y_k-1)+G(k)w(k)

低分辨率图像离散系统的量测方程：

Z(x_k,y_k)＝H(k)X(x_k,y_k)+v(k)

其中，Z(x_k,y_k)为k时刻的量测矩阵，w(k)和v(k)分别为状态噪声和量测噪声，均为互不相关的高斯白噪声向量序列，其协方差矩阵分别为Q(k)和R(k)；Φ(k)∈R^n×n，G(k)∈R^n×p和H(k)∈R^m×n分别为状态转移矩阵，输入转移矩阵和量测转移矩阵，这里的n和p指矩阵的维度，所述量测转移矩阵指测量得到的已知值；

已知状态方程X(x_k,y_k)＝Φ(k|k-1)X(x_k-1,y_k-1)+G(k)w(k)和所述量测方程Z(x_k,y_k)＝H(k)X(x_k,y_k)+v(k)，假设现在的系统状态是k，则根据系统的模型，

可以基于系统的上一状态k-1而预测出现在的状态：

X(x_k,y_k|x_k-1,y_k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)

其中X(x_k-1,y_k-1|x_k-1,y_k-1)代表上一状态最优结果。

则X(x_k,y_k|x_k-1,y_k-1)对应的协方差：

P(k|k-1)＝Φ(k|k-1)X(x_k-1,y_k-1|x_k-1,y_k-1)Φ^T(k|k-1)+G(k-1)Q(k-1)G^T(k-1)

得到现在状态k的预测结果，再收集现在状态的量测值，结合这二者，得到现在状态k的最优化估算值X(x_k,y_k|x_k,y_k)：

X(x_k,y_k|x_k,y_k)＝X(x_k,y_k|x_k-1,y_k-1)+K(k)[Z(x_k,y_k)-H(k)X(x_k,y_k|x_k-1,y_k-1)]

其中K(k)为卡尔曼增益：

K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1

再更新k状态下X(x_k,y_k|x_k,y_k)的协方差：

P(k|k)＝(I-K(k)H)P(k|k-1)

其中I为1的矩阵，对于单模型单测量，I＝1。

3.根据权利要求1所述的同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，其特征在于，所述步骤采用泰勒级数法对所述低分辨率图像序列I_L2运动估计结果进行图像配准，进一步为，

令参考图像为f，待配准的图像为g，假设相互之间的位移为a和b，将图像g通过泰勒级数展开，并忽略高次项得:

$g(x,y)=f(x+a,y+b)\&ap;f(x,y)+a\frac{\&part;f}{\&part;x}+b\frac{\&part;f}{\&part;b}$

定义误差函数E为:

$E=\mathrm{\&Sigma;}\&lsqb;f(x,y)+a\frac{\&part;f}{\&part;x}+b\frac{\&part;f}{\&part;b}-g(x,y)\&rsqb;$

则配准的过程为：

$a,b=\arg \underset{(a,b)}{min}E.$

4.根据权利要求1所述的同步轨道卫星姿态欠采样量测下超分辨率图像重建的方法，其特征在于，所述步骤采用时间信息标志的凸集投影POCS算法对低分辨率图像序列I_L2进行空域超分辨率重建，进一步为，

对于每一个低分辨率图像的像素y(m₁,m₂,t)定义如下凸约束集合：

$C_{(m_1,m_2,t)}=\{f(n_1,n_2,k):|r^{\left(f\right)}(m_1,m_2,t)|\&le;{\mathrm{\&delta;}}_0(m_1,m_2,t)\}$

其中y(m₁,m₂,t)表示t时间的低分辨率图像的(m₁,m₂)像素值，f(n₁,n₂,k)表示高分辨率图像(n₁,n₂)的像素值，r^(f)(m₁,m₂,t)为与y(m₁,m₂,t)有关的残差，δ₀(m₁,m₂,t)是阈值，

$r^{\left(f\right)}(m_1,m_2,t)\&ap;y(m_1,m_2,t)-\underset{n_1,n_2}{\&Sigma;}f(n_1,n_2,k)\&CenterDot;h(m_1,m_2;n_1,n_2,t)$

f(n₁,n₂,k)向凸集多次投影，最终可确定改点处的像素值，得到矩阵：

$\begin{array}{c}{P}_{(m_1,m_2,t)}f(n_1,n_2,k)=f(n_1,n_2,k)+\\ \left\{\begin{array}{c}\frac{r^{\left(f\right)}(m_1,m_2,k)-{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)>{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ 0,-{\mathrm{\&sigma;}}_0(m_1,m_2,t)<r^{\left(f\right)}(m_1,m_2,k)<{\mathrm{\&sigma;}}_0(m_1,m_2,t)\\ \frac{r^{\left(f\right)}(m_1,m_2,k)+{\mathrm{\&sigma;}}_0(m_1,m_2,t)}{\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}\underset{{\mathrm{\&sigma;}}_1}{\&Sigma;}{h}^2(m_1,m_2;{\mathrm{\&sigma;}}_1,{\mathrm{\&sigma;}}_2,t)}\&CenterDot;h(m_1,m_2;n_1,n_2,t),r^{\left(f\right)}(m_1,m_2,k)<-{\mathrm{\&sigma;}}_0(m_1,m_2,t)\end{array}\right.\end{array}.$