CN102129669B - 一种航空遥感影像的最小二乘区域网匀色方法 - Google Patents

Abstract

本发明公开了一种航空遥感影像的最小二乘区域网匀色方法，该方法借鉴光束法区域网空中三角测量的思想，首先统计单像正射纠正影像的相邻影像叠加区域的像素色调信息，再构建每幅影像像素色调信息补偿参数的方程组，利用最小二乘平差方法，解算每幅影像的像素色调信息补偿参数，最后根据解算得到的色调补偿参数调整每幅影像的色调信息。本发明方法能够对大测区以及测区内影像内容差异较大的影像进行色调差异消除处理，在消除色差的过程中还能够方便地指定测区内部任意位置、任意数量的影像作为色调控制影像来控制整个测区的基本色调。

Description

一种航空遥感影像的最小二乘区域网匀色方法

技术领域

本发明属于遥感影像处理领域，尤其涉及一种航空遥感影像的匀色方法。

背景技术

在测绘生产中，数字正射影像(DOM)制作的一个主要步骤是影像镶嵌，影像镶嵌就是将完成单像正射纠正的相邻影像镶嵌成标准图幅的正射影像成果，或者将完成单像正射纠正的相邻影像镶嵌成整个测区的正射影像然后输出为标准图幅。在影像镶嵌过程中不仅要处理好相邻影像间几何接边问题，还要保证镶嵌后的影像整体色调一致、色彩过渡均匀。影像的颜色处理效果直接关系到最终正射影像成果的质量，研究并解决影像镶嵌过程中颜色问题的自动处理，对于提高测绘生产效率，保证正射影像产品质量具有重要意义。

一般来说，影响影像颜色的因素主要有两个：一是单幅影像内部的亮度分布不均匀；二是相邻影像之间存在颜色差异。已有大量学者对单幅影像内部亮度分布不均匀问题的处理(也称匀光处理)进行了研究。对相邻影像之间颜色差异问题的处理(也称匀色处理)，较早提出的有基于Wallis滤波与色调基准影像相结合的方法^[1]，该方法处理影像内容比较一致的测区时效果较为理想，但对较大测区、大比例尺的航空摄影影像或低空摄影影像而言，测区内影像之间的内容往往差异较大(如图1)，此时，采用本方法进行匀色处理则无法取得理想效果。潘俊^[2]等提出根据影像间的邻接关系、采用最短路径算法确定影像间的传递路径、并利用线性关系对相邻影像进行匀色处理，由于航空影像内容的复杂性，传递路径上的影像不一定严格满足线性关系，所以在测区较大的情况下容易出现传递路径上的误差累积，从而影响影像整体的匀色处理效果。

文中文中涉及的参考文献如下：

[1]李治江.彩色影像色调重建的理论与实践[D].武汉：武汉大学：2005.4.

[2]潘俊.自动化的航空影像色彩一致性处理及接缝线网络生成方法研究[D].武汉：武汉大学，2008.11.

发明内容

针对现有技术存在的不足，本发明的目的是提供一种航空遥感影像的最小二乘区域网匀色方法，该方法能有效改善大测区、大比例尺的航空遥感影像的自动匀色效果。

为实现上述目的，本发明采用如下技术方案：

一种航空遥感影像的最小二乘区域网匀色方法，包括以下步骤：

步骤1，确认测区内所有待镶嵌的单像正射纠正影像在镶嵌过程中的接边关系，并获得q对接边影像对；

步骤2，对步骤1获得的q对接边影像对，按照地理坐标对每对接边影像对中的两张单像正射纠正影像进行叠加，叠加后排除叠加区域内单像正射纠正影像的背景像素和纹理差异较大的像素，然后，分别统计q对接边影像对的叠加区域内未被排除的像素数和包括像素灰度均值和像素灰度方差的像素色调信息；

步骤3，计算测区内所有待镶嵌的单像正射纠正影像的像素色调信息补偿参数，该步骤进一步包括以下子步骤：

3-1、利用步骤2统计得到的像素色调信息构建矩阵方程：A_BX＝L_B，

其中，

X＝[x₀ x₁ x₂ ... x_p-1]^T，x₀、x₁、x₂、x₃、...x_p-1分别为测区内所有待镶嵌的单像正射纠正影像的像素色调信息补偿参数，p为所有待镶嵌影像的张数；

A_B为q×p矩阵，且 $A_B=\left[\begin{array}{ccccc}1& -1& 0& ...& \\ 0& 1& -1& 0& ...\\ 0& ...& 1& -1& 0\\ ...& & & & \end{array}\right];$

$L_B={\left[\begin{array}{ccccc}{a}_0^1-a_0^0& a_1^1-a_1^0& a_2^1-a_2^0& ...& a_{q-1}^1-a_{q-1}^0\end{array}\right]}^T,$

为第i对接边影像对中的第一张影像x在叠加区域内的像素灰度均值或像素灰度方差，为第i对接边影像对中的第二张影像y在叠加区域内的像素灰度均值或像素灰度方差，且

和

同为像素灰度均值或像素灰度方差；

3-2、将p张单像正射纠正影像中色调较好的t张影像设置为色调控制影像，并构建矩阵方程：A_CX＝L_C，

其中，

X＝[x₀ x₁ x₂ ... x_p-1]^T，x₀、x₁、x₂、x₃、...x_p-1分别为p张单像正射纠正影像的像素色调信息补偿参数；

A_c[i，j]＝c_i，j，i∈[0，t-1]，j∈[0，p-1]，i表示在t张色调控制影像中的第i张影像，当j＝i时，c_i，j＝1，其它情况下，c_i，j＝0；

L_C为t行零常数的列向量，L_C＝[0 0…0]^T；

3-3、令 $A=\left[\begin{array}{c}{A}_B\\ A_C\end{array}\right],$ $L=\left[\begin{array}{c}{L}_B\\ L_C\end{array}\right],$ 构建矩阵方程AX＝L，则p张单像正射纠正影像的像素色调信息补偿参数X＝(A^TA)^-1A^TL；

3-4、构建权矩阵 $P=\left[\begin{array}{cc}\begin{array}{cccc}{p}_0& 0& \·\·\·& 0\\ 0& p_1& 0& 0\\ 0& \·\·\·& \·\·\·& 0\\ 0& \·\·\·& 0& p_{q-1}\end{array}& 0\\ 0& E\end{array}\right]$ 修正步骤3-3得到的像素色调信息补偿参数，则p张单像正射纠正影像的色调补偿参数X＝(A^TPA)^-1A^TPL，

其中，

r∈[0，q-1]，t_r为叠加区域内的有效像素数，E为(p-q)阶单位矩阵；

步骤4，根据步骤3所得的像素色调补偿参数，调整每张单向正射纠正影像的像素色调信息。

上述步骤1具体如下：

基于测区内所有待镶嵌的p张单像正射纠正影像的摄站坐标构建Voronoi图，利用Voronoi图的邻接关系来确认镶嵌过程中p张单像正射纠正影像之间的接边关系。

上述步骤2中排除叠加区域内纹理差异较大的像素采用如下的相关系数法：

利用预设尺寸的窗口对叠加区域进行逐个窗口内相关系数ρ的计算，如果窗口位置的相关系数值小于预先设定的阈值，则排除掉两幅影像上该窗口位置内的像素，相关系数ρ的计算公式如下：

$\mathrm{\ρ}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(f_{i,j}^1\·f_{i,j}^2)-\frac{1}{\mathrm{mn}}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^1\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^2\right)}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(f_{i,j}^1\right)}^2-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^1\right)}^2[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(f_{i,j}^2\right)}^2-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^2\right)}^2]}}$

其中，

m×n为窗口尺寸，

分别为叠加区域内的两张影像在窗口内(i，j)处的像素灰度值。

上述步骤4中是采用Wallis滤波算子调整p张单像正射纠正影像的像素色调信息。。

与现有技术相比，本发明具有以下优点和有益效果：

本发明方法能够对大测区以及测区内影像内容差异较大的影像进行色调差异消除处理，在消除色差的过程中还能够方便地指定测区内部任意位置、任意数量的影像作为色调控制影像来控制整个测区的基本色调。图7为利用本发明方法和常规方法对航空遥感影像进行匀色处理后的正射镶嵌结果对比，从图中可以看出，与常规匀色方法相比，采用本发明方法匀色处理后的镶嵌结果消除了影像间的色调差异。

附图说明

图1为165张航空遥感影像的正射镶嵌图；

图2为二维平面上一系列输入点的Voronoi图；

图3为原始航空遥感影像与其对应的单像正射纠正影像，(a)为原始航空遥感影像，(b)为图(a)的单像正射纠正影像；

图4为利用相关系数法对单像正射纠正影像进行纹理差异较大的像素排除结果，(a)为左影像正射纠正后的局部区域，(b)为右影像正射纠正后的局部区域，(c)为左右正射影像叠加图，(d)为纹理差异较大的像素的排除结果；

图5为利用本发明方法选用不同的色调控制影像进行匀色处理的结果，(a)为不同位置的色调控制方案，(b)为选用黑框位置的色调控制影像的匀色处理结果，(c)为选用白框位置的色调控制影像的匀色处理结果；

图6为利用本发明方法对253张原始洞窟壁画影像匀色处理前后的拼接效果，(a)为253张原始洞窟壁画影像的正射镶嵌效果，(b)为利用本发明方法进行匀色处理后的正射镶嵌效果；

图7为利用本发明方法和常规方法处理后的正射镶嵌结果对比，(a)为原始影像的正射镶嵌效果，(b)为利用常规方法匀色处理后的正射镶嵌效果，(c)为采用本发明方法匀色处理后的正射镶嵌效果。

具体实施方式

本发明采用基于最小二乘的区域网匀色方法解决航空影像间的颜色差异问题，该方法借鉴光束法区域网空中三角测量的思想，首先统计单像正射纠正影像的相邻影像叠加区域的像素色调信息，再构建每幅影像像素色调信息补偿参数的方程组，利用最小二乘平差方法，解算每幅影像的像素色调信息补偿参数，最后根据解算得到的色调补偿参数调整每幅影像的色调信息。

为了更好地理解本发明的技术方案，下面将结合附图对本发明做进一步的详细说明，本发明的具体步骤如下：

步骤1，

对测区内所有待镶嵌的p张单像正射纠正影像，构建基于它们摄站坐标的Voronoi图，利用Voronoi图的邻接关系确定正射影像镶嵌过程中单像正射纠正影像两两之间是否存在接边，并获得q对接边影像对。

Voronoi图是计算几何学中的一个基础数据结构，它的输入是分布于二维平面上的一系列点，输出是对应于每个点的一系列多边形区域，输出的多边形区域应满足多边形区域内部的任意位置的点到与多边形区域所对应的输入点的几何距离最短，具体如图2所示。在正射影像镶嵌过程中，往往需要选取靠近摄影镜头中心位置的影像拼接到大图象中，所以可以利用Voronoi图的邻接关系确定哪些影像存在接边关系。

步骤2，

本步骤包括以下两个子步骤：

2-1、对步骤1获得的q对接边影像对，按照地理坐标对每对接边影像对中的两张单像正射纠正影像进行叠加，叠加后排除叠加区域内单像正射纠正影像的背景像素和纹理差异较大的像素，

叠加区域内单像正射纠正影像的背景像素，即原始航空影像在正射纠正过程中边缘不能覆盖的区域中的像素(如图3(b)所示的黑色区域)。叠加区域内单像正射纠正影像的背景像素的排除可采用本领域公知的技术手段来实现，本实施例中排除叠加区域内单像正射纠正影像的背景像素采用如下方法：对叠加区域内的两张影像，若其中一张影像的像素灰度值与背景像素灰度值相同，则排除该像素，同时排除另一张影像中的对应像素。

单像正射纠正影像上纹理差异较大的叠加区域，比如由于投影差造成的房屋屋顶位置的不同(如图4(a)、(b)所示)。对存在纹理差异的重叠区域，可利用相关系数法进行鉴别排除，本实施例中采用的具体方法如下：将重叠区域按照一定的窗口尺寸(比如5×5)进行逐个窗口内相关系数ρ的计算，如果两幅影像上该窗口位置的相关系数ρ小于一定的阈值(比如0.7)，则将两幅影像上该窗口位置内的像素排除。相关系数ρ的计算公式为：

$\mathrm{\ρ}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(f_{i,j}^1\·f_{i,j}^2)-\frac{1}{\mathrm{mn}}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^1\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^2\right)}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(f_{i,j}^1\right)}^2-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^1\right)}^2[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(f_{i,j}^2\right)}^2-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^2\right)}^2]}}---\left(1\right)$

其中，

m×n为窗口尺寸，

分别为叠加区域内的两张影像在窗口内(i，j)处的像素灰度值。

2-2、统计q对接边影像对的叠加区域内未被排除的像素数t_r(即有效像素数)和像素色调信息，所述的像素色调信息包括像素灰度均值和像素灰度方差。

对第r对接边影像对中的第一张影像在叠加区域范围内的有效像素，利用下列公式分别计算其像素灰度均值

和像素灰度方差

$m_r^0=\frac{1}{t_r}\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{h}_i\·i---\left(2\right)$

$v_r^0=\sqrt{\frac{1}{t_r}\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{(i-m)}^2\·h_i}---\left(3\right)$

其中，

t_r为第r对接边影像对在叠加区域内的有效像素数；

k为图像的色阶数，对每个通道为8Bit的图像，k＝255；

h_i为叠加区域内第一张影像中像素灰度值为i且未被排除的像素个数，i∈[0，k]，h_i可由公式

表示，w、h分别为第r对单像正射纠正影像对的叠加区域的像素宽和高；

0≤x≤w，0≤y≤h；g⁰(x，y)表示第一张图像重叠区域内(x，y)坐标处的像素灰度值。

采用相同的方法计算得到第r对接边影像对的叠加区域内第二张影像的像素灰度均值和像素灰度方差

重复本步骤，逐一统计所有接边影像对在叠加区域内像素灰度均值

与像素灰度方差

其中，r∈[0，q-1]，s∈[0，1]。

步骤3，计算p张单像正射纠正影像的像素色调信息补偿参数，影像的像素色调信息包括像素灰度均值和像素灰度方差，则像素色调信息补偿参数也包括像素灰度均值补偿参数和像素灰度方差补偿参数，下面将以像素灰度均值补偿参数的计算为例来进行具体说明：

3-1、利用影像叠加区域的像素灰度均值信息构建矩阵方程A_BX＝L_B，矩阵方程A_BX＝L_B的构建具体如下：

首先，假设第r对接边影像对是由第x、y两张影像构成，根据影像镶嵌拼接过程中叠加区域应该具有相同色调的基本要求，针对第x、y张影像的叠加区域的像素灰度均值，得到如下方程：

$m_r^0+m_x=m_r^1+m_y---\left(4\right)$

其中，

x∈[0，p-1]，y∈[0，p-1]，且x≠y；

分别为第x张影像和第y张影像在叠加区域内的像素灰度均值；

m_x、m_y分别表示第x、y张影像的像素灰度均值的补偿参数。

将上述方程(4)写为： $m_x-m_y=m_r^1-m_r^0---\left(5\right)$

方程(5)等价为： $\left[\begin{array}{cccc}{b}_{r,0}& b_{r,1}& \·\·\·& b_{r,p-1}\end{array}\right]{\left[\begin{array}{cccc}{m}_0& m_1& \·\·\·& m_{p-1}\end{array}\right]}^T=[m_r^1-m_r^0]---\left(6\right)$

其中，

[b_r，0 b_r，1 … b_r，p-1]是由b_r，i组成的一维系数矩阵，i∈[0，p-1]；当i＝x时，b_r，i＝1；当i＝y时，b_r，i＝-1；其它情况下，b_r，i＝0；

然后，若按矩阵方程(6)写出q对接边影像对中所有单像正射纠正影像的矩阵方程，并将所得矩阵方程联立，即可构建矩阵方程：

$\left[\begin{array}{cccc}{b}_{\mathrm{0,0}}& b_{\mathrm{0,1}}& \·\·\·& b_{0,p-1}\\ b_{\mathrm{1,0}}& b_{\mathrm{1,1}}& \·\·\·& b_{1,p-1}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ b_{q-\mathrm{1,0}}& b_{q-\mathrm{1,1}}& \·\·\·& b_{q-1,p-1}\end{array}\right]\left[\begin{array}{c}{m}_0\\ m_1\\ \·\·\·\\ m_{p-1}\end{array}\right]=\left[\begin{array}{c}{m}_0^1-m_0^0\\ m_1^1-m_1^0\\ \·\·\·\\ m_{q-1}^1-m_{q-1}^0\end{array}\right]---\left(7\right)$

上述矩阵方程(7)可写为A_BX＝L_B形式， $A_B=\left[\begin{array}{cccc}{b}_{\mathrm{0,0}}& b_{\mathrm{0,1}}& \·\·\·& b_{0,p-1}\\ b_{\mathrm{1,0}}& b_{\mathrm{1,1}}& \·\·\·& b_{1,p-1}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ b_{q-\mathrm{1,0}}& b_{q-\mathrm{1,1}}& \·\·\·& b_{q-1,p-1}\end{array}\right],$ X＝[m₀ m₁ … m_p-1]^T， $L_B={\left[\begin{array}{cccc}{m}_0^1-m_0^0& m_1^1-m_1^0& \·\·\·& m_{q-1}^1-m_{q-1}^0\end{array}\right]}^T,$

其中，

X＝[m₀ m₁ … m_p-1]^T，m_i为第i张单像正射纠正影像的像素灰度均值的补偿参数，i∈[0，p-1]；

$A_B=\left[\begin{array}{cccc}{b}_{\mathrm{0,0}}& b_{\mathrm{0,1}}& \·\·\·& b_{0,p-1}\\ b_{\mathrm{1,0}}& b_{\mathrm{1,1}}& \·\·\·& b_{1,p-1}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ b_{q-\mathrm{1,0}}& b_{q-\mathrm{1,1}}& \·\·\·& b_{q-1,p-1}\end{array}\right],$ A_B可看作由b_i，j构成的矩阵，i∈[0，q-1]，j∈[0，p-1]，b_i，j中的i表示q对接边影像对中的第i对影像对，j表示p张单像正射纠正影像中的第j张影像。假设，第i对影像对由第x、y张影像构成，x∈[0，p-1]，y∈[0，p-1]，则当j＝x时，b_i，j＝1；当j＝y时，b_i，j＝-1，其它情况下，b_i，j＝0；

$L_B={\left[\begin{array}{cccc}{m}_0^1-m_0^0& m_1^1-m_1^0& \·\·\·& m_{q-1}^1-m_{q-1}^0\end{array}\right]}^T,$

为第i对接边影像对中的第一张影像在叠加区域内的像素灰度均值，

为第i对接边影像对中的第二张影像在叠加区域内的像素灰度均值。

3-2、将p张单像正射纠正影像中色调较好的t张影像设为色调控制影像，并构建矩阵方程：A_CX＝L_C；

将测区中一些色调较好的影像设为色调控制影像，使整个测区的最终色调与控制影像基本一致，色调控制影像的像素灰度均值补偿参数应该较小或者为0。假定第j张影像为色调控制影像(j∈[0，p-1])，设定其像素灰度均值补偿参数为0，可用如下方程表示：

[c_j，0 c_j，1 … c_j，p-1][m₀ m₁ … m_p-1]^T＝[0]    (8)

其中，[c_j，0 c_j，1 … c_j，p-1]为由c_j，i组成的一维系数矩阵，i∈[0，p-1]，当i＝j时，c_j，i＝1；其它情况下，c_j，i＝0。

假设选定t张影像设为色调控制影像，则可列矩阵方程：

$\left[\begin{array}{cccc}{c}_{\mathrm{0,0}}& c_{\mathrm{0,1}}& \·\·\·& c_{0,p-1}\\ c_{\mathrm{1,0}}& c_{\mathrm{1,1}}& \·\·\·& c_{1,p-1}\\ ...& ...& ...& ...\\ c_{t-\mathrm{1,0}}& c_{t-\mathrm{1,1}}& \·\·\·& c_{t-1,p-1}\end{array}\right]{\left[\begin{array}{cccc}{m}_0& m_1& \·\·\·& m_{p-1}\end{array}\right]}^T=\left[\begin{array}{c}0\\ 0\\ \·\·\·\\ 0\end{array}\right]---\left(9\right)$

上述矩阵方程(9)可写为A_CX＝L_C形式，

X＝[m₀ m₁ … m_p-1]^T，m₀、m₁ ... m_-p-1分别为p张单像正射纠正影像的像素灰度均值补偿参数；

$A_C=\left[\begin{array}{cccc}{c}_{\mathrm{0,0}}& c_{\mathrm{0,1}}& \·\·\·& c_{0,p-1}\\ c_{\mathrm{1,0}}& c_{\mathrm{1,1}}& \·\·\·& c_{1,p-1}\\ ...& ...& ...& ...\\ c_{t-\mathrm{1,0}}& c_{t-\mathrm{1,1}}& \·\·\·& c_{t-1,p-1}\end{array}\right],$ A_C可看作由c_i，j(i∈[0，t-1]，j∈[0，p-1])构成的矩阵，c_i，j中的i表示在t张色调控制影像中的第i张影像(其影像序号为z，有z∈[0，p-1])。当j＝z时，c_i，j＝1；其它情况下，c_i，j＝0；

L_C＝[0 0 … 0]^T。

3-3、将步骤3-1和步骤3-2得到的矩阵方程联立，得到 $\left[\begin{array}{c}{A}_B\\ A_C\end{array}\right]X=\left[\begin{array}{c}{L}_B\\ L_C\end{array}\right]---\left(10\right),$ 令 $A=\left[\begin{array}{c}{A}_B\\ A_C\end{array}\right],$ $L=\left[\begin{array}{c}{L}_B\\ L_C\end{array}\right],$ 则方程(10)可简化成AX＝L，则p张单像正射纠正影像的像素灰度均值补偿参数X＝(A^TA)^-1A^TL；

3-4、引入步骤2-2中统计的叠加区域的有效像素个数t_r，构建权矩阵P：

$P=\left[\begin{array}{cc}\begin{array}{cccc}{a}_0& 0& \·\·\·& 0\\ 0& a_1& 0& 0\\ 0& \·\·\·& \·\·\·& 0\\ 0& \·\·\·& 0& a_{q-1}\end{array}& 0\\ 0& E\end{array}\right]$

其中：

r∈[0，q-1]，t_r为叠加区域内的有效像素数，E为(p-q)阶单位矩阵。

采用权矩阵P对步骤3-3得到像素灰度均值补偿参数X＝(A^TA)^-1A^TL进行修正，即影像的最终像素灰度均值补偿参数为X＝(A^TPA)^-1A^TPL，即可得到每张影像的像素灰度均值补偿参数m_i。

同理，将矩阵X＝[m₀ m₁ … m_p-1]^T中的m_i改成影像的像素方差补偿参数v_i。i∈[0，p-1]；将矩阵 $L_B={\left[\begin{array}{cccc}{m}_0^1-m_0^0& m_1^1-m_1^0& \·\·\·& m_{q-1}^1-m_{q-1}^0\end{array}\right]}^T$ 中的

改成

i∈[0，q-1]，其中，

为第i对接边影像对中的第一张影像在叠加区域内的像素方差，

为第i对接边影像对中的第二张影像在叠加区域内的像素方差。然后，采用步骤3-1～3-4解算出每张影像的像素灰度方差补偿参数v_i。

步骤4，

根据步骤3所得的影像像素灰度均值补偿参数和像素灰度方差补偿参数，并利用Wallis等滤波算子调整p张单向正射纠正影像的像素色调信息。

如果单像正射纠正影像为彩色图像，则首先将单像正射纠正影像按R、G、B三个波段分解为灰度图像，再采用上述步骤分别对相同波段的图像进行匀色处理，并将匀色处理完毕后的灰度图像按照R、G、B波段重新组合成彩色图像。

本发明的方法的匀色效果可见图5所示。图7为利用本发明方法和常规方法匀色处理后的正射镶嵌结果对比，从图中可以看出，跟常规匀色方法相比，采用本发明方法匀色处理后的镶嵌结果消除了影像间的色调差异。

Claims (4)

1.一种航空遥感影像的最小二乘区域网匀色方法，其特征在于，依次包括以下步骤：

步骤1，确认测区内所有待镶嵌的单像正射纠正影像在镶嵌过程中的接边关系，并获得q对接边影像对；

步骤2，对步骤1获得的q对接边影像对，按照地理坐标对每对接边影像对中的两张单像正射纠正影像进行叠加，叠加后排除叠加区域内单像正射纠正影像的背景像素和纹理差异较大的像素，然后，分别统计q对接边影像对的叠加区域内未被排除的像素数和包括像素灰度均值和像素灰度方差的像素色调信息；

步骤3，计算测区内所有待镶嵌的单像正射纠正影像的像素色调信息补偿参数，该步骤进一步包括以下子步骤：

3-1、利用步骤2统计得到的像素色调信息构建矩阵方程：A_BX＝L_B，

其中，

X＝[x₀ x₁ x₂ ... x_p-1]^T，x₀、x₁、x₂、x₃、...x_p-1分别为测区内所有待镶嵌的单像正射纠正影像的像素色调信息补偿参数，p为所有待镶嵌影像的张数；

A_B为q×p矩阵，且 $A_B=\left[\begin{array}{ccccc}1& -1& 0& ...& \\ 0& 1& -1& 0& ...\\ 0& ...& 1& -1& 0\\ ...& & & & \end{array}\right];$

$L_B={\left[\begin{array}{ccccc}{a}_0^1-a_0^0& a_1^1-a_1^0& a_2^1-a_2^0& ...& a_{q-1}^1-a_{q-1}^0\end{array}\right]}^T,$

为第i对接边影像对中的第一张影像x在叠加区域内的像素灰度均值或像素灰度方差，为第i对接边影像对中的第二张影像y在叠加区域内的像素灰度均值或像素灰度方差，且

和

同为像素灰度均值或像素灰度方差；

3-2、将p张单像正射纠正影像中色调较好的t张影像设置为色调控制影像，并构建矩阵方程：A_CX＝L_C，

其中，

X＝[x₀ x₁ x₂ ... x_p-1]^T，x₀、x₁、x₂、x₃、...x_p-1分别为p张单像正射纠正影像的像素色调信息补偿参数；

A_c[i，j]＝c_i，j，i∈[0，t-1]，j∈[0，p-1]，i表示在t张色调控制影像中的第i张影像，当j＝i时，c_i，j＝1，其它情况下，c_i，j＝0；

L_C为t行零常数的列向量，L_C＝[0 0…0]^T；

3-3、令 $A=\left[\begin{array}{c}{A}_B\\ A_C\end{array}\right],$ $L=\left[\begin{array}{c}{L}_B\\ L_C\end{array}\right],$ 构建矩阵方程AX＝L，则p张单像正射纠正影像的像素色调信息补偿参数X＝(A^TA)^-1A^TL；

3-4、构建权矩阵 $P=\left[\begin{array}{cc}\begin{array}{cccc}{p}_0& 0& \·\·\·& 0\\ 0& p_1& 0& 0\\ 0& \·\·\·& \·\·\·& 0\\ 0& \·\·\·& 0& p_{q-1}\end{array}& 0\\ 0& E\end{array}\right]$ 修正步骤3-3得到的像素色调信息补偿参数，则p张单像正射纠正影像的色调补偿参数X＝(A^TPA)^-1A^TPL，

其中，

r∈[0，q-1]，t_r为叠加区域内的有效像素数，E为(p-q)阶单位矩阵；

步骤4，根据步骤3所得的像素色调补偿参数，调整每张单向正射纠正影像的像素色调信息。

2.根据权利要求1所述的航空遥感影像的最小二乘区域网匀色方法，其特征在于：

所述步骤1具体如下：

基于测区内所有待镶嵌的p张单像正射纠正影像的摄站坐标构建Voronoi图，利用Voronoi图的邻接关系来确认镶嵌过程中p张单像正射纠正影像之间的接边关系。

3.根据权利要求1或2所述的航空遥感影像的最小二乘区域网匀色方法，其特征在于：

所述步骤2中排除叠加区域内纹理差异较大的像素采用如下的相关系数法：

利用预设尺寸的窗口对叠加区域进行逐个窗口内相关系数ρ的计算，如果窗口位置的相关系数值小于预先设定的阈值，则排除掉两幅影像上该窗口位置内的像素，相关系数ρ的计算公式如下：

$\mathrm{\ρ}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(f_{i,j}^1\·f_{i,j}^2)-\frac{1}{\mathrm{mn}}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^1\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^2\right)}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(f_{i,j}^1\right)}^2-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^1\right)}^2[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(f_{i,j}^2\right)}^2-\frac{1}{\mathrm{mn}}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{f}_{i,j}^2\right)}^2]}}$

其中，

m×n为窗口尺寸，

分别为叠加区域内的两张影像在窗口内(i，j)处的像素灰度值。

4.根据权利要求1或2所述的航空遥感影像的最小二乘区域网匀色方法，其特征在于：

所述步骤4中是采用Wallis滤波算子调整p张单像正射纠正影像的像素色调信息。

Images