一种采用误差迭代补偿的高精度相位解缠方法

摘要

本发明公开一种采用误差迭代补偿的高精度相位解缠方法，它是基于传统的快速傅立叶变换(FFT)的最小二乘相位解缠算法，首先得到解缠相位的误差主值，通过将每次解缠相位与干涉相位差值的主值作为下次迭代的“缠绕干涉相位”，进行基于FFT的最小二乘解缠，并进行一次相位滤波，补偿解缠误差。然后对其进行基于FFT的最小二乘求解；最后反复迭代处理，将剩余相位误差进行相位滤波后补偿解缠相位，本发明与传统的快速傅立叶变换(FFT)的最小二乘相位解缠算法相比，在噪声较大、信噪比较低的情况下也能有效减小解缠误差，显著提高解缠精度具有高解缠精度的效果，为后续的InSAR高程反演精度提供保证。

说明书

技术领域

本发明属于雷达技术领域，它特别涉及干涉合成孔径雷达(InSAR)相位解缠技术领域。

背景技术

干涉合成孔径雷达(InSAR)在合成孔径雷达(SAR)成像技术的基础上发展而来，是获取数字高程的技术之一，其基本原理是以不同视角对同一地区的两幅SAR复图像的相位差(即干涉相位)与地形的几何关系来获取地形高度信息。InSAR具有全天时、全天候、高精度大面积测绘等特点，是目前提取大面积地表三维图像和地形高程变化信息的一项重要遥感技术，广泛运用于自然灾害监测、地形测绘和自然资源调查等领域。随着InSAR技术发展，地形测绘的精度要求越来越高，如何提高测绘精度是当今InSAR应用领域的迫切需求。由于InSAR处理过程中相位提取得到的干涉相位都是其位于[-π,π)之间的主值，称之为缠绕干涉相位，因此需要进行相位解缠，即由缠绕相位恢复其真实值。相位解缠在InSAR数据处理流程中至关重要，其精度的高低将直接影响后面高程反演结果的精度。实际上地形陡变、雷达阴影、相位噪声、图像失配等因素都会影响相位解缠精度和难度。因此，如何实现精确无误且高效快速的相位解缠仍然是一个难题。详见文献“星载合成孔径雷达干涉测量”，王超等编著，科学出版社。

基于快速傅立叶变换(FFT)的最小二乘相位解缠算法是一种全局算法，在最小二乘意义下寻找使缠绕相位的离散偏导数与解缠相位的偏导数整体偏差最小的解。它能得到最小二乘意义下的最优唯一解，具有稳定性强、不需识别残差点、处理简单和结果连续性好的优点，是目前实际应用中较常用的相位解缠方法。但由于用该方法解缠时没有绕过而是穿过包含残差点的相位不连续区域，因此会造成局部误差的全局传播,从而导致全局误差。详见文献“Least-squares two-dimensional phase unwrapping using FFTs”，Pritt MD&Shipman JS，IEEE Trans Geosci Remote Sense。

发明内容

本发明提出了一种采用误差迭代补偿的高精度相位解缠方法，基于传统的快速傅立叶变换(FFT)的最小二乘相位解缠算法，首先得到解缠相位的误差主值，然后对其进行基于FFT的最小二乘求解；最后反复迭代处理，将剩余相位误差进行相位滤波后补偿解缠相位。本发明与传统的快速傅立叶变换(FFT)的最小二乘相位解缠算法相比，具有高解缠精度的效果。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、干涉合成孔径雷达(InSAR)

干涉合成孔径雷达(InSAR)指利用同一观测场景不同观测角度的两组或者两组以上SAR数据进行干涉成像处理，然后结合雷达系统参数和雷达平台几何位置信息反演地形高度及高程变化信息的遥感技术，详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。

定义2、缠绕相位： 

缠绕相位是指在实际干涉处理中，经过三角运算，得到的干涉相位的主值。详见文献“星载合成孔径雷达干涉测量”，王超等编著，科学出版社。

定义3、解缠相位： 

解缠相位是指干涉处理中从缠绕相位经过相位解缠处理恢复的真实相位。详见文献“星载合成孔径雷达干涉测量”，王超等编著，科学出版社。

定义4、相位解缠： 

一切将相位由主值或相位差恢复为真实值的过程统称为相位解缠。详见文献“星载合成孔径雷达干涉测量”，王超等编著，科学出版社。

定义5、缠绕操作W[·]：

缠绕操作即取在(-π,π]的主值，详见文献“星载合成孔径雷达干涉测量”，王超等编著，科学出版 社。

定义6、基于FFT的最小二乘相位解缠：

一种最小范数类的相位解缠方法，在最小二乘意义下寻找使缠绕相位的离散偏导数与解缠相位的偏导数整体偏差最小的解。

假设缠绕相位为φ_i,j，解缠相位为i＝1,2,...,M,j＝1,2,...,N,M表示方位向点数，N表示距离向点数。则基于FFT的最小二乘相位解缠方法可大致分为以下几个步骤。

步骤1、对缠绕相位φ_i,j采用公式

${\tilde{\phi }}_{i,j}=\left(\begin{array}{cc}{\phi }_{i,j}& (1\le i\le M,1\le j\le N)\\ {\phi }_{2M-1-i,j}& (M<i\le 2M-\mathrm{1,1}\le j\le N)\\ {\phi }_{i,2N-1-j}& (1\le i\le M,N<j\le 2N-1)\\ {\phi }_{2M-1-i,2N-1-j}& (M<i\le 2M-1,N<j\le 2N-1)\end{array}\right)$

作镜像对称操作。

步骤2、采用公式 

${{\Delta }^a}_{i,j}=W[{\tilde{\phi }}_{i+1,j}-{\tilde{\phi }}_{i,j}]$

${{\Delta }^r}_{i,j}=W[{\tilde{\phi }}_{i,j+1}-{\tilde{\phi }}_{i,j}]$

计算距离向和方位向的一阶相位梯度。

步骤3、采用公式 

ρ_i,j＝[Δ^a_i+1,j-Δ^r_i,j]+[Δ^a_i,j+1-Δ^r_i,j]

计算二阶相位梯度和。

步骤4、对ρ_i,j作二维快速傅立叶变换，得到P_i,j。

步骤5、采用公式 

${\Phi }_{i,j}=P_{i,j}/(2\cos \frac{\mathrm{\pi i}}{M}+2\cos \frac{\mathrm{\pi j}}{N}-4)$

计算未缠绕相位的二维快速傅里叶变换。

步骤6、对Φ_i,j作二维逆傅立叶变换得到解缠相位

具体方法流程详见文献“星载合成孔径雷达干涉测量”，王超等编著，科学出版社。

定义7、Goldstein相位滤波： 

Goldstein等人于1998年提出的在频域上滤除干涉相位噪声的滤波方法。该滤波方法首先对干涉图进行分块,然后对每个小块干涉图进行傅立叶变换,得到其频谱，再采用经平滑处理的幅值对各小块干涉图进行处理。详见文献“基于信噪比的InSAR干涉图相位滤波方法研究”，孙倩，中南大学硕士学位论文。

定义8、快速傅里叶变换(FFT)：

快速傅里叶变换(FFT)，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。详见“数字信号处理”，程乾生等编著，北京大学出版社。

定义9、矩阵实验室(Matlab)软件：

一款常用的数学软件，randn函数是其函数库里产生高斯随机矩阵的函数。详见文献“MATLAB实用教程”，郑阿奇等编著，电子工业出版社。

本发明提供的一种采用误差迭代补偿的高精度相位解缠方法，它包括以下步骤：

步骤1、初始化采用误差迭代补偿的高精度相位解缠方法所需参数：

初始化采用误差迭代补偿的高精度相位解缠方法所需参数，包括：InSAR缠绕相位的方位向点数，记为N_a；InSAR缠绕相位的距离向点数，记为N_r，其中N_a、N_r为正整数；InSAR待解缠的缠绕相位，记为ω(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，其中a表示方位向第a个点，r表示距离向第r个点，a,r为正整数；误差迭代补偿次数Niter为自然数，迭代最大次数Nmax为正整数，相位加入的高斯噪声的标准差，记为σ。以上参数初始化后均为已知。

步骤2、将缠绕相位加入高斯噪声：

采用公式φ(a,r)＝σ·randn(N_a,N_r)+ω(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，计算得到加入高斯噪声的缠绕相位，记为φ(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，a表示方位向第a个点，r表示距离向第r个点；N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数；ω(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤1初始化的InSAR待解缠的缠绕相位；σ为步骤1初始化的相位加入的高斯噪声的标准差，randn(N_a,N_r)表示数学软件MATLAB库函数randn产生的均值为0、标准差为1的高斯噪声，并且为N_a×N_r的随机矩阵。

步骤3、对φ(a,r)进行传统的基于FFT的最小二乘相位解缠，得到初始解缠相位，记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中φ(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤2得到的缠绕相位；a表示方位向第a个点，r表示距离向第r个点，N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数。

步骤4、采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，计算初始解缠相位与缠绕相位的差的主值，记为Δε_w0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。其中，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤3得到的初始解缠相位；φ(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤2得到的缠绕相位；a表示方位向第a个点，r表示距离向第r个点；N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数；W[·]表示缠绕操作。

步骤5、对Δε_w0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，进行传统的基于FFT最小二乘法相位解 缠，得到结果记为Δε_u0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。其中Δε_w0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤4得到的初始解缠相位与缠绕相位的差的主值；a表示方位向第a个点，r表示距离向第r个点；N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数。

步骤6、相位误差补偿：

采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，对解缠相位进行误差补偿，结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤3得到的初始解缠相位；Δε_u0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤5得到的结果；a表示方位向第a个点，r表示距离向第r个点，N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数。

步骤7、迭代判定： 

当Niter＝1时，对采用与步骤4同样的方法，得到的结果记为Δε_w1(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_w1(a,r)采用与步骤5同样的方法，得到的结果记为Δε_u1(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_u1(a,r)采用与步骤6同样的方法，得到的结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。

当Niter＝k时，k＝2,3,4…,Nmax-1，对采用与步骤4同样的方法，得到的结果记为Δε_wk(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_wk(a,r)采用与步骤5同样的方法，得到的结果记为Δε_uk(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_uk(a,r)采用与步骤6同样的方法，得到的结果记为 a＝1,2,…,N_a，r＝1,2,…,N_r。

当Niter＝Nmax时，对采用与步骤4同样的方法，得到的结果记为Δε_wNmax(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_wNmax(a,r)采用与步骤5同样的方法，得到的结果记为Δε_uNmax(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_uNmax(a,r)采用与步骤6同样的方法，得到的结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。此时迭代结束。

其中，a表示方位向第a个点，r表示距离向第r个点，N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数；Niter为已经迭代的次数，Nmax为步骤 1初始化得到的最大迭代次数。

步骤8、相位滤波： 

采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，计算误差迭代补偿后的解缠相位与缠绕相位之差的主值，对结果Δε_e(a,r)进行Goldstein相位滤波，滤波结果记为Δε_fe(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。其中a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤7迭代结束后得到的结果；a表示方位向第a个点，r表示距离向第r个点；N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数；W[·]表示缠绕操作。

步骤9、对Δε_fe(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r采用与步骤3到步骤7相同的方法，得到的结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，补偿误差，结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中，Δε_fe(a,r)为步骤8得到的滤波后的结果；为步骤7迭代结束后的结果；a表示方位向第a个点，r表示距离向第r个点；N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数。

步骤10、采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，得到最终的解缠结果，记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中，为步骤9得到的结果，φ(a,r)为步骤1初始化的缠绕相位，a表示方位向第a个点，r表示距离向第r个点，N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数；W[·]表示缠绕操作。

经过以上步骤，得到最终高精度的解缠相位为后续的InSAR高程反演精度提供保证。

本发明的创新点在于提出了一种采用误差迭代补偿的高精度相位解缠方法，通过将每次解缠相位与干涉相位差值主值作为下次迭代的“缠绕干涉相位”进行基于FFT的最小二乘解缠，并进行一次相位滤波，补偿解缠误差。本发明的实质是通过误差的迭代补偿，使得每次迭代的解缠相位误差越来越小，最终得到高精度的解缠相位，为后续的InSAR高程反演精度提供保证。

本发明的特点是在噪声较大、信噪比较低的情况下也能有效减小解缠误差，显著提高解缠精度。

附图说明

图1给出了本发明所提供方法的流程示意图和基于FFT的最小二乘方法流程图(采用误差迭代补偿的相位解缠方法流程示意图)；图2是最小二乘相位解缠算法流程示意图。

具体实施方式

本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在MATLABR2013a软件上验证正确。具体实施步骤如下：

步骤1、初始化采用误差迭代补偿的高精度相位解缠方法所需参数：

初始化采用误差迭代补偿的高精度相位解缠方法所需参数，包括：InSAR缠绕相位的方位向点数，记 为N_a＝500；InSAR缠绕相位的距离向点数，记为N_r＝500；InSAR待解缠的缠绕相位，记为ω(a,r)，a＝1,2,…,500，r＝1,2,…,500，其中a表示方位向第a个点，r表示距离向第r个点；误差迭代补偿次数Niter＝0，迭代最大次数Nmax＝15，相位加入的高斯噪声的标准差，记为σ＝1。以上参数初始化后均为已知。

步骤2、将缠绕相位加入高斯噪声：

采用公式φ(a,r)＝σ·randn(N_a,N_r)+ω(a,r)，a＝1,2,…,500，r＝1,2,…,500，计算得到加入高斯噪声的缠绕相位，记为φ(a,r)，a＝1,2,…,500，r＝1,2,…,500，a表示方位向第a个点，r表示距离向第r个点；N_a＝500，N_r＝500；ω(a,r)，a＝1,2,…,500，r＝1,2,…,500，为步骤1初始化的InSAR待解缠的缠绕相位；σ＝1为步骤1初始化的相位加入的高斯噪声的标准差，randn(N_a,N_r)表示数学软件MATLAB库函数randn产生的均值为0、标准差σ为1的高斯噪声，并且为500×500的随机矩阵。

步骤3、对φ(a,r)进行传统的基于FFT的最小二乘相位解缠，得到初始解缠相位，记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中φ(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤2得到的缠绕相位；a表示方位向第a个点，r表示距离向第r个点，N_a＝500,N_r＝500。

步骤4、采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，计算初始解缠相位与缠绕相位的差的主值，记为Δε_w0(a,r)，a＝1,2,…,N_a，r＝1,2,…,Nr。其中，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤3得到的初始解缠相位；φ(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤2得到的缠绕相位；a表示方位向第a个点，r表示距离向第r个点；N_a为步骤1初始化的InSAR缠绕相位的方位向点数，N_r为步骤1初始化的InSAR缠绕相位的距离向点数；N_a＝500,N_r＝500；W[·]表示缠绕操作。

步骤5、对Δε_w0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，进行传统的基于FFT最小二乘法相位解缠，得到结果记为Δε_u0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。其中Δε_w0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤4得到的初始解缠相位与缠绕相位的差的主值；a表示方位向第a个点，r表示距离向第r个点；N_a＝500,N_r＝500。

步骤6、相位误差补偿：

采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，对解缠相位进行误差补偿，结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中，a＝1,2,…,N_a， r＝1,2,…,N_r，为步骤3得到的初始解缠相位；Δε_u0(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤5得到的结果；a表示方位向第a个点，r表示距离向第r个点；N_a＝500,N_r＝500。

步骤7、迭代判定： 

当Niter＝1时，对采用与步骤4同样的方法，得到的结果记为Δε_w1(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_w1(a,r)采用与步骤5同样的方法，得到的结果记为Δε_u1(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_u1(a,r)采用与步骤6同样的方法，得到的结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。

当Niter＝k时，k＝2,3,4…，Nmax-1，对采用与步骤4同样的方法，得到的结果记为Δε_wk(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_wk(a,r)采用与步骤5同样的方法，得到的结果记为Δε_uk(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_uk(a,r)采用与步骤6同样的方法，得到的结果记为 a＝1,2,…,N_a，r＝1,2,…,N_r。

当Niter＝Nmax时，对采用与步骤4同样的方法，得到的结果记为Δε_wNmax(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_wNmax(a,r)采用与步骤5同样的方法，得到的结果记为Δε_uNmax(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。对Δε_uNmax(a,r)采用与步骤6同样的方法，得到的结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。此时迭代结束。

其中，a表示方位向第a个点，r表示距离向第r个点，N_a＝500,N_r＝500；Niter为已经迭代的次数，Nmax＝15为步骤1初始化得到的最大迭代次数。

步骤8、相位滤波： 

采用公式a＝1,2,…,N_a，r＝1,2,…，N_r，计算误差迭代补偿后的解缠相位与缠绕相位之差的主值，对结果Δε_e(a,r)进行Goldstein相位滤波，滤波结果记为Δε_fe(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r。其中a＝1,2,…,N_a，r＝1,2,…,N_r，为步骤7迭代结束后得到的结果；a表示方位向第a个点，r表示距离向第r个点；N_a＝500,N_r＝500；W[·]表示缠绕操作。

步骤9、对Δε_fe(a,r)，a＝1,2,…,N_a，r＝1,2,…,N_r采用与步骤3到步骤7相同的方法，得到的结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。采用公式 a＝1,2,…,N_a，r＝1,2,…,N_r，补偿误差，结果记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中，Δε_fe(a,r)为步骤8得到的滤波后的结果；为步骤7迭代结束后的结果；a表示方位向第a个点，r表示距离向第r个点；N_a＝500,N_r＝500。

步骤10、采用公式a＝1,2,…,N_a，r＝1,2,…,N_r，得到最终的解缠结果，记为a＝1,2,…,N_a，r＝1,2,…,N_r。其中，为步骤9得到的结果，φ(a,r)为步骤1初始化的缠绕相位，a表示方位向第a个点，r表示距离向第r个点，N_a＝500,N_r＝500；W[·]表示缠绕操作。

经过以上步骤，得到最终的解缠相位a＝1,2,…,N_a，r＝1,2,…,N_r,N_a＝500,N_r＝500。

仿真结果表明，本发明提出的方法的最终解缠精度相比传统的基于FFT的最小二乘相位解缠方法能得到更高的解缠精度，特别是在相位噪声较大的情况下解缠优势更明显。