专利名称：一种图像复原方法
技术领域：
本发明属于遥感图像处理领域，涉及一种基于MTF的卫星影像复原方法。该方法结合 TDI CCD成像特点，在己知成像条件较少的情况下，只通过遥感图像本身的信息提取MTF 来有效地恢复图像，进一步改善图像的清晰度，增强图像视觉效果。具有简单实用、通用 性强的特点。
背景技术：
TDICCD是当前国内外新一代的高分辨率成像传感器的主要探测器类型，虽然在物理 硬件上，TDICCD器件的应用的确提高了高分辨率卫星影像的空间分辨率，但是由于成像 过程中各种因素的干扰，卫星影像的实际分辨率并不能达到传感器最初设计的地面分辨率。 为了改善成像的清晰度和成像质量，进而提高影像的分辨能力，需要对数据进行复原处理。 对退化的模糊图像的恢复一直是遥感图像处理中的热点课题之一，具有重要的理论意义和 应用价值。
影像的大部分信息存在于影像的边缘，因此边缘的清晰度对于影像质量至关重要。图 像退化的一种表现就是边缘的扩散，传统的边缘锐化方法通过提高边缘两边的灰度对比度 来突出边缘，提高边缘的清晰度，但会使整幅图像失真。针对边缘扩散的问题，有人釆用 缩小边缘宽度来达到提高边缘清晰度的目的，Leu提出了边缘宽度细化算法(Ramp Width Reduction)。 SKPatra讨论了利用高斯函数模拟点扩展函数的去巻积方法。以及以最大似然 估计为基础的盲去巻积法。
上述方法中成像系统的点扩展函数(PSF)都是未知的，还有一类方法是假定已知PSF， 如逆滤波，维纳滤波，约束最小二乘方滤波，Lucy-Richardson (L-R)算法等。
在很多情况下，我们只能获得遥感图像本身的信息，而点扩展函数(PSF)、成像系统平 台参数、外界成像环境(如大气条件)等是未知的，且精确模拟的计算量较大，或者很难通过 模拟来精确得到。盲去巻积方法虽然不需要点扩展函数，但是它需要迭代运算且收敛性比 较差。近年来，国内外学者对于基于调制传递函数的图像复原方法的研究越来越多。
调制传递函数(Modulation Transfer Function,以下简称MTF)是光学成像系统性能的一 个重要的综合评价指标，也是图像退化模型中最重要的部分，成像系统MTF的高低直接影 响到成像质量的好坏。发射前，传感器的MTF可以在实验室采用专门仪器精确测量，但由 于发射过程中的振动及从空气中进入真空中的变化会使传感器重新聚焦，另外又受大气MTF的影响，会使它的MTF发生衰减。因此，对于在轨运行的卫星，只能寻求其它的方 法来计算MTF 了。
利用布设的地面靶标和选取的地面标志物在遥感图像上的成像信息来计算MTF，对于 遥感卫星的在轨监测，较为方便和可行。其中，代表性的主要有两类做法 一类是美国、 法国等采用的用人工布设的地面靶标的方法。这种方法主要是适用于高分辨率的星载遥感 器和机载航拍仪器。另一类是自然地物法，是以美国为代表的直接利用地面标志物，如选 用桥梁、机场等大型的地面靶标，从含有这些目标的遥感图像中直接计算MTF的方法。
在国内，针对我国相继发射的在轨运行卫星，在其遥感器性能的地面监测方面也做了 大量的实质性的工作。比如针对CCD遥感器相机，国内主要采用了"样本对比法"和"地标 测量法"。样本对比法是用已知MTF的样本图像与卫星遥感图像进行比较和判读，从而确 定遥感卫星的MTF的方法；地标测量法借鉴了国外的方法，但它适用于较高分辨率的卫星， 因为要是分辨率低的话，所需靶标将非常庞大，实施起来不现实。
针对高分辨率卫星影像的特点，对基于MTF的图像复原方法进行了研究和改进，适用 于较高分辨率卫星影像的批量处理。
发明内容
本发明所要解决的问题是针对数据量大的遥感图像，提供一种有效的图像复原方法。
该方法基于MTF理论和影像质量退化理论，利用遥感图像本身的信息采用刀刃法提取光学 系统的MTF曲线来复原图像。其基本原理是通过计算成像系统的MTF曲线在不同空间频 率处的下降程度，以此来反推真实影像在高频部分的上升程度。该方法简单、易于实现， 具有较强的通用性。
本发明提供的技术方案是 一种基于MTF的图像复原方法，包括以下步骤
一、计算调制传递函数曲线
1、边缘检测
先对每一行像素值作简单差分，即计算每相邻两个像素的灰度差，找到差值最大的位 置即为边缘的像素级位置，采用该点附近4个点的值作三次多项式曲线拟合，确定其子像
素级位置，拟合多项式为
二 ，3 + a2w2 + + "4 (1)
其中n为像素点的位置(列号)，a(n)为该像素点的灰度值；将边缘附近四个点的位置 和灰度值代入求得多项式的系数。设边缘点的子像素级位置为m，对每一行，边缘点处的 二阶导为0，艮卩
"'，(m) = 6q附+ 2a2 = 0 (2)解得m = —2a2 / 6", = -"2 /
m即为所求的边缘点的子像素级位置；
2、边缘拟合
对边缘检测结果进行线性拟合，使得边缘点都位于一条直线上； 采用最小二乘法，假设直线方程为 少=ox + 6
其中，X为行号，y为该行边缘点的子像素级位置；由最小二乘法得:
,ff 、, —ff， 、广 m 、
》,乂 -》，&
(3)
附
附
乂 '=1
乂 '=1
(4)
6 =
》,2 》,_》,乂 》,
附
'〃 '(,/
IX -》,
(5)
其中，m为边缘点的个数，x'为刀刃边缘子图的第i行，^为第i行边缘点的子像素级 位置；将步骤(1)求出的边缘点的子像素级位置和行号代入式(4)和(5)得到拟合的边 缘直线；然后将行号代入(3)式重新计算该行边缘点的子像素位置；
3、 提取边缘扩展函数ESF
根据已有的采样点信息以每行边缘点子像素位置为起点分别向左和向右进行等间隔采 样，采样间隔为像元宽度，再采用三次样条插值方法以0.01 0.1 (优选0.05)像元间隔插 值采样，包括插值点在内的采样点灰度分布就是该行的边缘扩展函数ESF，对每一行的ESF 曲线以边缘点为基准求平均值可以获得最终的ESF曲线；
4、 计算线扩展函数LSF
对边缘扩展函数在每一采样点的子像素级位置进行简单差分，即计算每相邻两个采样
点的灰度差得到线扩展函数
= — -1), "=2,3'.../ (6)
其中n为第/个采样点；
5、 计算调制传递函数
对线扩展函数进行离散傅里叶变换，取变换之后各分量的模为各频率的调制传递函数 值，并以第一个调制传递函数值为基准，作归一化处理，得到一系列调制传递函数值；将频率点以截至频率为基准作归一化处理，则截止频率为l，取频率0 1处的调制传递函数 值构成光学系统的调制传递函数曲线；
二、 构建二维调制传递函数矩阵
将水平调制传递函数列向量乘以垂直调制传递函数列向量，艮[l-
MrF( v) = MTF x MTFV ⑦
式中，MTFu是在频率u处水平的调制传递函数值；MTFv是在频率v处垂直的调制传 递函数值；MTF(u,v)是二维频率坐标为(u,v)处的调制传递函数值；
取水平与垂直方向0.5频率处调制传递函数值的平均值再衰减90%作为45度方向0.5 频率处的调制传递函数值；再根据水平与垂直方向的调制传递函数向量之间的比例关系进 行插值得到二维插值调制传递函数矩阵；求出0 0.5频率处的调制传递函数值后，根据模 的对称性即可得到-0.5 0频率处的调制传递函数值；0.5频率即截止频率的一半；
三、 利用调制传递函数矩阵的图像复原
在频率域中，直接利用调制传递函数矩阵采用滤波器算子进行图像复原；其数学模型 表示为
/ (w， v) = /(w， v) xv)， (8)
其中，R(u，v)为复原图像的频谱，I(u,v)是原始图像频谱，P(u,v)是选取的滤波器算子； 采用维纳滤波法对图像进行复原，其滤波器算子为
MTPO,v)2一.
^ (9)
尸(",v)^ 1
MrF(w，v)
MTF(w，v)2 +、
其中，、=0.02是与图像信噪比有关的先验常数；MTF(u，v)为二维调制传递函数矩阵; 最后对复原图像的频谱进行反傅立叶变换即得到复原图像。
该方法的特点在于，只通过遥感图像本身的信息来提取调制传递函数，所求的调制传 递函数是光学系统、CCD器件、电子信号传输和大气扰动等各个过程的综合的调制传递函 数，并且在复原过程考虑了噪声的影响。实践表明，该方法能够能够有效地复原图像、改 善图像质量，算法效率也很高。目前，该方法已被成功应用于我国国产地面卫星预处理系 统中，实践证明了该方法的正确性、可行性和通用性。
图1为计算调制传递函数的流程示意图2为边缘扩展函数曲线三次样条插值的示意图；图3为二维调制传递函数矩阵示意图。
具体实施例方式
下面结合附图对本发明做进一歩详细描述。概括起来，本方法的实施可以分为四个步
第一步调制传递函数曲线的获取。
附图1为刀刃法计算MTF的流程图，具体步骤为
1、 边缘检测
先对每一行像素值作简单差分，即计算每相邻两个像素的灰度差，找到差值最大的位
置即为边缘的像素级位置，采用该点附近4个点的值作三次多项式曲线拟合，确定其子像
素级位置，拟合多项式为
<formula>formula see original document page 9</formula>
其中n为像素点的位置(列号)，a(n)为该像素点的灰度值；将边缘附近四个点的位置 和灰度值代入求得多项式的系数。设边缘点的子像素级位置为m，对每一行，边缘点处的 二阶导为0，艮P:<formula>formula see original document page 9</formula>
角军得<formula>formula see original document page 9</formula>m即为所求的边缘点的子像素级位置；
2、 边缘拟合
对边缘检测结果进行线性拟合，使得边缘点都位于一条直线上； 采用最小二乘法，假设直线方程为<formula>formula see original document page 9</formula>其中，X为行号，y为该行边缘点的子像素级位置；由最小二乘法得<formula>formula see original document page 9</formula>其中，m为边缘点的个数，《为刀刃边缘子图的第i行，X为第i行边缘点的子像素级 位置；将步骤(1)求出的边缘点的子像素级位置和行号代入式(4)和(5)得到拟合的边 缘直线；然后将行号代入(3)式重新计算该行边缘点的子像素位置；
3、 提取边缘扩展函数ESF
根据己有的采样点信息以每行边缘点子像素位置为起点分别向左和向右进行等间隔采 样，采样间隔为像元宽度，再采用三次样条插值方法以0.05像元间隔插值采样，如图2所 示；包括插值点在内的采样点灰度分布就是该行的边缘扩展函数ESF，对每一行的ESF曲 线以边缘点为基准求平均值可以获得最终的ESF曲线；
4、 计算线扩展函数LSF
对边缘扩展函数在每一采样点的子像素级位置进行简单差分，即计算每相邻两个采样 点的灰度差得到线扩展函数
<formula>formula see original document page 10</formula> (6)
其中w为第/个采样点；
5、 计算调制传递函数
对线扩展函数进行离散傅里叶变换，取变换之后各分量的模为各频率的调制传递函数 值，并以第一个调制传递函数值为基准，作归一化处理，得到一系列调制传递函数值；将 频率点以截至频率为基准作归一化处理，则截止频率为l，取频率0 1处的调制传递函数 值构成光学系统的调制传递函数曲线；
第二步构建二维调制传递函数矩阵
将水平调制传递函数列向量乘以垂直调制传递函数列向量，艮P:
<formula>formula see original document page 10</formula>
式中，MTFu是在频率u处水平的调制传递函数值；MTFv是在频率v处垂直的调制传 递函数值；MTF(u,v)是二维频率坐标为(u,v)处的调制传递函数值；
取水平与垂直方向0.5频率处调制传递函数值的平均值再衰减90%作为45度方向0.5 频率处的调制传递函数值；再根据水平与垂直方向的调制传递函数向量之间的比例关系进 行插值得到二维插值调制传递函数矩阵；求出0 0.5频率处的调制传递函数值后，根据模 的对称性即可得到-0.5 0频率处的调制传递函数值；此二维矩阵经过左右、上下对称变换 后，就形成了中心值为l，向四周递减的曲面。如图3所示第三步利用调制传递函数矩阵的图像复原
在频率域中，直接利用调制传递函数矩阵采用滤波器算子进行图像复原；其数学模型 表示为
i (w,v) =/(w,v)x尸(w，v)， (8)
其中，R(u,v)为复原图像的频谱，I(u，v)是原始图像频谱，P(u,v)是选取的滤波器算子； 采用维纳滤波法对图像进行复原，其滤波器算子为
1
MTF(m,v)2 MT尸O,v)2 +、
(9)
其中，~=0.02是与图像信噪比有关的先验常数；MTF(u,v)为二维调制传递函数矩阵；
最后对复原图像的频谱进行反傅立叶变换即得到复原图像-
权利要求
1.一种图像复原方法，包括以下步骤一、计算调制传递函数曲线(1)、边缘检测先对每一行像素值作简单差分，即计算每相邻两个像素的灰度差，找到差值最大的位置即为边缘的像素级位置，采用该点附近4个点的值作三次多项式曲线拟合，确定其子像素级位置，拟合多项式为a(n)＝a1n3+a2n2+a3n+a4(1)其中n为像素点的位置，a(n)为该像素点的灰度值；将边缘附近四个点的位置和灰度值代入求得多项式的系数；设边缘点的子像素级位置为m，对每一行，边缘点处的二阶导为0，即a″(m)＝6a1m+2a2＝0(2)解得m＝-2a2/6a1＝-a2/3a1m即为所求的边缘点的子像素级位置；(2)、边缘拟合对边缘检测结果进行线性拟合，使得边缘点都位于一条直线上；采用最小二乘法，假设直线方程为y＝ax+b(3)其中，x为行号，y为该行边缘点的子像素级位置；由最小二乘法得<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><mi>a</mi><mo>=</mo><mfrac> <mrow><mi>m</mi><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msub><mi>x</mi><mi>i</mi> </msub> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow><mo>-</mo><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msub><mi>x</mi><mi>i</mi> </msub> <mo>)</mo></mrow><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow> </mrow> <mrow><mi>m</mi><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msubsup><mi>x</mi><mi>i</mi><mn>2</mn> </msubsup> <mo>)</mo></mrow><mo>-</mo><msup> <mrow><mo>(</mo><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mi>m</mi></munderover><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mn>2</mn></msup> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo></mrow> </mrow>]]></math></maths><maths id="math0002" num="0002" ><math><![CDATA[ <mrow><mi>b</mi><mo>=</mo><mfrac> <mrow><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msubsup><mi>x</mi><mi>i</mi><mn>2</mn> </msubsup> <mo>)</mo></mrow><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow><mo>-</mo><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msub><mi>x</mi><mi>i</mi> </msub> <msub><mi>y</mi><mi>i</mi> </msub> <mo>)</mo></mrow><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msub><mi>x</mi><mi>i</mi> </msub> <mo>)</mo></mrow> </mrow> <mrow><mi>m</mi><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>m</mi> </munderover> <msubsup><mi>x</mi><mi>i</mi><mn>2</mn> </msubsup> <mo>)</mo></mrow><mo>-</mo><msup> <mrow><mo>(</mo><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mi>m</mi></munderover><msub> <mi>x</mi> <mi>i</mi></msub><mo>)</mo> </mrow> <mn>2</mn></msup> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中，m为边缘点的个数，xi为刀刃边缘子图的第i行，yi为第i行边缘点的子像素级位置；将步骤(1)求出的边缘点的子像素级位置和行号代入式(4)和(5)得到拟合的边缘直线；然后将行号代入(3)式重新计算该行边缘点的子像素位置；(3)、提取边缘扩展函数ESF根据已有的采样点信息以每行边缘点子像素位置为起点分别向左和向右进行等间隔采样，采样间隔为像元宽度，再采用三次样条插值方法以0.01～0.1像元间隔插值采样，包括插值点在内的采样点灰度分布就是该行的边缘扩展函数ESF，对每一行的ESF曲线以边缘点为基准求平均值可以获得最终的ESF曲线；(4)、计算线扩展函数LSF对边缘扩展函数在每一采样点的子像素级位置进行简单差分，即计算每相邻两个采样点的灰度差得到线扩展函数LSF(n)＝ESF(n)-ESF(n-1)，n＝2，3，…i(6)其中n为第i个采样点；(5)、计算调制传递函数对线扩展函数进行离散傅里叶变换，取变换之后各分量的模为各频率的调制传递函数值，并以第一个调制传递函数值为基准，作归一化处理，得到一系列调制传递函数值；将频率点以截至频率为基准作归一化处理，则截止频率为1，取频率0～1处的调制传递函数值构成光学系统的调制传递函数曲线；二、构建二维调制传递函数矩阵将水平调制传递函数列向量乘以垂直调制传递函数列向量，即MTF(u，v)＝MTFu×MTFv(7)式中，MTFu是在频率u处水平的调制传递函数值；MTFv是在频率v处垂直的调制传递函数值；MTF(u，v)是二维频率坐标为(u，v)处的调制传递函数值；取水平与垂直方向0.5频率处调制传递函数值的平均值再衰减90％作为45度方向0.5频率处的调制传递函数值；再根据水平与垂直方向的调制传递函数向量之间的比例关系进行插值得到二维插值调制传递函数矩阵；求出0～0.5频率处的调制传递函数值后，根据模的对称性即可得到-0.5～0频率处的调制传递函数值；0.5频率即截止频率的一半；三、利用调制传递函数矩阵的图像复原在频率域中，直接利用调制传递函数矩阵采用滤波器算子进行图像复原；其数学模型表示为R(u，v)＝I(u，v)×P(u，v)，(8)其中，R(u，v)为复原图像的频谱，I(u，v)是原始图像频谱，P(u，v)是选取的滤波器算子；采用维纳滤波法对图像进行复原，其滤波器算子为<maths id="math0003" num="0003" ><math><![CDATA[ <mrow><mi>P</mi><mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><mi>MTF</mi><mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo></mrow> </mrow></mfrac><mo>×</mo><mo>[</mo><mfrac> <mrow><mi>MTF</mi><msup> <mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo> </mrow> <mn>2</mn></msup> </mrow> <mrow><mi>MTF</mi><msup> <mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo> </mrow> <mn>2</mn></msup><mo>+</mo><msub> <mi>k</mi> <mi>w</mi></msub> </mrow></mfrac><mo>]</mo><mo>,</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中，kw＝0.02是与图像信噪比有关的先验常数；MTF(u，v)为二维调制传递函数矩阵；最后对复原图像的频谱进行反傅立叶变换即得到复原图像。
2.根据权利要求1所述的方法，其特征在于步骤一、(3)中采用三次样条插值方法以0.05像元间隔插值采样。
全文摘要
基于MTF的图像复原方法，包括以下步骤一、计算图像的MTF曲线；二、构建二维MTF矩阵；三、利用MTF矩阵在频率域中进行图像复原。该方法针对数据量大的遥感图像，基于MTF理论和影像质量退化理论，提供了一种有效的图像复原方法。实践表明该方法能够有效地复原图像、改善图像质量，算法效率也很高。目前，该方法已被成功应用于我国国产地面卫星预处理系统中，实践证明了该方法的正确性、可行性和通用性。
文档编号G06T5/10GK101635050SQ20091006286
公开日2010年1月27日 申请日期2009年6月26日 优先权日2009年6月26日
发明者俊 潘, 密 王, 苹 葛 申请人:武汉大学