CN111476723B - 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法 - Google Patents

一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法 Download PDF

Info

Publication number
CN111476723B
CN111476723B CN202010188037.7A CN202010188037A CN111476723B CN 111476723 B CN111476723 B CN 111476723B CN 202010188037 A CN202010188037 A CN 202010188037A CN 111476723 B CN111476723 B CN 111476723B
Authority
CN
China
Prior art keywords
pixel
pixels
strip
filling
gray
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010188037.7A
Other languages
English (en)
Other versions
CN111476723A (zh
Inventor
陈彦军
张玉红
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Normal University
Original Assignee
Harbin Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Normal University filed Critical Harbin Normal University
Priority to CN202010188037.7A priority Critical patent/CN111476723B/zh
Publication of CN111476723A publication Critical patent/CN111476723A/zh
Application granted granted Critical
Publication of CN111476723B publication Critical patent/CN111476723B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/90—Determination of colour characteristics
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20036—Morphological image processing

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Image Processing (AREA)

Abstract

一种Landsat‑7扫描线纠正器失效的遥感图像丢失像素恢复方法，它属于计算机图像处理技术领域。本发明解决了对有条带遥感图像进行丢失像素填充时引入图像中不存在的新灰度值导致填充后条带与条带周围图像之间区别明显，以及目前基于无参考图的丢失像素填充方法填充不精确、像素恢复效果不好的问题。本发明利用待填充图像中条带像素周围信息，用回归分析计算待填充像素周围各方向灰度变化趋势，在非监督分类准则下兼顾亮、暗纹理信息确定当前像素填充方向。在填充方向上使用三次样条插值，通过条带边缘以外两侧像素灰度值计算待填充像素灰度值。最后进行最大、最小值局部自适应滤波，去除突然出现的亮点和暗点。本发明可以应用于遥感图像丢失像素的恢复。

Description

一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法
技术领域
本发明属于计算机图像处理技术领域，具体涉及一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法。
背景技术
Landsat-7是美国的陆地卫星计划中的第七颗，于1999年4月15日在加利福尼亚范登堡空军基地用Delta II火箭发射。增强型专题制图仪(ETM+)是卫星携带的传感器。在卫星前进时，ETM+传感器对整个陆地卫星测绘带进行扫描。ETM+传感器的扫描线纠正器(SLC)是一个电化学装置(ECD)，可以对ETM+传感器正向和反向交叉轨迹扫描模式进行补偿以实现连续覆盖。2003年5月31日Landsat-7卫星ETM+传感器的SLC永久失效了，之后的扫描图像只有在图像中心的正下方，大约22公里宽的范围内才能看到有效的连续扫描。在此范围之外向两侧不再是平行排列，而是交替重叠。交叠扫描图像之间的条带宽度从0像素到14像素逐渐增加，从图像中间向两侧形成巨大的楔形间隙。图像中大约25％的数据丢失。SLC失效前采集的图像被称为“SLC-on”图像，而SLC失效后采集的有条带遥感图像被称为“SLC-off”图像。
自2003年6月以来，Landsat-7ETM+传感器一直采集并传输SLC失效导致的有条带数据。由于Landsat-7所采集的有条带数据仍有重要价值，学者们对Landsat-7SLC-off图像展开了各种填充丢失像素方法的研究。
遥感图像去条带方法主要分为两类，一类是无参考图方法直接恢复丢失像素。另外一类方法是有参考图恢复方法，采用多张参考图进行条带丢失像素的填充。多张参考图的方法又可分为基于Landsat-7多时段参考图方法和基于其他参考图如MODIS传感器图像的方法。而多张参考图方法中很大一部分都围绕着直方图匹配算法。此外，随着时间的推移，也出现一些其他方法，如多尺度分割、滤波或者图像修复等。
无参考图方法主要应用于MODIS传感器采集的图像中。MODIS传感器经常会产生丢失像素，主要是水平或者垂直的1个像素或者2个像素宽的条带，无参考图方法常通过建立一个滤波器进行条带去除。或者进行频率域滤波，如采用傅里叶变换，小波变换，此外也常采用直方图匹配方法。另外也常采用直接插值的方法进行填充。研究表明，基于FFT的方法经常劣于直接插值法。
在Landsat-7ETM+传感器的SLC-off图像去条带处理的无参考图方法中，多采用双线性插值、双三次插值和最近邻插值等传统插值方法。由于通常使用非条带区域来计算条带区域像素灰度，使得计算值并不是地面真实的反射率，经常会引入新的灰度值，导致填充后的条带与周围图像之间存在清晰的边界，效果不佳。
最早的多张参考图的方法来自美国联合地质调查局/美国国家航空航天局(USGS/NASA)Landsat团队，使用一个或多个SLC-off或SLC-on图像的直方图进行匹配的方法来填充条带区域。此方法应用局部线性直方图在一个移动窗口中匹配每个丢失的像素以获得重缩放功能。然后使用这个重新缩放函数来转换等效辐射测量中一个输入场景的辐射值填充空白的场景值，转换的数据用来填补场景的空白。
2018年，Liu等人提出一种基于自适应窗口线性直方图匹配(AWLHM)的方法。该算法从可用的多张参考图中通过相似性指标选择可靠的像素数据来填充目标图像的条带区域。Liu的方法除插值算法仍有改进余地外，辅助数据选择的可靠性以及相似性指标仍需利用概率指标进一步研究。2016年，Sadiq提出使用多线性回归模型在SLC-off目标图像和两个参考图像中的对应像素之间建立相关性。通过对基于多张参考图的局部线性直方图匹配方法和邻域相似像素插值方法的定量比较表明，该方法即使使用时间较远的填充图像，也能准确地估计出缺失值。
Zhang等人和Pringle等人分别于2007和2009年将地质统计学中的Kriging或Co-Kriging技术用于填充条带像素。虽然地质统计学方法对森林等同质区域，尚可获得满意的填充结果，但对于面积小于局部移动窗口尺寸的异质景观，直方图匹配的效果并不理想。另外，这些地质统计学方法计算量较大，受限于大规模的图像处理。
直方图匹配方法简单易实现，并可以解决许多丢失的数据问题。对参考图像质量要求比较高，当目标或填充图像中包含云或雪时，填充效果较差。当填充区域中小于移动窗口尺寸的区域发生重大变化时，填充效果差。根据Lark的理论，相邻像素空间上具有相关性，一般两个距离近的像素比距离远的像素灰度值更相似。使用一个或多个SLC-off或SLC-on图像的直方图匹配方法时没有充分利用图像的空间相关性信息。根据Jupp等人的理论，图像像素灰度的变化规律存在空间上的依赖性。
2011年，Jin Chen提出了称为邻域相似像素插值(NSPI)的方法，有可能准确地插值位于间隙中的像素值，特别是改善在异质景观区域的结果。2013年Chao Zeng等人在邻域相似像素插值(NSPI)方法的基础上，提出一种基于权重线性回归(WLR)的多张参考图像恢复为主的方法。目前认为是效果比较好的方法。
Roy等人2008年提出使用MODIS获取的图像为参考图来估计未扫描的条带像素值。而参考非Landsat卫星传感器的方法受制于光谱兼容性和空间分辨率。MODIS传感器具有类似光谱带，但空间分辨率要低于Landsat-7的ETM+传感器。2017年，Li提出一种使用中国环境减灾卫星(HJ-1A/1B)图像作为参考图的方法。该效果虽然低于Landsat-8卫星参考图，但由于HJ-1A/1B卫星与Landsat-7卫星采集图像可能只有数小时的间隔，可以在需要紧急恢复的时候提供一个选择。
学者们对于Landsat-7SLC-off图像去条带恢复方法的探索一直没有停止。Maxwell等人于2007年开发了一种利用多尺度分割填补Landsat-7ETM+SLC-off图像中空白的方法。该方法需要使用多幅SLC-on图像，属于多参考图方法。Bédard等人于2008年将该方法应用于一般的土地覆盖地图及视觉评估等领域。多尺度分割方法对于诸如溪流和城市道路等狭小对象具有较低的灰度预测精度。2015年，Sulong提出一种均值滤波器方法和一种距离权反比(IDW)插值方法。由于算法过于简单，条带与环境之间有明显边界，效果并不理想。2017年，Jain等人对形态学运算、均值滤波、中值滤波等各种滤波器进行了比较，认为形态学开运算滤波效果比较好，但其除处理速度比较快，填充效果并不理想。2017，Yin等人对常用的几种方法进行了比较。通过比较，Yin认为几种方法各有优缺点。对于同类区域，所有方法均能获得满意效果，对于异类区域，当使用的参考图像时间间隔较近时权重线性回归方法(WLR)与地质统计学邻域相似像素插值方法(GNSPI)方法效果类似。对于突变的场景或者缺少辅助数据的情况下，直接采样方法(DS)具有最佳性能。2019年，Miao提出一种基于参考图的图像修补(Inpainting)方法。该方法首先利用SLC-on图像建立并训练字典，进而使用一种基于低阶正则化的修复算法，并用交替方向乘子法对仿真图像进行修复。与其他几种方法比较，该方法可获得较好效果。其缺点是需要来自不同时段的大量参考图进行字典的训练。
综上，目前Landsat-7SLC-off图像条带填充以多张参考图像恢复方法为主。这些方法的优点是对于大多数条带像素，恢复计算更加准确可靠，但也存在一些很明显的缺点。对于使用Landsat-7多时段图像为参考的方法，要选取多张不同时间拍摄的卫星图像作为参考。Landsat-7卫星周期是几周，导致这些参考图像来自不同的季节、不同的观察条件。这些条件包括不同的阳光、大气、植被生长、土地湿度以及传感器状态等。这使得参考图像与待填充图像之间存在很大的差异。而且，类似人为或者地震等带来的突变导致填充图像和目标图像之间并不一致，其结果并非完全可靠。对于一些多张参考图片都无法获得参考点的待恢复像素点仍需要无参考图恢复方法。另外，搜索并选择同一地点的多时域图像也使研究者花费大量的时间，为获取最终恢复图像设置了障碍。
对于使用其他传感器获得的图片作为参考的方法，缺点除了不同传感器获得图像像素的反射率参考价值不大外，即使同一个传感器不同时间拍摄的图像也会因为光照角度、季节变化等原因影响参考价值。对于来自类似传感器的参考图像，若参考图分辨率较低，则目标图像待恢复区域的精度将被参考图分辨率拉低。而如果参考图分辨率较高，则恢复图像的任务变得没有意义。
无参考图方法进行Landsat-7卫星图像丢失像素填充时，经常会因为计算引入新的灰度值，导致填充后的条带与周围图像之间存在明显区别。同时，也因为算法不够细腻精确，恢复效果不够理想。
发明内容
本发明的目的是为解决目前基于无参考图方法对Landsat-7ETM+SLC-off遥感图像进行丢失像素填充时引入图像中不存在的新灰度值导致的填充后条带与条带周围图像之间区别明显，以及目前基于无参考图的丢失像素填充方法填充不精确、像素恢复效果不好的问题，而提出了一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法。
本发明为解决上述技术问题采取的技术方案是：一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，该方法包括以下步骤：
步骤一、从SLC-off遥感图像中初步提取出条带区域，再采用形态学闭运算对初步提取出的条带区域进行去噪处理，获得条带区域；
步骤二、对于条带区域内的任一像素(i,j)，分别从像素的位置(i,j)出发沿着各个选定方向进行搜索；
对于任意的一个选定方向，在选定方向的上方选取m个与条带上边界最相邻的非条带像素，在选定方向的下方选取m个与条带下边界最相邻的非条带像素；
将选取的与条带上边界最相邻的m个非条带像素记为Ak，k＝1,2,……,m，Ak的灰度值为ak；将选取的与条带下边界最相邻的m个非条带像素记为Bk，Bk的灰度值为bk；
以选定方向上选取的2m个像素的位置作为自变量，以选取的2m个非条带像素的灰度值为函数值进行回归分析，在该选定方向上进行回归分析计算出灰度变化的斜率值；
同理，在其他各个选定方向上进行回归分析计算出条带区域内各点的斜率值；
步骤三、在各个选定方向上，分别采用非监督分类准则，判断在选定方向上条带两边界外部的非条带像素是否为同类像素；
步骤四、计算条带区域各像素的局部平均灰度，根据步骤二计算出的在各个选定方向上的斜率值以及步骤三的判断结果确定条带区域内像素(i,j)的填充方向；
步骤五、根据条带区域内像素(i,j)是否为边界像素，采用不同的方式对像素(i,j)进行填充，获得像素填充后的条带区域；
步骤六、像素填充后，再分别对条带区域内包含的各像素进行滤波处理，滤波处理后，完成对Landsat-7扫描线纠正器失效的遥感图像丢失像素的恢复。
本发明的有益效果是：本发明提出了一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法。本发明方法不需要利用参考图，充分利用待填充图像中条带像素周围信息，用回归分析计算待填充像素周围各方向灰度变化趋势，在非监督分类准则下兼顾亮、暗纹理信息确定当前像素填充方向。在填充方向上使用三次样条插值，通过条带边缘以外两侧像素灰度值计算待填充像素灰度值。最后进行最大、最小值局部自适应滤波，去除突然出现的亮点和暗点。为避免出现原图中没有出现的颜色，插值和滤波过程中均从原图参与计算像素中选择与灰度计算值最接近的灰度代替当前像素灰度的计算值。
本发明方法的优点是无需参考图，填充灰度计算精确，能突出亮、暗纹理，而且不引入新的灰度值，避免了填充后条带与条带周围图像之间区别明显的问题。实验表明，在非监督分类时使用本发明方法得到的图像中有更少的条带恢复像素被误分类，本发明方法能够有效填充条带区域，本发明方法对条带像素的填充精确、图像条带像素恢复效果较好。
附图说明
图1是本发明的一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法的流程图；
图中，Gi,j为条带像素(i,j)填充后的灰度；
图2是形态学闭运算的结构元素图；
图3是在非监督分类准则下，判断填充方向上条带两侧是否属于同类像素的示意图；
图4是45°方向上参与线性回归分析求斜率计算的像素点示意图；
图5是含有置信区间不包括0点噪声点的残差图；
图6是当前待填充像素点三次样条插值图；
图7(a)是无条带噪声的原始遥感图像；
图7(b)是加入条带后的遥感图像；
图7(c)是对加入条带后遥感图像进行形态学闭运算的结果图；
图7(d)是对加入条带后遥感图像进行非监督分类的结果图；
图7(e)是采用LHM方法填充条带的效果图；
图7(f)是采用NSPI方法填充条带的效果图；
图7(g)是采用本发明方法填充条带的效果图；
图7(h)是采用LHM方法填充得到的图像与原始图像条带区域像素灰度的散点图；
图7(i)是采用NSPI方法填充得到的图像与原始图像条带区域像素灰度的散点图；
图7(j)是采用本发明方法填充得到的图像与原始图像条带区域像素灰度的散点图；
图7(k)是采用LHM方法得到的图像分类结果图；
图7(l)是采用NSPI方法得到的图像分类结果图；
图7(m)是采用本发明方法得到的图像分类结果图；
图8(a)是截取的400×400的部分有条带真彩色遥感图像；
图8(b)是对图8(a)中条带进行初步提取的结果图；
图8(c)是对图8(b)进行形态学闭运算后的结果图；
图8(d)是对图8(a)进行非监督分类结果图；
图8(e)是对图8(a)进行条带像素填充的结果图；
图8(f)是对图8(e)的滤波结果图；
图9(a)是图8(e)中矩形窗口1的放大图；
图9(b)是图8(e)中矩形窗口2的放大图；
图9(c)是图8(f)中矩形窗口3的放大图；
图9(d)是图8(f)中矩形窗口4的放大图；
图中，1代表矩形窗口1，2代表矩形窗口2，3代表矩形窗口3，4代表矩形窗口4；
图10(d),图10(a),图10(b),图10(c)分别是一幅Landsat-7SLC-off有条带原始图像从整体到局部逐渐放大的图；图10(e),图10(h),图10(g),图10(f)分别是图10(d),图10(a),图10(b),图10(c)条带像素填充后的结果。
具体实施方式
具体实施方式一：结合图1说明本实施方式。本实施方式所述的一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，该方法包括以下步骤：
步骤一、从SLC-off遥感图像中初步提取出条带区域，再采用形态学闭运算对初步提取出的条带区域进行去噪处理，获得条带区域；
以下步骤是针对于去噪后的条带区域进行的操作：
步骤二、对于条带区域内的任一像素(i,j)，分别从像素的位置(i,j)出发沿着各个选定方向进行搜索；
对于任意的一个选定方向，在选定方向的上方选取m个与条带上边界最相邻的非条带像素，在选定方向的下方选取m个与条带下边界最相邻的非条带像素；
将选取的与条带上边界最相邻的m个非条带像素记为Ak，k＝1,2,……,m，Ak的灰度值为ak；将选取的与条带下边界最相邻的m个非条带像素记为Bk，Bk的灰度值为bk；
以选定方向上选取的2m个像素的位置作为自变量，以选取的2m个非条带像素的灰度值为函数值进行回归分析，在该选定方向上进行回归分析计算出灰度变化的斜率值；
由于灰度变化趋势存在两种可能，斜率值有可能为负值，因此，采用斜率值的绝对值参与下述过程的比较。以下提到的斜率值均为其绝对值；
同理，在其他各个选定方向上进行回归分析计算出条带区域内各点的斜率值；
步骤三、在各个选定方向上，分别采用非监督分类准则，判断在选定方向上条带两边界外部的非条带像素是否为同类像素；
即在与条带垂直方向上，采用非监督分类准则，确定在与条带垂直方向上，与条带上边界最近邻的非条带像素和与条带下边界最近邻的非条带像素是否为同类像素，同理，进行其他选定方向上的判断；
步骤四、计算条带区域各像素的局部平均灰度，根据步骤二计算出的在各个选定方向上的斜率值以及步骤三的判断结果确定条带区域内像素(i,j)的填充方向；
步骤五、根据条带区域内像素(i,j)是否为边界像素，采用不同的方式对像素(i,j)进行填充，获得像素填充后的条带区域；
步骤六、像素填充后，再分别对条带区域内包含的各像素进行滤波处理，滤波处理后，完成对Landsat-7扫描线纠正器失效的遥感图像丢失像素的恢复。
具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤二中，分别从像素的位置(i,j)出发沿着各个选定方向进行搜索，各个选定方向分别为：垂直方向、与水平方向呈45°角方向以及与水平方向呈135°角方向。
具体实施方式三：本实施方式与具体实施方式二不同的是：所述步骤四中，计算条带区域的局部平均灰度，根据步骤二计算出的在各个选定方向上的斜率值以及步骤三的判断结果确定条带内像素(i,j)的填充方向，其具体过程为：
步骤四一、利用尺寸为17×17的窗口计算像素(i,j)的局部平均灰度G0；
步骤四二、若仅存在着一个选定方向满足条件：在选定方向上，条带两边界外部的非条带像素为同类像素，则将该选定方向作为填充方向；
步骤四三、若存在两个或两个以上的选定方向满足条件：在选定方向上，条带两边界外部的非条带像素为同类像素，则将满足条件的选定方向的斜率绝对值按照从小到大排序，并执行步骤四五；
步骤四四、若不存在选定方向满足：在选定方向上，条带两边界外部的非条带像素为同类像素，则将全部选定方向的斜率绝对值按照从小到大排序，并执行步骤四五；
步骤四五、选取出斜率绝对值最小和次最小的选定方向，根据步骤二找出的在每个选定方向上的2m个非条带像素的灰度值，分别将每个选定方向上找出的2m个非条带像素的灰度值取平均，获得斜率绝对值最小的选定方向对应的均值G1以及斜率次最小的选定方向对应的均值G2；
若满足：
|G1-G0|≥|G2-G0|
则将斜率绝对值最小的选定方向作为像素(i,j)的填充方向，否则将斜率绝对值次最小的选定方向作为像素(i,j)的填充方向。
具体实施方式四：本实施方式与具体实施方式三不同的是：所述步骤五中，根据条带区域内像素(i,j)是否为边界像素，采用不同的方式对像素(i,j)进行填充，其具体过程为：
步骤五一、若像素(i,j)所在位置是仅有1个像素宽的条带，则像素(i,j)为边界像素，在填充方向的上方或下方，随机选取出一个与像素(i,j)最相邻的非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五二、若像素(i,j)所在位置是有2个像素宽的条带，则像素(i,j)为边界像素；若像素(i,j)为上边界像素，则在填充方向的上方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；若像素(i,j)为下边界像素，则在填充方向的下方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五三、若像素(i,j)所在位置是有3个像素宽的条带，且像素(i,j)为上边界像素或下边界像素，则像素(i,j)的填充方法与步骤五二相同；否则，在像素(i,j)的填充方向的上方或者下方，随机选取出一个与像素(i,j)最相邻的非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五四、若像素(i,j)所在位置是有4个像素宽的条带，且像素(i,j)为填充方向的上边缘的2个条带像素中的一个，则在像素(i,j)的填充方向的上方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
若像素(i,j)所在位置是有4个像素宽的条带，且像素(i,j)为填充方向的下边缘的2个条带像素中的一个，则在像素(i,j)的填充方向的下方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五五、若像素(i,j)所在位置是有4个以上像素宽的条带，同时，像素(i,j)为填充方向的上边缘的2个条带像素中的一个或像素(i,j)为填充方向的下边缘的2个条带像素中的一个，则对像素(i,j)的填充方式同步骤五四；
若像素(i,j)不是填充方向的上边缘的2个条带像素中的一个且像素(i,j)也不是填充方向的下边缘的2个条带像素中的一个，则像素(i,j)为非边界像素；
根据步骤二选取出的在像素(i,j)填充方向上的2m个非条带像素进行三次样条插值，求出像素(i,j)的插值灰度P(i,j)；
Figure BDA0002414878950000091
其中：k0代表步骤二选取出的在像素(i,j)填充方向上的2m个非条带像素中的第k0个，k0＝1,2,…,2m，N(i,j,k0)为第k0个非条带像素的灰度值，N(i,j,m0)为步骤二选取出的在像素(i,j)填充方向上的2m个非条带像素的灰度值中，与P(i,j)最接近的灰度值；
将N(i,j,m0)作为像素(i,j)的灰度进行填充。
具体实施方式五：本实施方式与具体实施方式四不同的是：所述步骤二中，m的取值为3。
具体实施方式六：本实施方式与具体实施方式五不同的是：所述步骤六中，分别对条带区域内包含的各像素进行滤波处理，其具体过程为：
使用3×3的窗口对条带区域内包含的每个像素分别进行滤波处理，使窗口中心经过条带区域内的每个像素；
当窗口中心经过像素(i,j)时，则称像素(i,j)为中心像素，在3×3的窗口内，像素(i,j)的局部平均灰度为Pl；
若中心像素(i,j)的灰度(步骤五填充后的灰度)大于等于Pl(1+K)，其中K为经验值，则判断像素(i,j)为亮点，将在3×3窗口内的中心像素周围的8个像素的灰度按照从小到大排序，利用排序中的第5个值代替像素(i,j)的灰度；
若中心像素(i,j)的灰度小于等于Pl(1-K)，则判断像素(i,j)为暗点，将在3×3窗口内的中心像素周围的8个像素的灰度按照从小到大排序，利用排序中的第4个值代替像素(i,j)的灰度；
否则，不对像素(i,j)的灰度进行处理；
同理，分别对条带区域内包含的每个像素进行上述处理。
对条带区域内包含的每个像素进行上述处理，是指进行本实施方式的处理。
实施例
1.1条带填充流程
如图1所示，条带区域填充流程总体上分为预处理、回归分析、非监督分类准则下突出亮暗纹理的填充方向确定、条带边界像素点填充、条带非边界像素点插值以及对条带区域进行滤波处理五大部分。
首先，对图像进行预处理。为了填充条带区域丢失的像素，必须对条带进行精确定位。我们使用形态学图像处理方法去除条带图中的0值非条带像素，得到有效的条带区域。同时，在预处理阶段，我们还对有条带遥感图像进行非监督分类，包括条带区域共分成5类。将分类结果作为之后确定填充方向的准则。在预处理阶段，还对每个条带像素计算17×17窗口局部平均灰度，作为之后寻找局部亮纹理和暗纹理的依据。接下来对每个条带像素在垂直方向、45°方向、135°方向的条带外侧6像素灰度进行回归分析，计算斜率。随后，在非监督分类准则下，将条带像素根据三个方向条带两侧是否属于同类，来确定填充方向。此时兼顾突出亮、暗纹理。然后，在已经确定了当前像素的填充方向后，将条带像素按照是否为边界点分别处理。对于边界点在其非监督准则下的填充方向上就近随机取点填充。而对于非边界点，进行三次样条插值，计算其理想灰度，并选择填充方向上最接近的灰度代替插值灰度填充当前非边界点像素。最后使用最大值最小值滤波方法用局部窗口中的像素灰度代替灰度值过大或者过小的像素，既去除过亮和过暗的噪声点，又可以避免对图像整体平滑，从而保护图像纹理信息。下面详细介绍条带填充处理步骤。
1.2预处理阶段
1.2.1条带定位
要填充条带区域必须先对条带区域进行精确定位。观察有条带的遥感图像，可以看到条带呈楔形，从图像中部向两侧条带逐渐变宽。最窄处为1个像素宽，最宽处为14个像素宽，在遥感图像边缘处小范围内达到23个像素宽。
对于遥感图像5、4、3波段构成的真彩色图像，条带像素有三种。第一种是红、绿、蓝灰度都为0的黑色像素。第二种是纯红像素，这些像素的绿色和蓝色值为0，第三种是红色灰度值为0，蓝色和绿色不为0的蓝绿色像素，这些非地面真实反射率的条带像素需要利用特定算法进行填充处理，满足以上三种像素要求的像素被初步定位为条带区域内像素。
我们从待恢复图像中提取条带图，将条带区域像素灰度置0，非条带区域置255。从后面的实验可以看到，一些非条带区域被误当做条带区域提取到。这是因为某些非条带区域的灰度值恰好符合我们设定的条带灰度特征。观察发现，这些误提取的条带区域比较零散，明显与周期性出现的条带区域不同，需要采取某种策略去掉。由于图像数学形态学方法中的腐蚀和膨胀运算具有去掉小点和小桥等小区域噪声的特性，我们使用如图2所示的半径为2的圆盘形结构元素对二值条带图中条带宽度大于5像素的区域进行先膨胀后腐蚀的闭运算。若A为有条带遥感图像，B为图2中的结构元素，则A被B进行先膨胀后腐蚀的形态学闭运算可表示为：
Figure BDA0002414878950000101
这种闭运算的特性是，既能去除灰度值为0的微小区域的小桥和小点，又能保证大于结构元素的0值区域在闭运算前后不被破坏。实验表明，使用圆盘形结构元素进行闭运算，有效地去除了噪声，与原始图像相比，闭运算去噪处理前后条带的面积和形状没有发生变化，为准确进行条带像素的填充奠定基础。
利用形态学闭运算方法获得去噪处理后的条带区域位置，根据获得的条带区域位置，对原始图像中条带区域内的像素进行灰度填充。
1.2.2局部平均灰度
对所有的条带像素点计算它的局部平均灰度。由于条带宽度逐渐增加，这种计算局部平均灰度的窗口尺寸应随之逐渐增大。我们以17×17的窗口为例计算平均灰度值，与非监督分类一起用于之后突出亮暗纹理的填充方向的判断。
1.2.3非监督分类
利用ENVI软件4.5以上版本，选取非监督分类的参数，对有条带原始遥感图像进行分类初步分为5类。在这个部分的分类结果与回归分析计算得到的各方向的像素灰度变化斜率相配合，根据三个方向上条带两侧相邻像素是否属于同类，来判断确认填充方向。
如图3所示，中间浅灰色区域是条带区域，包括条带区域所有像素被分成5类。条带垂直方向上下外侧的邻近像素分别是2类和1类，不属于同一类型，这种情况意味着垂直方向不是连续的区域，不建议作为填充方向。45°方向上条带两侧右上和左下相邻像素都是1类，则此方向可以作为备选的填充方向。类似地，135°方向上，条带两侧左上和右下方向相邻像素都是2类，也可以作为备选的填充方向。
该非监督分类准则优先考虑条带两侧相邻像素类型，保证优先使用同类型区域方向为填充方向，在后面的插值计算中获得同类型区域的像素灰度作为填充灰度。这是因为，同一类型的区域往往是连续的，被条带遮挡部分是其他类的可能性要小，只有当发生突变阶跃的边界落在条带区时，两侧的类型才不同。
1.3回归分析
我们从齐齐哈尔扎龙湿地的有条带遥感图像中截取跨越条带两侧的一个10×10的区域，如图4所示。对条带上当前像素点在某个方向上进行线性回归分析时，从当前条带像素点的位置沿着选定方向进行搜索。如图4所示，以45°方向为例，当前像素点坐标为(i,j)，向右上方和左下方搜索，分别找到条带两侧外侧与条带相邻的非条带像素，两侧各取3点共6点，其灰度值分别为85、97、109、84、109、121。以像素位置为自变量，6点灰度为函数值进行回归分析。垂直方向和135°方向做法类似。
这里回归分析的置信度默认为小于0.05，若满足，则模型可用。如图5所示，在残差图上可以看到灰度为84的点其置信区间不包含零点，去掉该点重新进行回归分析得到斜率值为2.7429。去掉所有残差置信区间不包括0点的噪声点，重新计算得到斜率值。若置信度不小于0.05，则直接对其拟合求斜率。由于斜率代表灰度值变化趋势，其值也可以为负值，对于负的斜率值仅使用其绝对值参与斜率大小比较。
填充方向选择的主要策略就是将灰度变化最慢的方向，即斜率绝对值最小的方向作为插值填充的方向。而选择斜率最小方向为填充方向之前优先考虑条带两侧区域为同类型的各个方向。当条带两侧为不同类型时，极端情况下，三个方向条带两侧都不属于同一类型，此时三个方向的拟合斜率都比较大，比如在遥感图像中出现半岛状的区域。此时，我们优先选择斜率最小方向为填充方向来进行插值计算。
1.4确定填充方向
如图1的流程所示，在线性回归分析得到各个方向的斜率之后，使用非监督分类准则判断填充方向。
(1)若只有一个方向的条带两侧相邻像素为同类，则选择该方向为填充方向。
(2)若有两个及以上的方向上条带两侧相邻像素为同类，则将这些方向的斜率进行从小到大排序。
(3)若三个方向中没有一个方向上条带两侧相邻像素为同类，则同(2)的处理办法。
(4)在(2)和(3)的结果中选择最小和次最小斜率，比较这两个方向条带两侧相邻6个像素(两侧各取3个像素)的平均灰度与局部平均灰度之差的绝对值，如下面公式，
|G1-G0|≥|G2-G0|                           (2)
式中G1和G2分别为斜率最小和次最小方向上6像素平均灰度，G0为局部17×17窗口的平均灰度。若上式成立，则选取第一方向，即斜率最小方向为填充方向，反之选取第二方向，即斜率次最小方向为填充方向。
值得注意的是，对于(2)中条带两侧相邻像素类型相同的情况下，比较各个方向6像素平均灰度与局部平均灰度偏离程度，选取灰度偏离局部平均灰度更大的方向为填充方向的做法，可以突出亮纹理或者暗纹理。
1.5条带像素填充
1.5.1边界像素填充
对于条带像素，我们将其分为边界像素和非边界像素两类。对于边界像素分四种情况分别处理。对于只有1个像素宽的条带，该条带像素就是边界像素。在其填充方向上方或者下方随机取一点代替该点灰度。对于条带宽度只有2个像素的，对于上方的条带像素都认为是边界像素。上边界像素在当前像素填充方向上方的2点中随机取1点代替当前条带边界点，类似地，下边界点在填充方向上条带下方2点中随机取1点代替当前边界点。若条带宽度为3个像素，则上下两边界点填充方式参考宽度为2的条带边界点，而中间一个像素则随机取其填充方向上条带外侧上方或者下方一点来填充。对于4个及4个以上宽度的条带，分别定义上边缘的2个点及下边缘的2个点为边界点，在其填充方向上条带上侧的最近的2点中随机选取1点来代替该点灰度。而下侧的2个边界点则在其填充方向上条带下侧的最近的2点中随机选取1点来代替该点灰度。
根据Lark的理论，相邻像素空间上具有相关性，一般两个距离近的像素比距离远的像素灰度值更相似。据统计在卫星地图中平均至少有90％以上的像素与其邻近像素具有相同或相近灰度。对于条带边界的两像素，为了不破坏图像的纹理特性，我们随机选取条带外侧填充方向上相邻最近两像素中的一个来代替待填充边界像素灰度。
1.5.2非边界像素填充
对于条带的非边界点像素，是填充方向上条带两侧像素灰度的过渡。使用插值的方法恢复条带像素灰度值是很自然的选择。常用的插值算法有拉格朗日插值，牛顿插值以及三次样条插值等。由于三次样条插值可以避免高次插值可能出现的大幅度波动现象，同时能够克服拉格朗日插值的不收敛性，并且可以提高分段线性插值函数在节点处的光滑性。三次样条插值以其优异的特性成为插值领域的里程碑。
样条曲线S(x)是一个分段定义的公式。给定n+1个数据节点，共有n个区间，三次样条方程满足：
(a).在每个分段区间[xi,xi+1](i＝0,1,…,n-1)，S(x)＝Si(x)都是三次多项式。
(b).满足S(xi)＝yi(i＝0,1,…,n)。
(c).S(x)，导数S’(x)，二阶导数S”(x)在[a,b]区间都连续的，即曲线S(x)光滑。则n个三次多项式可分段写作：
Si(x)＝ai+bi(x-xi)+ci(x-xi)2+di(x-xi)3i＝0,…,n-1     (3)
其中ai，bi，ci，di代表4n个未知系数。对其进行求解。对于n+1个数据节点(x0,y0)，(x1,y1)，(x2,y2)，……，(xn,yn)，分别进行如下计算：
(a)计算步长hi＝xi+1-xi(i＝0,1,…,n-1)。
(b)将数据节点和首位的端点代入矩阵方程。
(c)解矩阵方程，求得二次微分值mi。我们给出样条曲线系数的表达式，
Figure BDA0002414878950000131
其中(i＝0,1,…,n-1)。
(d).在每个子区间xi≤x≤xi+1中创建方程Si(x)。对于其边界端点的数据节点不受任何约束为自由边界。
接下来我们使用三次样条插值算法计算条带非边界区域当前像素点的插值灰度。参与插值计算的节点为填充方向上条带两侧的6个点，这与图4所示的回归分析选点类似。如图6所示，将当前像素坐标作为插值点代入曲线方程，求出对应的插值灰度。其他非当前像素需要在自己的填充方向上单独计算各自的插值灰度。
计算出当前点的灰度值之后，不能直接使用，因为计算出的灰度有可能是遥感图像中没有的灰度，因此，在6个点中寻找与该计算灰度最接近的灰度，用来代替该点灰度。如下面公式所示，
N(i,j,m0)＝min{|P(i,j)-N(i,j,m)||m＝1…6}     (5)
其中，P(i,j)为当前条带像素点的插值计算灰度值，N(i,j,m)为参与插值计算的6个节点的灰度值，N(i,j,m0)为与当前像素点的插值计算结果最接近的节点的灰度值。用N(i,j,m0)代替当前像素点P(i,j)的值。这种从N1,N2,……,N6中寻找最接近灰度来代替待当前条带像素填充灰度值的方法，将最可能的灰度***插值点，同时不会引入附近区域没有出现过的新灰度。若引入新的灰度，在对遥感图像进行分类时会引入新的类型，在条带区域出现附近并不存在的新类型显然是不正确的。使用此种策略进行填充，提高了填充的准确性和真实性，同时保证了遥感图像的纹理特性。
1.6条带区域滤波
使用3×3的窗口进行滤波，其窗口中心经过条带的所有像素，且仅对条带像素进行滤波处理，并不破坏条带以外像素灰度。滤波过程中，首先判断中心像素是否为噪声点。我们借鉴双阈值的思想，若中心像素灰度高于3×3局部平均灰度Pl(1+K)(其中K为经验值)，则该点在3×3窗口中为亮点，需要去除。将周围8邻域像素灰度从小到大排序，取第5个值来代替当前亮点像素灰度。
若中心像素灰度低于3×3局部平均灰度Pl(1-K)，则该点在3×3窗口中为暗点，也需要去除。将周围8邻域像素灰度从小到大排序，取第4个值来代替当前亮点像素灰度。
在亮区域我们尽量选择稍亮的点，而在暗区域则相反选择稍暗的点来代替噪声点。若窗口中心像素灰度不是噪声点则不处理。其中的K为经验值。此方法选择性地对噪声点进行处理，既能去除亮点和暗点，又避免了滤波去除椒盐噪声时平滑图像的缺点，保留了遥感图像的大部分纹理特征。
1.7多波段遥感图像处理
对于彩色或者多波段的遥感图像，同时读取多幅多波段遥感图像，将其转化为灰度图像，根据灰度图像来确定条带像素的填充方向。在确定填充灰度之后，多个波段的多幅遥感图像选取同一相对位置的像素灰度来填充各自的待填充位置。就真彩色遥感图像而言，这样不会导致产生新的颜色，从而填充恢复更加客观真实。其他多波段遥感图像处理方法类似。
在实验部分，首先，对一幅无条带原始卫星图像添加条带进行仿真实验，使用几种方法及本发明所提出的算法对仿真图像进行处理，并比较各种方法的填充效果。之后，对齐齐哈尔扎龙国家级自然保护区(46°52′～47°32′N，123°47′～124°37′E)附近的多波段遥感图像中的5、4、3波段，使用所提出的算法进行条带像素恢复处理，并给出处理结果的真彩色图像。
2.1仿真实验
我们使用本发明所提出的算法对一幅无条带噪声的卫星遥感图像进行仿真实验。对图7(a)的原始图像加入有条带图像中可能出现的最大宽度14个像素宽的条带，得到图7(b)的仿真图像。在该仿真图像中条带像素占大约42.05％，若算法有好的效果，则算法对具有更小宽度条带图像填充效果会更好。
图7(c)给出了使用形态学闭运算对条带去噪的结果。除了下方少部分灰度为0的非条带像素没有被去掉外，绝大多数非条带像素被腐蚀掉了。而图7(d)是对有条带图像进行非监督分类的结果，包括条带在内所有像素被分成5类。图7(e)，7(f)，7(g)分别给出使用LHM，NSPI及本发明所提出的方法填充条带的结果。图7(h)，7(i)，7(j)分别是几种方法得到图像与原始图像条带区域像素灰度的散点图,可以看到使用所提出的策略恢复的条带像素散点图与其他方法得到的散点图差别不大，图中的大部分散点坐标都落在以45°直线为对称轴的椭圆形区域内，说明像素灰度恢复前后呈线性关系，证明了所提出方法恢复遥感卫星图像条带区域像素的有效性。
图7(k)，7(l)，7(m)是几种方法得到图像的分类结果。可以看到使用所提出方法对有条带图像的填充结果目视效果上不低于其他方法，大部分像素可以正确分类。比较黑色矩形框区域内像素的分类结果，使用所提出方法像素分类错误更少，而用LLHM方法获得的填充像素中有更多的数量被误分成其他类，效果最差。
为进一步定量评价几种方法对条带像素填充的效果，我们通过比较统计学方法中的均方根误差(RMSE)来定量比较几种算法对条带像素的恢复精度。
Figure BDA0002414878950000141
公式(6)中，Oi是原始图像中当前待恢复像素位置的灰度，Ri是对当前像素位置恢复之后的灰度。对图7(a)的原始卫星图像加入14个像素宽的条带得到图7(b)的仿真图像，对其使用几种方法填充。表1对条带区域像素使用几种方法进行恢复填充的结果进行了比较。
表1
Figure BDA0002414878950000142
2.2对实际卫星图像的处理
为方便展示实验结果，如图8(a)所示，截取400×400的部分真彩色遥感图像。接下来，按照流程对其进行处理。
为了避免误填充，首先定位条带像素。如图8(b)所示，可以看到很多非条带像素被误提取到，为了能够去除这些非像素，我们使用形态学图像处理的闭运算，先膨胀后腐蚀，既能够将黑色非条带像素去除，又保证了条带区域不受影响。选取尺寸为5×5的模板进行闭运算。从图8(c)中可以看出，大部分非条带像素都被去除了。
对图8(a)的有条带遥感图像进行非监督分类，得到图8(d)，这里我们将其自然分成5类。可以看到条带像素与一部分水域被自动分成一类。为了避免这种错误分类，同时也有助于目视解译遥感图像，我们考虑在非监督分类准则下来判断在回归分析计算得到的填充方向上的条带两侧区域是否为同类，我们优先以此原则来确定填充方向。
按照所提出的方法，在非监督分类准则下，对回归分析计算得到的填充方向进行排序，同时采用了突出明亮的和偏暗的纹理的策略，我们找到了每个条带像素的填充方向。在该方向上计算当前像素的三次样条插值，将计算得到的插值与参与插值计算的像素比较灰度。为了不引入新的颜色和灰度，用最接近的灰度代替当前像素灰度值，我们得到了图8(e)的条带像素填充处理结果。
图8(f)是对图8(e)使用提出的对亮区域的暗像素，或者暗区域的亮像素进行去除滤波的结果。对当前像素周围3×3窗口内的8像素进行排序，根据3×3区域平均灰度与局部小范围17×17窗口范围内平均灰度比较来确定使用排序后的哪个灰度来代替当前像素的滤波策略，可以看到滤波去除了很多噪声点。为了清晰起见，我们将图8(e)的矩形窗口和图8(f)的矩形窗口区域放大，得到图9(a)至图9(d)。可以看到，和中值滤波的所有像素灰度都被平滑相比，我们的滤波策略中，大部分噪声都被去掉，同时图像纹理并没有被破坏，非噪声像素没有被影响。
图10(a)至图10(f)给出了一幅完整有条带Landsat-7SLC-off卫星图像5，4，3波段的去条带噪声的处理结果。其中按照图10(d)，图10(a)，图10(b)，图10(c)顺序给出了有条带图由整体到局部的排列图。从图中可以看到条带的细节。而图10(e)，图10(h)，图10(g)，图10(f)的顺序给出了对应的条带图的像素恢复结果。目视非常理想，完全可以满足研究者的要求。
本发明针对因Landsat-7卫星ETM+传感器的扫描线纠正器(SLC)永久失效而产生的遥感卫星图像有条带问题，提出了一种无参考图条带丢失像素灰度恢复算法。该算法首先利用回归分析，有效计算并比较三个方向的像素灰度变化斜率，进而在非监督分类准则下，选择插值填充方向。随后，在填充方向上进行插值计算，并从原图中选择最接近插值的灰度值，保证了填充灰度的真实性与可靠性。此外，算法在选择填充方向时，为了最大限度保留并突出图像纹理，赋予图像亮暗纹理方向更高的优先权。与几种经典算法相比较，使用所提出的算法填充恢复条带区域，其灰度散点图与其他方法差别不大，而且，在目视效果和均方根(RMSE)统计学评价方面要优于其他两种算法。在非监督分类方面，目视某些区域所提出算法分类错误最少。实验结果证明，使用所提出的算法可有效恢复Landsat-7SLC-off多波段有条带图像。
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (5)

1.一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，其特征在于，该方法包括以下步骤：
步骤一、从SLC-off遥感图像中初步提取出条带区域，再采用形态学闭运算对初步提取出的条带区域进行去噪处理，获得条带区域；
步骤二、对于条带区域内的任一像素(i,j)，分别从像素的位置(i,j)出发沿着各个选定方向进行搜索；
对于任意的一个选定方向，在选定方向的上方选取m个与条带上边界最相邻的非条带像素，在选定方向的下方选取m个与条带下边界最相邻的非条带像素；
将选取的与条带上边界最相邻的m个非条带像素记为Ak，k＝1,2,……,m，Ak的灰度值为ak；将选取的与条带下边界最相邻的m个非条带像素记为Bk，Bk的灰度值为bk；
以选定方向上选取的2m个像素的位置作为自变量，以选取的2m个非条带像素的灰度值为函数值进行回归分析，在该选定方向上进行回归分析计算出灰度变化的斜率值；
同理，在其他各个选定方向上进行回归分析计算出条带区域内各点的斜率值；
步骤三、在各个选定方向上，分别采用非监督分类准则，判断在选定方向上条带两边界外部的非条带像素是否为同类像素；
步骤四、计算条带区域各像素的局部平均灰度，根据步骤二计算出的在各个选定方向上的斜率值以及步骤三的判断结果确定条带区域内像素(i,j)的填充方向；
步骤五、根据条带区域内像素(i,j)是否为边界像素，采用不同的方式对像素(i,j)进行填充，获得像素填充后的条带区域，具体过程为：
步骤五一、若像素(i,j)所在位置是仅有1个像素宽的条带，则像素(i,j)为边界像素，在填充方向的上方或下方，随机选取出一个与像素(i,j)最相邻的非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五二、若像素(i,j)所在位置是有2个像素宽的条带，则像素(i,j)为边界像素；若像素(i,j)为上边界像素，则在填充方向的上方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；若像素(i,j)为下边界像素，则在填充方向的下方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五三、若像素(i,j)所在位置是有3个像素宽的条带，且像素(i,j)为上边界像素或下边界像素，则像素(i,j)的填充方法与步骤五二相同；否则，在像素(i,j)的填充方向的上方或者下方，随机选取出一个与像素(i,j)最相邻的非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五四、若像素(i,j)所在位置是有4个像素宽的条带，且像素(i,j)为填充方向的上边缘的2个条带像素中的一个，则在像素(i,j)的填充方向的上方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
若像素(i,j)所在位置是有4个像素宽的条带，且像素(i,j)为填充方向的下边缘的2个条带像素中的一个，则在像素(i,j)的填充方向的下方，与像素(i,j)最相邻的2个非条带像素中随机选取出1个非条带像素，将选取出的非条带像素的灰度作为像素(i,j)的灰度进行填充；
步骤五五、若像素(i,j)所在位置是有4个以上像素宽的条带，同时，像素(i,j)为填充方向的上边缘的2个条带像素中的一个或像素(i,j)为填充方向的下边缘的2个条带像素中的一个，则对像素(i,j)的填充方式同步骤五四；
若像素(i,j)不是填充方向的上边缘的2个条带像素中的一个且像素(i,j)也不是填充方向的下边缘的2个条带像素中的一个，则像素(i,j)为非边界像素；
根据步骤二选取出的在像素(i,j)填充方向上的2m个非条带像素进行三次样条插值，求出像素(i,j)的插值灰度P(i,j)；
Figure FDA0004087745910000021
其中：k0代表步骤二选取出的在像素(i,j)填充方向上的2m个非条带像素中的第k0个，k0＝1,2,…,2m，N(i,j,k0)为第k0个非条带像素的灰度值，N(i,j,m0)为步骤二选取出的在像素(i,j)填充方向上的2m个非条带像素的灰度值中，与P(i,j)最接近的灰度值；
将N(i,j,m0)作为像素(i,j)的灰度进行填充；
步骤六、像素填充后，再分别对条带区域内包含的各像素进行滤波处理，滤波处理后，完成对Landsat-7扫描线纠正器失效的遥感图像丢失像素的恢复。
2.根据权利要求1所述的一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，其特征在于，所述步骤二中，分别从像素的位置(i,j)出发沿着各个选定方向进行搜索，各个选定方向分别为：垂直方向、与水平方向呈45°角方向以及与水平方向呈135°角方向。
3.根据权利要求2所述的一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，其特征在于，所述步骤四中，计算条带区域的局部平均灰度，根据步骤二计算出的在各个选定方向上的斜率值以及步骤三的判断结果确定条带内像素(i,j)的填充方向，其具体过程为：
步骤四一、利用尺寸为17×17的窗口计算像素(i,j)的局部平均灰度G0；
步骤四二、若仅存在着一个选定方向满足条件：在选定方向上，条带两边界外部的非条带像素为同类像素，则将该选定方向作为填充方向；
步骤四三、若存在两个或两个以上的选定方向满足条件：在选定方向上，条带两边界外部的非条带像素为同类像素，则将满足条件的选定方向的斜率绝对值按照从小到大排序，并执行步骤四五；
步骤四四、若不存在选定方向满足：在选定方向上，条带两边界外部的非条带像素为同类像素，则将全部选定方向的斜率绝对值按照从小到大排序，并执行步骤四五；
步骤四五、选取出斜率绝对值最小和次最小的选定方向，根据步骤二找出的在每个选定方向上的2m个非条带像素的灰度值，分别将每个选定方向上找出的2m个非条带像素的灰度值取平均，获得斜率绝对值最小的选定方向对应的均值G1以及斜率次最小的选定方向对应的均值G2；
若满足：
|G1-G0|≥|G2-G0|
则将斜率绝对值最小的选定方向作为像素(i,j)的填充方向，否则将斜率绝对值次最小的选定方向作为像素(i,j)的填充方向。
4.根据权利要求3所述的一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，其特征在于，所述步骤二中，m的取值为3。
5.根据权利要求4所述的一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法，其特征在于，所述步骤六中，分别对条带区域内包含的各像素进行滤波处理，其具体过程为：
使用3×3的窗口对条带区域内包含的每个像素分别进行滤波处理，使窗口中心经过条带区域内的每个像素；
当窗口中心经过像素(i,j)时，则称像素(i,j)为中心像素，在3×3的窗口内，像素(i,j)的局部平均灰度为Pl；
若中心像素(i,j)的灰度大于等于Pl(1+K)，其中K为经验值，则判断像素(i,j)为亮点，将在3×3窗口内的中心像素周围的8个像素的灰度按照从小到大排序，利用排序中的第5个值代替像素(i,j)的灰度；
若中心像素(i,j)的灰度小于等于Pl(1-K)，则判断像素(i,j)为暗点，将在3×3窗口内的中心像素周围的8个像素的灰度按照从小到大排序，利用排序中的第4个值代替像素(i,j)的灰度；
否则，不对像素(i,j)的灰度进行处理；
同理，分别对条带区域内包含的每个像素进行上述处理。
CN202010188037.7A 2020-03-17 2020-03-17 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法 Active CN111476723B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010188037.7A CN111476723B (zh) 2020-03-17 2020-03-17 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010188037.7A CN111476723B (zh) 2020-03-17 2020-03-17 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法

Publications (2)

Publication Number Publication Date
CN111476723A CN111476723A (zh) 2020-07-31
CN111476723B true CN111476723B (zh) 2023-04-18

Family

ID=71747498

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010188037.7A Active CN111476723B (zh) 2020-03-17 2020-03-17 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法

Country Status (1)

Country Link
CN (1) CN111476723B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112150385B (zh) * 2020-09-29 2023-05-12 四川虹美智能科技有限公司 红外图像滤波方法及装置
CN112669237B (zh) * 2020-12-31 2022-09-20 同济大学 一种Landsat 7 SLC-off条带修复方法
CN113034390B (zh) * 2021-03-17 2022-10-18 复旦大学 一种基于小波先验注意力的图像修复方法及***

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101694718A (zh) * 2009-10-13 2010-04-14 西安电子科技大学 基于感兴趣区域的遥感图像变化检测方法
CN104065946A (zh) * 2014-06-17 2014-09-24 四川虹微技术有限公司 基于图像序列的空洞填充方法
CN108109126A (zh) * 2018-01-12 2018-06-01 适普远景遥感信息技术（北京）有限公司 一种基于卫星遥感影像的目标区域填充与融合处理方法

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5592571A (en) * 1994-03-08 1997-01-07 The University Of Connecticut Digital pixel-accurate intensity processing method for image information enhancement
JPH10222688A (ja) * 1997-01-21 1998-08-21 Xerox Corp 画像処理方法
CN103034981B (zh) * 2012-12-18 2015-06-24 武汉大学 基于多时相数据的遥感影像加权回归恢复方法
US20170178297A1 (en) * 2014-02-19 2017-06-22 Yissum Research Development Company Of The Hebrew University Of Jerusalem Ltd. Method and system for dehazing natural images using color-lines
CN104200426A (zh) * 2014-08-25 2014-12-10 北京京东方视讯科技有限公司 图像插值方法和装置
WO2018022648A1 (en) * 2016-07-25 2018-02-01 Iteris, Inc. Image-based field boundary detection and identification
CN106780485B (zh) * 2017-01-12 2020-01-03 西安电子科技大学 基于超像素分割和特征学习的sar图像变化检测方法
JP6953183B2 (ja) * 2017-05-24 2021-10-27 キヤノン株式会社 画像処理装置、画像処理方法およびプログラム
CN107274361B (zh) * 2017-05-26 2019-11-19 深圳先进技术研究院 Landsat TM遥感影像数据除云方法及***
US10796145B2 (en) * 2017-06-29 2020-10-06 Samsung Electronics Co., Ltd. Method and apparatus for separating text and figures in document images

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101694718A (zh) * 2009-10-13 2010-04-14 西安电子科技大学 基于感兴趣区域的遥感图像变化检测方法
CN104065946A (zh) * 2014-06-17 2014-09-24 四川虹微技术有限公司 基于图像序列的空洞填充方法
CN108109126A (zh) * 2018-01-12 2018-06-01 适普远景遥感信息技术（北京）有限公司 一种基于卫星遥感影像的目标区域填充与融合处理方法

Also Published As

Publication number Publication date
CN111476723A (zh) 2020-07-31

Similar Documents

Publication Publication Date Title
CN111476723B (zh) 一种Landsat-7扫描线纠正器失效的遥感图像丢失像素恢复方法
Zhou et al. Object-based land cover classification of shaded areas in high spatial resolution imagery of urban areas: A comparison study
Zhang et al. Object-oriented shadow detection and removal from urban high-resolution remote sensing images
CN102132323B (zh) 用于自动图像矫直的***和方法
US8135174B2 (en) Automatic image object identification using threshold gradient magnitude based on terrain type
CN104361582B (zh) 一种对象级高分辨率sar影像洪水灾害变化检测方法
CN111696123A (zh) 超像素分类识别的遥感图像水域分割提取方法
Yang et al. River delineation from remotely sensed imagery using a multi-scale classification approach
CN106295562A (zh) 一种高分辨率遥感影像道路信息提取方法
Karami et al. Gully erosion mapping using object-based and pixel-based image classification methods
Ranjbarzadeh et al. LNPSS: SAR image despeckling based on local and non-local features using patch shape selection and edges linking
US20140081605A1 (en) Dtm estimation method, dtm estimation program, dtm estimation device, and method for creating 3-dimensional building model, and region extraction method, region extraction program, and region extraction device
Yang et al. Fully constrained linear spectral unmixing based global shadow compensation for high resolution satellite imagery of urban areas
CN111353371A (zh) 基于星载sar影像的海岸线提取方法
CN110889840A (zh) 面向地物目标的高分6号遥感卫星数据的有效性检测方法
CN111310771A (zh) 遥感影像的道路图像提取方法、装置、设备及存储介质
de Vries et al. Unmixing water and mud: Characterizing diffuse boundaries of subtidal mud banks from individual satellite observations
Basso et al. Combining morphological filtering, anisotropic diffusion and block-based data replication for automatically detecting and recovering unscanned gaps in remote sensing images
Huang et al. Multi-feature combined for building shadow detection in GF-2 Images
CN113408547A (zh) 一种多时相多极化sar滑坡提取方法
Aytekin et al. Automatic and unsupervised building extraction in complex urban environments from multi spectral satellite imagery
He et al. Narrow-linear and small-area forest disturbance detection and mapping from high spatial resolution imagery
CN113361532A (zh) 一种图像识别方法、***、存储介质、设备、终端及应用
CN107153823B (zh) 一种基于视觉关联双空间的车道线特征提取方法
Sghaier et al. A novel approach toward rapid road mapping based on beamlet transform

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant