本发明涉及一种遥感影像处理方法，尤其涉及一种基于NSCT的积分平差遥感影像超分辨率重建方法。
背景技术：
非下采样轮廓波变换(nonsubsampled Contourlet transform,NSCT)具有真正意义上的多尺度、多方向分析，目前已被广泛地应用到遥感影像融合当中。NSCT用于超分辨率重建目前已经有些研究，但目前的方法只是将NSCT分解后的高频信息进行简单的处理，并没用充分利用各影像的高频信息。
技术实现要素：
本发明要解决的技术问题是克服现有技术的不足，提供一种超分辨率重建能力强，并且鲁棒性好的基于NSCT的积分平差遥感影像超分辨率重建方法。
一种基于NSCT的积分平差遥感影像超分辨率重建方法，包括以下步骤1)设置NSCT变换级数N和方向数2N，至少对两幅原始遥感影像a和b进行NSCT分解获得各NSCT系数：
其中与分别为低分辨率影像a、b经N层NSCT分解后获得的低频子带与2N个方向上的高频子带；
2)对各子带影像进行插值放大，最终获得a和b两个影像各自子带的插值结果；
3)将插值后的各子带影像分别进行积分平差处理求取最优解，获得平差后的各子带影像；针对二维NSCT系数建立积分平差模型，模型表达式为：
其中Li为i点的像素值，Vi为改正数，vi为点i积分空间的范围，f(x,y)为客体表面函数，(x,y)表示i点坐标，ds对应dxdy；
以上对各NSCT系数积分平差的过程可表示：
其中分别表示子带iLa与iLb，iHa1与与的积分平差结果，In-Ad表示积分平差过程；
4)对平差后的各NSCT系数做NSCT逆变换，得最终积分平差法的重建结果；
其中SR为超分辨重建结果，iNSCT表示NSCT逆变换。
上述的一种基于NSCT的积分平差遥感影像超分辨率重建方法，优选的，所述步骤2)中采用双三次插值方法对NSCT分解后的各个频段进行双三次插值：
其中Bi-In表示双三次插值的过程，最终获得两影像各子带的插值结果。
与现有技术相比，本发明的优点在于：本发明的基于NSCT的积分平差遥感影像超分辨率重建方法是一种实用的超分辨率重建算法，重建能力强，并且鲁棒性好。
附图说明
图1为NSCT非下采样滤波器组的示意图。
图2为NSCT频域划分示意图。
图3为误差方程系数阵。
图4为第一组原始影像。
图5为对图4进行IBP法重建的结果。
图6为对图4进行POCS法重建的结果。
图7为对图4进行文献[1]重建方法的结果。
图8为对图4进行本发明的重建方法的重建结果。
图9为第二组原始影像。
图10为对图9进行IBP法重建的结果。
图11为对图9进行POCS法重建的结果。
图12为对图9进行文献[1]重建方法的结果。
图13为对图9进行本发明的重建方法的重建结果。
图14为对图4进行隔点降采样获得的tj01A影像图。
图15为对图4进行隔点降采样获得的tj01B影像图。
图16是混合噪声为0.02/0.001时对图4采用IBP法重建结果。
图17是混合噪声为0.02/0.001时对图4采用POCS法重建结果。
图18是混合噪声为0.02/0.001时对图4采用文献[1]法重建结果。
图19是混合噪声为0.02/0.001时对图4采用本方法重建的结果。
图20为表3积分平差法在NSCT分解不同级数与不同方向数时的结果图。
具体实施方式
为了便于理解本发明，下文将结合较佳的实施例对本发明作更全面、细致地描述，但本发明的保护范围并不限于以下具体的实施例。
除非另有定义，下文中所使用的所有专业术语与本领域技术人员通常理解的含义相同。本文中所使用的专业术语只是为了描述具体实施例的目的，并不是旨在限制本发明的保护范围。
实施例1
一种基于NSCT的积分平差遥感影像超分辨率重建方法，包括以下步骤1)设置NSCT变换级数N和方向数2N，至少对两幅原始遥感影像a和b进行NSCT分解获得各NSCT系数。
NSCT通过将非降采样方向滤波器、塔式滤波器相结合可实现真正意义上的多分辨率分析。利用滤波器的下采样代替数据的下采样过程，使得分解后影像的尺寸不会发生变化。该过程如图1和图2所示：
经NSCT将影像进行N层分解，得1个低频子图与2N个方向高频子图。通过图1可知方向滤波器的方向数由分解层数来确定，当进行N层分解时对应2N个方向。此分解过程体现了NSCT多方向、多尺度分解等特性，是影像多分辨率分析的有用工具。其分解过程可用式(1)表示。
(1)
其中L1，H1，H2...分别为低频子带与2N个方向上的高频子带，NSCTN为N层NSCT分解，i为原图。
其中与分别为低分辨率影像a、b经N层NSCT分解后获得的低频子带与2N个方向上的高频子带；
2)对各子带影像进行插值放大，最终获得a和b两个影像各自子带的插值结果。
重建过程涉及到影像图幅的放大，为使影像的图幅增大，过程中需要进行插值处理。本发明中采用双三次插值方法对NSCT分解后的各个频段进行双三次插值：
其中Bi-In表示双三次插值的过程，最终获得两影像各子带的插值结果。
3)将插值后的各子带影像分别进行积分平差处理求取最优解，获得平差后的各子带影像；
针对二维NSCT系数建立积分平差模型，模型表达式为：
其中Li为i点的像素值，Vi为改正数，vi为点i积分空间的范围，f(x,y)为客体表面函数，(x,y)表示i点坐标，ds对应dxdy。
通过以上模型可得观测值和待求量之间是一种积分表达。该模型适用于空间域，此时观测值表示的是影像的灰度值；同时模型也适用于频率域，通过小波、轮廓波、NSCT为工具变换获得同时具有空间域、频率域特性的各高低频子图，观测值便可以是这些高低频子图上的某一坐标位置处的对应值。
客体表面函数f(x,y)的形式可以是多项式、小波函数、傅里叶函数等，取f(x,y)为三次多项式时(4)式化为：
其中Li为观测值，Vi、vi同式(4)，XT＝(a0,a1,a2,a3,a4,a5,a6,a7,a8,a9)表示三次多项式的系数值，Bi为误差方程的系数矩阵，且Bi＝(Bi,1,Bi,2,Bi,3,Bi,4,Bi,5,Bi,6,Bi,7,Bi,8,Bi,9,Bi,10)，各项求取式：
当点i(m,n)为子带影像上确定的某点时如图3：
点i(m,n)处的参数Bi,2,Bi,3,Bi,4,Bi,5,Bi,6,Bi,7,Bi,8,Bi,9,Bi,10的求取过程可结合图3中Bi,1的积分过程与式(6)中的积分函数依次求得，积分空间vi为影像左上角原点至i点的空间范围，这样通过积分求得该点误差方程的系数矩阵Bi。对于该点与其周围8点均可以建立(4)式的积分模型，共产生9个观测方程，则对应另一观测影像的对应窗口范围内也可产生9个观测方程。对于上述10参数Bi模型只要两组以上对应NSCT系数就能进行平差求解。根据最小二乘原理，(4)式须满足VTPV＝min原则，此时法方程：
NBBX-W＝0
NBB＝BTPB,W＝BTPL
(7)
其中B为点i对应的系数阵Bi的集合，此时近似取权阵P为单位阵，L为点i对应的观测阵Li的集合。
根据求得的结果可以获得最终观测点的平差值：
这样通过建立积分平差模型，再进行平差解算，以求得NSCT系数的最优解值。以上对各NSCT系数积分平差的过程可表示：
其中分别表示子带iLa与iLb，iHa1与与的积分平差结果，In-Ad表示积分平差过程。
4)对平差后的各NSCT系数做NSCT逆变换，得最终积分平差法的重建结果；
其中SR为超分辨重建结果，iNSCT表示NSCT逆变换。
对比实验
不含噪声实验
采用两组遥感影像数据验证所提出方法的有效性。第一组为原影像tj01(256×256)进行隔点降采样获得的两幅影像tj01A、tj01B。第一组影像如图4所示，图14和图15分别为对图4进行隔点降采样获得的tj01A和tj01B的影响图。第二组为原影像tj02(512×512)进行隔点降采样获得的两幅影像tj02A、tj02B，图9为第二组原始影像图。隔点降采样是模拟遥感影像的成像过程，这样处理能使影像对之间存在信息的交互，符合实际遥感影像的成像原理。同时可将重建获得的影像与原图tj01、tj02来定量对比，以分析方法的优劣。
实验中本文算法先对影像(tj01A,tj01B)、(tj02A,tj02B)进行亚像素级配准，获得影像间各像素的配准信息。为能够很好地对影像细节进行平差，NSCT采用1级2向分解，方向滤波器与塔式滤波器分别选取效果较好的'dmaxflat7'与'maxflat'，然后与经典的迭代反投影法(Iterative Back Projection,IBP)(Gevrekci M,Gunturk B K.Super Resolution Approaches for Photometrically Diverse Image Sequences[C]//ICASSP(1).2007:757-760.)与凸集投影法(Projections Onto Convex Set,POCS)(Y.Zhou,"A POCS method for iterative deblending constrained by a blending mask,"Journal of Applied Geophysics,vol.138,no.pp.245-254,2017.)重建方法及文献[1](YANG Zhenping,NIU Hai-jun,DING Lei.A New Image Super-Resolution Restoration Method By NSCT[J].Journal of XidianUniversity,2012,39(1):34-37.)算法进行对比。
图5-图8分别为对图4进行IBP法、POCS法、文献[1]方法和本发明方法重建的结果。
图10-图13分别为对图9进行IBP法、POCS法、文献[1]方法和本发明方法重建的结果。
实验选用峰值信噪比(peak signal to noise ratio，PSNR)与均方根误差(root-mean-square error，RMSE)作为重建结果的定量评价指标，两参数均能客观反映重建影像与原参考影像间的差异，用来评价重建影像分辨率的变化与损失。实验结果见表1。
表1不含噪声时各超分辨方法重建结果
(a)tj01A-tj01B (b)tj02A-tj02B
含有噪声的实验
为验证本发明的算法在实际遥感影像重建过程中的鲁棒性，进行含噪声实验。对影像(tj01A，tj01B)加入遥感影像中常见的四种混合噪声类型：0.01/0.0005、0.02/0.001、0.05/0.001、0.1/0.01，每组数据的前一个值代表的是高斯噪声方差，后一个值代表的是脉冲噪声的密度。然后分别用POCS法、IBP法、文献[1]算法和本文算法进行对比实验。实验结果见表2。图14和图15分别为对图4进行隔点降采样获得的tj01A和tj01B的影响图。图16-图19分别为混合噪声为0.02/0.001时对图4采用IBP法、POCS法、文献【1法】和本发明方法重建的结果。
表2含混合噪声时各超分辨方法的重建结果
PSNR与RMSE是用来衡量两幅低分辨率影像重建得到的高分辨率影像与原影像分辨率的对比结果，是目前广泛被接受为超分辨率重建结果的客观评价指标，其中PSNR值越大且RMSE值越小表示重建效果越好，反之则重建效果较差。针对不含噪声的两对影像对(tj01A,tj01B)、(tj02A,tj02B)从图4中可以看出本文方法的重建效果无论是在建筑物线条的平滑度还是在影像总体清晰度上都较IBP法、POCS法、文献[1]方法等重建方法优越，且重建结果的PSNR指标要明显高于其他方法，RMSE这一项指标也是明显低于其他方法的，证明本文算法的有效可行。
从各方法鲁棒性对比实验中可以看出，针对遥感影像中常见的混合噪声情况，虽然在不存在噪声时IBP法的重建效果较POCS法要优越，但在各类混合噪声情况下能看出IBP法重建过程的鲁棒性明显不及POCS法。而且从图4的定性分析与表2中实验结果的定量对比可以看出本文重建算法的鲁棒性要明显优于IBP法与POCS法，而且在重建过程中并未用到影像降噪处理，这说明本文算法的鲁棒性较强，具有抗噪重建的优良性能。
通过不含噪声和含有常见混合类型噪声的遥感影像进行对比实验，从定性的视觉分析和定量的指标数据对比可得本文重建算法用于遥感影像的超分辨率重建其效果要明显优于经典的的IBP法与POCS法，同时也优于文献[1]方法，并且本文重建方法在各类噪声下的鲁棒性也明显较优，为此通过实验进一步证明了本文重建方法在常见遥感影像混合噪声情况下的可行性与有效性。
NSCT分解尺度实验
为了研究NSCT的分解级数对本文重建算法结果的影响，选取NSCT分解级数1至4，方向数1至4，共16种组合进行重建实验。实验选取的影像为(tj01A,tj01B)，结果见表3。选取PSNR与均方误差(mean square error，MSE)作为重建结果的定量评价指标。
表3积分平差法在NSCT分解不同级数与不同方向数时的结果
其中PSNR值越大MSE值越小，表示重建效果越好。从表3和对应图20中可以看出，随着NSCT分解级数和方向数的增加，积分平差法的效果逐渐降低，且每一级分解范围内趋势变化较缓慢，不同级数间的变化较大，说明级数对重建效果的影响要明显高于方向数对重建效果的影响，且积分平差法与NSCT结合时适用于低分解级数和低方向数的情况。