基于RNMU的多源星载SAR影像融合与土地覆盖分类

李 俐 陈琦琦 张 超 尤淑撑 魏 海 付 雪

(1.中国农业大学土地科学与技术学院，北京 100083；2.农业农村部农业灾害遥感重点实验室，北京 100083；3.中国国土勘测规划院，北京 100035；4.中国农业大学信息与电气工程学院，北京 100083)

0 引言

及时准确地获取土地覆盖信息，能够为土地资源的合理利用与监管提供有力的数据支持[1]。遥感技术由于具有覆盖范围广、获取频率高的优势，在大面积土地覆盖信息提取方面得到了广泛应用[2-3]。光学遥感影像凭借其与人眼视觉系统接近、直观易解译的特点首先被广泛应用于土地覆盖分类[4-6]。然而，光学数据的获取受到光照条件的限制，在云、雨、雪、雾等天气下，存在数据获取困难、品质不佳等问题。以主动方式工作的合成孔径雷达(Synthetic aperture radar，SAR)摆脱了光照的影响，能够为土地覆盖分类提供全天候、全天时微波遥感信息[7-9]。因此，SAR数据在土地覆盖分类中的应用越来越受到重视。自2014年以来，相继发射的Sentinel-1A/B和高分三号(GF-3)等SAR卫星，提供了越来越多的高分辨率、多极化、多波段的SAR数据，大大提高了SAR数据在土地覆盖监测应用中的数据获取能力。

然而，单星单极化SAR影像的分辨率与影像可解译能力限制了SAR数据在土地覆盖分类中的识别效果。多源遥感数据融合技术通过融合不同时相、不同来源的遥感数据，将目标的更多信息从不同角度展现出来，可有效提高土地覆盖分类精度[10-12]。已有研究将光学影像与SAR影像融合用于土地覆盖分类[13-15]。SAR与光学影像成像机理不同，简单地选取不同波段或不同成分进行融合，可能造成信息损失和噪声引入等问题，进而影响土地覆盖分类精度[16]。多源SAR影像融合可以较好地规避以上问题，为及时获取多波段、多极化、多角度的SAR数据提供了依据。目前，对于不同源SAR影像的融合，研究者进行了大量研究。针对多波段SAR影像，德国宇航局(Deutsches zentrum für luft- und raumfahrt，DLR)采用模糊专家系统法进行决策级融合[17]，该方法实时性好，并且具有一定的容错能力，但预处理代价较高，影像中的原始信息损失较多。为了在尽可能多保留原始信息的同时降低信息冗余，文献[18]将Gram-Schmidt(G-S)算法用于多时相多波段SAR图像，有效提高了检测性能。针对多极化SAR影像，文献[19]采用非下采样Contourlet变换(Nonsubsampled contourlet transform，NSCT)和脉冲耦合神经网络(Pulse coupled neural network，PCNN)相结合的方法进行了像素级融合，保留了NSCT的多尺度性、多方向性、平移不变性和PCNN全局耦合性，但融合模型中参数设置复杂，不当的参数设置易导致影像融合效果不佳。针对多时相SAR影像，文献[20]利用非负矩阵分解算法(Non-negative matrix factorization，NMF)进行特征级融合，有效保证了矩阵分解结果的可解释性，实现简单、且占用存储空间较小，在SAR影像变化检测应用中效果较好，但非负矩阵分解后，出现数据的稀疏性降低，因此可能导致融合后影像信噪比下降。

递归非负矩阵下近似(Recursive nonnegative matrix underapproximation，RNMU)算法在NMF算法的基础上进行了改进，有效避免了数据稀疏性的降低，成功应用于高光谱影像的特征提取、分类与光谱混合分析[21]。针对不同源多极化多时相微波影像，本文尝试将RNMU算法用于SAR影像融合，并将融合数据用于土地覆盖分类，以有效融合多源信息，并降低数据冗余。

1 研究区概况及数据获取

1.1 研究区概况

图1 研究区及采样点分布

选择吉林省大安市为研究区。大安市位于123°8′45″～124°21′56″E、44°57′～45°45′51″N，属于松嫩平原(图1)，总面积约4 879 km2，海拔为120～160 m，地势较为平坦。大安市属中温带季风气候，年平均降雨量413.7 mm[22]，降雨集中在6—9月，期间高品质、高分辨率的光学遥感数据获取困难。根据大安市政府网提供的2019年土地覆盖统计数据：现有耕地140 293 hm2；未利用地近217 607 hm2；大安市境内水资源丰富，水域面积占全省的1/7；湿地富集，种类繁多。大安市属于农、牧生态系统的交界过渡区域，由于土地的不合理开发利用，近年来出现土地退化、湿地减少等生态问题[23]，因此及时了解大安市土地覆盖现状，可为土地利用规划和生态保护提供基础数据。

1.2 SAR数据源

选取GF-3 SAR数据和Sentinel-1 SAR数据进行多源SAR数据融合与土地覆盖分类。GF-3卫星于2016年8月发射，是我国首颗最高分辨率可达1 m的C频段多极化SAR卫星，GF-3 SAR数据有精细条带2(FSⅡ)等12种成像模式，数据产品包括Level-0至Level-3级标准产品及Level-4级行业应用产品。Sentinel-1是由Sentinel-1A与Sentinel-1B两颗卫星组成的星座，Sentinel-1(以下简称S-1)可提供重复观测的 C 波段 SAR 数据[24]。S-1 SAR数据成像模式主要有干涉宽幅(IW)模式等4种，数据产品包括SAR Level-0 RAW数据、Level-1单视复数(Single look complex，SLC)、Level-1地距(Ground range detected，GRD)和Level-2海洋(Ocean，OCN)等产品。本文所用SAR数据具体参数见表1。

表1 SAR数据具体参数

1.3 样本数据获取

用于分类训练与精度验证的样本数据，由2018年8月实地采集，并结合Google Earth影像人工辨读选取得到。考虑到研究区为含有丰富湿地资源的粮食主产区，本文选取水体、建设用地、耕地、湿地、其他用地(未利用地、盐碱地和草地等)共5类作为待分土地覆盖类型。样本(图1)选取遵循平均分布于整个研究区域的原则，对于易分地类选取较少样本，易混淆地类选用较多样本，选取训练样本数为：水体27个、建设用地41个、耕地72个、湿地69个、其他用地70个。验证样本不同于训练样本，采用水体14个、建设用地21个、耕地46个、湿地35个、其他用地38个作为验证样本。

2 研究方法

2.1 技术流程

针对输入的GF-3双极化SAR数据和S-1双极化SAR数据，在预处理的基础上，进行影像配准、影像融合及土地覆盖分类，处理流程如图2所示。

图2 技术流程图

经预处理得到2颗卫星4种极化方式的后向散射系数图像后，采用互相关[25](Cross-Correlation, CC)算法对图像进行精确配准处理。然后，将RNMU方法引入多源SAR图像融合并用于研究区土地覆盖分类。作为对比，另一种同样具有较小冗余度的融合算法G-S被用于融合和土地覆盖分类。G-S融合方法基于多维正交线性变换, 通过正交化消除影像的冗余信息，并在第一正交分量中加入另一景影像信息，从而实现多源影像融合。

2.2 SAR数据预处理

SAR数据预处理主要包括辐射定标、斑点噪声滤波、地形校正等。

相干斑噪声是SAR遥感图像独有的乘性噪声，会对影像解译造成较大的干扰，斑点噪声滤波是为了降低相干斑噪声。由于GF-3和S-1卫星搭载不同的SAR传感器且成像模式不同，导致不同源SAR影像相干斑噪声不尽相同，因此本文在分别尝试了若干种滤波器的基础上，选择窗口大小为5×5的Gamma滤波降低GF-3影像噪声，在尽可能地保持影像空间细节的情况下提高信噪比[26]；选择Refined Lee滤波降低S-1影像斑点噪声，并将滤波窗口设为7×7[27]。

(1)

式中QV——该景影像量化前的最大值，可通过元数据获取[28]

KdB——定标参数，通过元数据获取

(2)

式中A——定标参数，通过元数据文件获取[29]

2.3 RNMU多源SAR影像融合算法

2.3.1RNMU基本原理

RNMU算法核心思想是：对于任意l×k(l、k∈R+)非负实数的输入矩阵W=[W1W2…Wk]，将W分解为非负的l×r的权重矩阵u和k×r的特征矩阵v，使得uvT≈W，求解u、v的优化过程即为W与uvT间距离的元素平方和最小化迭代求解的过程

uvT≤W

针对上述目标函数求解过程中可能出现的局部最优和病态优化问题，RNMU算法限制了权重矩阵和特征矩阵秩的大小，只保留r=1的矩阵，并增加约束条件uvT≤W。对于影像输入来说，低秩性有效地去除了影像中的噪声。每步仅计算秩为1的矩阵，最优解唯一，有效地避免了NMF算法应用于遥感影像时近似秩改变带来的计算重复、计算效率低下以及融合后影像的品质易受初始化矩阵影响等问题。待融合的雷达影像中，每景影像对应特定的时相和极化方式，增加的约束条件保障了所提取的信息为原始数据的根本特征，使结果的稀疏性和局部特征的可分离性更优[30]。

(3)

式中m、n——输入矩阵的行数和列数

式(3)的求解问题转换为u、vT和Λ迭代优化的过程：

(1)初始化Λ为零矩阵O，初始化u、v为秩为1的向量，且满足uvT≤W。

(2)更新u、vT

(3)用次梯度优化算法迭代更新拉格朗日乘子Λ

Λ=max(O,Λ+α(W-uvT)) (α→0)

(4)

(4)迭代步骤(2)、(3)，直到系统稳定，所得结果即为理想u、vT。

2.3.2基于RNMU的SAR影像融合

将GF-3单个极化影像和已配准的S-1单个极化影像分别作为输入影像a、b，按列展开方式排列为两个列向量Wa、Wb，构造出输入矩阵W

[a11a21…am1…a1n…amn]T=Wa

(5)

b⟹[b11b21…bm1…b1n…bmn]T=Wb

(6)

W=[WaWb]

(7)

v=[v11v21…vm1…v1n…vmn]⟹

(8)

2.3.3融合结果评价指标

针对多源SAR影像融合，本文选用标准差(Standard deviation，SD)、平均绝对误差(Mean absolute error，MAE)、信息熵(Entropy，H)[32]、信噪比(Signal-to-noise ratio，SNR)和峰值信噪比(Peak signal-to-noise ratio，PSNR)[33]等评价指标，分析评价RNMU融合影像的品质，其中，SD反映了融合前后影像的灰度分布情况；MAE反映融合后影像与原始影像之间的差异；H反映融合后影像信息丰富程度，SNR和PSNR反映影像噪声情况。

3 实验与结果分析

3.1 融合结果评价

为定量评价RNMU算法在多源SAR影像融合的效果，针对RNMU融合(R:HHVV, G:HVVH, B:HHVV-HVVH)影像与6月S-1(R:VV, G:VH, B:VV-VH)影像、6月GF-3(R:HH, G:HV, B:HH-HV)影像和G-S融合(R:HHVV, G:HVVH, B:HHVV-HVVH)影像，分别计算SD、MAE、H、SNR和PSNR等指标以比较融合前后和不同融合算法的性能差异。RNMU和G-S融合影像中，R：HHVV波段数据为GF-3 HH极化SAR数据和S-1 VV极化SAR数据分别经RNMU和G-S算法融合得到的双极化SAR产品，G：HVVH波段数据为GF-3 HV极化SAR数据和S-1VH极化SAR数据分别经RNMU和G-S算法融合得到的双极化SAR产品，融合后双极化产品的差值HHVV-HVVH作为B波段。如表2所示，G-S融合和RNMU融合后影像与S-1和GF-3影像相比，SD较大表明融合后影像的灰度级分布较分散，包含的信息量较多；H较大表明两种融合方法都增强了影像的细节表现能力。RNMU融合和G-S融合结果相比，RNMU融合后影像MAE较小，表明了RNMU融合算法具有更好的原始影像信息保留能力；而RNMU融合后的SNR和PSNR大于G-S融合后影像，体现了RNMU融合方法良好的影像噪声抑制能力。GF-3的HH影像和S-1的VV影像经RNMU融合后SD较大，影像具有较大的对比度，影像灰度级分布较为分散，有利于地物分类。而对于两个交叉极化的影像的融合，G-S算法得到的影像SD大于RNMU算法的结果。由于一般交叉极化信息弱于同极化，分类应用中同极化数据应用更为广泛，因此，RNMU融合影像更适于土地覆盖分类。

表2 不同融合算法定量评价

3.2 土地覆盖分类结果与分析

3.2.1典型地物类型散射特性分析

地物的微波后向散射特性是利用SAR影像进行土地覆盖分类的基础[34]。图3给出研究区6月GF-3双极化(HH、HV)影像和S-1双极化(VH、VV)影像以及10月S-1双极化(VH、VV)影像中典型地物的后向散射系数分布图。可以看出，无论是在两时相S-1的VH和VV影像还是在单时相的GF-3 HH和HV影像中，水体均具有最低的后向散射系数，建设用地则具有较高的后向散射系数，而耕地、其他用地和湿地的后向散射系数较为接近，特别是在HH极化的GF-3数据中，水体与建设用地具有和其他的土地覆盖类型良好的可分性。而耕地、其他用地和湿地的后向散射系数较为接近，耕地与湿地在S-1的VH极化下典型地物后向散射系数分离度略优于VV极化，但在GF-3的HH极化后向散射系数分离度优于HV极化，因此，单景影像中后向散射系数较为接近的耕地、湿地和其他用地需要利用SAR影像不同极化、不同时相的变化特性进行进一步区分。

图3 不同时相不同极化影像典型地物后向散射系数分布图

图4 典型地物后向散射系数的均值和标准差的比较

为了更直观看出不同地物的散射特性差异，图4给出了融合前GF-3 SAR影像和S-1 SAR影像与RNMU融合后SAR图像的均值和标准差。其中基于RNMU的融合，分别给出了研究区10月S-1 VV极化影像和6月GF-3 HH极化影像经RNMU算法融合得到的VVHH影像，10月S-1 VH极化影像和6月GF-3 HV极化影像经RNMU算法融合得到的VHHV影像(VVHH影像和VHHV影像叠加，简写为RNMU1)。6月S-1 VV极化影像、10月S-1 VV极化影像和6月GF-3 HH极化影像经RNMU算法融合得到VVVVHH影像，6月S-1 VH极化影像、10月S-1 VH极化影像和6月GF-3 HV极化影像经RNMU算法融合得到VHVHHV影像(VVVVHH影像和VHVHHV影像叠加，简写为RNMU2)。由图4可以看出，6月的S-1影像中湿地和耕地上的植被长势较好，因此比其他用地具有更高的后向散射系数，10月研究区的作物成熟或收获，耕地的后向散射系数降低明显，不同地物随时间的变化特性差异相对明显。而经RNMU融合后的影像，无论是RNMU1影像还是RNMU2影像各地类均值之间差异增大，因此增加了地类之间的可分性；而每种地类的标准差较S-1影像降低，进而提升了待融合影像(S-1)地类特征的一致性。因此利用不同时相、不同极化的多源SAR影像进行融合，能融合更多互补信息，得到更好的分类结果。

3.2.2分类结果分析

为分析融合影像在土地覆盖分类应用中的性能改进，分别对研究区 RNMU1影像，RNMU2影像，6月、10月S-1影像和6月GF-3影像经G-S融合后影像(简写为G-S)，应用最大似然分类器(Maximum likelihood classifier，MLC)进行土地覆盖分类，作为比较，同时给出6月单景S-1影像，6月单景GF-3影像，6月、10月S-1影像和6月GF-3经简单波段叠加后影像的土地覆盖分类结果。由图5可以看出，利用G-S和RNMU2影像识别的水体基本相同，建设用地、湿地和其他用地存在不同程度差别。选取土地覆盖类型多样的区域(图5矩形框内区域)作为验证区域，利用3 m分辨率的Google高清图像与S-1影像、G-S影像和RNMU2影像进行分类结果实验验证。如图6所示，与S-1影像分类结果相比，G-S影像和RNMU2影像分出的地类边界较清晰。RNMU2影像分出的建设用地(红色椭圆框圈出)更接近真实值，G-S影像分类结果次之，S-1影像分出的建设用地较为破碎。

图5 不同数据土地覆盖分类图

图6 验证区域的土地覆盖分类比较

利用验证样本对分类结果进行精度评价，各数据的分类总体精度(OA)、制图精度(PA)、用户精度(UA)和Kappa系数(Kc)如表3所示。由表3可知，GF-3数据的土地覆盖分类整体精度优于同时相S-1数据，不同时相不同源(S-1和GF-3)SAR数据经RNMU融合后数据分类OA分别由单时相SAR数据74.96%(S-1结果)和83.52%(GF-3结果)提高至86.08%(RMNU1影像结果)；Kc分别由0.55和0.71提高至0.74；耕地、建设用地和湿地的分类精度都有明显提高。由于研究区湿地种类繁多，同时间不同类别的湿地变化情况不一。虽然加入其他时相数据可以提高其识别精度，但其分类精度与其他地物相比仍较低。使用相同数据源的情况下，G-S融合后分类结果要低于所有SAR数据直接叠加或者RNMU融合后的结果，可能是因为G-S融合图像为了降低数据冗余，牺牲了一定的纹理信息，导致纹理特征明显的耕地UA较低，而耕地在研究区面积占比较大，导致OA和Kc整体下降。基于RNMU融合后影像(RNMU2)地类识别效果，整体优于基于G-S算法融合结果和简单波段叠加后影像，说明RNMU融合算法充分利用多时相、多极化SAR数据的互补信息，有效提高土地覆盖分类精度，且融合后影像数据量减少了1/3，实现了数据压缩，适用于较大区域的SAR数据。但对于易于分辨的水体地物，单时相S-1 SAR数据的水体分类PA高于RNMU2的结果，而融合后并没有达到较好的效果，说明对于后向散射特征足够明显的地物，融合后数据可能出现过拟合现象，反而可能引入噪声，降低分类精度。

表3 不同数据融合方式分类精度评价

注：表头第2行表示数据的月份、来源以及分类所用的极化影像。

4 结论

(1)为充分利用多时相、多极化SAR互补数据，将RNMU算法引入多源SAR融合。RNMU融合后，影像SD、H均优于原始S-1、GF-3双极化SAR影像；与G-S融合算法相比，MAE、H、SNR和PSNR均体现出优良的特性。说明RNMU融合算法有效降低了SAR数据噪声，提高了融合后影像的信噪比，表明RNMU融合算法在降低数据冗余的同时有效融合了多源数据之间更多的互补信息。

(2)RNMU融合后，SAR影像综合利用多时相、多极化SAR信息提高了土地覆盖分类精度。融合两时相SAR数据进行土地覆盖分类的精度明显高于融合前SAR数据，表明RNMU融合后多源SAR数据更适于土地覆盖分类应用。

(3)对于后向散射特征明显的水体地物，RNMU融合数据分类的PA较单时相S-1略有下降，这说明对于明显特征地类可能存在过拟合现象。在进一步融合研究中，将不同区域、不同地类设定不同的权重，以期在增加互补信息的同时，降低噪声干扰。