遥感卫星成像时，因大气状况、光照条件、传感器差异等因素的影响，同一地区不同影像之间可能存在明显的辐射差异，这种差异将直接影响影像的下一步应用，如镶嵌制作正射影像图^[1]、目视判读提取信息^[2]等。因此，不同时相影像镶嵌前都需要进行匀色处理，使其具有色彩一致性^[3]。

在多源卫星影像匀色方面，已有学者建立了多种算法，包括直方图匹配法、色彩空间变换法、MASK匀色法(马斯克匀光法)^[4]、小波变换法^[5]、辐射校正法^[6]，PhotoShop人工调色等^[7]，但这些方法或者操作复杂对人员要求高、或者计算复杂耗费时间长，都难以用于自动化、工程化的处理。本文基于多元变化检测方法，自动提取影像之间未发生变化的像元点，然后利用这些不变像元点建立线性回归模型，最后将模型应用于影像进行匀色处理。

1 数据与预处理

本文选择国内应用较为广泛的国产高分系列卫星数据开展匀色方法研究，数据包括高分一号(GF-1)和高分二号(GF-2)的PMS(全色多光谱相机)数据，数据覆盖范围主要为重庆市南川区，如图 1所示。

1) 高分一号卫星是我国首颗民用高分辨率对地观测卫星。GF-1传感器包括两台分辨率为2 m的全色(PAN)分辨率为8 m多光谱相机(PMS)和4台分辨率为16 m的宽幅多光谱相机(WFV)^[8]。本文所用高分一号数据为两景PMS数据，获取时间和覆盖范围如图 1所示。

2) 高分二号是我国自主研制的首颗空间分辨率优于1 m的民用光学遥感卫星。高分二号搭载了1 m(星下点为0.8 m)分辨率全色(PAN)和分辨率为4m多光谱相机(PMS)，具有亚m级空间分辨率、高定位精度和快速姿态机动能力等特点。本文所用高分二号数据为3景PMS数据，获取时间和覆盖范围如图 1所示。

3) 数据预处理。将研究区数据分成GF-1影像对、GF-2影像对、GF-1与GF-2混合影像对3组进行匀色，并分别进行正射校正和配准处理，每组影像配准时以基准影像为参考。1~3组基准影像依次为GF-1_20160726、GF-2_20160228和GF-2_20161002，匀色影像依次为GF-1_20160505、GF-2_20161205和GF-1_20160726。第3组数据中，由于GF-1影像和GF-2影像分辨率不一致，在配准时将GF-1影像分辨率采样为4 m。

2 研究方法

遥感影像匀色是为了消除不同时相影像之间，因成像条件和传感器差异等因素导致的辐射差异。其目标是使不同影像具有相同的辐射尺度，即在不同时相的多源遥感图像上，相同地物具有相同或相似的颜色。目前，常用的遥感影像匀色方法主要分为两大类：一是人工目视调整，如PhotoShop调色；二是利用两景影像的全部或部分像元值，然后以一景影像为参考，统计并建立算法模型对另一景影像进行匀色处理，如直方图匹配法。首先本文利用遥感影像多元变化检测方法，自动提取多源影像对之间未变化的像元点；然后，利用这些未变化像元点建立正交回归模型，对影像进行匀色处理。

2.1 多元变化检测

多元变化检测(multivariate alteration detection，MAD)由Nielsen等^[9]提出的，主要用于两景影像之间的变化检测^[10]。该算法基于典型相关分析，具有线性不变的特点。对于相同地区不同时相的两景影像，设随机向量X =[X₁…X_N]，Y =[Y₁…Y_N]分别对应两影像的像元值，N为波段数。对X、Y分别做线性变换U= a^TX，V= b^TY，其中a =[a₁…a_N]^T，b =[b₁…b_N]^T，则两景影像之间的变化信息可以集中到变化影像D=U－V中。不同的系数a和b得到不同的变化影像D，D的值表示变化的大小。要最大限度地突出变化，则D的方差应最大。如果对a和b乘上一个常数c，则方差会变c²倍，因此对a和b作出限制，令U和V具有单位方差，即

$ {\rm{var}}\{{\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{X}}\} = {\rm{var}}\{\mathit{\boldsymbol{b}}{^{\rm{T}}}\mathit{\boldsymbol{Y}}\} = 1 $

(1)

则MAD算法可描述求满足式(1)的a和b，使var{D}=var{a^TX － b^TY}达到最大。MAD算法假设a^TX和b^TY正相关，则D= a^TX － b^TY可表示变化信息，具体为:

$ \begin{array}{l} \;\;\;\;\;\;{\rm{var}}\left\{D \right\} = {\rm{var}}\{\mathit{\boldsymbol{a}}{^{\rm{T}}}\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{b}}{^{\rm{T}}}\mathit{\boldsymbol{Y}}\} = \\ {\rm{var}}\{\mathit{\boldsymbol{a}}{^{\rm{T}}}\mathit{\boldsymbol{X}}\} - {\rm{var}}\{{\rm{}}\mathit{\boldsymbol{b}}{^{\rm{T}}}\mathit{\boldsymbol{Y}}\} - 2{\rm{cov}}\{{\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{X}}, {\mathit{\boldsymbol{b}}^{\rm{T}}}\mathit{\boldsymbol{Y}}\} = \\ \;\;\;\;\;2(1 - {\rm{corr}}\{{\mathit{\boldsymbol{a}}^{\rm{T}}}\mathit{\boldsymbol{X}}, {\mathit{\boldsymbol{b}}^{\rm{T}}}\mathit{\boldsymbol{Y}}\} ) \end{array} $

(2)

当线性变换系数a、b使a^TX、b^TY的相关系数最小时，D的方差最大。通过典型相关分析理论(canonical correlation analysis, CCA)^[11]，可以求解满足条件的极值a、b。

由式(2)可知变化图像的方差和相关系数成反比，因此定义MAD变量为:

$ {\rm{MA}}{{\rm{D}}_i} = {U_{N + 1 - i}} - {V_{N + 1 - i}} = {\mathit{\boldsymbol{a}}_{N + 1 - i}}^{\rm{T}}\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{b}}_{N + 1 - i}}^{\rm{T}}\mathit{\boldsymbol{Y}} $

(3)

式中, i=1, …, N，则MAD变量的方差为:

$ {\rm{var}}({\rm{MA}}{{\rm{D}}_i}) = {\sigma _{{M_i}}}^2 = 2(1 - {\rho _{N + 1 - i}}) $

(4)

式中, ρ₁²≥ρ₂²≥…≥ρ_N²。根据典型变量的性质，MAD变量之间相互正交且var(MAD₁)>var(MAD₂)>…>var(MAD_N)，这表示不同的MAD变量代表不同的变化信息，且各MAD包含的变化信息量依次减少。

对于任一像元点，其所有的变化信息为$\sum\limits_{i = 1}^N {{\rm{MA}}{{\rm{D}}_i}} $，但$\sum\limits_{i = 1}^N {{\rm{MA}}{{\rm{D}}_i}} $的值的范围并不确定，因此无法利用一个合适的阈值来判断像元是否发生变化。令$Z = \sum\limits_{i = 1}^N {{{\left({\frac{{{\rm{MA}}{{\rm{D}}_i}}}{{{\sigma _{{M_i}}}}}} \right)}^2}} $将MAD变量标准化，则Z包含了所有的变化信息。对于未变化的像元，其MAD变量的值应近似为0，由U、V的相关性可知MAD变量互不相关，则Z近似服从N维自由度的卡方分布，因此，可用卡方分布检验像元是否发生变化:

$ {\rm{Pr}} = 1 - {p_{{\chi ^{2, N}}}}\left( Z \right) $

(5)

式中, Pr为未变化的置信度(本文取t=0.95);p_{χ^{2, N}}(Z)为Z的N维自由度卡方分布函数。

为了最大限度的突出变化像元与未变化像元的区别，Nielsen对MAD算法进一步改进，提出一种加权迭代的方法(iteratively re-weighted modification of the MAD transformation, IR-MAD)^[12]。相对于MAD算法，IR-MAD变化检测结果更准确，尤其对于变化较小的像元。其方法是对MAD变量进行多次加权迭代：对任一像元j，第k次迭代时，权重w_k=Pr_k－1, 其中w₁=1。权重w参与计算均值和协方差矩阵，如对于向量X，有：

$ \left\{\begin{array}{l} {\mathit{\overline X} _i} = \frac{{\sum\limits_{j = 1}^C {{w_j}{X_{ji}}}}}{{\sum\limits_{j = 1}^C {{w_j}}}}\\ {\rm{var}}\left( X \right) = \frac{{\sum\limits_{j = 1}^C {{{\left( {{X_{ji}} - \mathit{\overline X}} \right)}^2}}}}{{\left( {C - 1} \right)\sum\limits_{{\rm{}}j = 1}^C {{w_j}/C}}} \end{array} \right. $

(6)

式中, ${\overline X _i}$为波段i的均值; C为像元数。迭代终止条件为当前后两次相关系数的最大值之差的绝对值小于0.001。

2.2 回归分析

利用未变化的像元点建立线性回归模型，以一景影像为基准对另一景影像进行匀色处理。线性回归模型中，应用最广泛的是最小二乘法，它假设自变量x是精确的，而因变量y由于测量误差等原因在不确定性。本文中自变量x与因变量y为遥感影像像元值，均可能有误差，利用最小二乘法回归分析，由于只考虑了一个方向的扰动，拟合效果稳定性较差。而基于几何距离的正交回归法，则可以克服单方向最优带来的回归误差。

假设在二维平面中有n个数据点P_i(x_i, y_i), i=1, …, n，回归分析是通过某种原则寻找一条直线L来描述这些点的关系：

$ y = \alpha + \beta x $

(7)

其中, 最小二乘法是使$D = \sum\limits_{i = 1}^n {{{\left({{y_i} - \beta {x_i} - \alpha} \right)}^2}} $最小，而正交回归则是使$D = \sum\limits_{i = 1}^n {d{{\left({{P_i}, L} \right)}^2}} $最小，其中d(p_i, L)是点P_i到L的几何距离^[13]。可见，最小二乘法只考虑了Y方向的误差，而正交回归则考虑了两个方向的误差。

在相对辐射归一化中，由于匀色影像和基准影像的像元值都可能有误差，传统的最小二乘法的假设不成立。如有一组匀色影像和基准影像，线性回归所得方程为y=α+βx，将匀色影像和基准影像互换，对于正交回归，此时的斜率为1/β，而最小二乘回归的斜率不一定等于1/β，这与实际情况不符，所以正交回归的结果通常要优于最小二乘回归，Canty等^[14]的研究结果也证明了这一结论。因此本文采用线性正交回归进行匀色处理。

3 结果

在IDL+ENVI环境中编写代码模块，将上述多元变化检测和回归分析方法实现自动化处理。然后对表 1所示的3组数据进行自动匀色处理，分别检验本文方法对于GF-1影像对、GF-2影像对、GF-1 GF-2混合影像对匀色的适用性。

对试验区3组数据分别进行匀色处理后，首先采用目视效果对比的方法检验匀色效果。目视效果对比采用两景影像交叉镶嵌的方法，即交叉选取基准影像和匀色前影像进行镶嵌，交叉选取基准影像和匀色后影像进行镶嵌，并取波段321合成彩色图像，3组数据镶嵌后视觉效果如图 2所示。

由图 2可知，匀色处理前，3组影像对的交叉镶嵌结果中拼接处有明显的色差，图像整体呈现条纹现象(左图)，而经过匀色处理后的镶嵌结果中，拼接处的色差基本都被消除(云覆盖区域除外)，拼接结果的色彩一致性大大提高，图 3也证实了这一结果。图 3为图 2数据中选择一列像元值绘制成列剖面曲线(取匀色前后相同的列)。图 3显示，在匀色前的镶嵌影像列剖面曲线中，像元值在镶嵌行所在位置变化明显，曲线呈现阶跃状态；在匀色后的镶嵌影像列剖面曲线中，像元值在镶嵌行过度自然，与其他行无明显区别。

此外，本文还统计了基准影像、匀色影像处理前后各波段的均值，定量评价本文匀色方法对像元值的改变，结果如表 1所示。表 2为表 1数据中各影像对的基准影像波段均值分别减去匀色前后的均值，取绝对值并除以基准影像均值，得出的百分比。表 2中数值越大，表示两景影像差异越大。

表 1、表 2显示，匀色处理前，各影像对的大部分波段均值差的百分比都较大，其中，波段1最大，主要因为波段1(蓝波段)受大气影响最大。而不同数据组之间，GF-1、GF-2混合影像对中波段1的均值差的百分比最大，达到37.12%，主要是因为两景影像是不同传感器在不同时间获取的数据，除成像几何、光照条件不一致外，两传感器之间还存在系统辐射差异。经过匀色处理后，除GF-2影像对的波段4外，其余波段均值差异均明显减小，最大仅为2.99%，表明本文匀色方法较好地对3组试验影像进行匀色处理。对于GF-2影像对的波段4，经过匀色处理后，均值差百分比略有增大，可能的原因是GF-2影像对的成像时间分别是2月28日和12月5日，而2月28日部分植被已开始进入返青期。因此，在近红外波段有明显变化。

4 结束语

多时相遥感影像匀色处理是影像拼接镶嵌前的重要步骤，直接影响遥感影像判读解译、数字化测图等后续处理工作。本文基于多元变化检测算法，自动提取影像重叠部分未变化的像元点，基于未变化像元点建立正交回归模型对影像进行匀色处理。试验数据匀色处理结果表明，本文方法可有效消除GF-1、GF-2数据之间的颜色差异，匀色后镶嵌结果色彩一致性较好，有利于下一步的处理工作。

本文建立的匀色处理方法可自动对多波段影像进行匀色处理，无需人工干预。该方法的核心是基于光谱向量的逐像元变化检测，因此，基准影像和待匀色影像的波段数和空间分辨率需一致，且影像配准误差小于一个像元。