多基线/多频段干涉相位解缠频域快速算法

摘要

本发明提出一种多基线/多频段相位解缠的频域算法。技术方案的思路是：首先，利用干涉相位图的干涉相位观测值计算在水平方向及垂直方向上的相位梯度值，并根据时域无旋条件计算相位梯度图边界上的相位梯度值。然后，对相位梯度值进行傅立叶变换得到相位梯度值的频域函数，在相位梯度估计值的频域函数与各基线/各频段相位梯度值的频域函数加权和之差的平方和最小约束条件下计算出满足无旋条件的傅立叶系数的近似值。其后，对傅立叶系数近似值进行傅立叶反变换，得到满足相位无旋条件的相位梯度估计值。最后，对得到的满足无旋条件的相位梯度估计值进行积分，从而得到相位解缠值。本发明在不影响精度的情况下，减少了相位解缠算法的计算量。

说明书

技术领域

本发明属于遥感和信号处理的交叉技术领域，特别涉及一种利用干涉合成孔径雷达进行多基线/多频段条件下的干涉相位解缠的频域快速方法。 

背景技术

一切将相位由主值或相位差值恢复为真实值的过程统称为相位解缠，除了干涉合成孔径雷达应用外，相位解缠在合成孔径声纳、自适应光学、核磁共振、地震处理等方面都有重要应用。利用多部干涉合成孔径雷达形成多基线或多频段干涉可提高相位解缠的性能。 

本发明以干涉合成孔径雷达应用为例。传统的单基线或单频段干涉合成孔径雷达系统受干涉相位模糊和高程叠掩影响，在复杂地形区域相位解缠难度较大，极大地限制了单基线或单频段干涉合成孔径雷达系统的高精度全球测绘能力。多基线或多频段干涉合成孔径雷达系统的提出与实现则有效地提高了干涉合成孔径雷达对复杂地形的测量精度和测量覆盖能力。多基线或多频段干涉合成孔径雷达系统的最大优点就是可以充分利用其长短基线或高低频段获取疏密不同的干涉相位条纹来提高相位解缠的性能，短基线或低频段可以保证相位解缠的可靠性，长基线或高频段可以提高测量精度。因此多基线/多频段干涉合成孔径雷达系统更具吸引力，是未来发展的趋势。 

目前多基线/多频段相位解缠方法主要有：中国余数定律法、投影法以及线性组合法、迭代法、时域最小二乘法、Kalman滤波法、最大似然法、最大后验法、空-像域联合子空间正交投影法和网络流法等，其中多基线/多频段时域最小二乘法的基本思想是使相位梯度估计值与多个相位梯度值加权和之差的平方和最小，等效于求解具有牛曼边界的泊松方程，这种方法本质是对误差进行平均， 特点是十分稳健，但效率不高。 

发明内容

本发明提出一种计算量小的多基线/多频段相位解缠的频域算法，该算法的基本思想是在频域范围内使相位梯度估计值与多个相位梯度值加权和之差的平方和最小。 

本发明技术方案的思路是：首先，针对各基线/各频段的干涉相位图，利用干涉相位图的干涉相位观测值计算在水平方向及垂直方向上的相位梯度值，并根据时域无旋条件计算相位梯度图边界上的相位梯度值。然后，对两个方向上的相位梯度值进行傅立叶变换得到相位梯度值的频域函数，在相位梯度估计值的频域函数与各基线/各频段相位梯度值的频域函数加权和之差的平方和最小约束条件下计算出满足无旋条件的傅立叶系数的近似值。其后，对计算得到的满足无旋条件的傅立叶系数近似值，进行傅立叶反变换，得到满足相位无旋条件的相位梯度估计值。最后，对得到的满足无旋条件的相位梯度估计值进行沿着任意路径的积分，从而得到相位解缠值。 

本发明技术方案是： 

假定通过多基线/多频段观测得到K幅干涉相位图，其中第k幅干涉相位图对应第k条基线/第k个频段，第k条基线/第k个频段对应的缠绕相位函数为： k=1,2,糑,;m=0,1,2,糓,-1;n=0,1,2,糔,-1，M和N分别表示干涉相位图的方位向和距离向点数。记第k幅干涉相位图的垂直有效基线为b_k（或记第k幅干涉相位图的波长为l_k），记第k幅干涉相位图与第一幅干涉相位图的垂直有效基线之比为a_k=b₁/b_k（或记第k幅干涉相位图与第一幅干涉相位图的波长之比为a_k=l_k/l₁）。利用上述的观测信息，完成以下步骤： 

第一步：干涉相位图的相位梯度值计算。 

本步骤对多基线/多频段干涉相位进行相位梯度值计算，同时根据相位的时域无旋条件计算相位梯度图边界上的相位梯度值。对每一幅干涉相位图进行下述计算： 

第（1）步，计算非边界上的相位梯度值。 

对第k幅干涉相位图，用下式计算在水平方向即x方向的相位梯度 在垂直方向即y方向的相位梯度 

上式中W{×}表示取相位主值运算。 

第（2）步：计算边界上的相位梯度值。 

第二步，傅立叶系数近似值计算。 

计算第k幅干涉相位图x方向相位梯度值的傅立叶变换系数 和y方向相位梯度值的傅立叶变换系数 计算公式为： 

其中，p=0,1,K,M-1;q=0,1，K,N-1。 

计算x方向相位梯度的傅立叶变换系数近似值 和y方向相位梯度的傅立叶变换系数近似值 计算公式为： 

其中，$>C_1=\exp \left(j2p\left(\frac{q}{N}\right)\right)-1,$>$>C_2=\exp \left(j2p\left(\frac{p}{M}\right)\right)-1,$>分别是C₁和C₂的共轭，p=0,1,K,M-1;q=0,1,K,N-1。 

第三步：相位梯度估计值计算。 

用下式计算x方向相位梯度估计值 和在y方向的相位梯度估计值 

第四步：相位解缠值计算。 

用下式计算第一幅干涉相位图的相位解缠值： 

上式中 m ii＝0,1，2,K M-2;n ii＝0,1,2,K N-2， 表示第一幅干涉相位图在m=0,n=0时缠绕相位函数的值。 

利用上面计算的相位解缠值，可进行干涉合成孔径雷达系统的后续高程反演和地形形变估计等应用。 

采用本发明可取得以下技术效果： 

本发明的第二步计算两个方向相位梯度的傅立叶变换系数近似值，计算公式是在频域范围内使相位梯度估计值与多个相位梯度值加权和之差的平方和最小这一约束条件推导得到的，无需镜像对称操作，相对时域最小二乘法而言大大减少了计算量，同时保持了相位解缠值的计算精度。 

附图说明

图1为本发明提供的多基线/多频段相位解缠频域快速算法流程示意图； 

图2为基线长度为100米情况下的去平地后的干涉相位图； 

图3为基线长度为200米情况下的去平地后的干涉相位图； 

图4为基线长度为300米情况下的去平地后的干涉相位图； 

图5为利用多基线时域最小二乘法的解缠结果； 

图6为利用本发明得到的解缠结果； 

图7为图5和图6利用的两种解缠算法的差值图； 

图8为计算图5和图6的处理速度对比图。 

具体实施方式

图1为本发明提供的多基线相位解缠频域快速算法流程示意图。整个流程分为四步。第一步，干涉相位图的相位梯度值计算；这一步包括计算非边界上的相位梯度值和计算边界上的相位梯度值这两个步骤，可以得到边界满足无旋条件的相位梯度场。第二步，傅立叶系数近似值计算；基于第一步中建立的相位梯度场，利用最小二乘方法进行傅立叶系数计算，得到满足无旋条件的傅立叶系数近似值。第三步，相位梯度估计值计算；这一步利用第二步中的傅立叶系数估计值，进行傅立叶反变换，获得无旋条件下的相位梯度估计值。第四步，相位解缠值计算；对于第三步所得的满足无旋条件的相位梯度估计值进行沿着任意路径的积分，从而得到相位解缠值。 

图2～图4是利用实验室的天基雷达仿真系统仿真得到的干涉基线分别为100米、200米和300米条件下的干涉相位图，仿真过程中用到的地形是Etna火山口，图像大小为512×768像素。 

图5～图8是利用图2～图4提供的3幅干涉相位图进行仿真实验的结果，即K=3。 

图5为利用多基线时域最小二乘法的解缠结果；图6为利用本发明得到的解缠结果；图7为图5和图6利用的两种解缠算法的差值图。由图7可看出，两种方法的差异在10的负3次方量级，说明两种算法的精度很接近。 

图8为计算图5和图6的处理速度对比图，即利用本发明的解缠算法和传统的多基线时域最小二乘法的解缠方法处理速度对比表。 

从理论上对比：时域最小二乘方法等效于求解具有牛曼边界的泊松方程，为快速求解也采用傅立叶变换到频域，但为了满足牛曼边界条件，需要作镜像对称操作，其傅立叶变换的总运算量为： 

$>2·(\frac{4\mathrm{MN}}{2}·{\log }_{2}2N+4\mathrm{MN}·{\log }_{2}2N)+2·(\frac{4\mathrm{MN}}{2}{·\log }_{2}2M+4\mathrm{MN}·{\log }_{2}2M)=12\mathrm{MN}{\log }_{2}4\mathrm{MN}$>

本发明的方法直接在频域上满足最小二乘约束，避免了牛曼边界条件约束，无需镜像对称操作，其傅立叶变换总运算量为： 

$>2·(\frac{\mathrm{MN}}{2}·{\log }_{2}N+\mathrm{MN}·{\log }_{2}N)+2·(\frac{\mathrm{MN}}{2}·{\log }_{2}M+\mathrm{MN}·{\log }_{2}M)=3\mathrm{MN}{\log }_2\mathrm{MN}$>

利用上述两个公式，比较两种方法的傅立叶变换运算量，可知由于本发明不用做镜像对称操作，FFT数据处理点数小，因此可以更为有效快速地得到解缠相位。 

对于两种算法，实验环境均为Inter Core2 Quad CPU 2.33GHz，内存2GB。对于上述多基线干涉相位，图8中结果表明，本发明提出的算法处理时间小于多基线时域最小二乘解缠方法的1/3，速度更快。