CN103679675A - 一种面向水质定量遥感应用的遥感影像融合方法 - Google Patents

Abstract

一种面向水质定量遥感应用的遥感影像融合方法，最大限度保持水域提取精度以及水域光谱特性，结合PCA融合、SSVR融合设计了针对水质定量遥感应用的两级融合方法，融合过程引入决策级面向对象的地物分类解译方法，处理结果具备保持精准水体对象像素轮廓和水域光谱特性的特点。大量的实验结果表明。本发明得到水域解译精度达到90%，水域定量反演结果与直接利用高光谱/多光谱影像反演结果类似。

Description

一种面向水质定量遥感应用的遥感影像融合方法

技术领域

本发明涉及一种面向水质定量遥感应用遥感影像融合方法，属于卫星影像处理技术领域，用于水利遥感应用中多光谱/高光谱遥感影像和高分全色影像的融合。

背景技术

现代遥感技术为对地观测提供了多空间、多光谱、多时相分辨率的海量遥感影像数据，广泛用于地形测绘与地图更新、土地利用、城市识别、农业与森林调查、地质与洪涝灾害以及军事领域等。

多源遥感影像所提供的信息具有冗余性、互补性和合作性，冗余性表示它们对环境或目标的表示、描述或解译结果相同；互补性是指信息来自不同的自由度且相互独立；合作性是不同传感器在观测和处理信息时对其它信息有依赖关系。多源遥感影像数据融合是汇集多源遥感影像的最有效途径之一，通过同一地区不同数据源间的信息互补为多源遥感影像的处理、分析和应用提供最有效的应用。

通常，多光谱/高光谱影像光谱分辨率高、光谱信息丰富，但是空间分辨率低；高分全色影像具有较高的空间分辨率但缺乏光谱信息。如何将高光谱/多光谱影像和高分辨率全色影像结合起来开展应用，以实现优势互补和冗余控制，成为遥感要解决的关键问题，也是未来遥感面临的重要问题。

影像融合主要指在统一的地理坐标系中，将对同一目标检测的多幅遥感影像数据进行空间配准，再采用一定的算法，生成一幅新的、更能有效表达该目标的影像信息。通过融合能够消除多传感器之间的冗余性、增加数据间的互补性，从而提高影像的空间分辨率、清晰度等，以便能够最有效地利用各种遥感影像数据，获得更加有用的信息。

遥感影像融合从融合层次分为三级：像元级融合、特征级融合、决策级融合。像素级融合是在两幅或多幅影像实现空间配准的基础上，将各影像的像元的物理量加权求和，该值对应为同一坐标上的新影像的像元值，其作用是增加影像中的有用信息成分，以便改善如分割与特征提取等处理效果，在该层次的融合没有信息损失，因此这种处理提供了一种最优决策和识别性能；特征级融合是对于不同影像进行特征提取，按照各自影像相同类型的特征进行融合处理，其使得能够以高的置信度来提取有用的影像特征；决策级融合是高水平的融合，主要是先按照应用的要求对影像进行分类，确定各个类别中的特征影像，然后进行融合处理。

像元级融合方法处理信息量大、费时、实时性差，是最低级的融合技术，但其是信息损失最少的技术，精度较高，但决策级融合普遍不够成熟，人工干预高。对于以水质定量遥感应用目的而言，一方面需要保持尽可能真实的光谱信息满足定量遥感应用，另一方面还需保持获取水体的精确边缘为后续水利应用做好基础工作。因此，针对具体行业应用和定量遥感的融合方法采用多级融合方式，一方面重点保证了兴趣对象的光谱特性，另一方面可在融合时对兴趣对象的兴趣特征进行输出。

自80年代以来，我国在洪涝干旱灾害的监测与评估，水土流失调查与动态监测，水利工程前期规划，大型水库工程的地质调查，生态环境及水资源、水污染调查，干旱沙漠区的水资源调查以及有效灌溉面积调查等水利遥感应用方面开展了大量的应用与研究。从应用需求分析，水域提取精度和水域光谱信息是水利遥感应用需要重点关注的问题。

如何针对水利遥感应用需求，克服单一融合方法使用范围都比较狭窄的缺陷，建立统一的数据融合模型，保持精准水体对象像素轮廓和最大限度的保持水域的光谱特性是重点解决的问题。目前尚未有相关的文献报道。

发明内容

本发明的技术解决问题是：克服现有遥感影像融合技术的不足，提供一种面向遥感水质定量遥感的两级遥感影像融合方法，在精确保持水体对象像素轮廓的同时，最大限度的保持了高光谱/多光谱影像的光谱特征。

本发明的技术解决方案是：一种面向水质定量遥感应用的遥感影像融合方法，实现步骤如下：

一种面向水质定量遥感应用遥感影像融合方法，包含以下步骤：

（1）基于高分全色影像I₁的地理信息，对高光谱/多光谱遥感影像I₂进行配准和重采样，并对两幅影像的地理位置重叠区域进行裁剪，得到具有同一地理区域且相同大小的高分全色影像I₁和高光谱/多光谱影像I₂，两幅影像分辨率一致，都为高分全色影像的分辨率；所述高分影像为高空间分辨率遥感影像，高光谱影像为影像波段数一般多于100个的遥感影像，多光谱影像为影像波段数一般多于3个少于100个的遥感影像;

（2）计算高光谱/多光谱影像I₂的波段相关矩阵和矩阵的特征向量，进行主成分分析，得到主成分分量影像I_PCA；

（3）将高分全色影像I₁和主成分影像I_PCA的第一波段影像进行直方图匹配，并用直方图匹配的结果I′₁替换I_PCA的第一主成分，即替换掉I_PCA的第一波段为I′₁，得到替换后的主成分结果；

（4）将替换后的主成分结果进行主成分分析逆变换，得到保持高光谱/多光谱影像光谱信息和全色影像空间信息的初步融合影像I_M；

（5）将初步融合影像I_M进行均值漂移聚类，提取初步融合影像I_M的边缘影像I_edge，边缘影像为二值影像，其中边缘像素标记为1，非边缘像素标记为0；

（6）针对获取的边缘影像I_edge中可能存在的边缘轮廓不完全闭合的情况，采用3×3的圆形结构元素对地物边缘影像I_edge进行闭运算处理，得到边缘完全闭合的闭合边缘影像

（7）在闭合边缘影像

内逐一对每一个非边缘像素进行不同的标记N，N为取值大于等于1的整数，将这些像素作为种子点进行区域生长，在区域生长中已被其他像素区域生长过的像素不作为种子点；全图遍历完成后，每一个闭合边缘内部都被唯一标记给标识，不同闭合边缘内部像素的标识不同，边缘像素标识为0，得到标识影像g；

（8）基于标识影像g对初步融合影像I_M进行分水岭分割，得到同质影像对象；针对同质影像对象细碎的情况，计算影像对象及其邻接影像对象的光谱异质性指数和形状异质性指数，作为衡量影像对象性质相似程度，对细碎影像对象性质相似的进行反复迭代合并，当无法再合并时，此时同质影像对象合并得到具有明确地理意义的地物影像对象；

（9）对影像对象合并后的影像计算合并后影像对象的光谱特征、形状特征和纹理特征，基于人工解译信息训练粗判定贝叶斯分类器，并保存训练结果，相同类型影像只需要训练一次，对水体对象和非水体对象进行粗判定，得到判定后的水体对象和非水体对象的影像光谱特征；（粗判定贝叶斯分类器指判定对象为影像对象的贝叶斯分类器，影像对象包括水体对象和非水体对象等）

（10）对判定的水体对象和非水体对象的光谱特征作为训练对象，训练精分类贝叶斯分类器，对水体对象中每个像素进行精判定，并输出水域影像I_water；所述精分类贝叶斯分类器指分类对象为影像像素的贝叶斯分类器，即此分类器对每个像素进行分类；

（11）输入高光谱/多光谱影像每个波段宽度和全色影像波段宽度，计算水域影像每个像素的能量值，对水域影像I_water的水体像素进行像素级的SSVR融合，保证水域像素的光谱特性，形成水域融合影像I_m-water；SSVR为Simplified Synthetic Variable Ratio的缩写，即简易合成变量比；

（12）基于初步融合影像I_M和水域融合影像I_m-water进行镶嵌处理，得到最终的融合影像。

所述步骤（5）和（6）具体实现如下：

（2.1）将全图进行遍历，将未进行分割的像素选择为样本点的特征，设{x_i|i=1,2,...,N}为图像像素灰度值集合，其中x_i为图像第i个位置灰度值，N为图像像素个数，其分割步骤如下：

（2.1.1）选择存储的第一个像素灰度值为c₁，此时令聚类中心集合C={c₁}；

（2.1.2）计算x_i与c_j∈C的欧式距离，c_j为第j个聚类中心，如果||x_i-c_j||<T，T为距离阈值，则把像素分到S_j类中，S_j是第j类的灰度值集合，继续下一个像素的计算；

（2.1.3）如果||x_i-c_j||>T，则把x_i定义为新的聚类中心，继续下一个像素的计算；

（2.1.4）分类结束时，更新聚类中心

其中聚类中心的集合为C={c_j|j=1,2,...,m}，m为聚类中心个数，n_j为第j个聚类的像素个数；

（2.1.5）对c_j∈C，根据下式开始迭代，其中令

最后将迭代后的收敛点记为

$c_j^{(k+1)}=\frac{\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{n}_p{c}_{p}g\left({\left|\right|\frac{c_j^{\left(k\right)}-c_p}{h}\left|\right|}^2\right)}{\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{n}_{p}g\left({\left|\right|\frac{c_j^{\left(k\right)}-c_p}{h}\left|\right|}^2\right)}$

式中，h为用户指定的参数，

为第k次迭代第j个聚类中心的值，g(·)为用户指定的函数；

（2.1.6）最后对聚类收敛点进行合并，

(1≤a,b≤m,a≠b)，如果

则将a，b两类合并为一类，其中r为合并阈值；

（2.2）在（2.1）的过程之后，得到分割后的边缘线影像I_edge，对于整体边缘线采用闭运算处理。闭运算处理的过程如下式所示：

${\tilde{I}}_{\mathrm{edge}}=I_{\mathrm{edge}}\&CenterDot;\left[\begin{array}{ccc}0& 1& 0\\ 1& 1& 1\\ 0& 1& 0\end{array}\right].$

所述步骤（9）和（10）具体实现如下：

(3.1)基于相同影像只需要进行一次水体对象特征训练，基于贝叶斯分类器对水体影像对象进行粗提取。水体对象特征包括光谱特征、形状特征、纹理特征，纹理特征包括灰度-梯度共生矩阵提取结果，其中有能量、相关、梯度平均、混合熵、惯性矩和逆差矩；和Gabor小波滤波提取结果，包括Gabor小波均值和Gabor小波方差；

（3.2）基于获取的水体影像对象和非水体对象的光谱特征，对贝叶斯分类器进行训练，对水体内每个像素进行二次判定和识别，以精确提取水体像素，获取水体区域影像。

本发明与现有技术相比的优点在于：

(1)本发明在精确保持水体对象像素轮廓的同时，最大限度的保持了高光谱/多光谱影像的光谱特征。本发明采用像素级和决策级融合相结合的方法，引入面向对象的遥感影像解译方法，得到的融合结果一方面保证了水体的光谱特性，另一方面同时精确提取输出水体边缘轮廓，本发明可广泛适用于高光谱/多光谱传感器和单波段高分辨率遥感器影像的融合。

(2)本发明为最大限度保持水域提取精度以及水域光谱特性，结合PCA融合、SSVR融合设计了针对水质定量遥感应用的两级融合方法，融合过程引入决策级面向对象的地物分类解译方法，处理结果具备保持精准水体对象像素轮廓和水域光谱特性的特点。大量的实验结果表明。本方法得到水域解译精度达到90%，水域定量反演结果与直接利用高光谱/多光谱影像反演结果类似。

附图说明

图1为本发明的整体流程图；

图2为本发明主成分分析融合的流程图，对应具体实施方式中步骤2、步骤3、步骤4；

图3为本发明的水域分割提取的流程图，对应具体实施方式中步骤5、步骤6、步骤7、步骤8、步骤9；

图4为本发明的SSVR融合的流程图，对应具体实施方式中步骤10。

具体实施方式

下面结合附图1、图2、图3和图4对本发明的具体实施方式进行进一步的详细描述：

如图1所示，本发明具体实现如下：

1、基于高分全色影像I₁的地理信息，对高光谱/多光谱遥感影像I₂的重叠区域进行裁剪和配准，对裁剪和配准后高光谱/多光谱遥感影像进行三次样条插值重采样到高分全色影像的分辨率，重采样后的高光谱/多光谱遥感影像和高分全色影像具有相同的宽高。

2、计算高光谱/多光谱影像I₂的波段相关矩阵和特征向量，进行主成分分析，此时影像的主要光谱信息存储在前3个主成分分量影像I_PCA的波段中。

（1）影像主分量分析步骤：

（1.1）设有n幅影像，每幅影像观测p个分量，将原始数据标准化，得到

其中，X为数据矩阵，包含了所有图像的像素灰度值，X中的每个元素x_ij为第i幅影像第j个像素的灰度值，1≤i≤n，1≤j≤p。

（1.2）建立变量的协方差矩阵R=(r_ij)_p×p

（1.3）求R的特征值λ₁≥λ₂≥...≥λ_p及相应的单位特征向量

$A_1=\left[\begin{array}{c}{a}_{11}\\ a_{21}\\ .\\ .\\ .\\ a_{p1}\end{array}\right],A_2=\left[\begin{array}{c}{a}_{12}\\ a_{22}\\ .\\ .\\ .\\ a_{p2}\end{array}\right]...A_p=\left[\begin{array}{c}{a}_{1p}\\ a_{2p}\\ .\\ .\\ .\\ a_{\mathrm{pp}}\end{array}\right]$

其中，A_j为特征值λ_j对应的特征向量，即矩阵A的第j列，其中1≤j≤p。a_ij为特征向量矩阵A的的元素，1≤i≤p，1≤i≤p。

（1.5）写出主成分

F_i=A_1iX₁+A_2iX₂+...+A_piX_p，(i=1,2,…,p)

其中，F_i为第i主成分，A_pi为特征向量矩阵A的第p行i列的元素，X_p为数据矩阵X的第p列。

式中有协方差矩阵求特征值和特征向量，由于各波段影像方差不同，将导致各波段重要程度不一致。而如果采用相关矩阵求特征值和特征向量，由于相关矩阵中各波段的方差都归一化，从而使各波段具有同等的重要性。

（2）在本发明中，高光谱/多光谱主分量分析法的具体步骤为：

（2.1）计算参与融合的n波段TM影像的相关矩阵

（2.2）由相关矩阵计算特征值λ_i和特征向量A_i(i=1,...n)。

（2.3）将特征值按由大到小的次序排列，即λ₁≥λ₂≥...≥λ_n，特征向量A_i也要做相应变动。

（2.4）计算各主分量影像：式中k为主分量序数(k=1,...,n),PC_k为第k主分量，i为输入波段序数，n为总的TM影像波段数，d_i为i波段TM影像数据值，A_ik为特征向量矩阵在i行、k列的元素。经过上述主分量变换，第一主分量影像的方差最大，它包含原多光谱影像的大量信息（主要是空间信息），而原多光谱影像的光谱信息则保留在其他分量影像里（主要在第二、第三主分量中）。

3、将高分全色影像I₁和高光谱/多光谱影像的主成分影像I_PCA的第一波段影像进行直方图匹配，并用匹配后的高分全色影像I′₁替换高光谱/多光谱影像的第一主成分。

4、将替换后的主成分结果进行PCA逆变换，得到保持高光谱/多光谱影像光谱信息和全色影像空间信息初步融合影像I_M。

5、将初步融合影像I_M进行均值漂移聚类，提取初步融合影像I_M的边缘影像I_edge（此时边缘像素标记为1，非边缘像素标记为0）。采用3×3的圆形结构元素对地物边缘影像I_edge进行闭运算处理，得到边缘完全闭合的闭合边缘影像

（1）将全图进行遍历，将未进行分割的像素选择最为样本点的特征。设{x_i|i=1,2,...,N}为的像素灰度值构成的几何，其分割步骤如下：

（1.1）选择存储的第一个像素灰度值为c₁，此时C={c₁}。

（1.2）计算x_i与c_j∈C的欧式距离，如果||x_i-c_j||<T，则把像素分为S_j类，继续下一个像素的计算。

（1.3）如果||x_i-c_j||>T，则把x_i定义为新的聚类中心，继续下一个像素的计算。

（1.4）分类结束时，更新聚类中心其中聚类中心的集合为C={c_j|j=1,2,...,m}。

（1.5）对c_j∈C，根据下式开始迭代，其中

最后将迭代后的收敛点记为

$c_j^{(k+1)}=\frac{\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{n}_p{c}_{p}g\left({\left|\right|\frac{c_j^{\left(k\right)}-c_p}{h}\left|\right|}^2\right)}{\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{n}_{p}g\left({\left|\right|\frac{c_j^{\left(k\right)}-c_p}{h}\left|\right|}^2\right)}$

式中，h为用户指定的参数，

为第k次迭代第j个聚类中心的值，g(·)为用户指定的函数。

（1.6）最后对聚类收敛点进行合并，(1≤a,b≤m,a≠b)，如果则将a，b两类合并，其中r为合并阈值。

（2）在（1）的过程之后，得到分割后的边缘线影像I_edge，对于整体边缘线采用闭运算处理。闭运算处理的过程如下式所示：

${\tilde{I}}_{\mathrm{edge}}=I_{\mathrm{edge}}\&CenterDot;\left[\begin{array}{ccc}0& 1& 0\\ 1& 1& 1\\ 0& 1& 0\end{array}\right].$

式中

为闭运算处理后的边缘线影像。

6、在闭合边缘影像

内逐一对每一个非边缘像素进行不同的标记（1～N），作为种子点进行区域生长（已被区域生长过的像素不作为种子点）。全图遍历完成后，每一个闭合边缘内部都被唯一标记给标识，不同闭合边缘内部像素的标识不同，边缘像素标识为0，形成标识影像g。

7、基于标识影像g对初步融合影像I_M进行分水岭分割，形成同质影像对象，计算影像对象的光谱异质性指数和形状异质性指数，对性质相似的同质影像对象进行反复迭代合并。当无法再合并时，此时同质影像对象合并成为具备地理意义的地物影像对象。

（1）影像对象的光1谱特征包括光谱均值、光谱均方差、归一化差异水体指数（NDWI），计算公式如下（其中k为波段，N为k波段的像素总数）：

（1.1）光谱均值计算公式：

${\mathrm{Mean}}_k=\frac{\mathrm{\&Sigma;}{g}_k}{N_k}$

式中Mean_k为第k波段的光谱均值，∑g_k为标识影像g的第k个波段的像素总和。N_k为影像第k个波段的像素个数。

（1.2）光谱均方差计算公式：

${\mathrm{Sigma}}_k=\frac{\sqrt{\mathrm{\&Sigma;}{(g_k-{\mathrm{Mean}}_k)}^2}}{N_k}$

式中Sigma_k为第k波段的光谱均方差，g_k为标识影像g的第k个波段。

（1.3）归一化差异水体指数，其中Green为绿色波段灰度值，NIR为近红外波段灰度值：

$\mathrm{NDWI}=\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}$

式中NDWI归一化差异水体指数。

（2）影像对象的形状特征包括：面积（像素数目）、长宽比（最小外接矩形长宽比值）、长度（最小外接矩形长度）、宽度（最小外接矩形长度）、边界长度（边界像素的数目）、形状指数、密度和对称性：

（2.1）影像对象的形状指数S描述了对象边界的光滑性，当边界越破碎时，对象形状指数越大；当对象的形状越趋于圆时，对象的形状指数也就越小。计算公式如下，其中e是边界长度，A为面积：

$S=\frac{e}{4\&CenterDot;\sqrt{A}}$

（2.2）影像对象的密度d描述了一个对象的紧致性，若此对象越接近正方形，此数值越大。计算如下式所示，其中A是影像对象的面积，X和Y分辨是该影像对象所有像素x和y坐标组成的矢量矩阵，Var(x)和Var(y)分别是X和Y的方差：

$d=\frac{\sqrt{A}}{1+\sqrt{\mathrm{Var}\left(x\right)+\mathrm{Var}\left(y\right)}}$

（2.3）影像对象的对称性x定量描述了对象的对称特征，一个对象越长，并且在外接椭圆长半轴m和短半轴n两边信息越丰富，此数值越大，对称性的计算方法如下：

$x=1-\frac{n}{m}$

（3）本发明采用的异质性指数f是基于光谱差异性量度准则和形状差异性度量准则进行区域合并，计算公式如下：

f=ω_color·h_color+ω_shape·h_shape

（3.1）上式中，ω_color为光谱因子权重，ω_shape为形状因子权重，ω_color+ω_shape=1。

（3.2）光谱异质性指数计算公式如下：

$h_{\mathrm{color}}=\underset{c}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_c\&CenterDot;(n_{\mathrm{Merge}}\&CenterDot;{\mathrm{\&sigma;}}_c^{\mathrm{Merge}}-(n_{\mathrm{obj}1}\&CenterDot;{\mathrm{\&sigma;}}_c^{\mathrm{Obj}1}+n_{\mathrm{obj}2}\&CenterDot;{\mathrm{\&sigma;}}_c^{\mathrm{Obj}2}))$

在上式中，ω_c表示c波段所占的权重，n_Merge表示两对象obj1和obj2合并后的对象所含的数据个数，

表示合并后对象在c波段的标准差。

（3.3）形状异质性指数是由影像对象的光滑度指数h_smooth和紧凑度指数h_compact加权求和所得，计算公式如下：

h_shape=ω_smooth·h_smooth+ω_compact·h_compact

$h_{\mathrm{smooth}}=n_{\mathrm{Merge}}\&CenterDot;\frac{l_{\mathrm{Merge}}}{b_{\mathrm{Merge}}}-(n_{\mathrm{Obj}1}\&CenterDot;\frac{l_{\mathrm{Obj}1}}{b_{\mathrm{obj}1}}+n_{\mathrm{Obj}2}\&CenterDot;\frac{l_{\mathrm{Obj}2}}{b_{\mathrm{obj}2}})$

$h_{\mathrm{compact}}=n_{\mathrm{Merge}}\&CenterDot;\frac{l_{\mathrm{Merge}}}{4\sqrt{n_{\mathrm{Merge}}}}-(n_{\mathrm{Obj}1}\&CenterDot;\frac{l_{\mathrm{Obj}1}}{4\sqrt{n_{\mathrm{Obj}1}}}+n_{\mathrm{Obj}2}\&CenterDot;\frac{l_{\mathrm{Obj}2}}{4\sqrt{n_{\mathrm{Obj}2}}})$

其中，ω_smooth和ω_compact分别是紧凑度和光滑度所占的权重，ω_smooth+ω_compact=1。l为对象的周长，n为对象所包含像素的个数，b为最小外接矩形的周长。对象的光滑异质性即为对象周长和包含此对象的最小外接矩形周长的比值，紧凑度异质性即为对象的周长和对象所含像素总数的平方根之比。

8、计算完全合并后对象的光谱特征、形状特征和纹理特征，基于贝叶斯分类器和外部输入的训练结果，对水体对象进行粗判定。

（1）纹理特征的灰度-梯度共生矩阵（GGCM）计算：

（1.1）对原始影像的灰度进行归一化处理，归一化处理的方法为将原始影像灰度的最小值-最大值的区间内的灰度值线性归化到[0，N_g]。

（1.2）采用Sobel算子对原始影像计算梯度，获取梯度影像进行归一化处理，将梯度影像灰度的最小值-最大值的区间内的灰度值线性归化到[0，N_s]。

（1.3）灰度-梯度共生矩阵定义为：

{M(i,j),i=1,2,...,N_g,j=1,2,...,N_s}

矩阵M的矩阵元素m_i,j标识在归一化后的灰度和梯度影像中，具有第i级灰度和第j级梯度的像素点总是，或者取灰度-梯度的规划矩阵，其元素表示为：

P(i,j)=m_i,j/(N_xN_y)i=1,2,....,N_g,j=1,2,...,N_s

（1.4）用于描述影像纹理特征如下表所示：

（2）Gabor小波变换是小波基为Gabor函数的小波变换，本文取Gabor函数g(x,y)如下：

$g(x,y)=\left(\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}\right)\exp [-(\frac{x^2}{2{\mathrm{\&sigma;}}_x^2}+\frac{y^2}{2{\mathrm{\&sigma;}}_y^2})+2\mathrm{\&pi;\&omega;x}]$

上式的傅立叶变换G(u,v)为

$G(u,v)=\exp \{-\frac{1}{2}[\frac{{(u-\mathrm{\&omega;})}^2}{{\mathrm{\&sigma;}}_u^2}-\frac{v^2}{{\mathrm{\&sigma;}}_v^2}\left]\right\}$

${\mathrm{\&sigma;}}_u=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_x}$

${\mathrm{\&sigma;}}_v=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_y}$

（2.1）以g(x,y)为母小波，通过对g(x,y)进行旋转变换和尺度变换，就可以得到一组自相似的滤波器，即为Gabor小波。形式如下所示：

g_m,n(x,y)=a^-mg(x',y')，a>1，m,n∈Z

式中：

x'=a^-m(xcosθ+ysinθ)

y'=a^-m(-xsinθ+ycosθ)

其中，θ=nπ/k，n∈[0,k-1]，m∈[0,s-1]，k表示总的方向数目，a^-m为尺度因子，s为总的尺度数目。同时为了确保其总能量和m无关，再由傅立叶变换的线性特征可得：

$\left.\begin{array}{c}{u}^{\&prime;}=u\cos \mathrm{\&theta;}+v\sin \mathrm{\&theta;}\\ v^{\&prime;}=-u\sin \mathrm{\&theta;}+v\cos \mathrm{\&theta;}\end{array}\right\}$

（2.2）因此，通过改变m和n的值，便可以得到一组尺度和方向都不同的Gabor滤波器。Gabor小波集具有非正交性，在滤波后得到的影像具有大量的冗余信息，因此采取中心频率范围为[U₁,U_h]进行多通道滤波。Gabor小波多通滤波的参数计算方法如下式：

$a={\left[\frac{U_h}{U_1}\right]}^{\frac{1}{s-1}}$

${\mathrm{\&sigma;}}_u=\frac{(a-1)U_h}{(a+1)\sqrt{2\ln 2}}$

${\mathrm{\&sigma;}}_v=\tan \left(\frac{\mathrm{\&pi;}}{2k}\right)[U_h-2\ln \left(\frac{{\mathrm{\&sigma;}}_u^2}{U_h}\right)]\&times;[2\ln 2-\frac{{\left(2\ln 2\right)}^2}{U_h^2}]$

在本发明对影像对象进行纹理特征提取时，本文选择的总尺度数目为4，方向数为6，多通滤波器的频率范围是

即各个参数设置如下：

$\mathrm{\&omega;}{=U}_h,k=6,s=4,U_1=\frac{\sqrt{2}}{8},U_h=\frac{\sqrt{2}}{2}$

（2.3）在本发明中通过Gabor小波对影像进行滤波，将所有系数模的均值和方差值组合而成的48维纹理特征向量T=[μ₀,σ₀,μ₁,σ₁,...,μ₂₃,σ₂₃]作为影像纹理特征。在本发明中，采用均值方法降低特征向量的维度，则此时 $T^{\&prime;}=[\stackrel{\&OverBar;}{u},\stackrel{\&OverBar;}{\mathrm{\&sigma;}}]:$

$\stackrel{\&OverBar;}{\mathrm{\&mu;}}=[({\mathrm{\&mu;}}_{\max }+{\mathrm{\&mu;}}_{\min })+\frac{1}{M}\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i]/2$ $\stackrel{\&OverBar;}{\mathrm{\&sigma;}}=[({\mathrm{\&sigma;}}_{\max }+{\mathrm{\&sigma;}}_{\min })+\frac{1}{M}\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}{\mathrm{\&sigma;}}_i]/2$

其中，μ_max和μ_min是所有系数模的最大值和最小值，σ_max和σ_min是所有系数模方差的最大值和最小值。

（3）本发明中采用的分类器为朴素贝叶斯分类器。朴素贝叶斯分类以贝叶斯定理为基础，通过对大量样本的训练估算对象属性的后验概率。设有离散或连续的属性变量X₁，X₂，…，X_n，C是待分类对象，D是训练集。因此，由贝叶斯定理和公式可以得到：

$\begin{array}{c}P\left(C\right|X_1,X_2,...,X_n)=\frac{P\left(C\right)P(X_1,X_2,...,X_n|C)}{P(X_1,X_2,...,X_n)}\\ =\mathrm{\&alpha;P}\left(C\right)P(X_1,X_2,...,X_n|C)=\mathrm{\&alpha;P}\left(C\right)\mathrm{\&Pi;P}\left(X_i\right|C)\end{array}$

（3.1）

其中N(C_k)是C_k类中属性类别的个数。

（3.2）若X_i是离散变量，则：

$P(X_i=x_i|C_k)=\hat{P}(X_i=x_i|C_k)=\frac{N_{C_k}^{\left(x_i\right)}}{N_{C_k}}$

其中

是在第C_k类中，X_i=x_i时的情况数量。

（3.3）如果在存在X_i=x_i情况数量为0，则：

$\hat{P}(X_i=x_i|C_k)=\frac{\frac{1}{N}}{N_{C_k}+\frac{N_{x_i}}{N}}$

其中是属性变量X_i的取值个数。

（3.4）若属性变量X_i是连续变量，则：

$P(X_i=x_i|C_k)=g(x_i,{\mathrm{\&mu;}}_{c_k},{\mathrm{\&sigma;}}_{c_k})$

即P(X_i=x_i|C_k)等于第C_k类中的样本平均值：

$g(x_i,{\mathrm{\&mu;}}_{c_k},{\mathrm{\&sigma;}}_{c_k})=\frac{1}{\sqrt{2\mathrm{\&pi;\&sigma;}}}{e}^{\frac{{(x_i-{\mathrm{\&mu;}}_{C_k})}^2}{2{\mathrm{\&sigma;}}^2}}$

${\mathrm{\&mu;}}_{C_k}=\frac{x_{\mathrm{li}}\left(C_k\right)+...+x_{\mathrm{Ni}}\left(C_k\right)}{N_{C_k}}$

$x_{\mathrm{ji}}\left(C_k\right)=\left\{\begin{array}{c}{x}_{\mathrm{ji}},x_{\mathrm{ji}}\&Element;C_k\\ 0,x_{\mathrm{ji}}\&NotElement;C_k\end{array}\right.$

$N_{C_k}=\left\{\begin{array}{c}{N}_{C_k},x_{\mathrm{ji}}\&Element;C_k\\ 0,x_{\mathrm{ji}}\&NotElement;C_k\end{array}\right.$

${\mathrm{\&sigma;}}_{C_k}^2=\frac{{(x_{\mathrm{li}}\left(C_k\right)-{\mathrm{\&mu;}}_{C_k}\left(C_k\right))}^2+...+{(x_{\mathrm{Ni}}\left(C_k\right)-{\mathrm{\&mu;}}_{C_k}\left(C_k\right))}^2}{N_{C_k}}$

其中，为第C_k类样本元素个数

为第C_k类样本的方差，

为第C_k类样本的均值：

${\mathrm{\&mu;}}_{C_k}\left(C_k\right)=\left\{\begin{array}{c}{\mathrm{\&mu;}}_{C_k},x_{\mathrm{ji}}\&Element;C_k\\ 0,x_{\mathrm{ji}}\&NotElement;C_k\end{array}\right.$

（3.5）通过对已知训练集D属性值

进行训练，并得到待分类对象C的后验概率P(X₁|C)，P(X₂|C)，…，P(X_n|C)的估计值，并将分类对象C的类别落入概率最大的类别中。

（3.6）输入水体和非水体对象，通过光谱、形状和纹理特征对贝叶斯分类器进行训练，形成训练文件。

（3.7）通过训练文件判定影像对象是否为水体对象，将判定结果对训练文件进行更新。

9、对获取的水体对象和非水体对象的光谱特征作为训练对象，重新训练贝叶斯分类器，基于光谱特征，对水体对象中每个像素进行精判定，并输出水体区域影像I_water。

10、输入高光谱/多光谱影像每个波段宽度和全色影像波段宽度，计算水域影像每个像素的能量值，对水域影像I_water的水体像素进行像素级的SSVR融合，保证水域像素的光谱特性，形成水域融合影像I_m-water。

（1）SSVR融合算法流程：

（1.1）将遥感影像的灰度值转换为辐射亮度，转换公式为：

$L_k=\frac{{\mathrm{DN}}_k\&CenterDot;(C_k^{\max }-C_k^{\min })}{{\mathrm{DN}}_{\max }}+C_k^{\min }$

其中L_k是第k个波段的辐射亮度，和

是辐射亮度校正系数，通过影像头文件中获取。DN_k为第k波段影像的灰度值，DN_max为第k波段影像最大的灰度值。

（1.2）计算遥感影像各像元的能量值，计算公式为：

$E_{{\mathrm{XS}}_i}=L_i\&CenterDot;W_i$

其中，

是第i波段的能量值，L_i是第i波段的辐射亮度，W_i为第i波段的波段宽度。

（1.3）计算

计算出高光谱/多光谱影像每个波段和全色影像的能量值比例。

（1.4）计算SSVR，

XSP_i为融合后。为融合后的高分辨率影像中第i个波段的辐射亮度值，Pan_H是原始高分辨率全色波段影像的辐射亮度值，

是原始低分辨率多光谱影像中第i个波段的能量值，

是由原始高分辨率全色波段影像重采样到低分辨率后的能量值。

11、基于初步融合影像I_M和水域融合影像I_m-water进行镶嵌处理，得到最终的影像融合影像。

总之，本发明最大限度保持水域提取精度以及水域光谱特性，结合PCA融合、SSVR融合设计了针对水质定量遥感应用的两级融合方法，融合过程引入决策级面向对象的地物分类解译方法，处理结果具备保持精准水体对象像素轮廓和水域光谱特性的特点。大量的实验结果表明。本发明得到水域解译精度达到90%，水域定量反演结果与直接利用高光谱/多光谱影像反演结果类似。

本发明未详细阐述的部分属于本领域公知技术。

Claims (3)

1.一种面向水质定量遥感应用遥感影像融合方法，其特征在于包含以下步骤：

（1）基于高分全色影像I₁的地理信息，对高光谱/多光谱遥感影像I₂进行配准和重采样，并对两幅影像的地理位置重叠区域进行裁剪，得到具有同一地理区域且相同大小的高分全色影像I₁和高光谱/多光谱影像I₂，两幅影像分辨率一致，都为高分全色影像的分辨率；所述高分影像为高空间分辨率遥感影像，高光谱影像为影像波段数一般多于100个的遥感影像，多光谱影像为影像波段数一般多于3个少于100个的遥感影像;

（2）计算高光谱/多光谱影像I₂的波段相关矩阵和矩阵的特征向量，进行主成分分析，得到主成分分量影像I_PCA；

（3）将高分全色影像I₁和主成分影像I_PCA的第一波段影像进行直方图匹配，并用直方图匹配的结果I′₁替换I_PCA的第一主成分，即替换掉I_PCA的第一波段为I′₁，得到替换后的主成分结果；

（4）将替换后的主成分结果进行主成分分析逆变换，得到保持高光谱/多光谱影像光谱信息和全色影像空间信息的初步融合影像I_M；

（5）将初步融合影像I_M进行均值漂移聚类，提取初步融合影像I_M的边缘影像I_edge，边缘影像为二值影像，其中边缘像素标记为1，非边缘像素标记为0；

（6）针对获取的边缘影像I_edge中可能存在的边缘轮廓不完全闭合的情况，采用3×3的圆形结构元素对地物边缘影像I_edge进行闭运算处理，得到边缘完全闭合的闭合边缘影像

（7）在闭合边缘影像

内逐一对每一个非边缘像素进行不同的标记N，N为取值大于等于1的整数，将这些像素作为种子点进行区域生长，在区域生长中已被其他像素区域生长过的像素不作为种子点；全图遍历完成后，每一个闭合边缘内部都被唯一标记给标识，不同闭合边缘内部像素的标识不同，边缘像素标识为0，得到标识影像g；

（8）基于标识影像g对初步融合影像I_M进行分水岭分割，得到同质影像对象；针对同质影像对象细碎的情况，计算影像对象及其邻接影像对象的光谱异质性指数和形状异质性指数，作为衡量影像对象性质相似程度，对细碎影像对象性质相似的进行反复迭代合并，当无法再合并时，此时同质影像对象合并得到具有明确地理意义的地物影像对象；

（9）对影像对象合并后的影像计算合并后影像对象的光谱特征、形状特征和纹理特征，基于人工解译信息训练粗判定贝叶斯分类器，并保存训练结果，相同类型影像只需要训练一次，对水体对象和非水体对象进行粗判定，得到判定后的水体对象和非水体对象的影像光谱特征；（粗判定贝叶斯分类器指判定对象为影像对象的贝叶斯分类器，影像对象包括水体对象和非水体对象等）

（10）对判定的水体对象和非水体对象的光谱特征作为训练对象，训练精分类贝叶斯分类器，对水体对象中每个像素进行精判定，并输出水域影像I_water；所述精分类贝叶斯分类器指分类对象为影像像素的贝叶斯分类器，即此分类器对每个像素进行分类；

（11）输入高光谱/多光谱影像每个波段宽度和全色影像波段宽度，计算水域影像每个像素的能量值，对水域影像I_water的水体像素进行像素级的SSVR融合，保证水域像素的光谱特性，形成水域融合影像I_m-water；SSVR为Simplified Synthetic Variable Ratio的缩写，即简易合成变量比；

（12）基于初步融合影像I_M和水域融合影像I_m-water进行镶嵌处理，得到最终的融合影像。

2.根据权利要求1的面向水质定量遥感应用遥感影像融合方法，其特征在于：所述步骤（5）和（6）具体实现如下：

（2.1）将全图进行遍历，将未进行分割的像素选择为样本点的特征，设{x_i|i=1,2,...,N}为图像像素灰度值集合，其中x_i为图像第i个位置灰度值，N为图像像素个数，其分割步骤如下：

（2.1.1）选择存储的第一个像素灰度值为c₁，此时令聚类中心集合C={c₁}；

（2.1.2）计算x_i与c_j∈C的欧式距离，c_j为第j个聚类中心，如果||x_i-c_j||T，T为距离阈值，则把像素分到S_j类中，S_j是第j类的灰度值集合，继续下一个像素的计算；

（2.1.3）如果||x_i-c_j||T，则把x_i定义为新的聚类中心，继续下一个像素的计算；

（2.1.4）分类结束时，更新聚类中心

其中聚类中心的集合为C={c_j|j=1,2,...,m}，m为聚类中心个数，n_j为第j个聚类的像素个数；

（2.1.5）对c_j∈C，根据下式开始迭代，其中令

最后将迭代后的收敛点记为

$c_j^{(k+1)}=\frac{\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{n}_p{c}_{p}g\left({\left|\right|\frac{c_j^{\left(k\right)}-c_p}{h}\left|\right|}^2\right)}{\underset{p=1}{\overset{m}{\mathrm{\&Sigma;}}}{n}_{p}g\left({\left|\right|\frac{c_j^{\left(k\right)}-c_p}{h}\left|\right|}^2\right)}$

式中，h为用户指定的参数，为第k次迭代第j个聚类中心的值，g(·)为用户指定的函数；

（2.1.6）最后对聚类收敛点进行合并，(1≤a,b≤m,a≠b)，如果则将a，b两类合并为一类，其中r为合并阈值；

（2.2）在（2.1）的过程之后，得到分割后的边缘线影像I_edge，对于整体边缘线采用闭运算处理。闭运算处理的过程如下式所示：

${\tilde{I}}_{\mathrm{edge}}=I_{\mathrm{edge}}\&CenterDot;\left[\begin{array}{ccc}0& 1& 0\\ 1& 1& 1\\ 0& 1& 0\end{array}\right].$

3.根据权利要求1的面向水质定量遥感应用遥感影像融合方法，其特征在于：所述步骤（9）和（10）具体实现如下：

(3.1)基于相同影像只需要进行一次水体对象特征训练，基于贝叶斯分类器对水体影像对象进行粗提取。水体对象特征包括光谱特征、形状特征、纹理特征，纹理特征包括灰度-梯度共生矩阵提取结果，其中有能量、相关、梯度平均、混合熵、惯性矩和逆差矩；和Gabor小波滤波提取结果，包括Gabor小波均值和Gabor小波方差；

（3.2）基于获取的水体影像对象和非水体对象的光谱特征，对贝叶斯分类器进行训练，对水体内每个像素进行二次判定和识别，以精确提取水体像素，获取水体区域影像。

Images