2. 2. 浙江工业大学 计算机科学与技术学院,浙江 杭州 310023；
3. 中国科学院 新疆生态与地理研究所,新疆 乌鲁木齐 830011
2. College of Computer Science and Technology, Zhejiang University of Technology, Hangzhou 310023, China;
3. Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences，Urumqi 830011,China
1 引 言
影像镶嵌是在几何配准的基础上将两景或者多景具有重叠区域的影像合成包含整个研究区影像的过程，包含配准和合成两个步骤[1]。受获取时刻天气、光照以及时相等条件的影响，大区域的待镶嵌影像之间存在对比度、亮度分布等色彩不一致现象。色彩一致性处理就是通过相邻影像之间重叠区建立对应关系，消除影像间的色彩差异，使得待镶嵌各景影像具有相同或者相近的色彩特征。目前，常用的色彩一致性方法主要有色彩匹配法和线性回归法。
色彩匹配法是直接建立影像重叠区的灰度映射关系。文献[2]提出了一种将参考影像的均值和方差作为影像特征，并通过构建单调的映射方程，使得处理后的影像具有与参考影像相同的均值和方差；文献[3, 4]等采用重叠区域直方图匹配建立从输入影像到参考影像的灰度映射方程，从而使两景影像具有相似的亮度分布。色彩匹配法根据统计特征或者分布特征直接建立影像之间灰度值的变换关系，容易造成原始光谱特征的扭曲，不利于后续的应用；地物类型变化、云和云阴影都会造成均值、方差和直方图的改变，从而影响处理的精度。
线性回归方法从两幅影像重叠区内选取时相不变特征(pseudo invariant features,PIFs)，建立两者亮度变化的回归关系，对输入影像进行处理，并消除色差。针对PIFs的选取，文献[5]从近红外波段的散点图中，寻找大范围陆地和水体，以其中心作为样本点，构建线性变化方程；文献[6]认为经主成分变换后，PIFs分布在靠近第一主成分周围一定范围内，采用阈值分割的方法选取这些像素作为PIFs；文献[7]通过典型相关分析去除影像之间的相关性的多元变化检测(multivariate alteration detection，MAD)提取PIFs；随后又提出根据前一次确定的PIFs可能性大小赋予相应权重进行迭代加权的多元变化检测(iteratively reweighted MAD,IRMAD)。文献[9]将其应用于航空遥感影像的镶嵌之中，取得了较好的效果。简单线性回归、最小二乘回归、正交回归以及具有抗噪声能力的Theil-Sen回归[10, 11]都能够在正确选取PIFs的基础上得到准确的回归方程。因而，线性回归方法成功应用的关键在于重叠区域选取PIFs。
大区域影像镶嵌中由于空间跨度过大，造成相同时相的影像难于收集，且不同时段、不同区域受光照、地形条件的影响也更为明显，特别是植被在不同时相的影像上光谱差异较大，以上变化都会造成直方图形状、均值、亮度值的改变，限制了色彩匹配法所能达到的精度和PIFs选取的效果。而任何一景影像色彩一致性处理失真或者色彩一致性处理不彻底，都会造成后续处理影像的误差传递和累积。针对上述问题，本文考虑通过重叠区域所有的像素来建立变换关系，在获得较好的视觉效果的基础上，提高色彩一致性处理稳定性。
2 色彩一致性处理方法在不考虑地物类型发生变化的情况下，引起影像I1和I2重叠区域灰度值变化的因素可以概括为
式中，E1表示植被等自然地物随时相变化，该因素仅对植被等自然要素产生影响，并随着各种植被物候的差异性，而表现出不同的变化模式，变化相对缓慢；E2表示云/云阴影遮挡造成的色彩畸变，在空间分布上的随机性，可能遮挡任何类型的地物，表现出突变和变化幅度大的特点；E3表示光照条件及其成像系统畸变引起的变化，对影像上的所有像素均发生影响，可以认为这种因素在空间分布上是均匀的，并可以通过线性方程来表示。
将同一地区的两景影像分别作为两个随机变量，统计联合概率密度并构建概率密度图，所有仅受E3影响的像素(称为稳定特征)会形成一个“脊”(ridge),而通过脊的直线包含了两景影像之间亮度变化的线性关系[12]；受E1影响的地物会分布在脊的周围，受E2影响的像素分布在远离脊的位置，但是都不会影响脊的位置。因而，脊具有较高的稳定性，已经用于变化检测、分类等的相对辐射校正[13]。但是，目前脊提取需要目视判别，对于大范围影像镶嵌中的色彩一致性处理并不现实。这里提出一种脊自动提取的方法并将其应用到大范围影像镶嵌中的色彩一致性处理。
2.1 联合概率密度图与脊一般的，可以把不同时相获取的两景影像I1和I2看做描述地表辐射特征的二维随机变量，在获取过程中受E1、E2和E3等因素的影响。两景影像上对应波段之间的亮度值变化的关系可以用联合概率密度来表示。影像获取过程中的采样对地表辐射值进行离散化，该二维随机变量是一种离散性的随机变量，概率密度的计算方法可以表示为
式中，P(x,y)和N(x,y)表示影像I1和I2上亮度值分别为x和y的像素出现的概率和频数，概率密度函数构成了一个以亮度值取值范围为值域的二维空间。
图 1(a)和(b)分别为某一地区不同时相获取的TM 5两期影像，大小为1200像素×1200像素，影像之间存在上述3种因素引起的亮度值变化，红波段联合概率密度如图 1(c)所示，横坐标和纵坐标分别为参考影像和输入影像的亮度值，颜色由红、黄、蓝的渐变表示联合概率密度逐渐减小，图 1(d)是将概率密度作为Z轴三维显示的结果(图 1)。
|
| 图 1 联合概率密度图 Fig. 1 The joint probability density map |
从图中可以看出：联合概率密度图呈现出一种多层分布的结构，核心为一个不规则椭圆形的红色区域，主要为受E3影响和E1影响的像素，亮度值变化较小，呈现以直线l为对称轴的近似轴对称分布，且分布在直线l上点概率密度较高。红色区域周围的黄色过渡区域和远离直线l的蓝色的椒盐状散点，这些区域亮度值变化较大且概率密度较小，主要是受E2的影响。根据以上分析，直线l上的点具有亮度变化小，且概率密度较高的特点，受E1和E2地物分布在直线l的周围，因而可以认为直线l为待求的脊。
2.2 脊提取方法根据2.1节的分析，脊提取的问题就可以转换为求取椭圆的长半轴，这里认为脊具有两个特征：① 经过脊上点的概率密度大；② 分布在脊左右两侧的概率密度尽可能平衡。任意通过概率密度不为0区域确定的直线可以将概率密度图划分为3个部分：线的左侧、线的右侧和线上，分别统计3个部分概率分布P(L)、P(R)和P(ON)(图 2)。从概率密度所有可能的直线中，选择满足上述两个条件的直线作为待求方程。
|
| 图 2 脊提取方法 Fig. 2 The ridge extraction method |
通过概率密度不为0的区域可以确定无数条直线，这就给脊提取带来了困难。根据条件(1)的线上概率密度尽可能大，可以认为通过脊的直线至少经过两个概率密度不为0的点，并进行离散化。假设两景影像的灰度级均为m，可知概率密度图不为0的点的数目一定小于m×m，从而可以采用枚举法取其中最满足上述两个条件的点作为离散条件下的最优解。具体实现方法如下：
(1) 假设概率密度图上有n个P(x,y)≠0的点，从中任意选取两个点，构成一条直线，一共可以构成M条直线，记为
式中，M=n(n－1)/2，斜率和截距分别记为(为了说明问题的方便，这里不考虑斜率无穷大的情况，这种情况的直线对结果没有意义)
(2) 从M条直线中任取一条li，其直线方程为y=kix+bi，该直线将概率密度的二维空间分成了3个部分，直线的左侧、直线的右侧和直线上，分别记这3个部分的概率分布为P(L)、P(R)和P(ON)。
(3) 重复第(2)步，计算M条直线的P(L)、P(R)和P(ON)，并统计P(ON)的最大值P(ON)max；
(4) 计算M条直线的可选系数
式中，|P(L)－P(R)|表示左右两侧概率分布平衡性；|P(ON)－P(ON)max|表示和最大值接近程度。W值越小，表明是待求方程的可能性越高。
(5) 从M条直线中选取可选系数最小的直线作为校正方程。
计算影像上各个对应波段的变换方程，逐波段对色彩进行处理，并得到色彩一致的结果影像。
3 试验与分析选取覆盖新疆维吾尔自治区113景Landsat TM5影像进行试验，主要是2009年和2010年无云或者少云(小于30%)影像，时相分布上为 8月份到10月份(图 3)，该时相内的影像具有丰富的地表信息，利于目视解译和后续的处理，且天气晴朗，影像质量较高。影像为从USGS网站获取的1 T级产品，已经进行了标准地形校正(standard terrain correction)，几何定位精度达到0.5像素以上[14]。影像未进行辐射校正和大气校正，直接采用亮度值进行处理。大区域镶嵌首先将影像从UTM分带投影转为等积割圆锥投影(即阿尔伯斯投影，投影参数为标准纬线25°N和47°N，中央经线105°E)。
|
| 图 3 影像的时相分布 Fig. 3 he time and phase distribution of the images |
大区域影像镶嵌的问题可以定义为图
式中，设G表示影像镶嵌过程;V表示一景影像;E表示相邻的影像之间的边[15]。色彩处理优化的问题就转化为从图G中提取最短路径生成树(shortest path spanning tree,SPST)，从而使得传递的次数最少，文献[16]采用Dijkstra方法来确定影像镶嵌顺序；若E以色彩处理或者配准误差为权重，可以通过求取图G的最小代价生成树(minimum spanning tree,MST)使连接像对边总体误差最小来实现全局误差最小化，如文献[17]采用最小路由代价生成树对大规模显微影像进行镶嵌。本文的主要目标是提供一种大范围影像镶嵌中色彩一致性处理的稳健方法，对色彩一致性处理采用一种类似于SPST的策略，具体实现方法为：通过选择位于中心区域的且视觉特征良好影像作为参考，以减少传递次数，并对邻接的影像进行色彩一致性处理。重复上述过程，以进行色彩一致性处理结果影像作为参考对邻接的影像进行处理，直到所有影像都进行了处理，该过程结束。对于待处理的影像有多景相邻已经处理过的影像，则使用所有的已经进行处理过的影像作为参考。接缝线处理是另一个关键问题，文献[18, 19, 20]均有讨论，这里就不再赘述，本文试验直接将色彩一致性处理的结果利用商业软件ERDAS[21]进行处理和生成镶嵌结果。
原始的标准假彩色合成影像存在相邻轨道的 “条带性”以及同一轨道的“补丁式”的色彩差异(图 4(a))，按照本文方法进行色彩一致性处理之后镶嵌结果整体色彩基本一致，目视判读难以看出相对色彩差异，也不存在明显的接缝线；对于原始偏暗的影像也进行拉伸，对比度得到增强；靠近边缘的影像没有出现色彩畸变，说明误差的累积现象不明显，本文方法具有较好的稳定性。
|
| 图 4 中国新疆地区影像镶嵌结果 Fig. 4 The mosaic result of Xinjiang, China |
为了对算法定量评价和进行算法之间的比较，选取了图 4(a)中两组具有代表性的影像进行试验，并选取了一组SPOT 5的影像，记为试验1、2和3，元数据见表 1，原始数据叠加显示见图 5(a)、(d)和(g)。试验1中存在获取时刻天气条件引起的整体色彩差异以及植被随时相变化引起的色彩差异，重叠区域主要为植被和裸地，采用IRMAD[8]选取的样本主要为裸地，IRMAD方法处理结果(图 5(c))与本文所述方法(图 5(b))相比较，两侧的色彩仍然存在明显的整体差异，即校正结果仍存在较大的残差。试验2中输入影像部分地表被雪覆盖，导致参考影像和输入影像直方图形状发生变化，若强制进行直方图匹配，造成图 5(f)所示的雪呈现绿色的色彩畸变，而本文方法由于采用线性回归，保持了色彩之间的相对关系，雪的颜色未出现失真。试验3中包含植被随季节变化以及农作物收割引起的色彩变化，采用局部直方图匹配方法结果如图 5(i)所示，结果中出现了较为严重的色彩畸变，如蓝框内的河流出现红色，而本文方法(图 5(h))色彩未出现扭曲的现象。(图 4、图 5)
| 影像 | 获取日期 | path/row | 色彩差异原因 | |
| 试验1 | 参考 | 2009-08-08 | 144/29 | 天气条件和时相变化 |
| 输入 | 2009-10-02 | 145/29 | ||
| 试验2 | 参考 | 2009-10-21 | 149/32 | 部分地表被雪覆盖 |
| 输入 | 2007-09-23 | 150/32 | ||
| 试验3 | 参考 | 2008-09-07 | 214/260 | 植被时相变化、作物收割等 |
| 输入 | 2008-10-08 | 214/261 | ||
|
| 图 5 影像镶嵌效果对比 Fig. 5 Comparison of the image mosaic |
将输入影像划分为重叠区域(O)和非重叠区域(N)。统计本文方法、IRMAD方法、直方图匹配方法得到的结果影像红波段的均值(μ)、标准差(σ)、色彩一致性处理影像的相关系数(C)和灰度级(G)等特征。其中，均方根误差的计算方法为
式中，fref和fpre分别表示参考和预测的灰度值。重叠区域的残差说明影像的预测精度，残差越小说明预测越准确。统计结果如表 2所示。
| 输入影像 | 本文方法 | 对比方法 | ||||||||||||
| μ | σ | G | μ | σ | G | C | RMSE | μ | σ | G | C | RMSE | ||
| 试验1 | O | 61.5 | 36.2 | 125 | 71.3 | 42.5 | 145 | 1 | 3.6 | 69.2 | 39.3 | 125 | 1 | 5.6 |
| N | 63.7 | 38.9 | 145 | 74.6 | 45.7 | 160 | 1 | - | 72.6 | 42.1 | 147 | 1 | - | |
| 试验2 | O | 74.3 | 39.6 | 157 | 84.1 | 40.8 | 167 | 1 | 2.5 | 85.1 | 40.7 | 160 | 0.78 | 3.5 |
| N | 78.9 | 39.7 | 162 | 85.2 | 42.3 | 170 | 1 | - | 86.2 | 41.3 | 167 | 0.75 | - | |
| 试验3 | O | 100.6 | 79.1 | 207 | 109.2 | 83.9 | 219 | 1 | 4.5 | 111.3 | 81.2 | 208 | 0.83 | 6.7 |
| N | 115.9 | 80.1 | 218 | 117.9 | 83.9 | 225 | 1 | - | 116.1 | 82.3 | 217 | 0.79 | - | |
试验中本文方法的RMSE显著小于对比方法所得出的残差，说明具有较高的预测精度；试验2与试验3中，对比方法的结果影像与输入影像的相关系数较低，说明两者的线性相关性较低，亮度值出现非线性的扭曲。两组试验中本文的预测结果影像的方差(σ)和灰度级(G)均大于对比方法，说明影像的对比度得到了增强。
4 结 论针对不同时相遥感影像镶嵌中的色彩差异问题，提出了一种根据联合概率密度图中的脊进行色彩一致性处理的方法，采用中国新疆地区113景TM5影像镶嵌进行了试验，镶嵌结果具有较高的稳定性和较好的视觉效果，表明适合于大范围影像的镶嵌处理。由于通过考虑所有像素提取“脊”，因而方法稳定；采用线性方程进行色彩调整，避免色彩的非线性变化，因而具有保持原始色彩的相对关系。
但是本文方法也存在着以下问题：① 主要是采用形态学特征进行提取，而没有严格考虑构成样本点的具体地物类型，缺乏物理模型的支持；② 对大区域影像镶嵌在色彩一致性处理中的误差传递与累积的问题没有进行深入的研究，而是采取了一种较为简单的两两色彩一致性处理的策略。以上两个问题都需要作进一步研究。
| [1] | SOILLE P. Morphological Image Compositing[J]. IEEE Tran sactions on Pattern Analysis and Machine Intelligence, 2006,28(5):673-683. |
| [2] | XIAO Fu, WU Huizhong, XIAO Liang,et al. An Ambient Light Independent Image Mosaic Algorithm[J]. Journal of Image and Graphics, 2007,12(9):1671-1675.(肖甫,吴慧中,肖亮,等.一种光照鲁棒的图像拼接融合算法[J].中国图象图形学报,2007,12(9):1671-1675.) |
| [3] | CHEN Jianle, LIU Jilin, YE Jianhong, et al. Luminance and Chrominance Correction for Multi-view Video Using Overlapped Local Histogram Matching[J]. Journal of Image and Graphics, 2007,12(11):1992-1999.(陈建乐,刘济林,叶建洪,等.多视点视频中基于局部直方图匹配的亮度和色差校正[J].中国图象图形学报, 2007,12(11):1992-1999.) |
| [4] | HELMER E H, RUEFENACHT B. Cloud-free Satellite Image Mosaics with Regression Trees and Histogram Matching[J]. Photogrammetric Engineering and Remote Sensing, 2005,71(9): 1079-1089. |
| [5] | ELVIDGE C D, YUAN D, WEERACKOON R D, et al. Relative Radiometric Normalization of Landsat Multispectral Scanner(MSS) Data Using an Automatic Scatter Gram-controlled Regression [J].Photogrammetric Engineering and Remote Sensing, 1995,61(10): 1255-1260. |
| [6] | DU Yong. Radiometric Normalization, Compositing, and Quality Control for Satellite High Resolution Image Mosaics over Large Areas[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001,39(3):623-634. |
| [7] | CANTY M, NIELSEN A. Automatic Radiometric Normali-zation of Multitemporal Satellite Imagery[J].Remote Sensing of Environment, 2004,91(5): 441 - 451. |
| [8] | CANTY M. Automatic Radiometric Normalization of Multitemporal Satellite Imagery with the Iteratively Re-weighted MAD Transformation[J]. Remote Sensing of Environment, 2008,112(3):1025-1036. |
| [9] | PAN Jun, WANG Mi, LI Deren. A Network-based Radiometric Equalization Approach for Digital Aerial Ortho Image[J]. IEEE Geosciences and Remote Sensing Letters, 2011,7(2): 401-405. |
| [10] | OLTHOF I, POULIOT D, FERNANDES R, et al. Landsat -7 ETM + Radiometric Normalization Comparison for Northern Mapping Applications[J]. Remote Sensing of Environment, 2005,95(4):388-398. |
| [11] | HEALEY S, YANG Z, COHEN W. et al. Application of Two Regression-based Methods to Estimate the Effects of Partial Harvest on Forest Structure Using Landsat Data[J]. Remote Sensing of Environment, 2006,10(1):115-126. |
| [12] | SONG C, WOODCOCK C, SETO K. Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects?[J]. Remote Sensing of Environment, 2001, 75 (2): 230-244. |
| [13] | SCHROEDER T, COHEN W, SONG C.et al. Radiometric Correction of Multi-temporal Landsat Data for Characterization of Early Successional Forest Patterns in Western Oregon[J]. Remote Sensing of Environment, 2006, 103(1):16-26. |
| [14] | COLLECTIVE. U.S. Geological Survey[EB/OL].[2012-10-08].http://www.usgs.gov/. |
| [15] | ZHOU H. Graph-based Global Optimization for the Registration of a Set of Image[J].Lecture Notes in Computer Science, 2006,4319: 1206-1214. |
| [16] | PAN Jun. Research on Automated Color Consistency Processing and Generation of Seam Line Network for Aerial Images[D]. Wuhan:Wuhan University,2008.(潘俊.自动化的航空影像色彩一致性处理及接缝线网络生成方法研究[D].武汉:武汉大学,2008.) |
| [17] | GONG Yongxi, TIAN Yuan, XIE Yubao, et al. A Method for Large-scale Microscope Images Mosaicing Based on Minimum Routing Cost Spanning Tree [J]. Journal of Image and Graphics, 2009,14(6):1178-1187.(龚咏喜,田原,谢玉波,等. 基于最小路由代价树的大规模显微图像拼接方法[J].中国图象图形学报, 2009,14(6):1178-1187.) |
| [18] | PAN Jun, WANG Mi, LI Deren. Approach for Automatic Generation of Optimization of Seamline Network[J]. Acta Geodaetica et Cartographica Sinica,2010,39(3):289-294.(潘俊,王密,李德仁.接缝线网络的自动生成与优化方法[J].测绘学报,2010,39(3):289-294.) |
| [19] | YUAN Xiuxiao, ZHONG Can. An Improvement of Minimizing Local Maximum Algorithm on Searching Seam Line for Orthoimage Mosaicing[J]. Acta Geodaetica et Cartographica Sinica,2012,41(2):199-204.(袁修孝,钟灿.一种改进的正射影像镶嵌线最小化最大搜索算法[J].测绘学报,2012,41(2):199-204.) |
| [20] | ZUO Zhiquan, ZHANG Zuxun, ZHANG Jianqing, et al.Seam Line Intelligent Detection in Large Urban Orthoimage Mosaicing [J]. Acta Geodaetica et Cartographica Sinica, 2011,40(1):84-89.(左志权,张祖勋,张剑清,等.DSM辅助下城区大比例尺正射影像镶嵌线智能检测[J].测绘学报,2011,40(1):84-89.) |
| [21] | COLLECTIVE.Intergraph[EB/OL].[2012-10-08].http://geospatial.intergraph.com/Homepage.aspx. |