一种基于空谱径向基函数插值的条带修复方法与流程

本发明涉及遥感图像处理技术领域，尤其是涉及一种基于空谱径向基函数插值的条带修复方法。

背景技术：

landsat系列卫星包括landsat7为地表监测提供了十分有效的数据。2003年5月31日，landsat7卫星的etm+线性扫描校正器(thescan-linecorrector，slc)永久失效，使得其获取的影像上会出现未扫描的条形间隙。这些未扫描的部分约占整幅影像的22％，严重影响了landsat7卫星数据在各方面的应用。通常将slc失效前获取的影像称为slc-on影像，slc失效后获取的影像称为slc-off影像。为了解决slc-off问题，使etm+的应用更可靠，发展相应的理论与方法来估计未扫描像素的数据显得极为必要。slc-off数据修复的基本原理是：采用插值方法将数值填充到slc-off数据对应的位置，使整幅图像完整且色调尽可能一致。目前使用的条带修复方法主要包括以下几大类。第一大类是使用slc-off影像自身信息的插值方法，这类方法使用同一景影像中的有效像元来填补未扫描像元，例如envi中的gapfill插件就是基于这种原理开发出来的。该方法因不借助任何辅助数据且简单易实现而获得广泛应用。然而，其修复后的条带与周围地物通常有明显的差异，修复效果通常较差。为了得到更好的修复效果，研究者们开发出了借助其他传感器数据修复条带的方法，即第二类方法。这类方法使用的辅助影像(又称已知影像)尽管获取时间与slc-off影像不一致，但空间覆盖完整，提供了slc-off影像上缺失值对应地理位置在其他时间点的数据，因而此类数据的使用会得到更准确的预测结果。其中，美国地质勘探局(usgs)在2004年编制的一份官方报告中提到的局部线性直方图匹配(localizedlinearhistogrammatch，llhm)方法为使用最多的方法之一。该方法使用一景或多景已知影像，在每个缺失像元的移动窗口中根据有效像元和对应位置的已知影像数据，建立线性变换函数，然后使用该函数将已知影像中缺失像元对应位置的数据转换为slc-off缺失像元的数据。此类方法简单易行，在同质区域中具有较高的精度。此外，邻域相似像元插值方法(neighborhoodsimilarpixelinterpolater，nspi)也是一种较为经典和常用的插值方法。nspi方法用于异质性较强的区域也可以得到较好的效果，但对于获取时间距slc-off影像较远的已知影像，其预测结果精度会明显下降。

近年来，各类插值方法在修复slc-off影像条带上的应用，尤其是在地物覆盖均匀的地区往往都有较好的修复效果，而在异质性较强的区域往往会出现较大的预测偏差。此外，大多数插值方法需要使用与slc-off影像在时间上较为接近的影像作为已知影像才能获得较为理想的预测结果。然而，由于云雾的影响，在很多地区往往无法获得时间足够接近的可用已知影像，这使得前述方法的应用有时达不到较好的效果。为此，有必要开发出更稳定、更高精度的插值方法来预测缺失像元。

技术实现要素：

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于空谱径向基函数插值的条带修复方法，本发明方法能够有效综合已知图像与待修复图像的有用信息，可充分挖掘邻近像元光谱信息，具有更佳的修复效果。

本发明的目的可以通过以下技术方案来实现：

一种基于空谱径向基函数插值的条带修复方法，用于slc-off条带影像的修复包括如下步骤：

s1：根据待修复slc-off条带影像，对已知影像进行全局直方图匹配(glhm，globallinearhistogrammatch)，获取匹配后的已知影像；

s2：利用步骤s1得到的实施glhm后的已知影像，采用空谱径向基函数方法对待修复slc-off条带影像的缺失像元进行预测。

进一步地，步骤s1的具体内容为：

将已知影像和待修复slc-off条带影像的获取时刻分别记为tk和tp，对已知影像和待修复slc-off条带影像的获取时刻所对应的有效像元的集合进行线性回归建模，对每个波段的数据构造线性模型，利用构造的线性模型所得到的系数将已知影像整体变换为与待修复slc-off条带影像在特征空间内更接近的新已知影像。具体地：

11)根据tk时刻已知影像的有效像元和tp时刻待修复slc-off条带影像的有效像元的集合进行线性回归建模，构建tp时刻待修复slc-off条带影像与tk时刻已知影像之间的线性关系模型，获取线性系数；

对每个波段的数据构造的线性模型的表达式为：

l(x,y,b,tp)＝ab×l(x,y,b,tk)+bb

式中，l(x,y,b,tp)和l(x,y,b,tk)分别为待修复slc-off影像和已知影像在位置(x,y)处波段b的反射率，ab和bb分别为线性系数。线性系数ab和bb根据待修复slc-off影像和已知影像对应的有效数据集，采用最小二乘法拟合求解。

12)对tk时刻已知影像通过线性系数进行整体线性变换，获取tk时刻的新已知影像。

进一步地，步骤s2具体包括下列步骤：

21)选择已知影像和待修复slc-off条带影像的相似像元，对于待修复slc-off条带影像的任一缺失像元，基于已知影像中同一地理位置的像元，找出与其相似的邻域像元的位置，并对于每一个缺失像元选择n个相似像元参与后续计算；

22)获取每一个缺失像元选择的相似邻域像元之间的空间径向基函数，并定义已知影像窗口内中心像元与相似像元之间的空间径向基函数；中心像元与有效像元之间的空间径向基函数的表达式为：

式中，(x0,y0)为中心像元的坐标，δ1为空间径向基参数，(xi,yi)为窗口内第i个和第j个有效像元的坐标。

23)构建已知影像窗口中相似像元之间光谱相似度的光谱径向基函数以及中心像元与其邻近相似像元光谱相似度的光谱径向基函数；已知影像窗口中相似像元之间光谱相似度的光谱径向基函数及中心像元与其邻近相似像元光谱相似度的光谱径向基函数的表达式为：

式中，δ2为光谱径向基参数，rmsdi为中心像元与窗口中第i个相似像元之间的rmsd值，rmsdij为窗口中第i个像元与第j个相似像元之间的rmsd值。

24)将根据步骤22)获取的空间径向基函数与步骤23)获取的光谱径向基函数相乘，获取空谱径向基函数；空谱径向基函数的表达式为：

式中，为空间径向基函数，为中心像元与有效像元之间的空间径向基函数，为待修复slc-off条带影像窗口中相似像元之间光谱相似度的光谱径向基函数，为已知影像窗口的中心像元与其邻近相似像元光谱相似度的光谱径向基函数。

25)采用选定的n个相似像元，基于已知影像中的像元和待修复slc-off条带影像的像元之间的反射率变化，计算径向基权重；

具体地，使用选定的n个相似像元，将其关于已知影像与待修复slc-off条带影像之间的反射率变化记作δl(xi,yi,b)＝l(xi,yi,b,tp)-l′(xi,yi,b,tk)，其中l(x,y,b,tp)为待修复slc-off条带影像和已知影像在位置(x,y)处波段b的反射率，l′(x,y,b,tk)为步骤1)中将已知影像整体变换至与待修复slc-off条带影像在特征空间内更接近的新已知影像，tk和tp分别为已知影像和待修复slc-off条带影像的获取时刻，则将径向基权重记作ωbi(i＝1,2,...,n)，则径向基权重的计算式如下：

26)基于步骤25)获取每个波段的径向基权重，计算对应波段已知影像的像元与待修复slc-off条带影像的像元之间的反射率变化；

27)将步骤1)中实施glhm后的已知影像与步骤26)中的反射率变化叠加，获取最终预测结果。

与现有技术相比，本发明扩展了经典的空间径向基插值法，创新性地加入了光谱径向基并构建空谱径向基插值模型，充分挖掘了邻近像元光谱信息；同时，glhm预处理将已知影像整体变换到与待修复影像在特征空间内更接近的新已知影像，可在一定程度上减少因时间跨度大，地物发生较显著变化而导致的已知影像和待修复影像之间差异大的影响，使得此方法在已知影像时间间隔较大的情况下也能得到较好的修复效果。

附图说明

图1为实施例中基于空谱径向基函数插值的条带修复方法的流程示意图；

图2为实施例仿真实验中区域1的结果图，其中(a)为模拟的slc-off影像，(b)～(e)分别为gapfill插件、llhm、nspi和本发明ssrbf方法的slc-off影像修复结果，(f)为参考影像；

图3为实施例仿真实验中区域2的结果图，其中(a)为模拟的slc-off影像，(b)～(e)分别为gapfill插件、llhm、nspi和本发明ssrbf方法的slc-off影像修复结果，(f)为参考影像。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

实施例

本发明涉及一种基于空谱径向基函数插值的条带修复方法，该方法用于slc-off条带影像的修复，通过对已知辅助影像进行glhm(globallinearhistogrammatch)预处理，得到在特征空间内更接近待修复影像的已知影像；然后使用空谱径向基函数对待修复像元在时间上的反射率变化进行插值计算，得到各波段已知影像像元与待修复像元之间的反射率变化；最后将预测的变化值与已知影像上对应的缺失像元位置的数据相加，得到该缺失像元的最终预测结果。如图1所示，本发明修复方法具体包括如下步骤：

步骤1、对已知影像实施glhm。将已知影像和待修复的slc-off条带影像的获取时刻分别记为tk和tp。根据tk和tp时刻影像有效像元的集合进行线性回归建模，对每个波段的数据构造如式(1)所示的线性模型。

l(x,y,b,tp)＝ab×l(x,y,b,tk)+bb(1)

其中，l(x,y,b,tp)和l(x,y,b,tk)分别是条带中待修复影像和已知影像在位置(x,y)处波段b的反射率。ab和bb是两个线性系数，可根据待修复影像和已知影像对应的有效数据集，使用最小二乘法拟合求解。利用得到的系数可将已知影像整体变换到与待修复影像在特征空间内更接近的新已知影像(记为l′(x,y,b,tk))。

具体地，通过数学推导来证明实施glhm后的已知影像与待修复影像在特征空间内更接近。影像之间的差异大小通过同一波段之间差分图像的期望值来量化。为简化表达，以下公式均基于影像的单一波段表示。

根据tk时刻已知影像和tp时刻slc-off影像有效像元的集合进行线性回归建模，两幅影像之间的关系如式(a1)所示：

lp＝alk+b+γ(a1)

其中lp和lk分别是条带中待修复影像和已知影像的反射率矩阵。a和b是两个线性系数，根据待修复影像和已知影像对应的有效数据集，使用最小二乘法拟合求解。γ是回归模型的残差。由此，在未实施glhm时，tk和tp时刻影像有效像元之间的差异δl可由式(a2)表示：

在实施glhm后，新已知影像通过系数a和b可将已知影像进行整体线性变换得到l′k

l′k＝alk+b(a3)

故tp时刻影像和glhm变换后的tk时刻已知影像之间的差异δl′可由式(a4)表示

δl′＝lp-l′k＝γ(a4)

下面分别计算glhm前δl和glhm后δl′的平方的期望e(δl²)和e[(δl′)²]。

1)未实施glhm的e(δl²)：

根据期望和方差之间的关系：

e(δl²)＝var(δl)+e²(δl)(a5)

对于(a5)中的第一项，将(a2)代入，因其中b为常数，由方差的基本性质

公中“·”表示两矩阵之间的内积。对于(a6)中的cov(lk·γ)，根据期望的基本性质将其展开：

cov(lk·γ)＝e(lk·γ)-e(lk)e(γ)(a7)

对于经典的最小二乘线性回归模型，有两个重要的性质：1)残差的期望为0；2)自变量和残差的乘积的期望为0。如式(a8)所示：

根据式(a8)所述性质可知式(a7)等于0，故式(a6)可简化为：

var(δl)＝var[(a-1)lk]+var(γ)(a9)

即式(a5)可表示为：

e(δl²)＝var[(a-1)lk]+var(γ)+e²(δl)(a10)

2)实施glhm后的e[(δl′)²]

由式(a4)知e[(δl′)²]＝e(γ²)。对于e(γ²)，根据期望与方差之间的关系得：

e(γ²)＝var(γ)+e²(γ)(a11)

由式(a8)可知(a11)中第二项为0，故(a11)可简化为：

e(γ²)＝var(γ)(a12)

比较式(a10)和(a12)可知：e[(δl′)²]＜e(δl²)(式(a10)等式右侧第一项和第三项均大于0)。由此，证明了相较于未实施glhm的已知影像，实施glhm后的已知影像与slc-off影像之间在特征空间内更为接近。

步骤2、利用步骤1中得到的glhm后的已知影像，使用空谱径向基函数(ssrbf)方法来对缺失像元进行预测。具体步骤如下：

步骤2.1，选择相似像元。对于任一缺失像元，基于已知影像中同一地理位置的像元，找出相似的邻域像元的位置。邻域像元的相似度通过光谱的均方根偏差(rootmeansquaredeviation，rmsd)值来度量。均方根偏差的定义如式(2)所示：

其中，l(xi,yi,b,t)和l(xj,yj,b,t)为t时刻(slc-off影像或已知影像的获取时刻)影像窗口内第i和第j(i,j＝1,2,…,m，m为窗口内有效像元的个数)个像元，该影像窗口是成对的，同时在slc-off影像上和已知影像上存在。n是波段数。根据式(2)，基于已知影像，在一个窗口内可计算中心像元和其邻域像元之间的rmsd值，选出该窗口中前n个与中心像元最相似的像元。中心像元与窗口中第i个相似像元之间的rmsd值记为rmsdi。

步骤2.2，构建空间径向基函数。使用常用的高斯函数定义空间径向基函数如式(3)所示。

其中，(xi,yi)和(xj,yj)分别是影像窗口内第i个和第j个有效像元的坐标，δ1为空间径向基参数。参数δ1的值根据中心像元与其所在影像窗口内其他像元的最大距离的2倍的值来定义。类似地，可定义中心像元与有效像元之间的空间径向基函数如式(2)所示。

其中，(x0,y0)为中心像元的坐标。

步骤2.3，构建光谱径向基函数。定义基于窗口中相似像元之间光谱相似度的光谱径向基函数和基于中心像元与其邻近相似像元光谱相似度的光谱径向基函数如式(5)所示：

其中，δ2为光谱径向基参数。根据每个窗口中像元之间rmsd值的集合来确定参数δ2的值，即取集合中最大值的2倍。例如，若统计结果中绝大部分的rmsd值不超过0.05，则将δ2的值定为0.1。rmsdij为窗口中第i个像元与第j个相似像元之间的rmsd值。

步骤2.4，构建空谱径向基函数。在考虑像元之间空间距离的基础上，同时考虑邻近像元光谱相似性的影响。将空间径向基函数与光谱径向基函数相乘，则得到空谱径向基和如式(6)所示：

步骤2.5，计算径向基权重。利用ssrbf估计中心缺失像元在时间上的反射率变化。首先，使用选定的n个相似像元，记其关于已知影像与待修复影像之间的反射率变化为δl(xi,yi,b)＝l(xi,yi,b,tp)-l′(xi,yi,b,tk)，径向基权重为ωbi(i＝1,2,...,n)。权重的计算如式(7)所示：

对于一个窗口中的n个相似像元，利用上式可求解出唯一一组权重值ωbi。

步骤2.6，计算ssrbf插值。将每个波段计算得到的权重ωbi代入式(8)中，即得到对应波段已知影像像元(待填补像元位置对应已知影像位置上的像元)与条带中待修复像元(待填补像元)之间的反射率变化

步骤2.7，将已知影像中心像元值与预测的变化值相加，得到最终预测结果。如式(9)所示，

为了验证本发明的有效性，本实施例采用本发明方法修复使用landsatoli影像模拟生成的slc-off影像，并将修复结果与现有经典算法进行比较。对比方法为envigapfill插件中方法、usgs官方发布的局部线性直方图匹配法(localizedlinearhistogrammatch，llhm)和经典的邻域相似像元插值法(neighborhoodsimilarpixelinterpolator，nspi)进行修复。两个测试区域分别位于意大利的维罗纳市(区域1)和中国浙江省(区域2)。两个区域的修复结果分别如图2和图3所示。在两组图中每一个子图中(例如图2(a))，左侧大图为整个实验区域的影像，其对应的右侧小图为放大的两个子区域，以便更直观地对比使用各方法修复slc-off影像条带结果之间的差异。

从图2、图3可以看出，待修复影像所在区域总体上异质性较强。gapfill插件修复结果的色调与原始影像大体一致，但存在明显的条带边缘。llhm和nspi方法都能较好地预测出地物的反射率数据，但修复结果在异质性较强的区域均存在较大的预测误差，例如当条带位于成块完整地物中心时，对该地物的预测存在边缘不完整以及其内部像元光谱预测信息存在明显偏差的问题。由于在本发明方法中使用全局直方图匹配降低了已知数据和待修复数据获取时间差异带来的影响，并基于ssrbf插值模型高精度的优势充分挖掘数据自身更多信息，因此能更准确地预测出地物的完整轮廓和其内部缺失像元的光谱信息。

采用均方根误差(rootmeansquareerror，rmse)和相关系数(correlationcoefficient，cc)评价指标对各方法的修复结果进行精度评价，如表1所示。表中第一行自左向右分别为实验涉及到的6个波段，最后一栏是6个波段实验结果的平均值，6个波段依次为蓝色波段、绿色波段、红色波段、近红外波段、短波红外波段1、短波红外波段2。其中rmse度量修复信息与真实信息的差异性，其值越大表明预测信息越偏离真实信息；cc反映预测信息与真实信息之间的相关度，其值越大表示预测信息与真实信息越接近。

表1修复结果精度评价

从表1的客观评价结果可以看出，本发明方法优于其他3种传统方法，各项指标均表明本发明方法能够得到更接近真实情况的修复图像。

综上所述，本发明的条带修复方法从视觉和精度评价上都具有很明显的优势，得到的修复图像能较好地保持地物的信息，是一种可行的slc-off影像修复方法。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的工作人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。