在遥感领域中，大范围的镶嵌影像是影像分析和应用中的重要数据源。但是，受获取时刻光照、天气以及时相等条件的影响，待镶嵌影像之间存在着不同程度的色彩差异和对比度差异，使影像镶嵌更为复杂，也更为困难。因此，对区域范围内的多幅影像进行匀色，也称色彩一致性处理(吴炜 等，2012)，消除影像之间色彩差异，使测区内影像在色调上保持一致，降低影像镶嵌的难度，具有重要的现实意义。

现有方法可以分为两类，一类是非线性方法，如直方图匹配法(Helmer和Ruefenacht，2005；陈建乐 等，2007)和Gamma校正法(Zhou 等，2015)；另一类是线性变换法(Canty和Nielsen，2008；吴炜 等，2013)，也是最受关注的方法，这类方法倾向于从影像的重叠区域选取不变像素作为样本参与统计，然后利用线性模型进行相对辐射校正。对于像素样本的筛选，张鹏强等(2006)采用人工选取的方法，Chen等(2014)使用迭代加权的多元变化检测法，Zhang等(2014)利用迭代慢特征分析法，Lin等(2015)采用加权主成分分析法。这些方法虽然可以不同程度地解决色彩不一致问题，但是通常不能很好地消除对比度不一致现象。然而进行大范围的影像镶嵌时，影像之间常有色彩不一致和对比度不一致同时存在的情况。因此，一种基于均值和方差的特殊线性变换即Wallis变换，引入到色彩一致性处理中。李治江(2005)采用全局Wallis变换对相邻影像进行色调匹配。王密和潘俊(2006)提出了一种兼顾整体与局部的匀色方法，在利用Wallis变换对影像进行处理后，采用直方图匹配法对重叠区域进一步处理。孙明伟(2009)在Wallis变换基础上提出了基于最小二乘区域网平差的匀色处理方法，整体计算区域中每幅影像的色调调整参数。李烁等(2016)针对正射影像提出了一种基于重叠区域的Wallis匀色方法，根据相邻影像有效区域的重叠区域的均值和方差信息进行匀色处理。上述方法本质上都是基于Wallis变换且利用同一线性关系进行匀色处理，但是影像中的地物复杂多样，色彩差异程度也不相同，采用同一个线性关系不能很好地消除影像之间的色彩差异。

此外，随着遥感数据获取方式朝着多平台、多轨道、多传感器、多光谱、高分辨率、短周期访问等方向发展，使得获取的遥感数据呈爆炸性增长(闸旋，2013)，使影像匀色处理的效率面临挑战。传统基于CPU的串行处理模式运算能力已明显不足，无法满足在应急应用领域中对海量遥感数据快处理需求(方留杨 等，2013)。随着GPU(Graphic Processing Unit)通用计算能力的不断增强，为提高遥感数据处理效率提供了新的有效手段。基于GPU的并行技术在遥感领域中的应用日益广泛，越来越多的算法以并行结构被重新设计(王宗跃 等，2014；常方正 等，2016)，但是在影像匀色处理中应用相对较少。

针对上述问题，本文在Wallis变换的基础上，提出一种基于GPU的自适应分块加权Wallis并行匀色方法。一方面，综合利用线性模型和非线性模型的优势，从整体和局部两方面提高匀色处理的质量。另一方面，对匀色方法的并行性进行分析并设计相应的并行算法，提高匀色处理的效率。

Wallis匀色法是一种特殊的线性变换法，通过将待处理影像的均值和标准偏差映射到参考影像的均值和标准偏差，实现不同影像间的色彩均衡。线性数学模型如下

$f(x, y) = g(x, y){r_1} + {r_0}$

(1)

式中， $g(x, y)$ 为待处理影像的灰度值， $f(x, y)$ 为处理后的结果影像的灰度值。 ${r_0}$ 为加性系数， ${r_1}$ 为乘性系数，可以表示为

${r_0} = b{m_f} + (1 - b - {r_1}){\kern 1pt} {m_g}$

(2)

${r_1} = \frac{{c{s_f}}}{{c{s_g} + (1 - c){s_f}}}$

(3)

式中， ${m_g}$ 和 ${m_f}$ 分别为待处理影像和参考影像的灰度均值， ${s_g}$ 和 ${s_f}$ 分别为待处理影像和参考影像的标准偏差。 $b$ 为影像亮度系数， $b \in [0, 1]$ 。 $c$ 为影像方差扩展系数， $c \in [0, 1]$ 。

Wallis变换的目的是将处理影像的灰度均值和标准偏差分别被强制到 ${m_f}$ 和 ${s_f}$ ，通常取 $b{\rm{ = }}1$ 、 ${{c = }}1$ ，此时式(1)变为

$f(x, y) = \left( {g(x, y) - {m_g}} \right) \cdot \left({s_f}/{s_g}\right) + {m_f}$

(4)

由式(4)可知，Wallis匀色法是根据影像整体的均值和标准偏差，利用同一个线性关系对影像的每一个像素进行处理。但是影像中的地物复杂多样，颜色信息也各不相同，影像整体的均值和标准偏差无法准确地反映局部地物色彩特征，采用同一个线性关系显然是不合理的。

为了解决上述问题，本文将待处理影像进行分块并统计各块的均值和标准偏差，然后利用双线性插值获取每一个像素对应的线性变换参数，对不同像素采用不同的线性关系进行处理。

采用分块策略时，匀色处理的质量受分块数目影响。若分块数目太多，即各影像块太小时，容易过度校正造成地物失真偏色并且计算量较大；若分块数目太少，即各影像块太大时，统计的均值方差不能准确反映地物分布，不能很好地消除影像之间的色彩差异。

为了提高匀色处理的质量，本文利用变异系数(Zhang 等，2017)自适应确定分块大小。变异系数是标准偏差与均值之比，也称离散系数，可以描述影像内地物的丰富程度。当变异系数越大时，影像内地物种类越丰富，对应的分块数目应该越大即影像块越小，才能获得较好地匀色效果。因此，本文采用如下公式计算自适应确定最优分块数目

$W = r \times w, {\kern 1pt} {\kern 1pt} {\kern 1pt} H = r \times h$

(5)

$r = \frac{{CV}}{{C{V_{{\rm{Ref}}}}}}$

(6)

式中， ${{CV}}$ 为影像的变异系数， $C{V_{{\rm{Ref}}}}$ 为参考影像的变异系数。 ${{w}}$ 和 ${{h}}$ 分别为预设的行、列方向的参考分块数，对于航空遥感影像，均可设置为8。

在确定影像分块大小后，计算各影像块的灰度均值和标准偏差。为了避免出现“块效应”，在对各影像块处理时，考虑相邻影像块的均值和标准偏差，利用双线性插值计算每一个像素的线性变换参数。具体步骤如下：

(1) 将待处理影像分为互不重叠的 $W \times H$ 个影像块，计算各影像块的灰度均值和标准偏差。

(2) 计算各影像块角点对应的均值和标准偏差，若角点只属于一个影像块，则将该影像块的均值和标准偏差赋给该角点；若角点为多个相邻影像块之间的公共角点，则将所属多个影像块参数的平均值赋给该角点。

(3) 将所属影像块的四个角点按照距离赋予不同的权值相加，即利用双线性插值计算每一个像素对应的均值和标准偏差。

对于影像块 $B(w, h)$ 中的一点 $p(x, y)$ ，如图1所示，其均值和标准偏差由角点 ${P_{w, h}}$ 、 ${P_{w + 1, h}}$ 、 ${P_{w, h + 1}}$ 、 ${P_{w + 1, h + 1}}$ 的均值和标准偏差以及该点到影像块边缘的距离 $\Delta x{\text{、}}\Delta {\rm{y}}$ 确定，计算公式如下

$\begin{split} m(x, y) = & \left(1 - \frac{{\Delta x}}{X}\right)\left(1 - \frac{{\Delta y}}{{{Y}}}\right)m(w, h) + \\ & \frac{{\Delta x}}{X}\left(1 - \frac{{\Delta y}}{{{Y}}}\right)m(w + 1, h) + \\ & \left(1 - \frac{{\Delta x}}{X}\right)\frac{{\Delta y}}{{{Y}}}m(w, h + 1) + \\ & \frac{{\Delta x\Delta y}}{{XY}}m(w + 1, h + 1) \end{split}$

(7)

$\begin{split} s(x, y) =& \left(1 - \frac{{\Delta x}}{X}\right)\left(1 - \frac{{\Delta y}}{{{Y}}}\right)s(w, h) + \\ & \frac{{\Delta x}}{X}\left(1 - \frac{{\Delta y}}{{{Y}}}\right)s(w + 1, h)+ \\ & \left(1 - \frac{{\Delta x}}{X}\right)\frac{{\Delta y}}{{{Y}}}s(w, h + 1) + \\ & \frac{{\Delta x\Delta y}}{{XY}}s(w + 1, h + 1) \end{split} $

(8)

式中， $m(x, y)$ 和 $s(x, y)$ 分别为点 $p(x, y)$ 的均值和标准偏差， $m(w, h)$ 、 $m(w + 1, h)$ 、 $m(w, h + 1)$ 、 $m(w + 1, h + 1)$ 分别为角点 ${P_{w, h}}$ 、 ${P_{w + 1, h}}$ 、 ${P_{w, h + 1}}$ 、 ${P_{w + 1, h + 1}}$ 的均值， $s(w, h)$ 、 $s(w + 1, h)$ 、 $s(w, h + 1)$ 、 $s(w + 1, h + 1)$ 分别为角点 ${P_{w, h}}$ 、 ${P_{w + 1, h}}$ 、 ${P_{w, h + 1}}$ 、 ${P_{w + 1, h + 1}}$ 的标准偏差， $X$ 和 $Y$ 分别为影像块 $B(w, h)$ 的宽和高。

(4) 根据计算出每个像素的 $m(x, y)$ 和 $s(x, y)$ 进行Wallis变换处理，公式如下

$f(x, y) = \left( {g\left(x, y\right) - m\left(x, y\right)} \right) \cdot \left({s_f}/s(x, y)\right) + {m_f}$

(9)

利用双线性内插计算每个像素的线性变换参数，可以保证相邻影像块之间的平滑性。此外，利用影像的各角点而非中心点参与计算，可以避免影像边缘的分块出现锯齿现象。

在利用Wallis匀色时，如果待处理影像与参考影像的影像内容差距较大，处理结果容易偏色或退化。区域范围内的多幅影像之间的内容往往变化较大，单一的参考影像不适合区域内多幅影像之间的匀色处理。根据相邻影像之间具有相同或相似内容的性质，可采用多次调整参考影像两两处理相邻影像的方法。由于区域范围内多幅影像间的邻近关系复杂，空间分布不规则，影像之间的处理顺序即传递路径直接影响匀色的质量。

为了解决上述问题，本文采用Voronoi图与Dijkstra算法相结合的方法确定匀色处理顺序。根据影像中心点位置构造Voronoi图，对区域范围内的影像有效的组织，便于邻近关系的查询。利用Dijkstra算法计算最短传递路径(李烁 等，2016)，对应处理的相邻影像具有较大重叠区域的相邻影像，传递次数较少，可以减少传递过程中的误差(Pan 等，2010)，提高区域内影像之间的色彩一致性。具体步骤如下：

(1) 选取初始参考影像。初始参考影像的质量会影响区域内影像的整体质量，本文选取具有最大清晰度的影像作为初始参考影像。清晰度是利用影像的平均梯度来衡量影像在纹理和细节方面表达能力的质量评价参数，其值越大表示影像越清晰，影像质量越高，计算公式为

$D = \frac{1}{{(M - 1)(N - 1)}}\sum\limits_{x = 1}^{M - 1} {\sum\limits_{y = 1}^{N - 1} {\sqrt {\frac{{\Delta _x^2 + \Delta _y^2}}{2}} } } $

(10)

${\Delta _x} = f(x + 1, y) - f(x, y)$

(11)

${\Delta _y} = f(x, y + 1) - f(x, y)$

(12)

式中， $M$ 和 $N$ 分别为影像的宽和高， $f(x, y)$ 为影像在 $(x, y)$ 处的灰度值。

(2) 根据区域内影像的中心点构建Voronoi图，如图2中虚线所示，每个中心点对应一个Voronoi多边形。根据Voronoi多边形之间是否存在公共边，即可判断对应的影像是否相邻。由于Voronoi图是以距离中心点最近为原则划分的，可以合理地判断多度重叠影像之间的邻近关系，只有距离较近的影像才会判断为相邻，对应影像的重叠度较大，可以提高匀色效果。

(3) 利用Dijkstra算法计算初始参考影像中心点到其他影像中心点的最短路径，即图2中 ${v_1}$ 到其余各点 ${v_i}$ 的最短路径，并记录 ${v_i}$ 的前一点 ${v_j}$ ，影像 $j$ 的匀色结果为影像 $i$ 的参考影像。

(4) 利用分块Wallis匀色法两两处理。为了避免多次传递导致的色彩偏移，保持区域内影像整体的色彩一致性，利用初始参考影像的均值和标准偏差对新的参考影像的均值和标准进行约束。设影像 $j$ 的匀色结果为影像 $j'$ ，在对影像 $i$ 处理时，参考均值 ${m_f}$ 和标准偏差 ${s_f}$ 的计算公式为

${m_f} = w{m_{j'}} + (1 - w){m_1}$

(13)

${s_f} = w{s_{j'}} + (1 - w){s_1}$

(14)

式中， ${m_{j'}}$ 和 ${s_{j'}}$ 分别为影像 $j'$ 的均值和标准偏差， ${m_1}$ 和 ${s_1}$ 分别为初始参考影像的均值和标准偏差， $w \in [0, 1]$ 为权值常数，一般取0.5。

分块Wallis匀色法的主要计算任务包括3个步骤：计算各影像块的均值和标准偏差、逐像素双线性插值和线性变换计算新的灰度值。分析可知：

(1) 双线性插值和线性变换需要逐像素计算，是计算量最大的部分，各像素的计算是相互独立的，非常适合GPU并行处理，一个线程执行一个像素点的计算任务。

(2) 影像块均值和标准偏差的本质均为累计求和计算，耦合度低，不能直接进行并行设计。本文采用归约求和的方式进行优化，利用共享存储器在每一个线程块中进行并行运算。

双线性插值和线性变换为重复的密集计算任务，并行化比较简单，可以直接分配给各个线程同时计算，满负荷利用计算资源，可以大幅度缩短计算时间。

线程是GPU的最小执行单元，具体执行时，按照“线程格网—线程块—线程”的层次结构进行组织，如图3所示。假设待处理影像大小为 $M \times N$ ，分配相同大小的线程格网，设置线程块大小为 $l \times k$ ，则线程块数目为 $\dfrac{M}{l} \times \dfrac{N}{k}$ ，一个线程对应一个像素。各线程同时进行运算，并将结果按索引号赋给对应的像素。

归约求和是一种缩减计算方法，基于对数步长交替两两求和，如图4所示，可以将求和的时间复杂度由 $O(N)$ ( $N$ 为数据个数)降为 $O({\log _2}N)$ 。在每一个线程块内利用共享内存进行归约求和，可以达到并行加速的目的，并且交替策略可以避免存储片冲突并保持线程块的相邻线程处于活跃状态。

针对GPU的层次结构，本文采用两遍归约求和策略。第一阶段内核执行 $n$ 个并行规约，其中 $n$ 是指线程块数，得到一个中间结果数组；第二个阶段通过调用一个线程块对这个中间数组进行规约，从而得到最终结果，如图5所示。具体步骤如下：

(1) 对输入数组落入每个线程中的数据进行求和。每个线程把它得到的累计值写入共享内存，交替因子为 $n \times m$ ，并在执行对数步长的归约前进行同步操作。

(2) 对共享内存中的值进行对数步长的归约操作。共享内存中后半部分的值被加到前半部分，即 $a[i] = a[i] + a[m/2], (0 \leqslant i < m/2)$ ，参与的线程依次减半。在此操作执行 ${\log _2}m$ 次后，共享内存中第一个线程对应的值 $a[0]$ 即为该线程块的和。共享内存的大小等于线程块的线程数量 $m$ ，并且 $m$ 必须是2的幂次。

(3) 线程块的和写入全局内存。

利用并行规约求和方法计算影像像素值之和以及像素个数，即可求出影像均值，同理可求标准偏差。

为了达到尽可能高的并行加速比，本文结合算法的自身特点，从配置划分、存储器带宽和指令吞吐量等方面进行优化。具体如下：

(1) 合理组织线程。GPU中线程块独立被调度到流式多处理器中，来自同一个线程块的线程在同一个流式多处理器中执行。为了使流式多处理器的性能达到最优，流式多处理器中线程块的数量要小于8、线程数等于1536，线程块中的线程数要小于1024 (Kirk 等，2010)。分析可知，当线程块大小设置为256或512时，流式多处理器的性能最优，可以提高并行算法的运算速度。

(2) 存储器优化。由于常量存储器的访问速度明显优于全局存储器，合理利用常量存储器代替全局存储器可以有效地减少内存流量和内存带宽(闸旋，2013)。双线性插值中的四个角点的均值和标准偏差和线性变换中的参考均值和参考标准偏差数据量不大且每个线程都会读取，将这些参数分配为常量内存，可以大大加速访问数据速度和减少内存带宽，有效地提升程序的运行效率。

针对共享内存中的归约求和，为了减少不必要的线程同步，使用线程束同步优化。由于每个线程块中的线程束是按照锁步方式执行每条指令的，当线程块中的活动线程数低于硬件线程束的大小32时，无须再调用_syncthreads()内置函数。线程束同步可以避免每个Warp中出现分支导致效率低下，减小线程的闲置，提高并行度。

(3) 指令优化。为了使指令吞吐量最大化，在满足精度的前提下，尽可能使用单精度float类型代替双精度double类型，使用硬件函数代替常规函数，用最少的指令完成相同的运算。

为了验证本文方法的有效性和高效性，本文设计两组对比实验分别从匀色质量和匀色效率两个方面进行验证。

实验硬件环境为Intel(R) Core(TM) i5-6300HQ CPU 2.30 GHz、8 GB内存、NVIDIA GeForce GTX960M显卡和2 GB显存的便携式计算机。

软件环境为Windows 7操作系统、Microsoft Visual Studio 2010开发环境、C++编程语言和Cuda 7.5并行开发包。

为了验证本文方法的匀色质量，分别利用全局Wallis匀色法、商业软件Inpho5.6以及本文方法对两组航空正射影像进行对比实验，结果如图7和图8所示，图9和图10为图8中区域1和区域2的局部放大图。第一组实验数据(数据Ⅰ)为4幅存在明显的色彩不一致和对比度不一致的彩色影像，两两重叠度约为40%，如图7(a)所示。第二组实验数据(数据Ⅱ)为5条航带的234张影像，航向重叠度约为70%—80%，旁向重叠度约为40%—50%，如图8(a)所示。为了直观对比3种方法的匀色质量，图7—图10均为几何镶嵌的结果，未进行镶嵌缝羽化处理。

可以看出，全局Walllis匀色法只能在一定程度上减小影像间的色彩差异，影像间仍存在明显的色调差异，如图7(b)和图9(b)所示，这是由于影像整体的统计信息无法代表影像局部的地物特征。并且采用单一参考影像会产生过增强或偏色现象，如图10(b)所示。Inpho能够较好地消除影像间的色彩差异，使区域影像整体色调一致，但是影像整体对比度下降并且色调与原始影像产生了一定的偏差，如图8(c)所示。同时，当影像的色调与区域整体色调差异较大时，Inpho的处理结果不理想，如图7(c)中右下角影像。相比而言，本文方法对两组影像的处理结果都比较理想，在使区域影像整体色彩和对比度一致的同时消除或减小了相邻影像重叠区域的局部差异，并且没有产生偏色现象。同时，本文自动选取了清晰度最大的影像作为参考影像，区域内影像整体清晰度较高，具有较好的目视效果。

为了对匀光影像进行定量评价，统计区域范围内相邻影像重叠区域的灰度均值之差的平均值 $\overline {\Delta m} $ 和标准偏差之差的平均值 $\overline {\Delta s} $ ，计算公式如下

$\overline {\Delta m} = \frac{{\sum {\left| {{m_{ij}} - {m_{ji}}} \right|} }}{n}, (i \ne j)$

(15)

$\overline {\Delta s} = \frac{{\sum {\left| {{s_{ij}} - {s_{ji}}} \right|} }}{n}, (i \ne j)$

(16)

式中， ${m_{ij}}$ 和 ${m_{ji}}$ 分别为影像 $i$ 和影像 $j$ 在重叠区域中灰度均值， ${s_{ij}}$ 和 ${s_{ji}}$ 分别为影像 $i$ 和影像 $j$ 在重叠区域中标准偏差， $n$ 为重叠区域的个数。

$\overline {\Delta m} $ 可以反映影像间色彩差异， $\overline {\Delta s} $ 可以反映影像间对比度差异，其值越小表示影像间差异越小，处理结果越理想。数据Ⅰ和数据Ⅱ匀色结果的统计值如表1和表2所示。

由表1可知，数据Ⅰ的原始影像的 $\overline {\Delta m} $ 和 $\overline {\Delta s} $ 都较大，说明原始影像间的色调差异和对比度差异都较大。3种方法匀色影像的 $\overline {\Delta m} $ 都小于原始影像，说明3种方法不同程度地减小了影像间的色彩差异。本文方法的 $\overline {\Delta m} $ 最小，Inpho次之，全局Wallis方法最大，说明本文方法消除色彩差异的能力最强。同理，本文方法的 $\overline {\Delta s} $ 最小，说明本文方法消除对比度不一致的能力最好。Inpho的 $\overline {\Delta m} $ 虽然较小，但是 $\overline {\Delta s} $ 要大于原始影像，说明Inpho消除影像间对比度差异的能力较差。而全局Wallis则相反，可以较好地消除影像间对比度不一致，但是消除影像间色彩差异的能力较弱。对比表2中统计结果，可得到类似于表1的结论，与目视评价结果一致。

由上可知，本文方法在目视效果评价和定量评价两方面均优于全局Wallis匀色法和Inpho，可以同时有效地消除影像间的色彩差异和对比度差异，体现了明显的优势。

为了验证本文并行加速策略的有效性，分别统计CPU串行算法和GPU并行算法处理不同尺寸影像的运行时间和加速比，如表3和图11所示。实验影像为数据Ⅱ的缩放影像，加速比为CPU串行算法运行时间与GPU并行算法运行时间的比值。统计的运行时间不包括硬盘与内存之间数据传输消耗的时间，统计均值和标准偏差的运行时间为统计影像整体均值和标准偏差以及统计分块均值和标准偏差的运行时间之和，分块数目均设置为8×8。

分析表3和图11可知：

(1) 当影像尺寸较小如512×512时，统计影像均值和标准偏差的串行算法比并行算法快。这是由于统计每一个影像块的均值和标准偏差时都存在内存与显存之间的数据传输，而GPU并行的加速时间小于数据传输消耗的时间。但是，双线性插值和线性变换的并行加速效果较好，因此，整个过程GPU并行算法的运行时间仍优于CPU串行算法的运行时间。

(2) 随着影像尺寸的增加，GPU的计算资源得到更充分的利用，GPU并行算法的优势越来越明显，加速比越来越大。

(3) 当影像尺寸增加到一定程度时，GPU并行算法的加速比趋于稳定不再明显增加，最高加速比可以达到60倍以上。

针对区域范围内多幅影像之间的色彩不一致问题，本文提出一种基于GPU的自适应分块加权Wallis并行匀色方法。采用分块策略和双线性加权方法对Wallis变换进行改进，对不同像素采用不同的线性关系进行处理，可以在局部更好地消除相邻影像间的色彩差异和对比度差异。同时，本文采用最短传递路径确定匀色处理顺序，可以有效地保证区域内影像整体的色彩一致性。此外，本文结合GPU并行技术，设计了并行加速策略，大幅度提高了匀色处理的效率。与CPU串行算法相比，最高加速比可以达到60倍以上。本文方法有效地提高了匀色处理的质量和效率，能够降低影像镶嵌的难度，提高影像镶嵌的质量，同时可以满足海量遥感影像快速处理和应急响应的需求。

但是，在试验中也发现了本文方法的不足，本文未考虑遥感影像中的特殊区域(如水域和阴影等)，下一步应当对参与计算的像素进行检测和筛选，进一步提高影像匀色的质量。另外，随着多核GPU的出现，下一步可研究如何利用多核GPU实现影像匀色，进一步提高算法的效率。