跨轨宽幅扫描相机的条带噪声抑制方法及系统与流程

1.本发明涉及一种遥感影像去噪方法，具体地，涉及跨轨宽幅扫描相机的条带噪声抑制方法及系统。


背景技术：

2.跨轨宽幅扫描相机因可以实现快速获取全球影像，故而广泛应用于覆盖全球的遥感成像领域，典型的仪器包括terra卫星、aqua卫星上装载的中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer,modis)、npp卫星上装载的可见光红外成像辐射仪(visible infrared imaging radiometer,viirs)、fy-3卫星上装载的中分辨率光谱成像仪(medium resolut ion spectral imager)等。
3.跨轨宽幅扫描相机在进行定量应用时，需对获取的图像数据进行高精度辐射定标，然而由于线阵探测器的非均匀性以及在轨探测器响应的变化，导致获取的图像存在条带噪声，在经过绝对辐射定标系数转化为物理量后，图像仍然存在一定程度的条带噪声。
4.针对条带噪声，目前已有的噪声抑制方法包括：(1)基于统计的方法，如直方图匹配、矩匹配等：(2)基于滤波的方法，如带通滤波、小波域滤波、约束最优化等；(3)基于机器学习的方法。
5.针对跨轨宽幅扫描相机，现有文献[1](preesan rakwat in,wataru takeuchi,yoshifumi yasuoka,stripe noise reduction in modis data by combining histogram matching with facet filter,ieee transact ions on geoscience remote sensing,45(6),2007.)给出了一种结合直方图匹配与小面滤波的条带噪声抑制方法。文献[2](marouan bouali,said ladjal,toward optimal destriping of modis data using a unidirect ional variational model,49(8),2011.)给出了一种基于单向变分的去条纹方法。文献[3](gang zhou,houzhang fang,cen lu,robust destriping of modis and hyperspectral data using a hybrid unidirectional total variation model,optik,126(7-8),2015.)给出了一种基于单向变分以及权重矩阵的去条纹方法。文献[4](ranil basnayake,erik bollt,nicholas tufillaro,regularization destriping of remote sensing imagery.nonlinear processes in geophysics,24(3),2017.)给出了一种基于单向变分及约束阈值模板的去条纹方法。文献[5](pengfei xiao,yecai guo,peixian zhuang,removing stripe nose from infrared cloud images via deep convolutional networks,ieee photonics journal,10(4),2018.)给出了一种基于卷积神经网络的条带噪声去除方法。
[0006]
跨轨宽幅扫描相机幅宽一般可达2000km以上，由于地球曲率的影响，其图像具有一定的特殊性，一个扫描周期的图像对应的地球区域呈现“bowtie”(领结)形状，图像在偏离星下点较大的位置本就存在周期性的非连续性，这种非连续性应与条带噪声导致的非连续性进行区分，现有的条带噪声抑制方法均建立在图像相邻像元之间的连续性基础之上，未在去噪的同时考虑对“bowtie”效应引起的非连续性进行区别和保持。
[0007]
针对现有条带噪声抑制方法的不足，本发明采用结合空间位置重排与单向变分的最优化求解方法，对图像沿轨方向条带噪声进行抑制的同时尽量保证对原图细节的保留。


技术实现要素：

[0008]
针对现有技术中的缺陷，本发明的目的是提供一种跨轨宽幅扫描相机的条带噪声抑制方法及系统。
[0009]
根据本发明提供的一种跨轨宽幅扫描相机的条带噪声抑制方法，包括：
[0010]
步骤s1：基于原始影像数据矩阵计算每一像素在沿轨方向的位置，形成位置矩阵；
[0011]
步骤s2：根据位置矩阵计算排序索引矩阵，基于排序索引矩阵构建重排函数与还原函数；
[0012]
步骤s3：基于跨轨方向变分以及沿轨方向重排后的变分数据，构建最优化问题并进行迭代求解。
[0013]
优选地，所述步骤s1采用：基于尺寸为m行、n列的原始影像数据矩阵v计算每一像素v(i,j)在沿轨方向的位置y(i,j)，其中，i表示行号，i＝1,2,......,m；j表示列号，j＝1,2,......,n；形成尺寸为m行，n列的位置矩阵y。
[0014]
优选地，所述步骤s1采用：
[0015]
当跨轨宽幅扫描相机线阵共包含p个像元，则原始影像数据矩阵v以p行数据为一个扫描周期，位置矩阵y的计算方法按照扫描周期进行计算，若y(i,j)为第t个周期的第t像元，则行号i＝(t-1)p t；
[0016]
y(i,j)＝(t-0.5(p 1))α(j) 0.5(p 1) (t-1)p
ꢀꢀꢀ
(1)
[0017]
其中，系数α(j)为距离放大因子，计算方法为
[0018]
α(j)＝l(j)/l0ꢀꢀꢀ
(2)
[0019]
其中，l0为相机星下点观测时探测器中心对应的视线距离；l(j)为第j列数据成像时刻探测器中心对应的视线距离。
[0020]
优选地，所述步骤s2采用：对位置矩阵y的每一列[y(1,j),y(2,j),
……
,y(m,j)]
t
，按照由小到大排序，并存储排序索引[b(1,j),b(2,j),
……
,b(m,j)]
t
，y(b(1,j),j)≤y(b(2,j),j)≤
……
≤y(b(m,j),j)；其中，t表示转置；由每一列的排序索引，得到尺寸为m行、n列的排序索引矩阵b，按照排序索引矩阵b对数据矩阵v进行重新排序的操作记为重排函数f，原始数据矩阵v重新排序后的矩阵记为v
*
；
[0021][0022]
将重排后的矩阵v
*
还原至原始位置的操作记为还原函数f-1
；
[0023]
v＝f-1
(v
*
)
ꢀꢀꢀ
(4)
[0024]
v(b(i,j),j)＝v
*
(i,j)
ꢀꢀꢀ
(5)。
[0025]
优选地，所述步骤s3采用：跨轨与沿轨方向的变分计算采用如下方法：
[0026]
[0027][0028]
其中，分别为跨轨和沿轨方向的梯度算子；
[0029][0030][0031]
其中，系数ε1、系数ε2为恒定的正常数；可由公式8同理获得。
[0032]
优选地，所述最优化问题采用：
[0033][0034]
其中，表示最优解；u表示最优化问题的自变量；λ表示权重系数；函数g
x
表示跨轨方向的图像变分，函数gy表示沿轨方向的图像变分；u
*
表示矩阵u按照排序索引矩阵b进行重新排列后的矩阵。
[0035]
优选地，所述权重系数λ用于权衡图像在跨轨方向的保真性与在沿轨方向的平滑性；权重系数λ可通过如下方法进行计算：
[0036]
λ＝gy(v
*
)/g
x
(v)
ꢀꢀꢀ
(11)。
[0037]
优选地，最优化问题求解采用如下迭代计算方法，初值k＝0，u0＝v，
[0038][0039]
其中，k为迭代步数，uk为第(k-1)步迭代得到的结果，为对uk进行重排的图像，当|u
k 1-uk|小于设定阈值时终止迭代。
[0040]
根据本发明提供的一种跨轨宽幅扫描相机的条带噪声抑制系统，包括：
[0041]
模块m1：基于原始影像数据矩阵计算每一像素在沿轨方向的位置，形成位置矩阵；
[0042]
模块m2：根据位置矩阵计算排序索引矩阵，基于排序索引矩阵构建重排函数与还原函数；
[0043]
模块m3：基于跨轨方向变分以及沿轨方向重排后的变分数据，构建最优化问题并进行迭代求解。
[0044]
优选地，所述模块m1采用：基于尺寸为m行、n列的原始影像数据矩阵v计算每一像素v(i,j)在沿轨方向的位置y(i,j)，其中，i表示行号，i＝1,2,......,m；j表示列号，j＝1,2,......,n；形成尺寸为m行，n列的位置矩阵y；
[0045]
所述模块m1采用：
[0046]
当跨轨宽幅扫描相机线阵共包含p个像元，则原始影像数据矩阵v以p行数据为一个扫描周期，位置矩阵y的计算方法按照扫描周期进行计算，若y(i,j)为第t个周期的第t像元，则行号i＝(t-1)p t；
[0047]
y(i,j)＝(t-0.5(p 1))α(j) 0.5(p 1) (t-1)p
ꢀꢀꢀ
(1)
[0048]
其中，系数α(j)为距离放大因子，计算方法为
[0049]
α(j)＝l(j)/l0ꢀꢀꢀ
(2)
[0050]
其中，l0为相机星下点观测时探测器中心对应的视线距离；l(j)为第j列数据成像时刻探测器中心对应的视线距离；
[0051]
所述模块m2采用：对位置矩阵y的每一列[y(1,j),y(2,j),
……
,y(m,j)]
t
，按照由小到大排序，并存储排序索引[b(1,j),b(2,j),
……
,b(m,j)]
t
，y(b(1,j),j)≤y(b(2,j),j)≤
……
≤y(b(m,j),j)；其中，t表示转置；由每一列的排序索引，得到尺寸为m行、n列的排序索引矩阵b，按照排序索引矩阵b对数据矩阵v进行重新排序的操作记为重排函数f，原始数据矩阵v重新排序后的矩阵记为v
*
；
[0052][0053]
将重排后的矩阵v
*
还原至原始位置的操作记为还原函数f-1
；
[0054]
v＝f-1
(v
*
)
ꢀꢀꢀ
(4)
[0055]
v(b(i,j),j)＝v
*
(i,j)
ꢀꢀꢀ
(5)；
[0056]
所述模块m3采用：跨轨与沿轨方向的变分计算采用如下方法：
[0057][0058][0059]
其中，分别为跨轨和沿轨方向的梯度算子；
[0060][0061][0062]
其中，系数ε1、系数ε2为恒定的正常数；可由公式8同理获得；
[0063]
所述最优化问题采用：
[0064][0065]
其中，表示最优解；u表示最优化问题的自变量；λ表示权重系数；函数g
x
表示跨轨方向的图像变分，函数gy表示沿轨方向的图像变分；u
*
表示矩阵u按照排序索引矩阵b进行重新排列后的矩阵；
[0066]
所述权重系数λ用于权衡图像在跨轨方向的保真性与在沿轨方向的平滑性；权重系数λ可通过如下方法进行计算：
[0067]
λ＝gy(v
*
)/g
x
(v)
ꢀꢀꢀ
(11)；
[0068]
最优化问题求解采用如下迭代计算方法，初值k＝0，u0＝v，
[0069][0070]
其中，k为迭代步数，uk为第(k-1)步迭代得到的结果，为对uk进行重排的图像，当|u
k 1-uk|小于设定阈值时终止迭代。
[0071]
与现有技术相比，本发明具有如下的有益效果：
[0072]
1、本发明针对跨轨宽幅扫描相机的条带噪声，利用空间位置关系，在去除条带噪声的同时保留地球曲率造成的图像非连续性，可以对探测器响应不一致引起的条带噪声进行抑制，提升图像质量且不影响定量化应用；
[0073]
2、本发明方法合理、计算简单、实施简易，能够普遍应用于跨轨宽幅扫描相机的条带噪声抑制；
[0074]
3、本发明通过地理空间位置与沿轨平滑相结合的方式，实现了条带噪声的抑制，提升了跨轨宽幅扫描相机的图像质量，实现了为图像高精度定量化应用提供输入。
附图说明
[0075]
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
[0076]
图1为本发明的流程图。
[0077]
图2为跨轨宽幅扫描相机图像空间位置特性示意图。
[0078]
图3为某卫星跨轨宽幅扫描相机一级影像。
[0079]
图4为采用本发明方法对图3进行条带噪声抑制后的遥感影像。
[0080]
图5为条带噪声抑制前后图像边缘一列重排后的像素值对比。
具体实施方式
[0081]
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
[0082]
实施例1
[0083]
本发明所要解决的技术问题是提供一种跨轨宽幅扫描相机的条带噪声抑制方法，利用“bowtie”效应对图像沿轨方向进行平滑，并尽量保证对原图细节的保留。
[0084]
由于线阵探测器的非均匀性以及在轨探测器响应的变化，导致跨轨宽幅扫描相机获取的图像存在条带噪声，表现为噪声沿扫描方向分布，由于条带噪声沿跨轨方向分布，则图像跨轨方向的图像边缘远小于图像沿轨方向的图像边缘。利用该特性，本发明跨轨宽幅扫描相机的条带噪声抑制方法流程如图1所示，结合图像空间位置关系建立最优化问题。进行噪声抑制后的图像u为如下最优化问题的解：
[0085][0086]
其中，为最优解，u为最优化问题的自变量，λ为权重系数。函数g
x
为跨轨方向的图像变分，函数gy为沿轨方向的图像变分。u
*
为矩阵u按照排序索引矩阵b进行重新排列后的矩阵。
[0087]
以线阵进行跨轨宽幅扫描的相机图像受地球曲率影响，一个扫描周期的图像呈“bowtie”(领结)形状，多个周期排列后在偏离星下点较大的区域即存在周期性图像不连续，如图2所示。在针对探测器不一致性引起的条带噪声进行抑制时，需对“bowtie”效应引起的图像非连续性进行保留。
[0088]
若跨轨宽幅扫描相机线阵共包含p个像元，则图像v以p行数据为一个扫描周期，所述步骤1中位置矩阵y的计算方法按照扫描周期进行计算，若y(i,j)为第t个周期的第t像元，即行号i＝(t-1)p t，
[0089]
y(i,j)＝(t-0.5(p 1))α(j) 0.5(p 1) (t-1)p
ꢀꢀꢀ
(公式2)
[0090]
其中，系数α(j)为距离放大因子，计算方法为
[0091]
α(j)＝l(j)/l0ꢀꢀꢀ
(公式3)
[0092]
l0为相机星下点观测时探测器中心对应的视线距离，l(j)为第j列数据成像时刻探测器中心对应的视线距离。
[0093]
对位置矩阵y的每一列[y(1,j),y(2,j),
……
,y(m,j)]
t
，按照由小到大排序，并存储排序索引[b(1,j),b(2,j),
……
,b(m,j)]
t
，即有y(b(1,j),j)≤y(b(2,j),j)≤
……
≤y(b(m,j),j)。由每一列的排序索引，得到尺寸为m行、n列的排序索引矩阵b。将按照排序索引矩阵b对数据矩阵进行重新排列的操作记为重排函数f，原始数据矩阵v重新排列后的矩阵记为v
*
。
[0094][0095]
将重排后的矩阵v
*
还原至原始位置的操作记为还原函数f-1
，v＝f-1
(v
*
)，即有
[0096]
v(b(i,j),j)＝v
*
(i,j)
ꢀꢀꢀ
(公式5)
[0097]
原始图像在空间位置上具有跨轨方向连续而沿轨方向在远离星下点时逐渐不连续的特性，故而跨轨与沿轨方向的变分计算采用如下方法：
[0098][0099][0100]
分别为跨轨和沿轨方向的梯度算子。
[0101]
[0102][0103]
系数ε1、系数ε2为恒定的正常数，用于确保(公式1)所描述的最优化问题在定义域内均可微。可由(公式8)同理获得。不失一般性地，可以取ε1＝1
×
10-6
；ε2＝1
×
10-6
。
[0104]
(公式1)中的权重系数λ用于权衡图像在跨轨方向的保真性与在沿轨方向的平滑性。λ取值过大可能导致图像损失细节，λ取值过小可能导致图像条带噪声抑制不够。当图像中的条带噪声强度不一时，需要权重系数λ可以自适应地变化。故权重系数λ可通过如下方法进行计算：
[0105]
λ＝gy(v
*
)/g
x
(v)
ꢀꢀꢀ
(公式10)
[0106]gy
(v
*
)的计算方法可由(公式7)、(公式9)同理获得，g
x
(v)的计算方法可由(公式6)、(公式8)同理获得。
[0107]
(公式1)中的最优化问题求解采用梯度下降算法，第k次迭代过程如下，
[0108][0109]
初值k＝0，u0＝v，uk为第(k-1)步迭代得到的结果，为对uk进行重排的图像，每次迭代判断u
k 1
与uk的差异，当|u
k 1-uk|小于设定阈值时终止迭代。
[0110]
本发明提供的跨轨宽幅扫描相机的条带噪声抑制系统，可以通过本发明提供的跨轨宽幅扫描相机的条带噪声抑制方法中的步骤流程实现。本领域技术人员，可以将所述跨轨宽幅扫描相机的条带噪声抑制方法理解为跨轨宽幅扫描相机的条带噪声抑制系统的一个优选例。
[0111]
下面结合某卫星遥感图像对本发明方法进行验证，图3所示为某卫星跨轨宽幅扫描相机一级遥感数据影像，一级数据已经进行过绝对辐射标定，但由于非理想因素的存在，图像中仍存在条带噪声。图4为采用本发明方法对图3进行条带噪声抑制后的遥感影像。图4表明条带噪声已经得到抑制，且在图像左右两侧偏离星下点较远的区域，与扫描周期相关的非连续性得到了保留。为进一步检查图像左右两侧条带噪声的抑制情况，分别将图3与图4选取边缘一列(第5列)进行重排后比较，结果如图5所示。图5表明在沿轨方向，噪声抑制后的图像重排后连续性较好，且对原始图像像素值的修正较小，利于该一级遥感数据影像的定量化应用。
[0112]
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以，本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件，而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构；也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
[0113]
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述
特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本技术的实施例和实施例中的特征可以任意相互组合。