CN103439708B - 基于广义散射矢量的极化InSAR干涉图估计方法 - Google Patents

Abstract

本发明公开了一种基于广义散射矢量的极化InSAR干涉图估计方法，主要解决了现有的极化InSAR干涉图生成算法对配准误差不稳健的问题。本发明的实现步骤是：（1）输入图像数据；（2）图像粗配准；（3）构建广义散射矢量；（4）选取图像窗口；（5）估计广义相干矩阵；（6）估计广义干涉矩阵；（7）构建特征分解矩阵；（8）矩阵特征分解；（9）排列特征值；（10）生成干涉相位；（11）判断是否得到所有的干涉相位；（12）获得干涉相位图。本发明具有在极化通道和/或空间通道存在配准误差的情况下仍可获得高质量干涉图的优点，能有效减小配准误差对干涉图的影响，充分利用了每个像素点的所有极化信息。

Description

基于广义散射矢量的极化InSAR干涉图估计方法

技术领域

本发明属于通信技术领域，更进一步涉及雷达探测技术领域中的基于广义散射矢量的极化干涉合成孔径雷达（Synthetic Aperture Radar Interferomety，InSAR）干涉图估计方法。本发明对雷达图像中的像素点进行广义特征分解，再对原始干涉图加权来获得干涉相位图。本发明可用于极化通道和/或空间通道非完全配准情况下，极化InSAR干涉图的稳健估计。

背景技术

极化InSAR是现代雷达的研究热点，它既保持了InSAR的测量特点，又引入了极化合成孔径雷达（Synthetic Aperture Radar，SAR）的测量优势，所以能够同时获得地物的极化散射特性和干涉信息，并能够借助额外的信息来提高极化SAR与InSAR的应用性能。

干涉相位图是将雷达的主辅天线接分别收到信号的相位做干涉处理后，得到的一幅表示两个天线接收信号相位差的图像。极化相干最优分解是在极化InSAR的处理中，寻找一种极化状态组合，使相干系数在该极化状态组合下最大的过程。影响极化相干最优分解和干涉相位图生成的重要因素是图像配准误差，图像配准精度直接决定了干涉相位图结果的好坏。现有的基于Pauli基散射矢量的干涉相位图估计方法，当图像精确配准时，可以获得较好的干涉相位图估计；当图像之间存在配准误差时，干涉相位图估计性能下降，相位噪声增大；当存在较大的配准误差时，甚至不能获得正确的干涉图估计结果。

现有的极化InSAR的配准处理主要以一幅图像为参考图像，其余图像都与参考图像进行配准，并将配准结果应用到各个图像中。

熊涛在文章“极化干涉合成孔径雷达应用的关键技术研究”（清华大学博士论文，2009）中引入相似性参数对极化干涉数据进行配准，谈璐璐和杨汝良在文章“利用Cameron分解进行极化干涉图像对的配准研究”（系统工程与电子技术,2009,31(6):1284-1287）中也提出了利用Cameron分解来进行极化干涉图像的配准。上述方法相对于传统的配准方法配准精度有了较大提高。然而上述方法的不足之处是：在进行配准时，假设图像的极化通道之间是精确配准的，然而在通常情况下，由于受到不精确的天线相位中心位置和极化通道不完全正交等因素的影响，经过标定后的极化通道间仍然存在残余误差。

中国民航大学提出的专利申请“基于相关加权联合子空间投影的InSAR干涉相位估计方法”（申请日期：2012年01月09日，申请号：201210003415.5，公开号：CN102565798A）中公开了一种基于相关加权联合子空间投影的方法来稳健的估计干涉相位。该方法首先利用现有的相关法对主辅图像进行粗配准，然后构建出相关加权联合协方差矩阵和加权联合相关矩阵，最后利用噪声子空间投影构造出的代价函数来估计出干涉相位。该方法在配准精度很差的情况下可以得到准确的干涉相位结果。该方法存在的不足之处是：针对的只是单极化InSAR的数据处理，而在处理多极化InSAR数据时，因为极化通道和/或空间通道间存在残余误差，该方法会得到不理想的干涉相位图。

发明内容

本发明针对上述现有技术存在的极化InSAR干涉图生成技术对极化信息利用不充分不足的缺陷，提出了一种基于广义散射矢量的极化InSAR干涉图估计方法，通过建立广义散射模型，并基于该模型来降低极化通道和/或空间通道配准误差对极化相干最优分解及干涉相位图生成的影响，保证了在极化通道和/或空间通道间存在残余误差时，仍然可以得到质量较好的极化InSAR干涉相位图。

为实现上述目的，本发明的主要步骤如下：

（1）输入图像数据：

1a）将极化干涉合成孔径雷达主天线得到的主图像数据输入到系统；

1b）将极化干涉合成孔径雷达辅天线得到的辅图像数据输入到系统；

（2）图像粗配准：

利用InSAR粗配准方法，对主图像和辅图像进行粗配准处理；

（3）构建广义散射矢量：

3a）在主图像中任意选取一个像素点，将所选取的像素点与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出主图像中所选像素点的广义散射矢量；

3b）依次选取主图像中尚未构建广义散射矢量的像素点，将所有选取的像素点依次与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出主图像中所有选取的像素点的广义散射矢量；

3c）判断主图像中所有像素点的广义散射矢量是否全部构建完成，若全部构建完成，则执行步骤3d），否则，执行步骤3b）；

3d）在辅图像中任意选取一个像素点，将所选取的像素点与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出辅图像中所选像素点的广义散射矢量；

3e）依次选取辅图像中尚未构建广义散射矢量的像素点，将所有选取的像素点依次与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出辅图像中所有选取的像素点的广义散射矢量；

3f）判断辅图像中所有像素点的广义散射矢量是否全部构建完成，若全部构建完成，则执行步骤（4），否则，执行步骤3e）；

（4）选取图像窗口：

在主图像中，选取一个像素点作为中心，以固定长度为块半径，选取一个正方形的窗口，在辅图像中，选取与主图像中对应的像素点作为中心，获得同样大小的正方形窗口；

（5）估计广义相干矩阵：

5a）在主图像窗口中，将该窗口中所有像素点的广义散射矢量与自身的共轭转置相乘后取算术平均，得到主图像窗口中所选取像素点的广义相干矩阵；

5b）在辅图像窗口中，将该窗口中所有像素点的广义散射矢量与自身的共轭转置相乘后取算术平均，得到辅图像窗口中所选取像素点的广义相干矩阵；

（6）估计广义干涉矩阵：

6a）将辅图像窗口中所有像素点的广义散射矢量做共轭转置，获得共轭转置后的辅图像窗口中所有像素点的广义散射矢量；

6b）在共轭转置后的辅图像窗口中，将所有像素点的广义散射矢量与主图像窗口所对应的像素点的广义散射矢量相乘后取算术平均，得到广义干涉矩阵；

（7）构建特征分解矩阵：

采用下式，构建特征分解矩阵：

Χ=A^-1BC^-1B^H

其中，Χ表示特征分解矩阵，A表示主图像窗口中所选取像素点的广义相干矩阵，B表示辅图像窗口中所选取像素点的广义相干矩阵，C表示广义干涉矩阵，上标-1表示矩阵求逆，上标H表示对矩阵取共轭转置；

（8）矩阵特征分解：

8a）将特征分解矩阵进行特征分解操作，得到一列特征分解矩阵的乱序排列的特征值；

8b）将特征分解矩阵进行特征分解操作，得到一列与特征值一一对应的特征向量；

（9）排列特征值：

9a）将特征分解矩阵的乱序排列的特征值，按照从大到小进行排序，得到一列特征分解矩阵的顺序排列的特征值；

9b）将与特征值一一对应的特征向量，按照顺序排列的特征值的次序进行排序，得到一列特征分解矩阵的顺序排列的特征向量；

（10）生成干涉相位：

10a）利用顺序排列的特征值和特征向量，对广义干涉矩阵进行加权处理，生成与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位；

10b）利用顺序排列的特征值和特征向量，对广义干涉矩阵进行加权处理，生成与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位；

10c）利用顺序排列的特征值和特征向量，对广义干涉矩阵进行加权处理，生成与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位；

（11）判断是否得到所有的干涉相位：

遍历所有像素点，判断是否已得到所有的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位、与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位、与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，如果是，则执行步骤（12），否则，执行步骤（4）；

（12）获得干涉相位图：

12a）将所有像素点的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位图；

12b）将所有像素点的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位图；

12c）将所有像素点的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位图。

本发明与现有的技术相比具有以下优点：

第一，本发明构建了广义散射矢量，能够对当前估计像素点的所有极化信息自适应的恢复，克服了现有技术中基于相关加权联合子空间投影方法对极化信息利用不充分的不足，使得本发明具有极化信息利用充分、干涉相位图估计准确的优点。

第二，本发明对配准误差具有很好的鲁棒性，克服了现有技术利用Cameron分解进行极化干涉图像配准的方法中，极化通道和/或空间通道间存在配准误差时得到的干涉相位图质量差的不足，使得本发明具有在极化通道和/或空间通道存在配准误差的情况下仍可获得高质量干涉相位图的优点。

附图说明

图1为本发明的流程图；

图2为生成仿真数据所用的地形图；

图3为Pauli基散射矢量方法处理精确配准的仿真数据生成的干涉图；

图4为本发明处理精确配准的仿真数据生成的干涉图；

图5为Pauli基散射矢量方法处理存在1/10个像素误差的仿真数据生成的干涉图；

图6为本发明处理存在1/10个像素误差的仿真数据生成的干涉图；

图7为Pauli基散射矢量方法处理实测数据生成的干涉图；

图8为本发明处理实测数据生成的干涉图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照附图1，本发明的具体实施步骤如下：

步骤1，输入图像数据。

将极化干涉合成孔径雷达主辅天线分别得到的主图像数据和辅图像数据输入到系统中。输入的主图像数据和辅图像数据是对同一地区的成像结果，以此来保证主图像与辅图像之间的相干性。

步骤2，图像粗配准。

利用InSAR粗配准方法，对主图像和辅图像进行粗配准处理，在图像粗配准的处理过程中，配准精度只要求达到像素级，大大减轻了图像粗配准的难度。

InSAR粗配准具体步骤如下：

第一步，利用下式，计算主图像和辅图像的互相关矩阵：

L=IFFT 2(FFT 2(R)×conj(FFT 2(S)))

其中，L表示主图像和辅图像的互相关矩阵，IFFT2(·)表示做二维逆傅里叶变换，FFT2(·)表示做二维傅里叶变换，R表示主图像数据，S表示辅图像数据，conj(·)表示取共轭操作；

第二步，从互相关矩阵L的所有元素中寻找最大值，得到互相关矩阵L中最大值元素对应的坐标位置；

第三步，用互相关矩阵L中最大值元素的坐标位置的值，减去互相关矩阵L的中心坐标位置的值，得到配准时的偏移量；

第四步，对辅图像整体做与配准时的偏移量大小相等的平移，完成主辅图像的粗配准。

步骤3，构建广义散射矢量。

构建广义散射矢量是为了利用每一个像素点的所有极化信息，可以选择每一个像素与其周围相邻的1～8个像素点来构建广义散射矢量，本发明实施例中选取了每一个像素点与其周围相邻的8个像素点来构建广义散射矢量。

构建广义散射矢量的具体步骤为：

3a）在主图像中任意选取一个像素点，将所选取的像素点与其周围相邻的8个像素点的保利Pauli基散射矢量中元素排成一列，构建出主图像中所选像素点的广义散射矢量；

3b）依次选取主图像中尚未构建广义散射矢量的像素点，将所有选取的像素点依次与其周围相邻的8个像素点的保利Pauli基散射矢量中元素排成一列，构建出主图像中所有选取的像素点的广义散射矢量；

3c）判断主图像中所有像素点的广义散射矢量是否全部构建完成，若全部构建完成则执行步骤3d），否则，执行步骤3b）；

3d）在辅图像中任意选取一个像素点，将所选取的像素点与其周围相邻的8个像素点的保利Pauli基散射矢量中元素排成一列，构建出辅图像中所选像素点的广义散射矢量；

3e）依次选取辅图像中尚未构建广义散射矢量的像素点，将所有选取的像素点依次与其周围相邻的8个像素点的保利Pauli基散射矢量中元素排成一列，构建出辅图像中所有选取的像素点的广义散射矢量；

3f）判断辅图像中所有像素点的广义散射矢量是否全部构建完成，若全部构建完成则执行步骤4，否则，执行步骤3e）。

步骤4，选取图像窗口。

在主图像中，选取一个像素点t作为中心，以固定长度3～7个像素点为块半径，选取一个正方形的窗口，在辅图像中，选取与主图像中对应的像素点作为中心，获得同样大小的正方形窗口，本发明实施例中选取了大小为9×9的正方形窗口。

步骤5，估计广义相干矩阵。

在通常情况下，无法直接得到广义相干矩阵的统计平均值，所以利用广义相干矩阵的空间样本平均代替统计平均。

估计广义相干矩阵的具体步骤为：

第一步，在主图像窗口中，将该窗口中所有像素点的广义散射矢量与自身的共轭转置相乘后取算术平均，得到主图像窗口中所选取像素点的广义相干矩阵Τ；

第二步，在辅图像窗口中，将该窗口中所有像素点的广义散射矢量与自身的共轭转置相乘后取算术平均，得到辅图像窗口中所选取像素点的广义相干矩阵Υ。

步骤6，估计广义干涉矩阵。

在通常情况下，无法直接得到广义干涉矩阵的统计平均值，所以利用广义干涉矩阵的空间样本平均代替统计平均。

估计广义干涉矩阵的具体步骤为：

第一步，将辅图像窗口中所有像素点的广义散射矢量做共轭转置，获得共轭转置后的辅图像窗口中所有像素点的广义散射矢量；

第二步，在共轭转置后的辅图像窗口中，将所有像素点的广义散射矢量与主图像窗口所对应的像素点的广义散射矢量相乘后取算术平均，得到广义干涉矩阵Μ。

步骤7，构建特征分解矩阵。

采用下式，构建特征分解矩阵：

Χ=Τ^-1ΥΜ^-1Υ^H

其中，Χ表示特征分解矩阵，Τ表示主图像窗口所选取像素点的广义相干矩阵，Υ表示辅图像窗口所选取像素点的广义相干矩阵，Μ表示广义干涉矩阵，上标-1表示矩阵求逆，H表示对矩阵取共轭转置。

步骤8，矩阵特征分解。

将特征分解矩阵X进行特征分解操作，得到27个特征分解矩阵的乱序排列的特征值和27个与特征值一一对应的特征向量。

步骤9，排列特征值。

第一步，将特征分解矩阵的乱序排列的特征值，按照从大到小进行排序，得到27个特征分解矩阵的顺序排列的特征值λ₁,λ₂,…,λ₂₇；

第二步，将27个与特征值一一对应的特征向量，按照顺序排列的特征值的次序进行排序，得到27个特征分解矩阵的顺序排列的特征向量ω₁,ω₂,…,ω₂₇。

步骤10，生成干涉相位。

在实际处理中，由于受到噪声等因素的影响，无法直接得到精确的特征值与特征向量，即无法直接得到精确的干涉相位，故而本发明采用特征值和特征向量加权的方法来估计干涉相位。

生成干涉相位的具体步骤如下：

第一步，采用下式，对像素点t，生成与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位：

$\mathrm{\Φ}=\arg (\underset{i=1}{\overset{9}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\mathrm{\ω}}_i^{H}M{\mathrm{\ω}}_i/\underset{i=1}{\overset{9}{\mathrm{\Σ}}}{\mathrm{\λ}}_i)$

其中，Φ表示像素点t的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位，λ_i分别表示特征分解矩阵的顺序排列的特征值中第1个到第9个特征值，ω_i分别表示特征分解矩阵的顺序排列的特征向量中第1个到第9个特征向量，i=1,2,…9，Μ表示广义干涉矩阵，arg(·)表示取相位操作，H表示对矩阵取共轭转置；

第二步，采用下式，对像素点t，生成与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位：

$\mathrm{\Ψ}=\arg (\underset{i=10}{\overset{18}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\mathrm{\ω}}_i^{H}M{\mathrm{\ω}}_i/\underset{i=10}{\overset{18}{\mathrm{\Σ}}}{\mathrm{\λ}}_i)$

其中，Ψ表示像素点t的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位，λ_i分别表示特征分解矩阵的顺序排列的特征值中第10个到第18个特征值，ω_i分别表示特征分解矩阵的顺序排列的特征向量中第10个到第18个特征向量，i=10,11,…18，Μ表示广义干涉矩阵，arg(·)表示取相位操作，上标H表示对矩阵取共轭转置；

第三步，采用下式，对像素点t，生成与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位：

$\mathrm{\Ω}=\arg (\underset{i=19}{\overset{27}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\mathrm{\ω}}_i^{H}M{\mathrm{\ω}}_i/\underset{i=19}{\overset{27}{\mathrm{\Σ}}}{\mathrm{\λ}}_i)$

其中，Ω表示像素点t的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，λ_i分别表示特征分解矩阵的顺序排列的特征值中第19个到第27个特征值，ω_i分别表示特征分解矩阵的顺序排列的特征向量中第19个到第27个特征向量，i=19,20,…27，Μ表示广义干涉矩阵，arg(·)表示取相位操作，H表示对矩阵取共轭转置。

步骤11，判断是否得到所有的干涉相位。

遍历所有像素点，判断是否已得到所有的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位Φ、与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位Ψ、与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位Ω，如果是，则执行步骤12，否则，执行步骤4。

步骤12，获得干涉相位图。

第一步，将所有像素点的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位图；

第二步，将所有像素点的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位图；

第三步，将所有像素点的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位图。

下面结合仿真实验对本发明的效果做进一步的说明。

1、仿真条件：

参照附图2，利用软件PolSARPro，对附图2中的场景区域进行极化InSAR仿真成像，分别得到精确配准的主辅图像数据和存在1/10像素配准误差的主辅图像数据，仿真参数设置如下：

参数	参数值
		天线个数	2
主天线斜距	4242.6406m
		主天线下视角	45°
辅天线斜距	4250.4236m
		辅天线下视角	45.0858°
中心频率	1.3GHz
		方位向分辨率	1.5m
斜距分辨率	1.0607m
		平均植被高度	10.0m
每公顷植被密度	900
		方位向地面坡度	0.02
距离向地面坡度	0.01

2、仿真结果分析：

图3（a）表示Pauli基散射矢量方法处理精确配准的仿真数据生成的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位，横坐标轴表示距离向，单位为像素，纵坐标轴表示方位向，单位为像素，左侧的图像中像素的灰度值表示地面散射对应像素点的干涉相位的大小，图中干涉相位波动的大小，反应出了相位噪声的多少，右侧的图像条表示干涉相位值与左侧图像中灰度值的对应关系，右侧图像的数值表示干涉相位的大小。图3（b）表示Pauli基散射矢量方法处理精确配准的仿真数据生成的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位。图3（c）表示Pauli基散射矢量方法处理精确配准的仿真数据生成的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位。图4（a）表示本发明处理精确配准的仿真数据生成的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位。图4（b）表示本发明处理精确配准的仿真数据生成的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位。图4（c）表示本发明处理精确配准的仿真数据生成的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位。图5（a）表示Pauli基散射矢量方法处理存在1/10个像素误差的仿真数据生成的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位。图5（b）表示Pauli基散射矢量方法处理存在1/10个像素误差的仿真数据生成的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位。图5（c）表示Pauli基散射矢量方法处理存在1/10个像素误差的仿真数据生成的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位。图6（a）表示本发明处理存在1/10个像素误差的仿真数据生成的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位。图6（b）表示本发明处理存在1/10个像素误差的仿真数据生成的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位。图6（c）表示本发明处理存在1/10个像素误差的仿真数据生成的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位。

对比附图3和附图4，附图4中的干涉相位的波动要比附图3的小，尤其是在有植被覆盖的区域显现的尤为明显，从而可以得到，当主辅图像精确配准时，本发明得到的干涉相位图质量优于Pauli基散射矢量方法。对比附图5和附图6，附图5中的相位值的波动很大，附图6中的相位波动较为平稳，从而可以得到，当主辅图像的配准误差为1/10个像素时，本发明得到的干涉相位图质量明显高于Pauli基散射矢量方法。对比附图4和附图6，附图6中的干涉相位值的波动与附图4中的相位波动相近，从而可以得到，本发明对配准误差有很好的鲁棒性，在主辅图像间存在配准误差的情况下，仍然可以得到质量较好的干涉相位图。综合上述分析，可以得到本发明在处理极化InSAR数据时，得到的结果要优于现有Pauli基散射矢量方法，并且对配准误差有很好的鲁棒性。

3、实测数据处理

处理实测数据时，使用的数据是日本ALOS卫星PALSAR录取的Amazon区域得到的极化InSAR数据。

图7（a）表示Pauli基散射矢量方法处理实测数据生成的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位。图7（b）表示Pauli基散射矢量方法处理实测数据生成的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位。图7（c）表示Pauli基散射矢量方法处理实测数据生成的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位。图8（a）表示本发明处理实测数据生成的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位。图8（b）表示本发明处理实测数据生成的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位。图8（c）表示本发明处理实测数据生成的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位。

对比附图7和附图8，附图8（a）中的相位波动小于附图7（a），附图8（b）中的相位波动明显小于附图7（b）中的相位波动，附图7（c）中的相位波动非常大，含有非常多的相位噪声，而附图8（c）中的相位波动要比附图7（c）中相位波动小得多，可以看到较为清晰的干涉条纹。由上述分析可知，本发明处理配准误差未知的极化InSAR实测数据时，仍然可以得到比现有处理方法质量更好的干涉相位图，进一步验证了本发明对配准误差有很好的鲁棒性，并且说明了本发明的实用性和有效性。

Claims (7)

1.基于广义散射矢量的极化InSAR干涉图估计方法，包括如下步骤：

(1)输入图像数据：

1a)将极化干涉合成孔径雷达主天线得到的主图像数据输入到系统；

1b)将极化干涉合成孔径雷达辅天线得到的辅图像数据输入到系统；

(2)图像粗配准：

利用InSAR粗配准方法，对主图像和辅图像进行粗配准处理；

(3)构建广义散射矢量：

3a)在主图像中任意选取一个像素点，将所选取的像素点与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出主图像中所选像素点的广义散射矢量；

3b)依次选取主图像中尚未构建广义散射矢量的像素点，将所有选取的像素点依次与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出主图像中所有选取的像素点的广义散射矢量；

3c)判断主图像中所有像素点的广义散射矢量是否全部构建完成，若全部构建完成，则执行步骤3d)，否则，执行步骤3b)；

3d)在辅图像中任意选取一个像素点，将所选取的像素点与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出辅图像中所选像素点的广义散射矢量；

3e)依次选取辅图像中尚未构建广义散射矢量的像素点，将所有选取的像素点依次与其周围相邻像素点的保利Pauli基散射矢量中元素排成一列，构建出辅图像中所有选取的像素点的广义散射矢量；

3f)判断辅图像中所有像素点的广义散射矢量是否全部构建完成，若全部构建完成，则执行步骤(4)，否则，执行步骤3e)；

(4)选取图像窗口：

在主图像中，选取一个像素点作为中心，以固定长度为块半径，选取一个正方形的窗口，在辅图像中，选取与主图像中对应的像素点作为中心，获得同样大小的正方形窗口；

(5)估计广义相干矩阵：

5a)在主图像窗口中，将该窗口中所有像素点的广义散射矢量与自身的共轭转置相乘后取算术平均，得到主图像窗口中所选取像素点的广义相干矩阵；

5b)在辅图像窗口中，将该窗口中所有像素点的广义散射矢量与自身的共轭转置相乘后取算术平均，得到辅图像窗口中所选取像素点的广义相干矩阵；

(6)估计广义干涉矩阵：

6a)将辅图像窗口中所有像素点的广义散射矢量做共轭转置，获得共轭转置后的辅图像窗口中所有像素点的广义散射矢量；

6b)在共轭转置后的辅图像窗口中，将所有像素点的广义散射矢量与主图像窗口所对应的像素点的广义散射矢量相乘后取算术平均，得到广义干涉矩阵；

(7)构建特征分解矩阵：

采用下式，构建特征分解矩阵：

Χ＝A^-1BC^-1B^H

其中，Χ表示特征分解矩阵，A表示主图像窗口中所选取像素点的广义相干矩阵，B表示辅图像窗口中所选取像素点的广义相干矩阵，C表示广义干涉矩阵，上标-1表示矩阵求逆，上标H表示对矩阵取共轭转置；

(8)矩阵特征分解：

8a)将特征分解矩阵进行特征分解操作，得到一列特征分解矩阵的乱序排列的特征值；

8b)将特征分解矩阵进行特征分解操作，得到一列与特征值一一对应的特征向量；

(9)排列特征值：

9a)将特征分解矩阵的乱序排列的特征值，按照从大到小进行排序，得到一列特征分解矩阵的顺序排列的特征值；

9b)将与特征值一一对应的特征向量，按照顺序排列的特征值的次序进行排序，得到一列特征分解矩阵的顺序排列的特征向量；

(10)生成干涉相位：

10a)利用顺序排列的特征值和特征向量，对广义干涉矩阵进行加权处理，生成与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位；

10b)利用顺序排列的特征值和特征向量，对广义干涉矩阵进行加权处理，生成与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位；

10c)利用顺序排列的特征值和特征向量，对广义干涉矩阵进行加权处理，生成与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位；

(11)判断是否得到所有的干涉相位：

遍历所有像素点，判断是否已得到所有的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位、与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位、与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，如果所有像素点的三种干涉相位都已得到，则执行步骤(12)，否则，执行步骤(4)；

(12)获得干涉相位图：

12a)将所有像素点的与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位图；

12b)将所有像素点的与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位图；

12c)将所有像素点的与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，保存到与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位图中相应的位置，得到与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位图。

2.根据权利要求1所述的基于广义散射矢量的极化InSAR干涉图估计方法，其特征在于：步骤(2)中所述的InSAR粗配准方法的具体步骤如下：

第一步，利用下式，计算主图像和辅图像的互相关矩阵：

L＝IFFT2(FFT2(R)×conj(FFT2(S)))

其中，L表示主图像和辅图像的互相关矩阵，IFFT2(·)表示做二维逆傅里叶变换，FFT2(·)表示做二维傅里叶变换，R表示主图像数据，S表示辅图像数据，conj(·)表示取共轭操作；

第二步，从互相关矩阵L的所有元素中寻找最大值，得到互相关矩阵L中最大值元素对应的坐标位置；

第三步，用互相关矩阵L中最大值元素的坐标位置的值，减去互相关矩阵L的中心坐标位置的值，得到配准时的偏移量；

第四步，对辅图像整体做与配准时的偏移量大小相等的平移，完成主辅图像的粗配准。

3.根据权利要求1所述的基于广义散射矢量的极化InSAR干涉图估计方法，其特征在于：步骤3a)中用来构建广义散射矢量的相邻像素点的个数为1～8个。

4.根据权利要求1所述的基于广义散射矢量的极化InSAR干涉图估计方法，其特征在于：步骤(4)中的固定长度为3～7个像素点。

5.根据权利要求1所述的基于广义散射矢量的极化InSAR干涉图估计方法，其特征在于：步骤10a)中利用下式，生成与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位：

$\mathrm{\Φ}=\arg (\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\mathrm{\ω}}_i^{H}B{\mathrm{\ω}}_i/\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i)$

其中，Φ表示与极化干涉合成孔径雷达照射区域中地面目标对应的干涉相位，m表示构建广义散射矢量时所使用的像素点个数，λ_i分别表示特征分解矩阵的顺序排列的特征值中第一个到第m个特征值，ω_i分别表示特征分解矩阵的顺序排列的特征向量中第一个到第m个特征向量，i＝1,2,…m，B表示广义干涉矩阵，arg(·)表示取相位操作，H表示对矩阵取共轭转置。

6.根据权利要求1所述的基于广义散射矢量的极化InSAR干涉图估计方法，其特征在于：步骤10b)中利用下式，生成与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位：

$\mathrm{\Ψ}=\arg (\underset{i=m+1}{\overset{2m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\mathrm{\ω}}_i^{H}B{\mathrm{\ω}}_i/\underset{i=m+1}{\overset{2m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i)$

其中，Ψ表示与极化干涉合成孔径雷达照射区域中介于地面和树冠之间目标对应的干涉相位，m表示构建广义散射矢量时所使用的像素点个数，λ_i分别表示特征分解矩阵的顺序排列的特征值中第m+1个到第2m个特征值，ω_i分别表示特征分解矩阵的顺序排列的特征向量中第m+1个到第2m个特征向量，i＝1,2,…m，B表示广义干涉矩阵，arg(·)表示取相位操作，上标H表示对矩阵取共轭转置。

7.根据权利要求1所述的基于广义散射矢量的极化InSAR干涉图估计方法，其特征在于：步骤10c)中利用下式，生成与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位：

$\mathrm{\Ω}=\arg (\underset{i=m+1}{\overset{2m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{\mathrm{\ω}}_i^{H}B{\mathrm{\ω}}_i/\underset{i=m+1}{\overset{2m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i)$

其中，Ω表示与极化干涉合成孔径雷达照射区域中树冠目标对应的干涉相位，m表示构建广义散射矢量时所使用的像素点个数，λ_i分别表示特征分解矩阵的顺序排列的特征值中第2m+1个到第3m个特征值，ω_i分别表示特征分解矩阵的顺序排列的特征向量中第2m+1个到第3m个特征向量，i＝1,2,…m，B表示广义干涉矩阵，arg(·)表示取相位操作，H表示对矩阵取共轭转置。