基于正则化IR-MAD的GF-1影像辐射归一化

黄莉婷，焦伟利，龙腾飞，康传利

(1.桂林理工大学 测绘地理信息学院，广西 桂林 541006；2.中国科学院空天信息创新研究院，北京 100094)

0 引言

高分一号(GF-1)卫星开启了我国对地观测的新时代，该卫星具备高空间分辨率、高时间分辨率、宽覆盖和多源遥感数据的综合数据获取能力，能满足各行业的应用需求。同时，国内外各类遥感卫星的种类日益增多，使得遥感数据资源呈现多源性，降低了单颗卫星在进行长时间序列研究中的困难。Landsat-8陆地成像仪(operational land imager，OLI)和Sentinel-2多光谱仪(multi spectral instrument，MSI)数据已成为多源影像来源的研究热点，其数据的免费获取极大地推进了中分辨率陆地表面成像的虚拟星座计划。可见，GF-1、Landsat-8 OLI和Sentinel-2 MSI这3颗卫星数据在当前以及今后相当长时间内，是实现卫星长时间序列数据集构建的重要互补性数据源。

地表反射率(bottom-of-atmosphere corrected reflectance，BOA)是遥感定量反演的一个关键参数。目前可公开获取的Landsat-8 OLI Level 2级BOA产品和Sentinel-2 L2A级BOA产品都具有较高的精度，但是GF-1地表反射率产品精度问题依然没有得到足够重视。对GF-1影像大气校正的研究，目前主要是用基于大气辐射传输模型的绝对辐射校正方法。该方法对输入参数精度要求严格，并且精确的大气参数难以获取。而采用基于像元值经验模型的相对辐射校正方法来进行研究的还比较少，在现有的辐射归一化研究中，也多集中应用在多时序遥感影像的变化检测上。如果缺少相关大气参数(气溶胶、臭氧浓度等)，又要得到较为精确的BOA，那么相对辐射校正是一种较好的方法。此外，Landsat-8和Sentinel-2地表反射率产品为辐射校正提供了较好的参考数据。

相对辐射归一化方法假设同一区域的多时相或同时相影像相同波段的灰度值存在线性关系[1]，目前常用的线性相对辐射归一化方法包括暗集-亮集法(dark set-bright set normalization，DB)[2]、伪不变特征法(pseudo-invariant feature，PIF)[3]、自动散点控制回归算法(automatic scattergram-controlled regression，ASCR)[4]、多元变化检测法(multivariate alteration detection，MAD)[5]和迭代加权多元变化检测法(iteratively reweighted MAD，IR-MAD)[6]。不少学者基于这些方法进行了相关的辐射校正研究[7-13]。在这些方法的基础上，都采取了一些方法对2幅影像上没有发生变化且光谱特性稳定的同一地物点(即PIF点)进行高质量的提取，实现了多时相遥感影像的辐射归一化。但是这些研究大多数都是针对Landsat数据的，随着越来越多的不同传感器数据的综合应用，现有的相对辐射归一化方法均存在一定的限制，不能直接应用于GF-1数据[14]。如DB方法要求有足够的暗点和亮点，这在有些情况下满足不了。PIF方法选取的PIF点代表影像的大部分区域时，才能得到有效的辐射校正结果。通过水陆集群中心确定未变化集(no changed set，NC)可以排除掉一些辐射特性发生较大改变的像元，但是NC的确定需依据经验确定，当研究区水体较少时，不能确定NC。虽然IR-MAD方法能够有效地集中变化信息，减少背景信息对变化信息的影响，对影像的要求比较少，但仍然存在多个阈值需要依据经验确定的缺点，且在多次迭代的过程中可能会出现近奇异矩阵。

为使GF-1影像具有较高精度的BOA，降低大气校正的复杂度，本文提出一种简便实用的方法，对大气校正后的GF-1影像，以Landsat-8和Sentinel-2的BOA影像为基准，通过相对辐射归一化提高GF-1 BOA精度，同时也可减少GF-1和Landsat-8、Sentinel-2之间的辐射差异，使得多源影像数据可以作为互补性数据联合使用。

1 方法

本文提出一种基于正则化IR-MAD的GF-1影像辐射归一化方法。该方法不要求影像同时具备暗点和亮点特征，并能消除大样本集中样本噪音对辐射校正结果的影响。其主要技术流程为：首先根据目标影像及基准影像的红波段和近红外波段散点图确定NC；然后运用正则化IR-MAD规则对NC进一步筛选，得到不变特征点；最后由不变特征点建立线性回归模型求解系数，得到辐射归一化方程，进行辐射归一化(图1)。

图1 本文方法流程图

该辐射归一化方法主要包含确定NC和提取不变特征点2个关键步骤，下面分别对其进行介绍。

1.1 确定未变化集

与ASCR方法不同的是，本文方法不用散点图上的水陆集群中心确定NC，而是由正交回归法确定初始回归基线，再求解NC。首先，由目标影像及基准影像的红波段和近红外波段分别生成散点图，散点图横坐标是目标影像，纵坐标是基准影像，如图2所示。由线性变换确定的NC能初步表示“伪变化”的像元，即辐射特性几乎没有改变的像元。

图2 从散点图确定NC

图中HVW是NC选择的限制条件[4]，由HPW控制，HPW是不变区域一侧的垂直宽度。HVW和HPW的关系如式(1)所示。

(1)

式中：a是由正交回归法确定的初始回归线斜率。分别确定了红波段和近红外波段散点图的初始回归线及HVW后，由公式(2)确定最终的NC[4]。

NC=(X，Y)：|y1-b1-a1x1|≤
HVWNC1&|y2-b2-a2x2|≤HVWNC2

(2)

式中：y1、y2分别是基准影像红波段和近红外波段的BOA值；x1、x2分别是目标影像红波段和近红外波段的BOA值。

1.2 提取不变特征点

在NC中，对不变特征点的选择做进一步的优化。采用以下规则来寻找可以使用的拟合样本点数据[6]：对影像的每个像元设置初始权值为1，每一次迭代过程中计算均值向量与方差矩阵，用典型相关分析的方法计算每个MAD变量，根据新计算的MAD变量更新权值，值域区间为[0，1]；未发生变化的像元具有较大的权值，经过若干次迭代计算之后，每个像元的权值会趋于稳定，通过权值与阈值的比较，便可判定每个像元点是否属于不变特征点。更新权值的公式如式(3)～式(5)所示。

(3)

(4)

Pr(nochange)=1-Pχ2；N(Z)

(5)

式中：γi是第i波段的方差；Z是迭代赋值引入的随机变量，表示标准化MAD变量平方和，其服从自由度为N(波段数)的卡方检验，通过卡方分布概率密度函数对Z进行加权；Pr(nochange)表示迭代权值，可以设定一个固定的阈值t(如t=0.95，便能得到很好的结果)，当Pr(nochange)>t时，判断为不变特征点。

在多次迭代的过程中可能会出现近奇异矩阵，当矩阵中某一元素发生很小的变动就会产生很大的误差。所以，Nielsen[15]在迭代的过程中引入正则化(也称为“惩罚项”)，避免在计算中产生较大权重，影响不变特征点的判断。本文在求解协方差矩阵时，也加入了正则化参数，如式(6)～式(7)所示。

Var11=∑11+λ1Ω1

(6)

Var22=∑22+λ2Ω2

(7)

式中：λi是规则化参数；Ωi是N×N的对角线矩阵。

2 数据

2.1 实验数据

Sentinel-2A MSI其L2A级产品为经过辐射定标、大气校正的BOA反射率影像，可由欧空局(European Space Agency，ESA)发布的专门生产L2A级数据的插件Sen2cor根据需求自行生产[16]，或者在ESA官网(https：//scihub.copernicus.eu/dhus/#/home)进行下载，目前可下载到2018年12月后的中国区域L2A级数据。

Landsat-8 OLI的Level 2级BOA产品，由美国地质调查局(United States Geological Survey，USGS)专门为Landsat-8数据设计的大气校正程序LaSRC(landsat-8 surface reflectance code)生成。该程序利用海岸气溶胶波段(band1)进行气溶胶反演测试，运用MODIS辅助气候数据并使用专门的辐射传输模型，被认为是最精确的Landsat-8 OLI大气校正程序，其结果精度被广泛认可和使用[17]。Level 2级数据可在USGS的网站(https：//earthexplorer.usgs.gov/)下载。

为保证传感器在获取遥感数据时受大气影响的差异较小，应该尽量选取同一日期的太阳高度角和方位角接近的研究区，进行多源影像的相对辐射归一化[18]。但是，由于不同卫星的重访周期和幅宽不相同，因此难以获得同一地区的同一日期、同时过境的不同卫星影像。本研究的影像选取遵循影像前后日期在7 d左右的规则。

本文分别选取分辨率相近的GF-1 PMS2和Sentinel-2A MSI，以及GF-1 WFV4和Landsat-8 OLI影像进行2组实验。第1组选择了2017年12月19日的GF-1 PMS2影像和2017年12月11日Sentinel-2A MSI影像，土地覆被类型以丘陵、水体为主；第2组选择了2018年6月14日的GF-1 WFV4影像和2018年6月11日Landsat-8 OLI影像，土地覆被类型以耕地、裸土和不透水面(建筑物、道路等)为主。

3颗卫星影像的对应波段信息对比如表1及图3所示。由表1及图3可知，对于GF-1 PMS2和Sentinel-2A MSI而言，除了蓝波段，GF-1 PMS2的波谱带宽明显比Sentinel-2A宽，但中心波长差距较小；GF-1 WFV4和Landsat-8 OLI 2类传感器在可见光波段波谱带宽设置较相近，在近红外波段波谱带宽设置差异显著，GF-1在近红外波段处有较宽的波谱通道，并且在近红外波段二者中心波长相距较大。

表1 所选影像波段信息

图3 传感器光谱响应曲线

2.2 数据处理

获取BOA需要对GF-1的PMS2和WFV4影像进行绝对辐射校正，在ENVI5.3中进行。首先把影像灰度值(digital number，DN)值转换为辐射亮度，再将辐射亮度转换成BOA。在本文实验中具体处理步骤包括：辐射定标、大气校正、正射校正和影像重采样。Landsat-8 OLI和Sentinel-2A MSI影像经过波段组合后，需要和处理过的GF-1 BOA影像进行精确的几何配准，配准前GF-1需要重采样，配准精度要求达到均方根误差小于0.5个像元，影像配准的质量会影响辐射归一化的效果。

3 实验与分析

3.1 不变特征点提取

相对辐射归一化主要由不变特征点提取、求解回归系数和辐射变换3个步骤组成，在本文中选取影像大小为1 000像素×1 000像素进行实验。为减弱大气辐射的影响，首先将目标影像和基准影像的近红外和红波段分别生成散点密度图，确定初始回归线；再由初始回归线及式(1)、式(2)确定每个波段的NC。初始回归线及NC如图4(c2)、图4(d2)、图5(c2)和图5(d2)所示。接着运用正则化IR-MAD规则对NC做进一步筛选。NC中每个像元根据自由度N为4(波段数)以及影像差值的卡方变量Z，通过式(5)得到每个像元的权值。再将权值与阈值t进行比较，判断该像元是否属于变化和未变化像元，一般阈值t需要根据经验多次求取。为分析阈值t对求解辐射归一化方程的影响，本文在t∈[0.91，0.95]等间隔取5个阈值：0.91、0.92、0.93、0.94、0.95，分别计算各波段的回归判定系数R2(式(8))和均方根误差RMSE(式(9))。

图4 GF-1 PMS2和Sentinel-2A MSI 不变特征点的3次选择过程

(8)

(9)

R2和RMSE计算结果如表2所示。分析表2可知，2组实验中各波段在不同阈值t的情况下R2和RMSE均相差不大，但RMSE值会随着阈值t的升高略有降低。分析其原因，是因为随着阈值t的升高，样本数量逐渐降低，导致RMSE值也逐渐减小。再分别对5个阈值辐射归一化后的影像以相同的样本数量计算R2和RMSE(表3)。由表3可知，2组实验中各波段的RMSE的差距进一步减小，由于使用了相同的样本集，所以各波段的R2相一致。该结果说明，所有阈值都能减小目标影像与基准影像之间的辐射差异，也进一步说明了该方法的有效性和稳定性。

表3 相同样本集不同阈值t的R2和RMSE

由于阈值t的稳定性，2组实验均取相同的阈值t=0.95进行下一步分析。

图4和图5分别显示了2组实验影像各波段不变特征点的3次选择过程，图中绿色直线代表初始回归线，红色直线为最终求解出来的拟合线。各实验影像的BOA反射率值均被扩大了10 000倍。

图5 GF-1 WFV4和Landsat-8 OLI 不变特征点的3次选择过程

第1次图(a1)、图(b1)、图(c1)、图(d1)是目标影像和基准影像对应波段生成的散点图；第2次图(a2)、图(b2)、图(c2)、图(d2)显示的是由式(1)和式(2)确定的NC；第3次图(a3)、图(b3)、图(c3)、图(d3)是由正则化的迭代多元变化检测规则提取出来的不变特征点。这些散点图有一些共同的特征：不变特征点逐渐收敛，最后一次散点图中的点数量减少，沿直线分布；主轴的斜率从第2次到第3次有一定的变化；散点的变化范围逐渐变小，且各点分布较均匀。这些特征表明了该方法在不变特征点选择中的有效性和稳定性。

传感器波段t0.91t0.92t0.93t0.94t0.95R2RMSER2RMSER2RMSER2RMSER2RMSE蓝 0.9807.9160.9807.6480.9827.3380.9827.3490.9846.839PMS2和MSI绿 0.98610.7840.98810.3420.98810.1350.9889.3680.9908.736红 0.97610.9940.97810.6120.97810.7460.97810.0840.9809.832近红外0.91525.8500.93524.4750.99424.8560.92724.5230.92723.759蓝 0.96614.3880.98314.8600.97013.6920.96814.0210.97014.071WFV4和OLI绿 0.97417.5430.98717.2460.97616.7540.97415.8920.97615.743红 0.98422.7090.99122.6460.98421.5470.98420.9390.98419.911近红外0.96042.2900.97842.6770.96641.2950.96640.7430.96438.917

3.2 求解辐射归一化方程和辐射变换

将上述选出的不变特征点集运用正交回归法分别计算得到PMS2和MSI以及WFV4和OLI在蓝、绿、红、近红外波段的反射率归一化方程(表4)。最后把各波段辐射归一化方程应用在目标影像GF-1 PMS2和GF-1 WFV4相应波段中，得到辐射变换后的影像。为了提供能够衡量总体差异性的指标，1/3不变特征点集作为验证数据用来进行光谱辐射归一化结果的验证。根据式(8)计算回归判定系数R2，式(10)计算相关系数r，式(11)计算平均值差异Δmean，计算结果如表4所示。

表4 辐射归一化方程与R2、r、Δmean

(10)

(11)

在2组实验中，各波段的R2均在0.92以上，说明所建立的回归模型拟合效果较好。Δmean都较小，说明辐射变换后的GF影像和基准影像差异程度较小。

结合图4、图5(a1)、图5(b1)、图5(c1)、图5(d1)分析，GF-1 PMS2和Sentinel-2A MSI以及GF-1 WFV4和Landsat-8 OLI的BOA值在可见光和近红外波段都不同程度地存在一定偏差。对比表4中各传感器对应波段的光谱值的相关系数r，均大于0.96，这表明对应波段反射率都具有很强的线性相关性，基准影像能够显著提升GF-1卫星影像的辐射精度。

3.3 结果分析

1)目视效果分析。采用目视的方法对辐射归一化后的影像进行定性评价，将归一化后的影像与基准影像进行显示对比，如果2个影像颜色、亮度等十分相似，则说明归一化的效果较好[19]。图6(a)和图6(b)是GF-1 PMS2和Sentinel-2A MSI辐射归一化前后结果的对比；图6(c)和图6(d)是GF-1 WFV4和Landsat-8 OLI辐射归一化前后结果的对比，影像均在ENVI中选择相同的拉伸方式显示。观察图6(a)和图6(b)，在PMS2(右)和MSI(左)影像接边处，水体、植被、耕地等土地覆被类型的边界镶嵌处差异明显减小，2个影像各土地覆被类型的颜色、亮度等十分相似。在图6(c)和图6(d)中，WFV4(右)和OLI(左)影像接边处裸地、植被、耕地等土地覆被类型的边界镶嵌处差异较小，但在城市建筑用地区域，2个影像的接边处还存在着一定亮度的差异。整体来看，在2组实验中，目标影像与基准影像边界镶嵌处差异减小，可见数据辐射一致性明显改善。经过辐射归一化后的影像，BOA值的分布会有所改变，从辐射归一化前后4个波段的核密度估计(kernel density estimation，KDE)密度分布图(图7、图8)可以看出，辐射归一化后的影像BOA数值分布，更接近于基准影像BOA数值的分布特征。其中，图7(d)有2个峰谷，分别代表的是影像上水、陆集群中心。图8(d)中由于影像上没有水像元，所以只有一个峰谷。本文方法不仅适用于有水陆分布的影像，对于缺少水体分布的影像也能很好地进行辐射归一化。

图6 辐射归一化前后对比

图7 PMS2影像、MSI影像和归一化后影像BOA密度分布比较

图8 WFV4影像、OLI影像和归一化后影像BOA密度分布比较

2)典型地物波谱曲线分析。因受辐射畸变和不同传感器的特性影响，导致地物波谱曲线有很大的差异，相对辐射归一化的目的就是消除或减弱这些差异[20-21]。地物的波谱曲线越相似，说明地物的反射特性越相近，地物波谱曲线的变化在一定程度上可以说明相对辐射归一化的效果[20]。本文分别从目标影像、辐射归一化后影像与基准影像中选取几种典型地物进行波谱曲线分析。从图9和图10中可以看出，经辐射归一化后的影像，其典型地物波谱曲线向基准影像拉近，说明辐射归一化是有效果的。结合图9、图10及表4分析，在可见光波段(蓝、绿、红)对地物光谱特征的辐射归一化效果相对较好，而近红外波段辐射归一化效果不如可见光波段。分析其原因，主要是由于PMS2和WFV4相较于MSI和OLI传感器，在近红外波段的光谱响应函数差异明显大于其他3个可见光波段的差异，尤其是光谱范围在770～900 nm之间的光谱响应曲线(图3)。

3)与其他方法的比较。将本文方法和原ASCR方法、IR-MAD方法的归一化效果进行比较，以均方根误差RMSE作为归一化影像与基准影像的相似性度量指标，RMSE的值越小，说明辐射归一化效果越好。

以图4和图5最后一次筛选出来的不变特征点作为样本集，分别计算目标影像和基准影像的RMSE，以及ASCR方法、IR-MAD方法和本文方法辐射归一化前后影像的RMSE，计算结果如表5所示。由表5可知，经辐射归一化后的影像RMSE均小于未辐射校正前的，说明这3种方法都减小了目标影像和基准影像的辐射值差异，3种方法同样在近红外波段辐射归一化效果不及可见光波段，其中本文方法与ASCR方法、IR-MAD方法相比，整体精度都较高。

表5 不同方法辐射归一化的RMSE

ASCR方法采用整景影像比较多的像元信息进行辐射归一化，因此易受到离初始回归线较远的像元影响。那些离初始回归线较远的像元不仅线性作用较差，反射率还发生了变化，会给辐射归一化效果带来一定的影响。IR-MAD方法基于辐射特征的变化进行不变特征点的选取，其MAD变量中的像元值反映了变化信息的情况，对未变化像元的选取更加准确。而在本文方法中，用初始回归线先确定了NC，保证未变化像元的线性关系，再对NC采用IR-MAD的规则进一步筛选，再次控制了不变特征点的选取质量，所建立的辐射归一化方程的转换效果较好。

4 结束语

在GF-1大气校正后，本文基于IR-MAD，采用控制NC的线性关系、逐步筛选不变特征点的方法进行GF-1影像与其他影像的光谱辐射归一化。实验表明，GF-1 PMS2和Sentinel-2A MSI以及GF-1 WFV4和Landsat-8 OLI在BOA上有较好的线性关系，不变特征点的相关系数达到了0.96以上，说明通过辐射归一化提高GF-1影像的辐射精度是可行的；与其他辐射归一化方法相比，本文方法具有较好的精度。分别从目视效果、典型地物波谱曲线、均方根误差等方面分析评价了辐射归一化结果，取得了较为理想的效果，为提高GF-1影像BOA的辐射精度提供了一种解决思路。

本文主要针对大气校正后的GF-1影像进行研究。由于地表具有二向反射分布特性，反射率随太阳入射角及观测角的变化而不同，导致不同时期和不同传感器获取的影像其成像几何角度有所差异，所以应当进行BRDF校正以减小由成像几何角度差异引起的反射率差异[22-23]，这在今后的研究中需要进一步分析。