CN104808203A - 一种迭代最大似然估计多基线InSAR相位解缠方法 - Google Patents

Abstract

本发明提出了一种迭代最大似然估计多基线InSAR相位解缠方法，该方法首先采用最大似然函数估计方法，得到多基线InSAR初始的解缠相位，再求出几组较大似然函数值对应的解缠相位，然后通过多组较大似然函数对应的解缠相位迭代优化初始解缠相位，减小相位解缠过程中噪声带来的误差，与传统的最大似然估计多基线相位解缠方法相比，本发明方法在噪声较大情况下，依然能得到高精度的解缠相位，为后续的高精度InSAR高程反演提供保证。

Description

一种迭代最大似然估计多基线InSAR相位解缠方法

技术领域

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

背景技术

干涉合成孔径雷达(InSAR)是合成孔径雷达(SAR)成像技术的扩展，可以测量得到目标的高度信息，从而实现对目标场景的三维测绘，该技术主要利用两个或两个以上的天线对同一观测场景进行不同角度观测得到的两幅以上SAR复图像，通过两幅SAR复图像信息可以得到干涉相位信息，经过干涉处理后的相位信息通过高程反演后得到场景的地形高程信息。InSAR技术不但可以全天候、全天时成像，并且成像区域大、成像精度高，已经成为当今世界获取地形三维图像和地形高程变化信息的一项重要遥感技术。其中多基线InSAR有着测量精度更高、能更有效地处理复杂地形等优点，多基线InSAR也得到了越来越多的关注，也是当今InSAR研究的重要方向之一。

在多基线InSAR中，每个副天线的复图像和主天线复图像做干涉处理可以得到一组干涉相位，由于每个副天线的位置不同，所以导致了每个副天线和主天线得到的干涉相位存在一定的差异，而基线长度与干涉相位的条纹个数有关，并且长基线对应的副天线得到的干涉相位具有更多的干涉条纹，所以能够反演得到精度更高的地形信息。

相位解缠是InSAR处理的关键步骤，InSAR处理过程中相位提取得到的干涉相位都是缠绕的，进行地形高程重构前必须对缠绕的干涉相位进行解缠，相位解缠精度的高低，直接影响后面高程反演的精度。为了对复杂地形的缠绕相位和低信噪比条件下的缠绕相位进行相位解缠，高精度的多基线InSAR相位解缠方法有了深入的研究。

最大似然估计多基线相位解缠方法利用多个通道的缠绕相位进行相位解缠，根据最大似然函数值求得的解缠相位，解缠精度较高，然而，在有噪声存在时，最大似然函数值对应的解缠相位不一定是准确的，有可能是另外一个较大的似然函数值对应的解缠相位，这时得到的解缠相位就有较大误差。详见文献“G.Fornaro.Maximum likelihood multi-baseline SAR interferometry.2005.”。

发明内容

本发明提出了一种迭代最大似然估计多基线InSAR相位解缠方法，该方法首先采用最大似然函数估计方法，得到多基线InSAR初始的解缠相位，再求出几组较大似然函数值对应的解缠相位，然后通过多组较大似然函数对应的解缠相位迭代优化初始解缠相位，减小相位解缠过程中噪声带来的误差，为后续的高精度InSAR高程反演提供保证。为了方便描述本发明的内容，首先作以下术语定义：

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

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

定义2、相位解缠

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

定义3、似然函数：

设总体X服从分布P(x；θ)(当X是连续型随机变量时为概率密度，当X为离散型随机变量时为概率分布)，θ为待估参数，X₁,X₂,….X_n是来自于总体X的样本，x₁,x₂,…,x_n为样本X₁,X₂,….X_n的一个观察值，则样本的联合分布(当X是连续型随机变量时为概率密度，当X为离散型随机变量时为概率分布)称为似然函数。详见文献“概率论与数理统计”，徐全智等编著，电子科技大学出版社。

定义4、缠绕相位：

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

定义5、解缠相位：

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

定义6、高斯分布函数randn(m,n)

高斯分布函数randn(m,n)是MATLAB中产生均值为0，标准差为1的高斯分布，并且为m×n的随机矩阵。

本发明提供的一种迭代最大似然估计多基线InSAR相位解缠方法，它包括以下几个步骤：

步骤1、初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数：

初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数，包括：多基线InSAR方位向点数，记为N_x；多基线InSAR距离向点数，记为N_y；多基线InSAR待解缠的缠绕相位，记为ω(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点；多基线InSAR非待解缠的缠绕相位组数，记为t；多基线InSAR非待解缠的缠绕相位，记为β(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，其中i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位；多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，记为ε_i，i＝1,2,…,t；多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，记为γ；迭代最大次数，记为MI；相位加入的高斯噪声的标准差，记为λ。

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

采用公式ψ(a,r)＝λ·randn(N_x,N_y)+ω(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，计算得到加入高斯噪声的待解缠的缠绕相位，记为ψ(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，λ为步骤1初始化的相位加入的高斯噪声的标准差，ω(a,r)为步骤1初始化的多基线InSAR待解缠的缠绕相位，randn(N_x,N_y)表示MATLAB产生的均值为0、标准差为1的高斯噪声，并且为N_x×N_y的随机矩阵。

采用公式φ(i；a,r)＝λ·randn(N_x,N_y)+β(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，计算得到加入高斯噪声的非待解缠的缠绕相位，记为φ(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，其中i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，β(i；a,r)为步骤1初始化的多基线InSAR非待解缠的缠绕相位。

步骤3、利用最大似然函数估计多基线InSAR解缠相位：

多基线InSAR待解缠的缠绕相位的缠绕次数，记为k。采用公式

       $\begin{array}{c}{\hat{k}}_0=\arg \max \underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_a，r＝1,2,…,N_r，i＝1,2,…,t)，计算得到最大似然函数估计的缠绕次数，记为其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率，argmax(·)表示求

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。

采用公式a＝1,2,…,N_x，r＝1,2,…,N_y，计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，记为a＝1,2,…,N_x，r＝1,2,…,N_y。

将取代公式中的k：

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

得到：

       $\begin{array}{c}{A}_0=\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t)，令：A₀＝0。

步骤4、计算较大的似然函数值估计的解缠相位：

采用公式

       $\begin{array}{c}{\hat{k}}_m=\arg \max \underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，m＝1,2,…,MI)，计算得到第m次迭代的缠绕次数，记为m＝1,2,…,MI，其中m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数，a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率，argmax(·)表示求

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。

采用公式a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，计算得到方位向第a个、距离向第r个点的第m次迭代的多基线InSAR解缠相位，记为r＝1,2,…,N_y，m＝1,2,…,MI。

采用公式

      

a＝1,2,…,N_x，r＝1,2,…,N_y，计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，记为SUM₀，其中a-1表示方位向第a-1个点，a+1表示方位向第a+1个点，r-1表示距离向第r-1个点，r+1表示距离向第r+1个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，为步骤3得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，|·|表示求绝对值。

采用公式

      

a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，计算得到方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，记为SUM_m，m＝1,2,…,MI。

步骤5、优化最大似然估计的多基线InSAR解缠相位：

比较SUM₀和SUM_m的大小：

如果m＝1,2,…,MI，则令：a＝1,2,…,N_x，r＝1,2,…,N_y，优化方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数，SUM₀为步骤4得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，SUM_m为步骤4得到的方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和。

如果SUM_m≥SUM₀，m＝1,2,…,MI，的值不变。

步骤6、判断算法迭代条件：

如果m＜MI，迭代次数m加1，将代入

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

得到

       $\begin{array}{c}{A}_m=\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，m＝1,2,…,MI)，令：A_m＝0，重复步骤4至步骤6进行迭代计算，其中m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数，a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，Ny为步骤1初始化的多基线InSAR距离向点数，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，为步骤4得到的第m次迭代的缠绕次数，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率。

如果m＝MI，则终止迭代。

经过以上步骤得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位即为最终的解缠相位。

本发明的创新点在于提出了一种迭代最大似然估计多基线InSAR相位解缠方法，本发明利用最大似然估计多基线InSAR解缠相位，并利用几组较大似然函数值估计得到的解缠相位迭代优化最大似然估计的解缠相位，使得在有较大噪声时，多基线InSAR相位解缠精度很高，。

本发明的优点是通过几组较大似然函数值估计的解缠相位迭代优化最大似然估计的多基线InSAR解缠相位，减小相位解缠过程中噪声带来的误差，使得在信噪比较低的情况下，多基线InSAR相位解缠的精度依然很高。

附图说明

图1为本发明所提供方法的流程示意图；

图2为多基线InSAR示意图

其中A₀为多基线InSAR主天线，A₁，A₂，A₃为多基线InSAR副天线，p为场景中心点目标，H为平台高度，θ为主天线雷达入射角，α为基线与水平方向夹角。

图3为本发明所提供方法的多基线InSAR相位解缠参数值

具体实施方式

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

步骤1、初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数：

初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数，包括：多基线InSAR方位向点数N_x＝500；多基线InSAR距离向点数，记为N_y＝500；多基线InSAR待解缠的缠绕相位，记为ω(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点；多基线InSAR非待解缠的缠绕相位组数t＝6；多基线InSAR非待解缠的缠绕相位，记为β(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，其中i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位；多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数 ${\mathrm{\ϵ}}_1=\frac{3}{5},{\mathrm{\ϵ}}_2=\frac{5}{8},{\mathrm{\ϵ}}_3=\frac{29}{41},{\mathrm{\ϵ}}_4=\frac{67}{100},{\mathrm{\ϵ}}_5=\frac{37}{97},{\mathrm{\ϵ}}_6=\frac{59}{101};$ 多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ＝0.9；迭代最大次数MI＝6；相位加入的高斯噪声的标准差λ＝0.8。

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

采用公式ψ(a,r)＝λ·randn(N_x,N_y)+ω(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，N_x＝500，N_y＝500，计算得到加入高斯噪声的待解缠的缠绕相位，记为ψ(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，N_x＝500，N_y＝500，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数N_x＝500，N_y为步骤1初始化的多基线InSAR距离向点数N_y＝500，λ为步骤1初始化的相位加入的高斯噪声的标准差λ＝0.8，ω(a,r)为步骤1初始化的多基线InSAR待解缠的缠绕相位，randn(N_x,N_y)表示MATLAB产生的均值为0、标准差为1的高斯噪声，并且为N_x×N_y的随机矩阵。

采用公式φ(i；a,r)＝λ·randn(N_x,N_y)+β(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，N_x＝500，N_y＝500，t＝6，计算得到加入高斯噪声的非待解缠的缠绕相位，记为φ(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，N_x＝500，N_y＝500，t＝6，其中i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t＝6，β(i；a,r)为步骤1初始化的多基线InSAR非待解缠的缠绕相位。

步骤3、利用最大似然函数估计多基线InSAR解缠相位：

多基线InSAR待解缠的缠绕相位的缠绕次数，记为k。采用公式

       $\begin{array}{c}{\hat{k}}_0=\arg \max \underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_a，r＝1,2,…,N_r，i＝1,2,…,t，N_x＝500，N_y＝500，t＝6)，计算得到最大似然函数估计的缠绕次数，记为其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数N_x＝500，N_y为步骤1初始化的多基线InSAR距离向点数N_y＝500，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t＝6，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数 ${\mathrm{\ϵ}}_1=\frac{3}{5},{\mathrm{\ϵ}}_2=\frac{5}{8},{\mathrm{\ϵ}}_3=\frac{29}{41},{\mathrm{\ϵ}}_4=\frac{67}{100},{\mathrm{\ϵ}}_5=\frac{37}{97},{\mathrm{\ϵ}}_6=\frac{59}{101},$ γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ＝0.9，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率，argmax(·)表示求

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。

采用公式a＝1,2,…,N_x，r＝1,2,…,N_y，N_x＝500，N_y＝500，计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，记为a＝1,2,…,N_x，r＝1,2,…,N_y，N_x＝500，N_y＝500。

将代入公式

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

得到：

       $\begin{array}{c}{A}_0=\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，N_x＝500，N_y＝500，t＝6)，令：A₀＝0。

步骤4、计算迭代所需较大的似然函数值估计的解缠相位：

采用公式

       $\begin{array}{c}{\hat{k}}_0=\arg \max \underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，m＝1,2,…,MI，N_x＝500，N_y＝500，t＝6，

MI＝6)，计算得到第m次迭代的缠绕次数，记为m＝1,2,…,MI，MI＝6，其中m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数MI＝6，a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数N_x＝500，N_y为步骤1初始化的多基线InSAR距离向点数N_y＝500，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t＝6，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数 ${\mathrm{\ϵ}}_1=\frac{3}{5},{\mathrm{\ϵ}}_2=\frac{5}{8},{\mathrm{\ϵ}}_3=\frac{29}{41},{\mathrm{\ϵ}}_4=\frac{67}{100},{\mathrm{\ϵ}}_5=\frac{37}{97},{\mathrm{\ϵ}}_6=\frac{59}{101},$ γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ＝0.9，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率，argmax(·)表示求

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。

采用公式a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，N_x＝500，N_y＝500，MI＝6，计算得到方位向第a个、距离向第r个点的第m次迭代的多基线InSAR解缠相位，记为a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，N_x＝500，N_y＝500，MI＝6。

采用公式

      

a＝1,2,…,N_x，r＝1,2,…,N_y，N_x＝500，N_y＝500，计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，记为SUM₀，其中a-1表示方位向第a-1个点，a+1表示方位向第a+1个点，r-1表示距离向第r-1个点，r+1表示距离向第r+1个点，N_x为步骤1初始化的多基线InSAR方位向点数N_x＝500，N_y为步骤1初始化的多基线InSAR距离向点数N_y＝500，为步骤3得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，|·|表示求绝对值。

采用公式

      

a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，N_x＝500，N_y＝500，MI＝6，计算得到方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，记为SUM_m，m＝1,2,…,MI，MI＝6。

步骤5、优化最大似然估计的多基线InSAR解缠相位：

比较SUM₀和SUM_m的大小：

如果SUM_m＜SUM₀，m＝1,2,…,MI，MI＝6，则令：a＝1,2,…,N_x，r＝1,2,…,N_y，N_x＝500，N_y＝500，优化方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数N_x＝500，N_y为步骤1初始化的多基线InSAR距离向点数N_y＝500，m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数MI＝6，SUM₀为步骤4得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，SUM_m为步骤4得到的方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和。

如果SUM_m≥SUM₀，m＝1,2,…,MI，MI＝6，则不进行任何计算，的值不变。

步骤6、判断算法迭代条件：

如果m＜MI，MI＝6，迭代次数m加1，将代入

       $\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

得到：

       $\begin{array}{c}{A}_m=\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a.r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，m＝1,2,…,MI，N_x＝500，N_y＝500，t＝6，MI＝6)，令：A_m＝0，重复步骤4至步骤6进行迭代计算，其中m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数MI＝6，a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数N_x＝500，N_y为步骤1初始化的多基线InSAR距离向点数N_y＝500，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t＝6，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数 ${\mathrm{\ϵ}}_1=\frac{3}{5},{\mathrm{\ϵ}}_2=\frac{5}{8},{\mathrm{\ϵ}}_3=\frac{29}{41},{\mathrm{\ϵ}}_4=\frac{67}{100},{\mathrm{\ϵ}}_5=\frac{37}{97},{\mathrm{\ϵ}}_6=\frac{59}{101},$ γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ＝0.9_，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，为步骤4得到的第m次迭代的缠绕次数，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率。

如果m＝MI，MI＝6，则终止迭代。

最后得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位即为最终的解缠相位。

Claims (1)

1.一种迭代最大似然估计多基线InSAR相位解缠方法，其特征是它包括以下步骤：

步骤1、初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数：

初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数，包括：多基线InSAR方位向点数，记为N_x；多基线InSAR距离向点数，记为N_y；多基线InSAR待解缠的缠绕相位，记为ω(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点；多基线InSAR非待解缠的缠绕相位组数，记为t；多基线InSAR非待解缠的缠绕相位，记为β(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，其中i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位；多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，记为ε_i，i＝1,2,…,t；多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，记为γ；迭代最大次数，记为MI；相位加入的高斯噪声的标准差，记为λ；

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

采用公式ψ(a,r)＝λ·randn(N_x,N_y)+ω(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，计算得到加入高斯噪声的待解缠的缠绕相位，记为ψ(a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，λ为步骤1初始化的相位加入的高斯噪声的标准差，ω(a,r)为步骤1初始化的多基线InSAR待解缠的缠绕相位，randn(N_x,N_y)表示MATLAB产生的均值为0、标准差为1的高斯噪声，并且为N_x×N_y的随机矩阵；

采用公式φ(i；a,r)＝λ·randn(N_x,N_y)+β(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，计算得到加入高斯噪声的非待解缠的缠绕相位，记为φ(i；a,r)，a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，其中i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，β(i；a,r)为步骤1初始化的多基线InSAR非待解缠的缠绕相位；

步骤3、利用最大似然函数估计多基线InSAR解缠相位：

多基线InSAR待解缠的缠绕相位的缠绕次数，记为k；采用公式

$\begin{array}{c}{\hat{k}}_0=\arg \max \underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_a，r＝1,2,…,N_r，i＝1,2,…,t)，计算得到最大似然函数估计的缠绕次数，记为其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率，argmax(·)表示求

$\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值；

采用公式a＝1,2,…,N_x，r＝1,2,…,N_y，计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，记为a＝1,2,…,N_x，r＝1,2,…,N_y；

将取代公式中的k：

$\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

得到：

$\begin{array}{c}{A}_0=\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t)，令：A₀＝0；

步骤4、计算较大的似然函数值估计的解缠相位：

采用公式

$\begin{array}{c}{\hat{k}}_m=\arg \max \underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，m＝1,2,…,MI)，计算得到第m次迭代的缠绕次数，记为m＝1,2,…,MI，其中m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数，a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率，argmax(·)表示求

$\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值；

采用公式a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，计算得到方位向第a个、距离向第r个点的第m次迭代的多基线InSAR解缠相位，记为r＝1,2,…,N_y，m＝1,2,…,MI；

采用公式

a＝1,2,…,N_x，r＝1,2,…,N_y，计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，记为SUM₀，其中a-1表示方位向第a-1个点，a+1表示方位向第a+1个点，r-1表示距离向第r-1个点，r+1表示距离向第r+1个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，为步骤3得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，|·|表示求绝对值；

采用公式

a＝1,2,…,N_x，r＝1,2,…,N_y，m＝1,2,…,MI，计算得到方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，记为SUMm，m＝1,2,…,MI；

步骤5、优化最大似然估计的多基线InSAR解缠相位：

比较SUM₀和SUM_m的大小：

如果m＝1,2,…,MI，则令：a＝1,2,…,N_x，r＝1,2,…,N_y，优化方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位，其中a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数，SUM₀为步骤4得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和，SUM_m为步骤4得到的方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和；

如果SUM_m≥SUM₀，m＝1,2,…,MI，的值不变；

步骤6、判断算法迭代条件：

如果m＜MI，迭代次数m加1，将代入

$\begin{array}{c}\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

得到

$\begin{array}{c}{A}_m=\underset{i=1}{\overset{t}{\mathrm{\Π}}}(\frac{1}{2\mathrm{\π}}\frac{1-{\left|\mathrm{\γ}\right|}^2}{1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\πk}+\mathrm{\ψ}(a,r)))\right)}^2}\·\\ (1+\frac{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}{\hat{k}}_m+\mathrm{\ψ}(a,r)))\right){\cos }^{-1}(-\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}{\hat{k}}_m+\mathrm{\ψ}(a,r)))\right))}{{(1-{\left(\right|\mathrm{\γ}\left|\cos (\mathrm{\φ}(i;a,r)-{\mathrm{\ϵ}}_i(2\mathrm{\π}{\hat{k}}_m+\mathrm{\ψ}(a,r)))\right)}^2)}^{1/2}}))\end{array}$

(a＝1,2,…,N_x，r＝1,2,…,N_y，i＝1,2,…,t，m＝1,2,…,MI)，令：A_m＝0，重复步骤4至步骤6进行迭代计算，其中m为正整数，表示第m次迭代，MI为步骤1初始化的迭代最大次数，a和r为正整数，a表示方位向第a个点，r表示距离向第r个点，N_x为步骤1初始化的多基线InSAR方位向点数，N_y为步骤1初始化的多基线InSAR距离向点数，i为正整数，表示多基线InSAR第i组非待解缠的缠绕相位，t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数，ε_i，i＝1,2,…,t，为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数，γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数，φ(i；a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位，ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位，为步骤4得到的第m次迭代的缠绕次数，Π(·)为累乘符号，|·|表示求绝对值，cos(·)为余弦函数，cos^-1(·)为反余弦函数，π为圆周率；

如果m＝MI，则终止迭代；

经过以上步骤得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位即为最终的解缠相位。