基于超复数的极化SAR图像多通道联合匹配方法

摘要

本发明公开了一种基于超复数的极化SAR图像多通道联合匹配方法，首先对极化SAR图像进行Pauli分解，构造极化SAR图像的Pauli矢量超复数形式，然后对Pauli矢量超复数进行偶对分解，并基于分解结果进行快速傅里叶变换，最后通过计算超复数相位相关，获取待匹配极化SAR图像的匹配位置。实验结果表明，本发明公开的方法即使在极化SAR图像之间存在局部差异及包含乘性噪声的情况下，仍能够获得精确的匹配结果，具有较好的匹配性能。

说明书

技术领域

本发明涉及一种极化SAR图像匹配方法，尤其涉及一种基于超复数的极化SAR图像多通道联合匹配方法。

背景技术

近年来，极化合成孔径雷达(Synthetic Aperture Radar，SAR)图像匹配技术快速发展，RADARSAT-2、Terra SAR-X、ALOS-2等星载极化SAR陆续升空，其在军事、民事应用中起到越来越重要的作用，如图像匹配、目标检测、图像融合、飞行器定位与导航等领域。其中，基于极化SAR的很多应用需要综合利用不同成像条件下(不同时间、不同波段、不同下视角、不同入射方向或不同传感器等)的极化SAR图像的信息实现，如地形矫正、环境监测、飞行器定位与导航、精制导等。对于这些应用，极化SAR图像匹配是不可或缺的步骤，图像匹配是在两幅图像中寻找同一地理位置的对应像素的坐标关系。

与传统的单极化SAR图像相比，极化SAR图像包含4个极化通道的SAR图像，每个极化通道的图像对匹配效果的贡献不同，尤其是4个极化通道之间的相对关系，包含了能够反映地物目标散射结构的极化信息，若能有效地用于图像匹配，可以大大提高匹配效果。

国内外有关单极化SAR图像匹配方法的研究已经持续有30年之久，相关文献也很多，但有关极化SAR图像匹配方法的文献却很少。目前最常用的匹配方法是通过某种融合规则，如极化总功率(SPAN)、极化白化滤波器(PWF)等，将多个通道的图像融合为一个通道，然后使用单极化SAR图像匹配方法完成融合后的图像的匹配。其中，最具代表性的是Chureesampant等人在2014年的IEEE Transactions on Geoscience and Remote Sensing杂志中发表的文章“Automatic GCP extraction of fully polarimetric SAR images”中提出的SPAN+改进的SIFT-OCT(尺度旋转不变变换)方法，并通过实测数据对该方法的匹配效果进行了验证，该方法中，极化信息的贡献仅体现在图像融合过程，并没有在图像匹配过程中利用极化信息改善匹配效果。

其他方法，如Chou等人在2010年的论文“Multitemporal polarimetric SAR data fusion for land cover mapping”中提出的方法和Xie等人在2010年的论文“New approach for improving Mutual Information in POLSAR image registration”中提出的方法，也均未在图像匹配过程中利用极化信息。

综上，目前已有的极化SAR图像匹配方法，均是在单通道图像中实现匹配的，致使极化SAR图像中所包含的独特的极化信息并未能在图像匹配中发挥更好的效果，匹配效果有待提高。

另外，极化SAR图像中含有乘性噪声，导致图像的强度越大，噪声越强。国际上多篇论文中指出，例如文献“Applicability of the SIFT operator to geometric SAR image registration.International Journal of Remote Sensing，2010，31(8)：1959-1980”，乘性噪声大大降低了图像匹配的效果，是极化SAR图像匹配亟待解决的一大难题。

发明内容

发明目的：针对以上问题，本发明提出一种基于超复数的极化SAR图像多通道联合匹配方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于超复数的极化SAR图像多通道联合匹配方法，包括以下步骤：

(1)极化SAR图像Pauli分解，具体是将待匹配的两幅极化SAR图像进行Pauli分解，分别得到含有三个元素的Pauli矢量；

(2)Pauli矢量超复数化，具体是用Pauli矢量的三个元素对超复数的三个虚数通道赋值，实数通道置0，得到Pauli矢量的超复数形式；其中，超复数包含一个实数通道和三个虚数通道；

(3)超复数偶对分解，具体是将Pauli矢量形式的超复数分解成为包含单部和复部的广义复数形式；

(4)快速傅里叶变换，具体是对偶对分解完成的两幅图像进行左右傅里叶变换及逆变换；

(5)计算两幅图像的超复数相位相关，具体是利用两幅图像的左右傅里叶变换及其逆变换，通过计算得到两幅极化SAR图像的相位相关；

(6)计算极化SAR图像的最佳匹配位置，具体是在相位相关平面中寻找最大值，根据该最大值对应的横、纵坐标，计算获得两幅极化SAR图像的最佳匹配位置。

步骤(1)中极化SAR图像的Sinclair矩阵表示为：

其中，H、V分别表示水平、垂直极化方式，S_ij中i表示发射天线的极化方式，j表示接受天线的极化方式，S_ij表示在该发射接受极化状态下获得的目标复散射系数

Pauli分解后的矢量k为：

步骤(2)中超复数定义为：

q＝a+bi+cj+dk

其中，a、b、c、d为任意实数，i、j、k为复数算子；当a＝0时，q为纯超复数。

Pauli矢量的超复数形式为：

步骤(3)中将超复数偶对分解为：

q＝A′+B′μ₂

其中，μ₂为任意单位纯超复数，A′＝a′+b′μ₁为单部，B′＝c′+d′μ₁为复部，μ₁为垂直于μ₂的任意单位纯超复数，a′、b′、c′、d′为实数系数；因此：

q＝(a′+b′μ₁)+(c′+d′μ₁)μ₂

展开可得：

q＝a′+b′μ₁+c′μ₂+d′μ₃

其中，μ₃＝μ₁μ₂。

步骤(4)中设f(m,n)、g(m,n)为两幅极化SAR图像Pauli矢量的超复数矩阵，超复数的傅里叶变换分为左傅里叶变换与右傅里叶变换，以f(m,n)为例：

其中，M、N为矩阵的行数和列数，μ为傅里叶变换轴，L为左(逆)傅里叶变换标识，R为右(逆)傅里叶变换标识，+为傅里叶变换标识，－为逆傅里叶变换标识，v、u是傅里叶变换后频域矩阵的行、列。

以超复数左傅里叶变换为例，具体包含以下步骤：

步骤4.1：将f(m,n)偶对分解：

f(m,n)＝f₁(m,n)+f₂(m,n)μ₂

其中，f_i为偶对分解的分量，i＝1～2，f_i∈(1,μ₁)平面；μ₁为傅里叶变换轴，与步骤3中的复数算子相同。

步骤4.2：展开偶对分量：

f_i(m,n)＝f_i,1(m,n)+f_i,2(m,n)μ₁

其中，f_i,1(m,n)和f_i,2(m,n)为实数。

步骤4.3：构造等价复矩阵：

f_i′(m,n)＝f_i,1(m,n)+f_i,2(m,n)j

其中，j为传统复数算子。

步骤4.4：对f_i′(m,n)进行FFT：

F_i′(v,u)＝F′_i,1(v,u)+F′_i,2(v,u)j

步骤4.5：用μ₁代替j，形成傅里叶变换的单部与复部：

F_i′(v,u)＝F′_i,1(v,u)+F′_i,2(v,u)μ₁

步骤4.6：组合单部与复部，完成傅里叶变换。

步骤(5)中两幅图像的相位相关ph(m,n)：

其中，R_R(v,u)为两幅极化SAR图像的互功率谱。

其中，^——为共轭运算；G_R||(v,u)、G_R⊥(v,u)分别为G_R(v,u)平行、垂直于傅里叶变换轴的分量；F_L(v,u)、F_-L(v,u)和G_R(v,u)分别为f(m,n)的左傅里叶变换、f(m,n)的逆左傅里叶变换和g(m,n)的右傅里叶变换。

步骤(6)中假设相位相关ph(m,n)幅值最大的位置坐标为(m₀,n₀)，待匹配的两幅极化SAR图像f(m,n)、g(m,n)的匹配关系为：

f(m,n)＝g(m-m₀,n-n₀)

有益效果：与现有技术相比，本发明将超复数引入极化SAR图像处理领域，并提出了极化SAR图像的Pauli矢量超复数形式，为极化SAR图像匹配及其他应用的多通道联合处理提供了途径；与现有方法仅在融合后的单通道图像实施匹配不同，本发明在图像匹配的同时完成多通道图像的融合，是一种多通道联合匹配方法，从而在匹配过程中能够更充分地利用极化SAR图像的极化信息，对于提高图像匹配效果更有优势；本发明的方法具有良好的抗噪声性能，在频域对相位相关进行了幅值归一化，归一化后仅保留了相位信息，去除了大部分噪声的影响，因而得到的相位相关的相关峰非常明显，匹配结果受乘性噪声的影响比较小。

附图说明

图1是本发明的基于超复数的极化SAR图像多通道联合匹配方法的流程图；

图2是待匹配的极化SAR图像，图2a是L波段的图像，图2b是C波段的图像；

图3是待匹配极化SAR图像的相位相关图，图3a是相位相关平面图，图3b是相位相关峰值放大图，图3c是相位相关立体图；

图4是极化SAR图像匹配拼接图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

如图1所示，本发明所述的基于超复数的极化SAR图像多通道联合匹配方法，具体通过以下步骤实现：

(1)极化SAR图像的Pauli分解。

将待匹配的两幅极化SAR图像进行Pauli分解，分别得到含有三个元素的Pauli矢量。极化SAR图像通常使用2×2的复数矩阵Sinclair矩阵表示，水平垂直极化基时如式1所示：

其中，H、V分别表示水平、垂直极化方式，S_ij中i表示发射天线的极化方式，j表示接受天线的极化方式，S_ij表示在该种发射接受极化状态下获得的目标复散射系数。

若雷达为单站工作方式，且介质各向同性时，互易性条件成立，此时S_IIV＝S_VII。Pauli分解后的矢量形式k为：

Pauli矢量是极化SAR图像应用中一种常用的矢量形式，Pauli矢量的三个元素的强度分别反映了目标的平面散射成分、二面角散射成分和交叉通道散射成分。

(2)Pauli矢量的超复数化。

超复数的定义如式3所示：

q＝a+bi+cj+dk3

其中，a、b、c、d为任意实数，当a＝0时，q为纯超复数，i、j、k为复数算子，其运算关系如式4所示：

i²＝j²＝k²＝ijk＝-1

ij＝jk＝ki＝1

ji＝kj＝ik＝-14

如式3所示，超复数包含一个实数通道和三个虚数通道，四个通道的值均为实数。而Pauli矢量包含三个元素，且均为复数，无法直接构成超复数。本发明在Pauli矢量超复数化时，使用Pauli矢量三个元素的强度对超复数的三个虚数通道赋值，实数通道置0，从而得到Pauli矢量的纯超复数形式。

通过步骤1和步骤2，可将多通道的极化SAR图像用一个超复数表示，本发明通过这种间接升维的方式实现极化SAR图像的多通道联合匹配。

(3)超复数的偶对分解。

将Pauli矢量形式的超复数分解为包含“单部”和“复部”的广义复数形式，其中“单部”和“复部”又均是“复数算子”(定义为任意单位纯超复数)垂直空间内的虚数。需注意的是，偶对分解中的复数算子与步骤4中的傅里叶变换轴相同，通过偶对分解可实现步骤4中的快速傅里叶变换。

因此，任意一个超复数q均可写成：

q＝A′+B′μ₂6

其中，μ₂为任意单位纯超复数，A′＝a′+b′μ₁为单部，B′＝c′+d′μ₁为复部，μ₁为垂直于μ₂的任意单位纯超复数，a′、b′、c′、d′为实数系数，因此：

q＝(a′+b′μ₁)+(c′+d′μ₁)μ₂7

将式7展开：

q＝a′+b′μ₁+c′μ₂+d′μ₃8

其中，μ₃＝μ₁μ₂，可知μ₃⊥μ₁、μ₃⊥μ₂。

因此，式8可以理解为超复数q在空间基(μ₁，μ₂，μ₃)下的表达形式。从式3到式8的变换，可以理解为从空间基(i，j，k)到空间基(μ₁，μ₂，μ₃)的坐标变换。

需注意的是，偶对分解中的复数算子与步骤4中的快速傅里叶变换轴要完全相同。

(4)快速傅里叶变换。

在实际应用中，极化SAR图像的图幅通常非常庞大，导致极化SAR图像利用超复数升维后，傅里叶变换的计算量急剧升高。本发明基于偶对分解完成两幅图像的左右傅里叶变换及其逆变换的快速计算。设f(m,n)、g(m,n)分别为待匹配的两幅极化SAR图像Pauli矢量的二维超复数矩阵。

由于超复数不满足交换率，因此超复数的傅里叶变换分为左傅里叶变换与右傅里叶变换。以f(m,n)为例，分别如式9和式10所示。

其中，f(m,n)为二维超复数矩阵，M、N为矩阵的行数和列数，μ为傅里叶变换轴，L为左(逆)傅里叶变换标识，R为右(逆)傅里叶变换标识，+为傅里叶变换标识，－为逆傅里叶变换标识，v、u是傅里叶变换后频域矩阵的行、列。

以超复数左傅里叶变换快速计算为例，步骤如下：

步骤4.1：将f(m,n)偶对分解：

f(m,n)＝f₁(m,n)+f₂(m,n)μ₂11

其中，f_i为偶对分解的分量，i＝1～2，f_i∈(1,μ₁)平面；μ₁为傅里叶变换轴，与步骤3中的复数算子相同。

步骤4.2：展开偶对分量：

f_i(m,n)＝f_i,1(m,n)+f_i,2(m,n)μ₁12

其中，f_i,1(m,n)和f_i,2(m,n)为实数。

步骤4.3：构造等价复矩阵：

f_i′(m,n)＝f_i,1(m,n)+f_i,2(m,n)j13

其中，j传统复数算子。实际上，式12与式13在数值上是等价的，因为计算时仅存储f_i,1(m,n)与f_i,2(m,n)的值，μ₁与j的值是不存储的，且二者的运算规则相同。

步骤4.4：对f_i′(m,n)进行二维复FFT：

F_i′(v,u)＝F′_i,1(v,u)+F′_i,2(v,u)j14

步骤4.5：用μ₁代替j，形成傅里叶变换的简部与复部：

F_i′(v,u)＝F′_i,1(v,u)+F′_i,2(v,u)μ₁15

步骤4.6：组合简部与复部，完成傅里叶变换。

F(v,u)＝F₁(v,u)+F₂(v,u)μ₂16

同理，超复数的右傅里叶变换及逆左、右傅里叶变换均可使用以上步骤快速计算。

(5)计算两幅图像的相位相关ph(m,n)：

其中，R_R(v,u)为两幅极化SAR图像的互功率谱。

其中，^——为共轭运算；G_R||(v,u)、G_R⊥(v,u)分别为G_R(v,u)平行、垂直于傅里叶变换轴的分量；F_L(v,u)、F_-L(v,u)和G_R(v,u)分别为f(m,n)的左傅里叶变换、f(m,n)的逆左傅里叶变换和g(m,n)的右傅里叶变换。

(6)计算极化SAR图像的最佳匹配位置。

在步骤5得到的相位相关平面中寻找最大值，根据该最大值所对应的横、纵坐标计算得到两幅极化SAR图像的最佳匹配位置。

假设相位相关ph(m,n)幅值最大的位置坐标为(m₀,n₀)，那么待匹配的两幅极化SAR图像f(m,n)、g(m,n)的匹配关系为：

f(m,n)＝g(m-m₀,n-n₀)19

下面结合具体实验数据进行说明。

本实验中使用的数据为1993年8月12日，由AirSAR得到的加拿大北纬54°西经104°附近的图像，距离向分辨率为6.662m，方位向分辨率为12.1m。该数据包括一次飞行得到的L波段和C波段的极化SAR图像，待匹配的两幅极化SAR图像的大小均为400×400像素，二者的行、列位移实际为(100，100)。

步骤1：Pauli分解。

将待匹配的两幅极化SAR图像Pauli分解，并使用蓝色显示平面分量(|S_IIII+S_VV|²)，对应的典型地物为水域；使用红色显示二面角分量(|S_IIII-S_VV|²)，通常人造目标具有较强的二面角分量；使用绿色显示交叉通道分量(2|S_IIV|²)，典型地物为植物，与人类的视觉相符合。本实施例中将其灰度化显示，如图2a和图2b所示，图2a为L波段的图像，图2b为C波段的图像。

从图中可以看出，由于两幅极化SAR图像除波段不同外，其余的成像条件均相同，所以局部几何失真较小；但不同波段的电磁波穿透力不同，得到的目标散射特性不同，所以Pauli矢量三个元素的强度存在差异，导致图2a和图2b之间存在不同。另外也可以看出，两幅极化SAR图像中的斑点噪声非常严重。

步骤2：根据式5将Pauli矢量超复数化。

步骤3：超复数的偶对分解。

偶对分解中的复数算子取该复数算子对应的分量为极化SAR图像的强度分量。另外

偶对分解后有：

f(m,n)＝f₁(m,n)+f₂(m,n)μ₂20

其中，

步骤4：基于上述偶对分解实现快速傅里叶变换。

其中，傅里叶变换轴为

步骤5：根据式17～18计算两幅图像的相位相关ph(m,n)。

为了更清楚的展示相位相关的数值变化，本实施例中采用平面图和立体图两种方式展示，如图3所示。图3a为图2中两幅极化SAR图像的相位相关平面图，图3b为其峰值处的放大图，其中横、纵坐标分别为两幅图像的行、列平移量。相关平面中峰值位置(白色圆内)的坐标即为两幅图像的匹配位置。图3c为相位相关平面的立体图，其中平面坐标为两幅图像的行、列平移量，纵坐标为相位相关值。从图3a、3b、3c中可以看出，即使在斑点噪声非常严重的情况下，通过本发明方法获得的相关峰仍然非常陡峭，旁瓣较小，可见本发明的方法受噪声的影响较小。

另外，根据文献(New approach for improving Mutual Information in POLSAR image registration，Proceedings of IEEE International Conference on Signal Processing，2010：853-856)所述，相似性度量的峰越陡峭，匹配结果越精确，可知，本实施例中获得的相位相关平面对应的匹配结果非常精确。

步骤6：计算两幅极化SAR图像的最佳匹配位置(m₀,n₀)。

匹配位置和误差如表1所示，匹配拼接图如图4所示。图4为匹配透视图，图中使用绿色通道、红色通道分别独立显示两幅极化SAR图像。从匹配透视图中，可以清晰的看到两幅图像中每个像素的匹配情况。若匹配精确，匹配透视图将呈现清晰的细节纹理，否则将呈现混乱的纹理。本实施中将该匹配透视图灰度化显示。

表1

从表1可以看出，行和列的匹配误差均为0像素，匹配结果精确。图4中两幅极化SAR图像的重叠区域非常吻合，湖的边缘以及地形的纹理非常清晰，两幅图像得到了精确地匹配。另外，虽然两幅图像中一些局部区域差异较大(白色方框内的部分)，但并未对整体匹配结果造成影响，可见本发明的方法对局部几何失真具有较好的稳健性。