CN101685155A - 基于极化干涉合成孔径雷达数据优化干涉相干系数的方法 - Google Patents
基于极化干涉合成孔径雷达数据优化干涉相干系数的方法 Download PDFInfo
- Publication number
- CN101685155A CN101685155A CN200810223609A CN200810223609A CN101685155A CN 101685155 A CN101685155 A CN 101685155A CN 200810223609 A CN200810223609 A CN 200810223609A CN 200810223609 A CN200810223609 A CN 200810223609A CN 101685155 A CN101685155 A CN 101685155A
- Authority
- CN
- China
- Prior art keywords
- coherence
- phase
- omega
- coefficient
- optimum
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提供一种基于极化干涉合成孔径雷达数据优化干涉相干系数的方法,涉及合成孔径雷达技术,该方法先利用快速傅里叶变换计算各极化状态的平地相位;去除各极化状态的平地相位对极化干涉数据的影响;对极化干涉数据平滑滤波;利用特征分解提取最优投影矢量对;计算最优干涉相干系数与干涉相位,检验最优估计的有效性;将平地相位补偿至最优干涉相位中,按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位。本发明提供方法不受极化干涉数据矩阵可逆的条件限制,无需补偿最优投影矢量间的绝对相位差,并可检验优化结果的有效性。本发明的方法原理简单、计算方便,可用于极化干涉分类、参数反演等数据处理中。
Description
技术领域
本发明涉及合成孔径雷达遥感信息处理技术领域,是一种基于单基线极化干涉合成孔径雷达数据的优化干涉相干系数的方法,它能简单、快速实现干涉相干系数的优化并提取对应于散射中心的干涉相位。
背景技术
极化干涉合成孔径雷达技术将全极化数据引入干涉处理,不但可利用极化信息对散射体形状、指向、介电系数等参数的敏感性[1],而且还可基于干涉相位获得高度、运动、形变信息,因此它比单独的干涉或极化SAR在信息提取、分类等应用中具有更突出的优势,现已成为SAR领域的研究热点之一。
参考文献
[1]Papathanassiou K.P.and Cloude S.R.(2001),Single-Baseline Polarimetric SAR Interferometry,IEEE Trans.onGeoscience and Remote Sensing,39(11):2352-2363.
[2]Cloude S.R.and Papathanassiou K.P.(1998),PolarimetricSAR Interferometry[J],IEEE Trans.on Geoscience and RemoteSensing,36(5):1551-1565.
[3]Maxim Neumann,et al.(2006),PolInSAR Coherence SetTheory and Application,EUSAR 2006,Dresden:Germany.
[4]Papathanassiou K.P.(1999),Polarimetric SARInterferometry,Doctor Dissertation,Institut fürHochfrequenztechnik des DLR,Oberpfaffenhofen.
[5]Colin C.,et al.(2003),Investigation on DifferentCoherence Optimization Methods,the 1st Internaltional Workshop onScience and Applications of SAR Polarimetry and Polarimetric SARInterferometry,Frascati:ESA.
[6]Sagues L.,et al.(2000),Indoor Experiments onPolarimetric SAR Interferometry,IEEE Trans.on Geoscience andRemote Sensing,38(2):671-684.
[7]Colin C.,et al.(2006),An interferometric coherenceopmtimization method in radar polarimetry for high-resolutionimagery,IEEE Trans.on Geoscience and Remote Sensing,44(1):167-175.
[8]Lu Bai,.et al.(2006),An Improved Coherence OptimizationMethod in Polarimetric SAR Interferometry,InternationalConference on Radar,Shanghai:China.
[9]Tao Xiong,et al.(2006),A New Approach for DEM GenerationBased on Polarimetric SAR Interferometry,International Conferenceon Radar,Shanghai:China.
作为极化干涉SAR数据处理中一项关键技术,优化相干性的方法利用干涉相干系数对极化状态的依赖性[2],通过改变极化状态优化干涉相干系数、提高干涉相位的质量。从数学的角度讲,优化相干性是在所有复矢量或者极化状态的空间中寻找使干涉相干系数达到最大的一组投影矢量或极化状态的过程[3]:
优化相干性在改善干涉性能的同时,还可提取出最接近于点状散射体的多个确定性散射机理,从而将高度上分布的有效散射相位中心等效为点处理,抑制体散射去相干。该方法通过选择与散射中心散射机理有关的最优投影矢量与补偿时间去相干、优化信噪比、改善噪声去相干[4]。因此,基于单基线极化干涉数据的优化相干性的方法在分类、参数反演应用中具有巨大的应用潜力。
Papathanassiou于1997年引入Lagrange乘数法求解最优投影矢量,该方法可采用两次特征分解或一次奇异值分解即可求解除最优投影矢量,然而该算法应用前提为三维相干矩阵T11、T22可逆[2,4]。由于特征矢量绝对相位不唯一,因而还需引入额外的限定条件才能避免主、副天线选用最优投影矢量对干涉相位估计的影响。由于该方法未对投影矢量作严格的限定,因而造成临近的位置的最优矢量间变化程度较大,导致干涉相位不连续,说明其提取的干涉相位不一定准确[5]。
由于主要散射中心的散射机理通常是稳定不变的,因此应选用相同的极化状态或投影矢量计算干涉相干系数,从而出现了约束投影矢量对内两矢量相同的优化相干性的方法。与之前的极化全局空间的优化方法相比,约束类方法解决的是极化子空间的优化问题。Sagues从不同极化基中搜索出具有最高干涉相干性的极化状态作为优化结果,该方法计算量大且其估计精度受搜索步长的影响[6]。
Colin从Lagrange函数求解出发,通过引入完全相干假设将最优投影矢量的求解过程简化为一次特征矢量的求解过程。该算法要求矩阵(T11+T22)可逆,且受到完全相干假设的限制,其适用于处理基线去相关、时间去相关影响较小的数据[5]。另一种约束类相干性优化方法,利用近似关系T11≈T22≈(T11+T22)/2在矩阵(T11+T22)可逆条件下将优化问题转换为数值距离的计算;该算法同样受到(T11+T22)可逆的限制,且由于该算法利用数值逼近求解,因而其估计结果的准确性与初始估计、搜索步长密切相关[7]。此外,Lu Bai尝试利用特征分解结果作为初始估计,并利用干涉相干性能否改善作为结果收敛条件,但是该方法仍要求矩阵(T11+T22)可逆[8]。Tao Xiong比较HH、HV、VV三种极化状态的干涉结果中具有最大干涉相干性的干涉相位作为最优结果[9],虽然该算法运算效率高,但其干涉相干性的改善程度有限。
本发明方法通过分析一个极化干涉数据矩阵提取最优投影矢量对,且不要求该极化干涉数据矩阵可逆,因而普遍适用于各类极化干涉数据。本发明方法无需补偿最优投影矢量间的绝对相位差,简化了处理;并通过检验优化结果的有效性,可确保优化后干涉相干系数以及干涉相位结果的可靠性。
发明内容
本发明的目的在公开一种基于极化干涉合成孔径雷达数据优化干涉相干系数的方法,该方法计算速度快、相干性改善程度大、适用于对不同类型的极化干涉数据的优化相干系数。
为达到上述目的,本发明的技术解决方案是:
基于极化干涉合成孔径雷达数据优化干涉相干系数的方法,其包括步骤如下:
A)利用快速傅里叶变换计算各极化状态的平地相位;
B)去除各极化状态的平地相位对极化干涉数据的影响;
C)对极化干涉数据平滑滤波;
D)利用特征分解提取最优投影矢量对;
E)计算最优干涉相干系数与干涉相位,检验最优估计的有效性;
F)将平地相位补偿至最优干涉相位中,按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位。
所述的方法,其所述A)步先利用快速傅里叶变换计算各极化状态的平地相位,其各极化状态的平地相位计算包含如下步骤:
A1、计算指定范围内m1≤m≤m2,各极化状态下的干涉相位φXX、φXY、φYX、φYY:
φXX(m,n)=Arg(s1XX(m,n)·s* 2XX(m,n))
φXY(m,n)=Arg(s1XY(m,n)·s* 2XY(m,n))
φYX(m,n)=Arg(s1YX(m,n)·s* 2YX(m,n))
φYY(m,n)=Arg(s1YY(m,n)·s* 2YY(m,n))
其中s1表示主天线接收到的复图像数据;s2表示副天线接收到的复图像数据;(m,n)指示像素在图像中的行、列位置;m1、m2表示用于计算平地相位的数据所在行的起始位置、终止位置;XX、XY、YX、YY表示复图像数据的极化方式,符号.表示复数乘运算,*表示共轭运算,Arg表示复数的辐角主值。
A2、利用快速傅立叶变换计算指定范围内m1≤m≤m2,各极化状态的复干涉相位序列的频谱Xm,XX、Xm,XY、Xm,YX、Xm,YY:
其中φXX(m,..)表示由φXX(m,1)、φXX(m,2)、...φXX(m,N)构成的序列;
A3、计算指定范围内m1≤m≤m2各极化状态复相位序列频谱能量的平均值|XXX|、|XXY|、|XYX|、|XYY|:
A4、确定各极化状态复相位序列频谱能量的平均值的最大值位置KPeak,XX、KPeak,XY、KPeak,YX、KPeak,YY;
A5、计算各极化状态平地相位序列φflat,XX、φflat,XY、φflat,YX、φflat,YX:
其中N2表示快速傅里叶变换长度。
所述的方法,其所述B)步去除各极化状态的平地相位对极化干涉数据的影响,是去除平地相位的影响可减小相邻位置干涉相位梯度,增强相邻位置的极化干涉数据的一致性,具体以下式计算:
。
所述的方法,其所述C)步极化干涉数据平滑滤波,包含以下步骤:
C1、以待滤波位置为窗口中心,大小为wm×wn的滤波窗口内,生成包含待滤波位置的,沿行、列以及两对角线方向的,中心位置不同的滤波子窗口;
C2、计算各个子窗口内极化干涉数据总能量的方差σ2 L(m,n):
σ2 L(m,n)=<|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2>(m,m),L 2
-<(|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2)2>(m,n),L
其中σ2 L(m,n)表示以(m,n)中心的第L个子窗口内极化干涉数据总能量的方差;符号||表示复数的模;符号<>(m,n),L以(m,n)中心的第L个子窗口内极化干涉数据的算术平均;
C3、以极化干涉数据总能量方差最小为原则确定滤波子窗口LF;
C4、计算选定的滤波子窗口内极化干涉数据的算术平均T11(m,n)、T22(m,n)、Ω12(m,n):
D1、计算矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)的单位特征矢量:
其中,λ1、λ2、λ2是矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)的三个特征矢量; 分别是矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)对应于特征值λi的单位化特征矢量;
所述的方法,其所述E)步计算最优干涉相干系数与干涉相位,检验最优估计的有效性,包含以下步骤:
E1、利用最优投影矢量对计算最优干涉相干系数γi(m,n)与干涉相位φ′i(m,n):
;
E2、检验干涉相干系数值的合理性,设置错误估计标记:
所述的方法,其所述F)步将平地相位补偿至最优干涉相位中,是按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位,具体包含如下步骤:
F1、将平地相位φflat,XX(n)补偿至最优干涉相位φ′i(m,n)中,获得干涉相位φi(m,n):
F2、按照γ1(m,n)≥γ2(m,n)≥γ3(m,n)的顺序,排列生成三组最优干涉相干系数以及干涉相位结果{γ1(m,n),φ1(m,n)}、{γ2(m,n),φ2(m,n)}、{γ3(m,n),φ3(m,n)};其中第一最优结果的干涉相干系数γ1(m,n)最高,明显优于优化前的XX、XY、YX、YY极化的干涉相干系数。
本发明的方法计算速度快、相干性改善程度大、适用于对不同类型的极化干涉数据的优化相干系数。
附图说明
图1本发明中所使用的极化干涉数据示意图;
图2本发明方法的处理流程图;
图3a设定5×5窗口时统一滤波窗口示意图;
图3b设定5×5窗口时Lee自适应滤波时采用的滤波子窗口示意图;
图3c设定5×5窗口时本发明中在所选用的滤波子窗口示意图;
图4a SIR-C系统L波段数据去除平地相位影响前HH极化干涉相位量化图;
图4b SIR-C系统L波段数据去除平地相位影响前HV极化干涉相位量化图;
图4c SIR-C系统L波段数据去除平地相位影响前VH极化干涉相位量化图;
图4d SIR-C系统L波段数据去除平地相位影响前VV极化干涉相位量化图;
图4e图4a-图4d中干涉相位的量化方式示意图;
图5a SIR-C系统L波段数据计算的HH极化平地相位量化图;
图5b SIR-C系统L波段数据计算的HV极化平地相位量化图;
图5c SIR-C系统L波段数据计算的VH极化平地相位量化图;
图5d SIR-C系统L波段数据计算的VV极化平地相位量化图;
图5e图5a-图5d中干涉相位的量化方式示意图;
图6a SIR-C系统L波段数据去除平地相位影响后HH极化干涉相位量化图;
图6b SIR-C系统L波段数据去除平地相位影响后HV极化干涉相位量化图;
图6c SIR-C系统L波段数据去除平地相位影响后VH极化干涉相位量化图;
图6d SIR-C系统L波段数据去除平地相位影响后VV极化干涉相位量化图;
图6e图6a-图6d中干涉相位的量化方式示意图;
图7a SIR-C系统L波段数据经本发明方法处理获得第一最优干涉相干系数的量化图;
图7b SIR-C系统L波段数据经本发明方法处理获得第二最优干涉相干系数的量化图;
图7c SIR-C系统L波段数据经本发明方法处理获得第三最优干涉相干系数的量化图;
图7d图7a-图7c中干涉相干系数的量化方式示意图;
图8a SIR-C系统L波段数据经本发明方法处理获得第一最优干涉相干系数的量化图;
图8b SIR-C系统L波段数据经本发明方法处理获得第二最优干涉相干系数的量化图;
图8c SIR-C系统L波段数据经本发明方法处理获得第三最优干涉相干系数的量化图;
图8d图8a-图8c中干涉相干系数的量化方式示意图;
图9a SIR-C系统L波段数据用5×5滤波窗口计算的HH极化干涉相干系数的量化图;
图9b SIR-C系统L波段数据采用5×5滤波窗口计算的HV极化干涉相干系数的量化图;
图9c SIR-C系统L波段数据采用5×5滤波窗口计算的VH极化干涉相干系数的量化图;
图9d SIR-C系统L波段数据采用5×5滤波窗口计算的VV极化干涉相干系数的量化图;
图9e图9a-图9b中干涉相干系数的量化方式示意图;
图9d图9a-图9b最优干涉相位量化方式示意图;
图10SIR-C系统L波段数据的经本方法优化干涉相干系数与未经本方法优化的各极化干涉相干系数的对比图。
具体实施方式
如图2所示,基于极化干涉合成孔径雷达数据优化干涉相干系数的方法包括步骤如下:
A)利用快速傅里叶变换计算各极化状态的平地相位;
B)去除各极化状态的平地相位对极化干涉数据的影响;
C)对极化干涉数据平滑滤波;
D)利用特征分解提取最优投影矢量对;
E)计算最优干涉相干系数与干涉相位,检验最优估计的有效性;
F)将平地相位补偿至最优干涉相位中,按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位。
步骤A)利用快速傅里叶变换计算各极化状态的平地相位,是为去除平地相位造成的相邻位置上数据的波动。考虑到不同极化状态下干涉相位的差异,因此分别基于各极化下干涉相位估计平地相位;为提高计算效率,仅选用干涉相位变化规律的部分干涉相位数据计算即可。具体步骤如下:
A1、计算指定范围内m1≤m≤m2,各极化状态下的干涉相位φXX、φXY、φYX、φYY:
φXX(m,n)=Arg(s1XX(m,n)·s* 2XX(m,n))
φXY(m,n)=Arg(s1XY(m,n)·s* 2XY(m,n))
φYX(m,n)=Arg(s1YX(m,n)·s* 2YX(m,n))
φYY(m,n)=Arg(s1YY(m,n)·s* 2YY(m,n))
其中s1表示主天线接收到的复图像数据;s2表示副天线接收到的复图像数据;(m,n)指示像素在图像中的行、列位置;m1、m2表示用于计算平地相位的数据所在行的起始位置、终止位置;XX、XY、YX、YY表示复图像数据的极化方式,符号.表示复数乘运算,*表示共轭运算,Arg表示复数的辐角主值。
A2、利用快速傅立叶变换计算指定范围内m1≤m≤m2,各极化状态的复干涉相位序列的频谱Xm,XX、Xm,XY、Xm,YX、Xm,YY:
其中φXX(m,..)表示由φXX(m,1)、φXX(m,2)、...φXX(m,N)构成的序列。
A3、计算指定范围内m1≤m≤m2各极化状态复相位序列频谱能量的平均值|XXX|、|XXY|、|XYX|、|XYY|:
A4、确定各极化状态复相位序列频谱能量的平均值的最大值位置KPeak,XX、KPeak,XY、KPeak,YX、KPeak,YY。
A5、计算各极化状态平地相位序列φflat,XX、φflat,XY、φflat,YX、φflat,YX:
其中N2表示快速傅里叶变换长度。
步骤B)去除各极化状态的平地相位对极化干涉数据的影响,是为减小相邻位置干涉相位梯度,增强相邻位置的极化干涉数据的一致性。为了避免平地相位去除对主天线极化数据的影响,采用去除平地相位对副天线的全极化数据的处理方法。由于不同极化下平地相位存在差异,因而需对各极化状态的数据分别处理。具体步骤如下:
步骤C)极化干涉数据平滑滤波,是为抑制噪声影响。若采用如图3所示的统一矩形窗口滤波,能造成非同质、不均匀数据之间的平滑。为保证用于平滑滤波的极化干涉数据是均匀的,因而需根据极化干涉数据的波动性自适应选择具有方向性的如图3b所示的滤波子窗口。考虑到相邻位置数据的一致性高的特点,使得滤波子窗口能包含尽量多的与滤波位置临近的数据,本方法通过平移图3b所示的子窗口,提供了更丰富的滤波子窗口。由于在极化干涉数据滤波结果中补偿平地相位比较繁琐,因而暂不补偿平地相位。极化干涉数据平滑滤波包含步骤如下:
C1、以待滤波位置为窗口中心大小为wm×wn滤波窗口内,生成包含待滤波位置的,沿行、列以及两对角线方向的,中心位置不同的如图3b所示的滤波子窗口;
C2、计算各个子窗口内极化干涉数据总能量的方差σ2 L(m,n):
σ2 L(m,n)=<|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2>(m,n),L 2
-<(|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2)2>(m,n),L
其中σ2 L(m,n)表示以(m,n)中心的第L个子窗口内极化干涉数据总能量的方差;符号||表示复数的模;符号<>(m,n),L以(m,n)中心的第L个子窗口内极化干涉数据的算术平均。
C3、以极化干涉数据总能量方差最小为原则确定滤波子窗口LF。
C4、计算选定的滤波子窗口内极化干涉数据的算术平均T11(m,n)、T22(m,n)、Ω12(m,n):
步骤D)利用特征分解提取最优投影矢量对。本方法基于包含极化、干涉信息的矩阵Ω12,利用特征分解即可提取最优投影矢量对 本方法不要求矩阵T11、T22可逆,可用于各类极化干涉数据的处。利用特征分解提取最优投影矢量对的具体步骤为:
D1、计算矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)的单位特征矢量:
其中,λ1、λ2、λ2是矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)的三个特征矢量; 分别是矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)对应于特征值λi的单位化特征矢量;
步骤E)计算最优干涉相干系数与干涉相位,检验最优估计的有效性。已有的优化干涉相干系数的方法在矩阵T11、T22不可逆时,无法获得正确估计,却并未对优化结果进行有效性检验。本方法添加了对最优估计的有效性检验,其具体步骤如下:
E1、利用最优投影矢量对计算最优干涉相干系数γi(m,n)与干涉相位φ′i(m,n):
E2、检验干涉相干系数值的合理性,设置错误估计标记:
步骤F)将平地相位补偿至最优干涉相位中,按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位。为避免将平地相位补偿至极化干涉数据滤波结果的复杂过程,将平地相位补偿至最优干涉相位中;该平地相位补偿方法可保证最优干涉相位结果的连续性。具体包含如下步骤:
F1、将平地相位φflat,XX(n)补偿至最优干涉相位φ′i(m,n)中,获得干涉相位φi(m,n):
F2、按照γ1(m,n)≥γ2(m,n)≥γ3(m,n)的顺序,排列生成三组最优干涉相干系数以及干涉相位结果{γ1(m,n),φ1(m,n)}、{γ2(m,n),φ2(m,n)}、{γ3(m,n),φ3(m,n)}。
下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是,所描述的实施例仅旨在便于对本发明的理解,而对其不起任何限定作用。
首先介绍一下极化干涉数据及其表示方法:
如图1所示,全极化合成孔径雷达数据是由采用某两种正交极化状态X、Y发送、接收雷达信号的天线系统,采集的到的XX、XY、YX、YY四种极化状态下的雷达信号分别经过成像处理获得的这四种极化状态下的复数据矩阵sXX、sXY、sYX、sYY。而单基线极化干涉合成孔径雷达数据是由不同时间采集的或者同时在不同位置采集的两组全极化合成径雷达数据构成,且每组的各个极化状态下的雷达信号要分别经过成像处理后,两组全极化合成孔径雷达数据还需要经过极化定标处理、图像配准和重采样处理才能获得两组全极化复数据矩阵s1,XX、s1,XY、s1,YX、s1,YY与s2,XX、s2,XY、s2,YX、s2,YY。其中,标号1表示参考图像数据,本发明中叫做主天线全极化数据,标号2表示非参考图像数据,本发明中叫做副天线全极化数据。
本发明实例中使用的数据是美国喷气推进实验室的SIR-C系统于1994年10月7日、8日采集自中国天山山麓的L波段全极化数据构成的单基线极化干涉数据,这两组全极化数据的编号为122.20和138.20。图像尺寸为6000×900(行:方位位置,列:斜距位置,单位:像素)。SIR-C的全极化系统采用H、V极化状态。
如图2所示,基于极化干涉合成孔径雷达数据优化干涉相干系数的方法包含如下步骤:
A)利用快速傅里叶变换计算各极化状态的平地相位;
B)去除各极化状态的平地相位对极化干涉数据的影响;
C)对极化干涉数据平滑滤波;
D)利用特征分解提取最优投影矢量对;
E)计算最优干涉相干系数与干涉相位,检验最优估计的有效性;
F)将平地相位补偿至最优干涉相位中,按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位。
下面介绍各个步骤具体的实施过程:
A1、利用行标在3601-4600范围内的数据,计算HH、HV、VH、VV极化下干涉相位:
φHH(m,n)=Arg(s1HH(m,n)·s* 2HH(m,n))
φHV(m,n)=Arg(s1HV(m,n)·s* 2HV(m,n))
。 (1)
φVH(m,n)=Arg(s1VH(m,n)·s* 2VH(m,n))
φVV(m,n)=Arg(s1VV(m,n)·s* 2VV(m,n))
将各位置的干涉相位φHH(m,n)、φHV(m,n)、φVH(m,n)、φVV(m,n)分别按照如图4d所示的方式量化为灰度值,生成如图4a、图4b、图4c、图4d所示未去去平地相位前HH、HV、VH、VV极化干涉相位量化。观察图4a-图4d发现,干涉相位梯度较大、干涉条纹密集。
A2、将行标在3601-4600范围内的数据补零为长度1024的复数序列,利用快速傅立叶变换计算频谱:
K={0,2...,1023}
其中Xm,XX(K)表示行标号为m的复干涉相位序列快速傅立叶变换结果的第K个位置的值。
A3、计算指定范围内m1≤m≤m2各极化状态复相位序列频谱能量的平均值:
, (3)
A4、确定各极化状态复相位序列频谱能量的平均值的最大值位置
, (4)
通常情况下,不同方位位置的干涉条纹距离向频谱除了受到平地相位的影响还受到地形的影响。为了能准确估计平地相位,本发明方法采用一定范围内的干涉相位距离向平均频谱估计平地相位。由于各个范围位置的地形变化不同,因而在平均过程中将地形变化引起的频谱当作噪声抑制了,而强化了该方位范围内共同的平地相位引起的频谱,因而平均频谱获得的平地相位估计更准确。
A4、由于平均频谱主要受到平地相位的影响,因此可直接利用能量峰值位置计算出平地相位:
n=1,2,...,900
其中φflat,XX(n)表示XX极化状态下第n个距离向像素位置的平地相位。其中平地相位表示当观测区域为平坦地面时干涉相位,平地相位随着距离向位置而变化。图5a、图5b、图5c、图5d是HH、HV、VH、VV极化的平地相位结果的量化图。
步骤B:将副天线全极化图像中各个相位位置的复数据乘以对应平地相位:
。 (6)
为了对比说明去除平地相位前后各个极化状态下干涉相位的情况,利用去除平地相位后的副天线全极化数据与主天线全极化干涉数据计算
HH、HV、VH、VV极化状态下的干涉相位:
φ′HH(m,n)=Arg(s1HH(m,n)·s′* 2HH(m,n))
φ′HV(m,n)=Arg(s1HV(m,n)·s′* 2HV(m,n))
。 (7)
φ′VH(m,n)=Arg(s1VH(m,n)·s′* 2VH(m,n))
φ′VV(m,n)=Arg(s1VV(m,n)·s′* 2VV(m,n))
图6a、图6b、图6c、图6d是去平地相位影响后的HH、HV、VH、VV极化的干涉相位量化图。与图4a-图4d相比较,平地相位去处后干涉相位梯度明显减小,而且相邻像素之间干涉相位的差异也明显减小,从而保证了后续滤波窗口内极化干涉数据的一致性。
步骤C的具体实施过程为:
C1、以待滤波位置为窗口中心大小为5×5滤波窗口内,生成包含待滤波位置的,沿行、列以及两对角线方向的,中心位置不同的如图3c所示的滤波子窗口;
C2、计算各个子窗口内极化干涉数据总能量的方差σ2 L(m,n):
σ2 L(m,n)=<|s1XX|2+|s1XY|2+|s1XY|2+|x1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2>(m,n),L 2
(8)
-<(|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2)2>(m,n),L
其中σ2 L(m,n)表示以(m,n)中心的第L个子窗口内极化干涉数据总能量的方差;符号||表示复数的模;符号<>(m,n),L以(m,n)中心的第L个子窗口内极化干涉数据的算术平均。
C3、以极化干涉数据总能量方差最小为原则确定滤波子窗口LF。
C4、计算选定的滤波子窗口内极化干涉数据的算术平均T11(m,n)、T22(m,n)、Ω12(m,n):
该处理的目的是在适合方向的子窗口内进行平滑滤波,尽量保证滤波区域内的极化干涉数据是同质均匀的。该自适应滤波处理在抑制噪声影响的同时,尽可能的保留各个极化的幅度和相对相位关系,从而保证后续提取到的主要散射中心信息的准确性。
步骤D的具体实施过程为:
D1、计算矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)的单位特征矢量:
(10)
其中,λ1、λ2、λ3是矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)的三个特征矢量; 分别是矩阵Ω12(m,n)Ω12 H(m,n)和Ω12 H(m,n)Ω12(m,n)对应于特征值λi的单位化特征矢量;
步骤E的具体实施过程为:
E1、利用最优投影矢量对计算最优干涉相干系数γi(m,n)与干涉相位φ′i(m,n):
E2、检验干涉相干系数值的合理性,设置错误估计标记:
步骤F的具体实施过程为:
F1、将平地相位φflat,XX(n)补偿至最优干涉相位φ′i(m,n)中,获得干涉相位φi(m,n):
F2、按照γ1(m,n)≥γ2(m,n)≥γ3(m,n)的顺序,排列生成三组最优干涉相干系数以及干涉相位结果{γ1(m,n),φ1(m,n)}、{γ2(m,n),φ2(m,n)}、{γ3(m,n),φ3(m,n)}。
F3、重复执行步骤C-步骤F2,求得各个位置的最优干涉相干系数以及干涉相位,并构成最终优化结果:
(13)
图7a、图7b、图7c分别是第一、第二、第三最优干涉相干系数Γ1、Γ2、Γ3的量化图。图8a、图8b、图8c分别是第一、第二、第三最优干涉相位Φ1、Φ2、Φ3的量化图。
本发明的一个实例,是用美国喷气推进实验室的SIR-C系统重复飞行采集的中国天山山麓的L波段单基线极化干涉数据的优化干涉相干系数的结果。
图4a-图4d是步骤A1中利用未去除平地相位影响的行标号在3601-4000的极化干涉数据计算的HH、HV、VH、VV极化下干涉相位量化图。图4e指明了生成图4a-图4d的量化方式。从图4a-图4d中可发现,干涉相位在列方向形成明暗相间的条纹,且整体的干涉条纹十分密集,这说明干涉相位在距离向上变化很快、相位梯度大,有必要去除平地相位的影响。如图4a-图4b所示,干涉相位形成的条纹十分规律,适用于计算平地相位。
本实例采用方位向标号从3601-4600的干涉相位的距离向频谱估计平地相位。图5a-图5d为计算出的HH、HV、VH、VV极化的平地干涉相位主值的量化图。图5e指明了生成图5a-图5d的量化方式。从图5a-图5b可发现,平地相位在列方向形成平行的明暗相间的条纹。
本实例数据经去除平地相位处理后,行标号从3601-4600数据计算HH、HV、VH、VV极化下的干涉相位。图6a-图6d为计算出的HH、HV、VH、VV极化的平地干涉相位主值的量化图。图6e指明了生成图6a-图6d的量化方式。与图4a-图4d比较发现,平地相位去除后干涉相位梯度明显降低,相邻像素之间的干涉相位差异明显减小。
本实例数据经本发明方法处理,最终获得了三组最优干涉相干系数以及最优干涉相位结果。图7a、图7b、图7c分别是经过本方明的优化相干系数方法处理得到的第一、第二、第三最优干涉相干系数结果量化图;图7d指示了图7a-图7c量化方式。图8a、图8b、图8c分别是利用本方明优化干涉相干系数的方法计算得到的第一、第二、第三最优干涉相位的量化图。图8d指示了图8a-图8c的量化方式。从图8a-图8c发现,三组最优的干涉相位具有较好的连续性。
为对比应用本发明方法后干涉相干系数的改善程度,计算5×5矩形窗口内HH、HV、VH、VV极化的干涉相干系数γHH(m,n)、γHV(m,n)、γVH(m,n)、γVV(m,n):
(14)
图9a、图9b、图9c、图9d分别是HH、HV、VV极化状态下的干涉相干系数均匀量化图。图9e是图9a-图9d中干涉相干系数的量化方式。从图9a-图9d发现,HH、VV极化下的干涉相干系数图比HV、VH极化下干涉相干系数图亮,这说明HH、VV极化下干涉相干系数整体优于HV、VH极化下干涉相干系数,并且还证实了干涉相干系数对极化状态的依赖性以及利用极化状态优化干涉相干系数的可行性。将图9a-图9d与图7a-图7c对比后发现,第一最优干涉相干系数图明显亮于HH、VV极化的干涉相干系数图,这证实了本发明方法确实可以明显优化干涉相干系数。
为了更加清晰的描述本发明提出相干性优化方法提高干涉相干性的能力,图10展示了HH、HV、VV的干涉相干系数与三组最优干涉相干系数的统计直方图。具体的统计方法为将[0,1]区间划分为均匀50小区间,统计图像中干涉相干系数落入各个小区间的频数,从而计算各区间内干涉相干系数的出现频率,从而绘制出统计直方图。从图10发现,第一最优干涉相干系数总体分布更接近于1,其总体干涉相干性明显优于HH、VV的结果;验证了本发明方法的有效性。
通过对本实例结果的分析可证实,本发明提出的基于单基线极化干涉数据优化干涉相干系数的方法确实可以明显提高干涉相干系数、改善干涉相位质量。因而,本发明方法在极化干涉的分类、参数反演的数据处理中有巨大的应用潜力。
Claims (7)
1、基于极化干涉合成孔径雷达数据优化干涉相干系数的方法,其特征在于:包括步骤如下:
A)利用快速傅里叶变换计算各极化状态的平地相位;
B)去除各极化状态的平地相位对极化干涉数据的影响;
C)对极化干涉数据平滑滤波;
D)利用特征分解提取最优投影矢量对;
E)计算最优干涉相干系数与干涉相位,检验最优估计的有效性;
F)将平地相位补偿至最优干涉相位中,按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位。
2、如权利要求1所述的方法,其特征在于,所述A)步先利用快速傅里叶变换计算各极化状态的平地相位,其各极化状态的平地相位计算包含如下步骤:
A1、计算指定范围内m1≤m≤m2,各极化状态下的干涉相位φXX、φXY、φYX、φYY:
φXX(m,n)=Arg(s1XX(m,n)·s* 2XX(m,n))
φXY(m,n)=Arg(s1XY(m,n)·s* 2XY(m,n))
φYX(m,n)=Arg(s1YX(m,n)·s* 2YX(m,n))
φYY(m,n)=Arg(s1YY(m,n)·s* 2YY(m,n))
其中s1表示主天线接收到的复图像数据;s2表示副天线接收到的复图像数据;(m,n)指示像素在图像中的行、列位置;m1、m2表示用于计算平地相位的数据所在行的起始位置、终止位置;XX、XY、YX、YY表示复图像数据的极化方式,符号·表示复数乘运算,*表示共轭运算,Arg表示复数的辐角主值。
A2、利用快速傅立叶变换计算指定范围内m1≤m≤m2,各极化状态的复干涉相位序列的频谱Xm,XX、Xm,XY、Xm,YX、Xm,YY:
其中φXX(m,..)表示由φXX(m,1)、φXX(m,2)、...φXX(m,N)构成的序列;
A3、计算指定范围内m1≤m≤m2各极化状态复相位序列频谱能量的平均值|XXX|、|XXY|、|XYX|、|XYY|:
A4、确定各极化状态复相位序列频谱能量的平均值的最大值位置KPeak,XX、KPeak,XY、KPeak,YX、KPeak,YY;
A5、计算各极化状态平地相位序列φflat,XX、φflat,XY、φflat,YX、φflat,YX:
其中N2表示快速傅里叶变换长度。
3、如权利要求1所述的方法,其特征在于:所述B)步去除各极化状态的平地相位对极化干涉数据的影响,是去除平地相位的影响可减小相邻位置干涉相位梯度,增强相邻位置的极化干涉数据的一致性,具体以下式计算:
。
4、如权利要求1所述的方法,其特征在于:所述C)步极化干涉数据平滑滤波,包含以下步骤:
C1、以待滤波位置为窗口中心,大小为Wm×Wn的滤波窗口内,生成包含待滤波位置的,沿行、列以及两对角线方向的,中心位置不同的滤波子窗口;
C2、计算各个子窗口内极化干涉数据总能量的方差σ2 L(m,n):
σ2 L(m,n)=<|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2>(m,n),L 2-<(|s1XX|2+|s1XY|2+|s1XY|2+|s1YX|2+|s′2XX|2+|s′2XY|2+|s′2XY|2+|s′2YX|2)2>(m,n),L
其中σ2 L(m,n)表示以(m,n)中心的第L个子窗口内极化干涉数据总能量的方差;符号||表示复数的模;符号<>(m,n),L以(m,n)中心的第L个子窗口内极化干涉数据的算术平均;
C3、以极化干涉数据总能量方差最小为原则确定滤波子窗口LF;
C4、计算选定的滤波子窗口内极化干涉数据的算术平均T11(m,n)、T22(m,n)、Ω12(m,n):
7、如权利要求1所述的方法,其特征在于:所述F)步将平地相位补偿至最优干涉相位中,是按干涉相干系数大小顺序排列最优干涉相干系数以及最优干涉相位,具体包含如下步骤:
F1、将平地相位φflat,XX(n)补偿至最优干涉相位φ′i(m,n)中,获得干涉相位φi(m,n):
F2、按照γ1(m,n)≥γ2(m,n)≥γ3(m,n)的顺序,排列生成三组最优干涉相干系数以及干涉相位结果{γ1(m,n),φ1(m,n)}、{γ2(m,n),φ2(m,n)}、{γ3(m,n),φ3(m,n)};其中第一最优结果的干涉相干系数γ1(m,n)最高,明显优于优化前的XX、XY、YX、YY极化的干涉相干系数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN200810223609A CN101685155B (zh) | 2008-09-27 | 2008-09-27 | 基于极化干涉合成孔径雷达数据优化干涉相干系数的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN200810223609A CN101685155B (zh) | 2008-09-27 | 2008-09-27 | 基于极化干涉合成孔径雷达数据优化干涉相干系数的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN101685155A true CN101685155A (zh) | 2010-03-31 |
CN101685155B CN101685155B (zh) | 2012-08-29 |
Family
ID=42048406
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN200810223609A Expired - Fee Related CN101685155B (zh) | 2008-09-27 | 2008-09-27 | 基于极化干涉合成孔径雷达数据优化干涉相干系数的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN101685155B (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102253377A (zh) * | 2011-04-22 | 2011-11-23 | 哈尔滨工业大学 | 基于特征值分析的极化干涉合成孔径雷达目标检测方法 |
CN103323818A (zh) * | 2013-02-25 | 2013-09-25 | 中国科学院电子学研究所 | 多通道合成孔径雷达系统非均匀采样奇异点的方法和装置 |
CN103940834A (zh) * | 2014-05-09 | 2014-07-23 | 中国科学院电子学研究所 | 采用合成孔径雷达技术测量土壤湿度的方法 |
CN104035082A (zh) * | 2014-06-04 | 2014-09-10 | 中国科学院电子学研究所 | 一种合成孔径雷达图像检测方法及装置 |
CN103996175B (zh) * | 2014-05-13 | 2017-02-15 | 西安电子科技大学 | 森林或城市区域高分辨干涉相位滤波方法 |
CN103744079B (zh) * | 2013-12-12 | 2017-02-15 | 中国科学院深圳先进技术研究院 | 一种甘蔗植期的确定方法及系统 |
CN106772311A (zh) * | 2016-11-23 | 2017-05-31 | 中国科学院电子学研究所 | 混合极化合成孔径雷达、发射误差修正系统及修正方法 |
CN107179537A (zh) * | 2017-06-09 | 2017-09-19 | 中国人民解放军61540部队 | 一种基于子空间正交原理的分布式阵列sar相位中心定标方法 |
CN108051810A (zh) * | 2017-12-01 | 2018-05-18 | 南京市测绘勘察研究院股份有限公司 | 一种InSAR分布式散射体相位优化方法 |
CN108205134A (zh) * | 2016-12-16 | 2018-06-26 | 核工业北京地质研究院 | 一种极化合成孔径雷达影像的次地表信息增强方法 |
CN109711446A (zh) * | 2018-12-18 | 2019-05-03 | 中国科学院深圳先进技术研究院 | 一种基于多光谱影像和sar影像的地物分类方法及装置 |
CN109765554A (zh) * | 2018-11-14 | 2019-05-17 | 北京遥感设备研究所 | 一种雷达前视成像系统及方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN100410683C (zh) * | 2005-11-10 | 2008-08-13 | 复旦大学 | 基于全极化合成孔径雷达数据的地表分类方法 |
US8212717B2 (en) * | 2006-10-26 | 2012-07-03 | Raytheon Company | Radar imaging system and method using second moment spatial variance |
-
2008
- 2008-09-27 CN CN200810223609A patent/CN101685155B/zh not_active Expired - Fee Related
Cited By (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102253377A (zh) * | 2011-04-22 | 2011-11-23 | 哈尔滨工业大学 | 基于特征值分析的极化干涉合成孔径雷达目标检测方法 |
CN103323818A (zh) * | 2013-02-25 | 2013-09-25 | 中国科学院电子学研究所 | 多通道合成孔径雷达系统非均匀采样奇异点的方法和装置 |
CN103323818B (zh) * | 2013-02-25 | 2015-06-10 | 中国科学院电子学研究所 | 多通道合成孔径雷达系统非均匀采样奇异点的方法和装置 |
CN103744079B (zh) * | 2013-12-12 | 2017-02-15 | 中国科学院深圳先进技术研究院 | 一种甘蔗植期的确定方法及系统 |
CN103940834A (zh) * | 2014-05-09 | 2014-07-23 | 中国科学院电子学研究所 | 采用合成孔径雷达技术测量土壤湿度的方法 |
CN103940834B (zh) * | 2014-05-09 | 2016-04-27 | 中国科学院电子学研究所 | 采用合成孔径雷达技术测量土壤湿度的方法 |
CN103996175B (zh) * | 2014-05-13 | 2017-02-15 | 西安电子科技大学 | 森林或城市区域高分辨干涉相位滤波方法 |
CN104035082A (zh) * | 2014-06-04 | 2014-09-10 | 中国科学院电子学研究所 | 一种合成孔径雷达图像检测方法及装置 |
CN106772311A (zh) * | 2016-11-23 | 2017-05-31 | 中国科学院电子学研究所 | 混合极化合成孔径雷达、发射误差修正系统及修正方法 |
CN108205134A (zh) * | 2016-12-16 | 2018-06-26 | 核工业北京地质研究院 | 一种极化合成孔径雷达影像的次地表信息增强方法 |
CN108205134B (zh) * | 2016-12-16 | 2021-09-17 | 核工业北京地质研究院 | 一种极化合成孔径雷达影像的次地表信息增强方法 |
CN107179537A (zh) * | 2017-06-09 | 2017-09-19 | 中国人民解放军61540部队 | 一种基于子空间正交原理的分布式阵列sar相位中心定标方法 |
CN108051810A (zh) * | 2017-12-01 | 2018-05-18 | 南京市测绘勘察研究院股份有限公司 | 一种InSAR分布式散射体相位优化方法 |
CN108051810B (zh) * | 2017-12-01 | 2020-06-09 | 南京市测绘勘察研究院股份有限公司 | 一种InSAR分布式散射体相位优化方法 |
CN109765554A (zh) * | 2018-11-14 | 2019-05-17 | 北京遥感设备研究所 | 一种雷达前视成像系统及方法 |
CN109711446A (zh) * | 2018-12-18 | 2019-05-03 | 中国科学院深圳先进技术研究院 | 一种基于多光谱影像和sar影像的地物分类方法及装置 |
CN109711446B (zh) * | 2018-12-18 | 2021-01-19 | 中国科学院深圳先进技术研究院 | 一种基于多光谱影像和sar影像的地物分类方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
CN101685155B (zh) | 2012-08-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN101685155B (zh) | 基于极化干涉合成孔径雷达数据优化干涉相干系数的方法 | |
Stompor et al. | Forecasting performance of CMB experiments in the presence of complex foreground contaminations | |
Cui et al. | On complete model-based decomposition of polarimetric SAR coherency matrix data | |
CN108051810B (zh) | 一种InSAR分布式散射体相位优化方法 | |
CN102270341B (zh) | 一种自适应的高精度干涉sar相位估计方法 | |
Millea et al. | Sampling-based inference of the primordial CMB and gravitational lensing | |
CN103439708B (zh) | 基于广义散射矢量的极化InSAR干涉图估计方法 | |
CN103654789A (zh) | 磁共振快速参数成像方法和系统 | |
CN105388476B (zh) | 一种基于联合稀疏模型的层析sar成像方法 | |
CN103413292B (zh) | 基于约束最小二乘的高光谱图像非线性丰度估计方法 | |
CN102662171A (zh) | 一种sar层析三维成像方法 | |
Moradikia et al. | Joint SAR imaging and multi-feature decomposition from 2-D under-sampled data via low-rankness plus sparsity priors | |
CN108765313B (zh) | 基于类内低秩结构表示的高光谱图像去噪方法 | |
Girard et al. | Sparse representations and convex optimization as tools for LOFAR radio interferometric imaging | |
CN107423543B (zh) | 一种超复数磁共振波谱的快速重建方法 | |
CN103824287A (zh) | 一种基于稳健平面拟合的相位相关亚像素匹配方法 | |
CN110940638B (zh) | 一种高光谱图像亚像元级水体边界探测方法及探测系统 | |
CN106646470A (zh) | 基于广义正交匹配追踪的层析sar三维点云生成方法 | |
Shen et al. | A novel polarimetric PSI method using trace moment-based statistical properties and total power interferogram construction | |
CN111025294A (zh) | 基于均方容积卡尔曼滤波器的InSAR相位解缠绕方法 | |
CN105242236B (zh) | 宽带信号超分辨测向中的阵元位置误差校正方法 | |
WO2021068494A1 (zh) | 基于平面互质阵列虚拟域张量空间谱搜索的高分辨精确二维波达方向估计方法 | |
EP4249944A1 (en) | Direction of arrival (doa) estimation using circular convolutional network | |
Yashtini | Euler’s elastica-based algorithm for parallel MRI reconstruction using sensitivity encoding | |
Shen et al. | JPPL: A joint-polarization phase linking algorithm for phase optimization of TSPolInSAR data |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20120829 Termination date: 20160927 |