1 引言
近年来，光谱相机已经成为重要的光电遥感载荷，其中滤光片阵列多光谱相机因其较低的平台需求，被广泛应用到无人机多光谱遥感中[1]。滤光片阵列多光谱相机同一时刻可获取多个波段条带图像数据，再对不同时刻获取的多波段条带图像数据进行裁剪拼接处理，得到各单波段图像。由于相机在不同时刻成像曝光量存在一定差异，导致裁剪拼接后单波段图像中不同时刻条带间存在明显灰度差异，影响后续图像处理及应用需求。此外，地物灰度差异使得图像光谱特征提取变得困难、特征匹配时误匹配率增加[2]。为满足滤光片阵列多光谱数据的应用要求，必须进行灰度处理。
目前，图像灰度处理主要用于图像拼接领域[3-5]，较为常用的是对图像进行匀光处理，例如：处理单幅图像的MASK算法[6]、同态滤波算法[7-8]、Retinex算法[9]等，处理多幅图像的Wallis算法[10]、直方图匹配法[11-12]等。但上述算法在处理滤光片阵列多光谱图像的灰度差异时，效果并不理想。针对滤光片阵列多光谱图像灰度差异，方秀秀等[13]提出了图像灰度线性变换方法，通过图像灰度线性变换模型调整各条带图像灰度。但该方法的图像灰度调整效果过于依赖于图像重叠区域匹配点的数量及分布，当重叠区域提取的匹配点较少、无匹配点或分布较为集中时，无法准确建立图像之间的灰度变换关系，或变换关系存在较大偏差。
此外，该方法也是仅考虑了单波段图像间的灰度差异，却忽略了各单波段图像之间的灰度变化关系，灰度调整后多光谱图像地物光谱特征容易出现偏差，无法准确反映地物不同波段特性。
因此，文中提出一种基于重叠区域的滤光片阵列多光谱图像灰度调整算法，用于解决滤光片阵列多光谱图像地物灰度不一致问题，调整后的单波段图像整体灰度一致，且能较好保持地物光谱信息。
2 滤光片阵列多光谱成像特点
滤光片阵列多光谱相机在光路系统中加入了一个滤光板，滤光板放置在CCD探测器前面，在滤光板上沿着垂直于飞行方向镀制了几个窄带的带通滤光膜，每个带通滤光膜只能通过一个谱段的图像，如图1所示(图1中显示为8条滤光片)。
图1 滤光片阵列多光谱相机成像原理示意图
Fig.1 The imaging principle of filter array multispectral camera
多光谱相机的多波段图像均是在同一阵面上成像，需要通过图像拼接的方法才能够形成一幅单波段图像，条带图像裁切拼接如图2所示。
在t1、t2、t3、t4时刻分别获得4幅含有8个波段的条带图像，经过裁切拼接获得单波段图像，如图2中以波段2和波段3为例显示拼接过程。
图2 滤光片阵列多光谱数据阵面拼接处理示意图
Fig.2 Area array processing of filter array multispectral dataset
大部分多光谱相机单波段图像是同时成像[13]，其曝光时间相同，各区域地物灰度一致。但滤光片阵列多光谱相机曝光时间由快门自适应控制。因此，各时刻曝光量可能存在一定差异，裁切拼接后单波段图像可能存在条带状的灰度差异。
例如，图2中t3时刻与t1、t2、t4时刻存在明显的曝光量差异，t3时刻成像总体上偏暗(或偏亮)，裁切拼接后的单波段图像中均出现较暗(或较亮)的条带。
各单波段多光谱图像中出现较暗(或较亮)的条带，会影响到地物光谱信息的准确度，如图3所示，图像采用一维表示。
图3 滤光片阵列多光谱图像拼接示意图
Fig.3 Area array processing of filter array multispectral dataset
由图3可以看出，拼接后图像A处，t2时刻成像时图像曝光量相对偏小，引起λ2波段偏暗，提取的地物波谱特征向量(fλ1, fλ2, fλ3,…, fλ8)明显失真，如图4所示。图4中λ2处偏低(如图3所示)，λ6处偏高(图4中未显示)。
图4 地物光谱特征曲线
Fig.4 Spectral curves of ground objects
3 滤光片阵列多光谱图像灰度调整
滤光片阵列多光谱图像灰度调整技术过程与其他图像的拼接流程基本一致，主要包括空间位置确定、灰度级调整和条带图像裁剪拼接3个步骤。
3.1 图像空间位置确定
图像空间位置确定主要包括条带图像投影变换和平移量确定两部分。
3.1.1 图像投影变换
滤光片阵列多光谱成像时除了记录图像数据外，还记录成像时的平台状态参数信息，例如平台的位置信息(经度L、纬度B、高度H)和姿态信息(俯仰角φ、翻滚角ω、航向角κ等)。因此，选取t1、t2、…、tn时刻n幅多条带图像。根据式(1)和式(2)对各时刻图像进行投影变换[14]。
(1)
(2)
其中: X、Y为投影后像点坐标，单位像素；H为航高；x、y为条带图像上像点坐标，单位毫米；f为相机焦距，单位毫米。a1、a2、a3、b1、b2、b3、c1、c2、c3由式(2)计算得到。式(2)中， φ、ω、κ分别为相机的俯仰角、翻滚角和航向角。
3.1.2 平移量确定
一般情况下，经投影变换后图像只存在平移关系，考虑到经纬度采样误差影响平移量计算精度，故本文采用SIFT算法确定相邻图像间平移关系。首先利用SIFT算法依次提取相邻序列中第ti幅图像和第tj幅图像的特征点[15-16]；然后，利用匹配点对坐标(rik,cik)和(rjk,cjk)差均值，计算ti和tj时刻条带图像的行列平行量(Rij,Cij)，如式(3)所示：
(3)
例如，相邻两幅多光谱条带图像局部SIFT匹配点如图5所示。计算得到的条带图像行列平行量如图3中的R12、C12、R23、C23、R34、C34、…，这样就确定了条带图像空间上的位置关系。
图5 相邻图像的匹配点示意图
Fig.5 Matching points of adjacent images
3.2 各条带图像灰度调整
在对多幅图像进行拼接时，为了减小累计误差，一般选择中间图像为基准[17-18]。而滤光片阵列多光谱图像也是对多幅图像进行灰度调整、图像拼接处理，因此，选择中间序列的多光谱图像作为基准。
图像灰度调整主要利用条带图像重叠区域的灰度均值，通过计算相邻多光谱图像中各个波段重叠区域灰度均值，各波段灰度均值比例的平均值用于多光谱图像灰度调整系数。对于不相邻图像，通过相邻图像间灰度传递，依次进行灰度调整。涉及到关键技术主要包括确定各条带图像重叠区域、确定图像灰度调整系数、图像灰度调整3个方面。
3.2.1 确定各条带图像重叠区域
以8谱段滤光片阵列多光谱相机为例，其多光谱图像数据如图6(a)所示。波段图像过渡区域约为50行，见图6(a)中黄色框线区域。因此，需要计算各波段条带有效区域的4个顶点坐标，并构建条带图像有效区域模板，如图6(b)所示。
图6 多光谱图像
Fig.6 Multi-spectral image
以两幅相邻多光谱图像中第1波段图像为例，按照图像之间平移关系分别将模板1、模板2扩展，如图7(a)、图7(b)所示。再将扩展后模板1和模板2相乘，如图7(c)所示；得到重叠区域模板如图7(d)所示；再将扩展后模板拆分为第1、2幅第1波段原模板1，如图7(e)、图7(f)所示。
图7 图像重叠区域示意图
Fig.7 The image overlap area
然后，模板1分别与2幅多光谱图像中的第1波段图像相乘得到重叠区域图像，如图8所示。图8(a)、图8(b)分别为第1、2幅中的第1波段图像，图8(c)、图8(d)分别为第1、2幅中第1波段的重叠区域图像。
图8 第1波段重叠区域图像
Fig.8 The overlapping area in 1 band image
3.2.2 计算图像重叠区域灰度均值
基于单波段图像重叠区域，可以计算出图像灰度平均值以及相对于相邻影像的灰度系数。但使用单波段图像灰度系数调整图像灰度，容易破坏地物在不同波段下的光谱特性。因此，先计算出相邻图像间各个波段重叠区域灰度平均值的比例系数，然后，对各波段的比例系数求平均值，得到整张图像(各波段)的灰度调整系数。
以8谱段滤光片阵列多光谱序列图像中相邻的第i幅、第j幅图像为例进行调整，如图9所示。
图9 重叠区域图像灰度比调整示意图
Fig.9 Grayscale ratio of the overlapping area image
首先计算第i、 j幅图像中第1、2、3、4、5、6、7、8波段重叠区域的平均灰度值ai1、aj1、ai2、aj2、ai3、aj3、ai4、aj4、ai5、aj5、ai6、aj6、ai7、aj7、ai8、aj8，然后计算8个波段重叠区域的灰度比c1、c2、c3、c4、c5、c6、c7、c8，则灰度调整系数gij为
(4)
gij作为第j幅图像的灰度调整系数，调整后图像重叠区域的灰度与第i幅图像灰度基本一致。
3.2.3 图像灰度调整
求出图像灰度调整系数gij后，就可以对图像进行灰度调整。例如，以9幅多光谱条带图像为例，图像平移关系与重叠区域灰度值示意图如图10。
图10 序列图像平移与重叠区域平均灰度示意图
Fig.10 Sequence image translation with overlapping area average gray scale schematic
选择中间的第5幅图像为基准，依次将各灰度调整系数与对应序列图像中的每个像素灰度相乘，各灰度调整系数计算公式如下：
(5)
3.3 图像裁剪拼接
通过对灰度调整后序列图像的各单波段图像进行裁剪，并按照图像间平移关系进行拼接，就得到了灰度一致的滤光片阵列多光谱单波段图像，如图11所示(为了方便表示，将条带图像按照飞行方向投影形成一维数据进行显示)。
图11 多光谱图像拼接示意图
Fig.11 Multi-spectral image stitching
4 实验与分析
为了验证本文提出算法的可行性，采用了某航拍无人机滤光片阵列多光谱相机拍摄某地区的多光谱图像数据进行验证实验。其中，前6张序列多光谱图像如图12。
图12 实验图像
Fig.12 Experimental images
裁剪拼接和本文灰度调整后第1波段图像如图13。从图13(a)中可以明显看出，直接对多光谱图像数据进行拼接的单波段图像整体灰度不均匀，同一地物存在着较大的灰度差异。图13(b)所示的图像灰度整体比较均匀，说明本文方法较好地解决条带灰度不均的问题。
图13 灰度调整图像
Fig.13 Grayscale adjustment comparison chart
为了更好地说明本文算法的性能，选择了滤光片型多光谱数据灰度一致性校正算法(文献[12])与本文提出的灰度调整算法进行对比实验。文献12与本文算法处理后的单波段图像如图14，对比图像分别为1、3、6、7和8波段图像，其中图14(a)、图14(b)和图14(c)分别为未进行灰度调整、文献12处理和本文算法处理的图像。
由图14可知，文献12处理后的第1波段和第8波段图像较好的解决了图像灰度不均匀问题，但相邻时刻条带图像中相同波段重叠区域较小，提取的同名像素数量较少，容易造成灰度调整系数不准确，使得第3、6、7波段图像中仍存在地物灰度不均匀问题，而且图像中部分条带的灰度差异要比未调整图像还要明显。其局部放大图如图15，其中，图15(a1)、图15(b1)、图15(c1)，图15(a2)、图15(b2)、图15(c2)和图15(a3)、图15(b3)、图15(c3)分别为未进行灰度调整、文献12处理、本文算法处理后的3、6、7波段图像。
图14 灰度调整后图像
Fig.14 Comparison of the effect after grayscale adjustment
图15 灰度调整后局部放大图
Fig.15 Grayscale adjusted partial comparison
由图15可知，文献12处理后的图像灰度不一致现象更加明显。另外，第4与第5波段由于在重叠区域未提取同名像点，则无法进行灰度调整。所以，当重叠区域的匹配点数量较少时，基于匹配点灰度计算出的灰度线性变换模型具有较大偏差，会增加单波段相邻条带图像间的灰度差异。甚至提不出同名像点出现灰度无法调整的局面。
因此，当相邻条带图像重叠率低、重叠区域匹配点数量少以及数量分布不均匀时，文献12的算法无法解决大图像灰度不一致问题，甚至会加大图像中同一区域中地物的灰度差异。而本文算法调整后各单波段图像灰度分布均匀，较好解决了多光谱单波段图像灰度不一致问题。
为进一步说明本文算法优越性，分别计算调整后各单波段条带图像重叠区域灰度平均值的差值来说明，如图16所示(仅显示1、3、7波段图像)。
图16纵坐标为重叠区域平均灰度差，横坐标为幅序号。如图16(a)中幅序号为1时本文算法灰度差为0.001 7，是通过计算第一、二幅第一个波段重叠区域平均灰度差得到的，其他波段方法相同。由图16可知，文献12处理后单波段图像重叠区域灰度均值的差值整体变化较大，说明调整后图像部分区域仍存在较大灰度差异。而本文算法中，重叠区域灰度均值差值差异较小，近似一条直线，说明相邻条带图像灰度变化较小，图像灰度近似一致。
图16 重叠区域灰度均值差值曲线
Fig.16 Difference in mean grey value of overlapping areas
由此可知，本文基于相邻图像的重叠区域灰度均值计算的调整系数，充分利用了重叠区域图像的灰度信息，调整后各单波段图像灰度均匀，有效解决了单波段图像中灰度不一致问题。
5 结论
针对滤光片阵列多光谱图像灰度不一致的问题，提出了一种基于重叠区域的多光谱图像灰度调整算法，即根据各波段图像重叠区域灰度均值比的平均值，对整张图像各波段进行灰度调整，调整后的拼接图像同一地物灰度一致，降低了对后期图像配准过程中特征提取与匹配的影响。通过实验对比，表明本文针对滤光片阵列多光谱图像提出的灰度调整算法较目前算法有更好的灰度处理效果，较好解决了滤光片阵列式多光谱图像灰度不一致问题，可以满足影像的后期处理及应用需求。
[1] 赵宝玮.机载多光谱相机数据的高质量获取与处理技术研究[D].西安：西安电子科技大学， 2014.
Zhao B W.Study of high quality data acquisition and processing technology for airborne multi-spectral camera[D].Xi’an:Xidian University,2014.
[2] Banerjee B P,Raval S A,Cullen P J.Alignment of uav-hyperspectral bands using keypoint descriptors in a spectrally complex environment[J].Remote Sensing Letters,2018,9(4/6):524-533.
[3] 吕楠，苟永刚，龙川，等.多相机图像拼接匀色算法[J].测绘通报， 2016(07):44-47.
Lv N,Gou Y G,Long C,et al.An image dodging algorithm based on multiple cameras[J].Bulletin of Surveying and Mapping,2016(07):44-47.
[4] 杨国鹏，周欣，韦红波，等.面阵摆扫航空相机序列图像的大区域无缝拼接[J].测绘科学， 2020,45(03):46-52,60.
Yang G P,Zhou X,Wei H B,et al,Large-area seamless mosaic for the sequence images from a frame-sweep aerial cameras with plane array cameras with plane array[J].Science of Surveying and Mapping,2020,45(03):46-52,60.
[5] 么鸿原，王海鹏，林雪原，等.遥感图像拼接中改进的图像预处理算法研究[J].计算机与数字工程， 2020,48(02):428-432.
Yao H Y,Wang H P,Lin X Y,et al,Research on improved image preprocessing algorithm in remote sensing image mosaicing[J].Computer & Digital Engineering,2020,48(02):428-432.
[6] 李烁，王慧，王利勇，等.遥感影像变分Mask自适应匀光算法[J].遥感学报，2018,22(03):450-457.
Li S Wang H,Wang L Y,et al.Adaptive dodging method based on variational mask for remote sensing images[J].National Remote Sensing Bulletin,2018,22(03):450-457.
[7] Kaur K,Jindal N,Singh K.Improved homomorphic filtering using fractional derivatives for enhancement of low contrast and non-uniformly illuminated images[J].Multimedia Tools and Applications,2019,78(19):27891-27914.
[8] Zaheer Ud Din S,Suganthi K.Image contrast enhancement by homomorphic filtering based parametric fuzzy transform[J].Procedia Computer Science,2019,165(14):166-172.
[9] 李烁.遥感影像色彩一致性处理和镶嵌线优化技术研究[D].战略支援部队信息工程大学，2018.
Li S.Research on the optimization of color consistency processing and seamline determination of remote sensing image[D].PLA Strategic Support Force Information Engineering University,2018.
[10] 范铀，陈旭帅，王丹丹，等.超分辨率重建影像Wallis匀光拼接算法[J].测绘通报，2018(02):136-141.
Fan Y,Chen X S,Wang D D,et al.Wallis dodging and mosaic algorithm for super resolution reconstruction images[J].Bulletin of Surveying and Mapping,2018(02):136-141.
[11] Ashiba M I,Tolba M S,El-Fishawy A S,et al.Gamma correction enhancement of infrared night vision images using histogram processing[J].Multimedia Tools and Applications,2019,78(19):27771-27783.
[12] 方秀秀，黄旻，王德志，等.基于高程和地物光谱约束的多光谱图像预处理算法[J].半导体光电， 2020,41(02):264-267,272.
Fang X X,Huang W,Wang D Z,et al.Multispectral image preprocessing based on elevation and surface feature spectrum constraints[J].Semiconductor Optoelectronics,2020,41(02):264-267,272.
[13] 朱豪男，胡孟晗，张健，等.低成本便携式多光谱成像系统的研发及优化[J].中国图象图形学报， 2021,26(08):1796-1808.
Zhu H N,Hu M H,Zhang J,et a.Development and optimization of a low-cost and portable multispectral imaging system[J].Journal of Image and Graphics,2021,26(08):1796-1808.
[14] 冯相辉.一种改进的同态滤波图像增强算法[J].重庆邮电大学学报(自然科学版)， 2020,32(01):138-145.
Feng X H.An improved homomorphic filtering image enhancement algorithm[J].Journal of Chongqing University of Posts and Telecommunications(Natural Science Edition),2020,32(01):138-145.
[15] Jia Y J,Su Z B,Zhang Q,et al.Research on uav remote sensing image mosaic method based on sift[J].International Journal of Signal Processing Image Processing & Pattern Recognition,2015,8(11):365-374.
[16] Zhao J Q,Zhang X H,Gao C X,et al.Rapid mosaicking of unmanned aerial vehicle (uav) images for crop growth monitoring using the sift algorithm[J].Remote Sensing,2019,11(10):1226-1245.
[17] Ghosh D,Kaabouch N.A survey on image mosaicing techniques[J].Journal of Visual Communication and Image Representation,2016,34(03):1-11.
[18] 朱标.基于优化ORB算法的图像角点特征匹配方法[J].重庆工商大学学报(自然科学版)，2021,38(04):42-51.
Zhu B.An image corner features matching method based on optimized orb algorithm[J].Journal of Chongqing Technology and Business University(Natural Science Edition),2021,38(04):42-51.