结合波段相关性的FY-3C MERSI异常条带检测方法

阎宏波，樊宪磊，瞿瑛

(东北师范大学 地理科学学院，长春 130024)

0 引言

风云三号是我国最新一代极轨气象卫星，可以实现对全球大气、海洋和陆表特征参量全天候、多光谱和定量长期观测。风云三号C星(FY-3C)作为02批的首发卫星，其核心遥感传感器的技术指标和性能相对于01批的A/B星有了进一步的提升。FY-3C卫星搭载的中分辨率光谱成像仪(medium resolution spectral imager，MERSI)共有20个光谱波段，包括从可见光到短波红外范围(0.4～2.1 μm)的19个波段和1个热红外范围(10～12.5 μm)波段，其中有5个250 m空间分辨率的波段和15个1 000 m空间分辨率的波段[1-2]。其光谱波段及性能指标见表1，光谱响应曲线如图1所示。

图1 FY-3C MERSI光谱响应函数

表1 FY-3C MERSI光谱波段及性能指标

MERSI反射率数据集是进行高精度定量反演大气、陆地和海洋特征参量产品的基础数据，MERSI反射率数据集存在的异常条带会降低数据的可用性，影响基于影像数据分析结果和结论的可靠性。因此，在定量处理与分析前，需要首先对FY-3C MERSI图像进行异常条带检测。通过图像判断和解译发现，FY-3C MERSI 反射率数据集中存在大量的异常像元，其中异常条带是最主要的、频率最高的异常像元出现形式，包括异常行条带和左侧边缘异常列条带2种(图2)。

图2 异常条带示意图

目前遥感图像异常条带检测与去除大致可分为空间域滤波和频率域滤波2种方法[3]。前者主要基于空间域模板和滤波算子对遥感影像的像元值进行处理[4-5]，主要包括直方图匹配法和矩匹配法；后者将遥感影像从空间域转换到频率域，在频率域选择合适的滤波算子进行处理，最后把处理结果再转换回到空间域，完成遥感影像的滤波处理，主要包括傅里叶变换和小波分析等方法[6-7]。这些方法虽然都可以对遥感图像异常条带检测，但是其处理效果通常受到遥感影像数据和研究区特性的制约。如直方图匹配法对遥感影像的均一性要求较高[8]，而传统的矩匹配法对于范围较小且地物覆盖类型较复杂的遥感影像，处理后会产生“带状效应”[9-11]。基于傅里叶变换的条带检测方法需要在频率域中找到合适的滤波算子对噪声信号进行剥离。如果滤波算子不合适，则条带噪声检测去除效果较差[12]。小波分析方法可以很好地剥离噪声信号和高频信息，但处理后的遥感影像容易出现边缘模糊和振铃等现象，造成遥感影像的视觉失真[13]。

在MERSI数据异常条带检测与去除研究中，通常使用固定阈值法和数学形态法2种方法。其中，固定阈值法是假设反射率的取值范围在0～1之间，因此将0和1作为阈值进行异常条带检测，这种方法比较简单但准确率相对较低。而数学形态法是通过构建相邻行差矩阵，设置阈值来检测异常条带。该方法已经被应用于FY-3A MERSI植被指数产品异常条带的检测和去除[14]，取得了比较好的检测结果。但该方法在阈值选取和定位异常行条带位置方面缺乏普适性，并且假设异常条带间隔均为一个扫描带，而该假设在FY-3C MERSI数据中并不成立。因此，本研究提出了一种结合反射光谱波段相关性的FY-3C MERSI数据异常条带检测识别方法，基于反射率影像数据进行了验证，并将检测结果与数学形态法的检测结果进行了对比分析。

1 异常条带检测方法

1.1 异常条带出现规律分析

异常条带是卫星传感器在扫描记录数据时，受到了信号干扰或者系统错误的影响，使某些遥感波段数值出现异常或缺失而导致的。通过大量影像数据解译和判读分析，发现FY-3C MERSI数据的异常条带具有很强的规律性。

1)并非所有的光谱波段均存在异常条带，而是左侧边缘异常条带仅出现在1 000 m空间分辨率的地表反射率数据集(EV_250_Aggr.1KM_RefSB)中的第2和第4波段，且在这2个波段的异常列条带之前的所有列的值均为0，其出现的位置为数据的前10列 (如图3所示，其中，红色箭头所在列为异常列，黑色列的值为0)，而异常行条带主要集中在数据集中的第1～4波段，其余波段虽然也存在异常行条带，但异常像元值均为1，可以通过阈值检测直接去除，因此在本研究中此类条带不作为研究对象。如2014年8月15日的1景FY-3C MERSI影像，在EV_1KM_RefSB数据集中，在第6波段的第1 202、1 205、1 206、1 208和1 210行出现异常值。

图3 部分影像第2、4波段左侧边缘列条带

2)无论是行条带还是列条带均与正常行、正常列之间存在着非常明显的分割界线，且分割界线均为水平或垂直的直线 (图4)。

图4 异常条带检测流程图

3)与FY-3A/B MERSI数据比较，发现FY-3C MERSI 数据中不存在“Z”字形式和斜线形式的异常条带，这主要是因为FY-3C卫星搭载的MERSI传感器在各项性能上比FY3-A/B上搭载的MERSI有较大的提高[15]。

1.2 基于波段相关性的异常条带检测方法

1)波段相关性统计分析。本研究基于波段相关性，对FY-3C MERSI反射率数据集第1～4波段中存在的异常条带进行检测。具体思路是：分别寻找与第1～4波段反射率值相关性最好的匹配波段。通过对FY-3C MERSI图像进行分析，发现第 1、2、4波段分别与第10、12、16波段具有很强的相关性，而第3波段与第 13和第14波段的相关性都很强，因此采用第13和第14波段的均值((波段13+波段14)/2)与第4波段匹配来建立相关关系 (表2)。本研究基于波段相关性原理，计算待检测波段与匹配波段反射率值的差值，然后通过对行/列差值的和设置阈值，识别异常行/列条带出现的准确位置。该算法的实施步骤是：首先对影像进行辐射定标处理；然后结合波段相关性和异常列条带规律，分别对异常行列条带进行检测；最后识别出所有异常行列所在位置(图4)。

表2 波段匹配统计结果

2)辐射定标处理。为了便于后续异常条带的检测分析，首先对FY-3C MERSI 数据进行辐射定标，将原始DN值定标转换为反射率值。在定标过程中，可将值为1的异常条带直接掩膜。其中，辐射定标可以通过式(1)、式(2)实现。

DN**=slope×(DN*-intercept)

(1)

ρ=(k0+k1×DN**+k2×DN**∧2)/cosθ

(2)

式中：DN**为辐射定标值；DN*为对原始值DN经过多探元归一化处理后的计数值；ρ表示反射率；k0、k1和k2为相应的定标系数；θ为太阳天顶角；slope和intercept为对应科学数据集中的定标斜率和截距。

3)异常行条带检测。对FY-3C MERSI影像中的异常行条带的检测算法具体步骤如下。

①定义波段差值矩阵DRow(式(3))。

DRow=Bandi-Bandj，
(i=10，12，0，16；j=1，2，3，4)

(3)

通过此矩阵构建匹配波段与对应待检测波段的反射率之差。其中，该矩阵是一个大小与波段反射率矩阵相同的二维矩阵。

②定义一个行数为N的列向量S(式(4))。

(4)

此矩阵计算DRow矩阵按行求和的结果，是一个N行的列向量，可以基于向量在下一步中确定分割阈值来识别异常行条带。

③定义一个0，1标记矩阵FA，对不同的波段组合分别确定典型分割阈值Ti1进行二值化，见式(5)。

(5)

FA矩阵即为异常行的位置标记矩阵。当满足FAx=1时的x0就是异常行所在的行号。

4)异常列条带检测。基于FY-3C MERSI异常行列出现规律分析可知，异常列条带规律性更好，因此其检测也相对简单，其具体检测算法步骤如下。

①定义列差值矩阵DCol，见式(6)。

DCol=|Bandi(x，y)-Bandj(x，y)|
(0≤x≤N，0≤y≤10；x，y∈N+；
i=2，4，j=12，16)

(6)

此矩阵是对波段i的反射率矩阵的前10列逐列做差并取绝对值，是一个大小为10×N二维矩阵。

②定义标记矩阵FB，见式(7)。

(7)

该矩阵用来判断列差值与匹配波段对应列值是否相等，进而确定异常列条带的列号。当FBx，y-1=0且FBx，y=1时，y值即为异常列的列号，第y列为异常列条带。

1.3 数学形态法

目前用于检测MERSI数据异常行条带方法主要是数学形态法识别方法，该方法已经被广泛应用于图像处理。本研究将其作为对照算法，对基于波段相关性的条带检测方法进行评价。设为N×M的反射率影像数据矩阵，其中x，y分别代表行下标和列下标，则数学形态法的异常行条带检测算法实现步骤如下。

1)计算影像数据的“相邻行差值”矩阵D，见式(8)。

D=|Rx，y+Rx+1，y|，
(0≤x≤N-1；0≤y≤M；x，y∈N+)

(8)

此矩阵的意义是，Rx，y相邻2行的反射率值的差是一个行数为N-1、列数为M的二维矩阵，其中N+代表正整数。

2)定义一个行数为N-1、列数为M的掩膜矩阵FC，且设置阈值T2，根据式(9)进行二值化，标记每个像元的像元值是否与其同列下行的像元值存在显著差异。

(9)

定义一个行数为N-1、列数为M的求和矩阵N，同时固定行标x=x0，0≤x≤N-1，对任意一行的有效像元值进行累加求和，如式(10)所示。

N|x=x0=∑FC|x=x0，
(Rx0，y×Rx0，y+1>0，0

(10)

3)定义一个行数为N-1、列数为1的掩膜矩阵FD，且设置阈值T3，通过式(11)进行二值化，提取出异常行的分割线标记矩阵。

(11)

在影像矩阵中，利用得到的分割线标记矩阵，将在分隔线之间的像元判断为异常像元。

1.4 条带检测效果评价与比较

选取2014年FY-3C MERSI数据的67景影像(数据获取于风云卫星遥感数据服务网http：//satellite.nsmc.org.cn)作为实验数据，分别对影像数据采用本文提出的基于波段相关性的行列条带检测方法和数学形态法，进行异常行列条带检测。将自动检测结果与目视判读进行比较，统计检测出的异常行数、准确率和漏检率，对算法的检测效果进行综合评价。

2 结果与讨论

基于67景FY-3C MESRI影像数据进行条带检测实验，其中图像数据共有155条异常行条带，134条异常列条带。基于本文算法可以检测出152条异常行条带，检测率为98.1%，误检率为0，漏检率为1.9%，并且可以检测出所有的列异常条带(列检测率为100.0%)；而采用数学形态法仅能检测出100条异常行条带，检测率为64.5%，误检率为45.2%，漏检率为35.5%，如表3、表4所示。图5展示了基于本文算法和数学形态法对FY-3C MERSI影像进行异常条带检测的一些典型案例，其中图5(a)、图5 (b)显示对于第1波段，2种方法均能很好地将异常行条带检测出来；图5(c)、图5 (d)显示对于第2波段，数学形态法漏检了1行，其中黄色箭头所在行即为未能正确检测出来的异常行条带；图5(e)、图5 (f)显示对于第3波段，数学形态法不仅未能正确检测出异常行条带，反而误检了1行，其中红色箭头所在位置为误检行，黄色与绿色箭头指示了漏检的异常行条带；图5(g)、图5 (h)显示对于第4波段，数学形态法误检了1行，其中黄色箭头为误检的异常行条带。

表3 行条带检测结果

表4 列条带检测结果

通过表3、表4和图5可以发现，基于数学形态法进行FY-3C MERSI数据条带检测，存在比较严重的误检和漏检等问题，其中误检率高达45.2%，漏检率也有35.5%。造成这些问题的原因主要有2个：首先数学形态法是在同一波段数据中逐行做差，因此基于这种方法不能使异常条带处差值与正常行处的差值形成明显的差距，特别是在异常行与相邻行差异不明显的地方更是如此，这就很难确定一个具有普适性的分割阈值来批量处理影像数据，并且第2个阈值还受到第1个阈值的影响，使其普适性大大降低；其次数学形态法最终检测出的是正常行和异常行之间的分割界线，理论上分割界线之间的就是异常行，但是如何确定分割界线的起止位置是一个难点。研究发现，相邻分割界线之前的行数并不都是一个扫描带的宽度，因此尽管找到了异常行与正常行之间的分割界线，但是要准确确定异常行位置仍然相对困难。

图5 行检测结果对比图

本文提出的算法本质上只需要确定一个阈值，并且正常值与异常值区分度较大，容易确定普适性强的阈值。此外，本文算法最终直接标记的是异常行、列所在的行、列号，而非分割界线，这就避开了确定的分割界线起止顺序的难点。但是在实验中发现，无论是数学形态法还是本文算法，对于极少出现的断续条带(即在影像一行中，只有若干异常像元，而剩余像元均为正常像元的情况)，均无法正确进行检测，这需要在研究中根据图像特点设置更为复杂的阈值判断规则才能够加以解决。

3 结束语

本研究提出了一种结合波段相关性的FY-3C MERSI数据异常条带检测方法，基于反射光谱波段之间的相关关系，构建波段差矩阵，最终检测出异常行列条带所在位置。基于67景FY-3C MERSI数据的条带检测实验表明，本文算法适用于大批量检测FY-3C MERSI数据异常条带，显著高于数学形态法的检测率，无误检条带，大大降低了误检率和漏检率，为风云系列卫星反射率数据集进一步的应用研究奠定良好的基础。