CN104808203B - 一种迭代最大似然估计多基线InSAR相位解缠方法 - Google Patents
一种迭代最大似然估计多基线InSAR相位解缠方法 Download PDFInfo
- Publication number
- CN104808203B CN104808203B CN201510094461.4A CN201510094461A CN104808203B CN 104808203 B CN104808203 B CN 104808203B CN 201510094461 A CN201510094461 A CN 201510094461A CN 104808203 B CN104808203 B CN 104808203B
- Authority
- CN
- China
- Prior art keywords
- phase place
- insar
- solution
- many
- baseline insar
- 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.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Transmission And Conversion Of Sensor Element Output (AREA)
- Measuring Phase Differences (AREA)
Abstract
本发明提出了一种迭代最大似然估计多基线InSAR相位解缠方法,该方法首先采用最大似然函数估计方法,得到多基线InSAR初始的解缠相位,再求出几组较大似然函数值对应的解缠相位,然后通过多组较大似然函数对应的解缠相位迭代优化初始解缠相位,减小相位解缠过程中噪声带来的误差,与传统的最大似然估计多基线相位解缠方法相比,本发明方法在噪声较大情况下,依然能得到高精度的解缠相位,为后续的高精度InSAR高程反演提供保证。
Description
技术领域
本发明属于雷达技术领域,它特别涉及干涉合成孔径雷达(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为离散型随机变量时为概率分布),θ为待估参数,X1,X2,….Xn是来自于总体X的样本,x1,x2,…,xn为样本X1,X2,….Xn的一个观察值,则样本的联合分布(当X是连续型随机变量时为概率密度,当X为离散型随机变量时为概率分布)称为似然函数。详见文献“概率论与数理统计”,徐全智等编著,电子科技大学出版社。
定义4、缠绕相位:
缠绕相位是指在实际干涉处理中,经过三角运算,得到的干涉相位的主值。详见文献“星载合成孔径雷达干涉测量”,王超等编著,科学出版社。
定义5、解缠相位:
解缠相位是指干涉处理中从缠绕相位经过相位解缠处理恢复的真实相位。详见文献“星载合成孔径雷达干涉测量”,王超等编著,科学出版社。
定义6、高斯分布函数randn(m,n)
高斯分布函数randn(m,n)是MATLAB中产生均值为0,标准差为1的高斯分布,并且为m×n的随机矩阵。
本发明提供的一种迭代最大似然估计多基线InSAR相位解缠方法,它包括以下几个步骤:
步骤1、初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数:
初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数,包括:多基线InSAR方位向点数,记为Nx;多基线InSAR距离向点数,记为Ny;多基线InSAR待解缠的缠绕相位,记为ω(a,r),a=1,2,…,Nx,r=1,2,…,Ny,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点;多基线InSAR非待解缠的缠绕相位组数,记为t;多基线InSAR非待解缠的缠绕相位,记为β(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,其中i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位;多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数,记为εi,i=1,2,…,t;多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数,记为γ;迭代最大次数,记为MI;相位加入的高斯噪声的标准差,记为λ。
步骤2、将缠绕相位加入高斯噪声:
采用公式ψ(a,r)=λ·randn(Nx,Ny)+ω(a,r),a=1,2,…,Nx,r=1,2,…,Ny,计算得到加入高斯噪声的待解缠的缠绕相位,记为ψ(a,r),a=1,2,…,Nx,r=1,2,…,Ny,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数,Ny为步骤1初始化的多基线InSAR距离向点数,λ为步骤1初始化的相位加入的高斯噪声的标准差,ω(a,r)为步骤1初始化的多基线InSAR待解缠的缠绕相位,randn(Nx,Ny)表示MATLAB产生的均值为0、标准差为1的高斯噪声,并且为Nx×Ny的随机矩阵。
采用公式φ(i;a,r)=λ·randn(Nx,Ny)+β(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,计算得到加入高斯噪声的非待解缠的缠绕相位,记为φ(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,其中i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位,t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数,β(i;a,r)为步骤1初始化的多基线InSAR非待解缠的缠绕相位。
步骤3、利用最大似然函数估计多基线InSAR解缠相位:
多基线InSAR待解缠的缠绕相位的缠绕次数,记为k;采用公式
(a=1,2,…,Na,r=1,2,…,Nr,i=1,2,…,t),计算得到最大似然函数估计的缠绕次数,记为其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤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得到的加入高斯噪声的待解缠的缠绕相位,∏(·)为累乘符号,|·|表示求绝对值,cos(·)为余弦函数,cos-1(·)为反余弦函数,π为圆周率,argmax(·)表示求
取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。
采用公式a=1,2,…,Nx,r=1,2,…,Ny,计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,记为a=1,2,…,Nx,r=1,2,…,Ny;
将取代公式中的k:
得到:
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t),令:A0=0。
步骤4、计算较大的似然函数值估计的解缠相位:
采用公式
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,m=1,2,…,MI),计算得到第m次迭代的缠绕次数,记为m=1,2,…,MI,其中m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数,a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤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得到的加入高斯噪声的待解缠的缠绕相位,Π(·)为累乘符号,|·|表示求绝对值,cos(·)为余弦函数,cos-1(·)为反余弦函数,π为圆周率,argmax(·)表示求
取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。
采用公式a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,计算得到方位向第a个、距离向第r个点的第m次迭代的多基线InSAR解缠相位,记为r=1,2,…,Ny,m=1,2,…,MI。
采用公式
a=1,2,…,Nx,r=1,2,…,Ny,计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,记为SUM0,其中a-1表示方位向第a-1个点,a+1表示方位向第a+1个点,r-1表示距离向第r-1个点,r+1表示距离向第r+1个点,Nx为步骤1初始化的多基线InSAR方位向点数,Ny为步骤1初始化的多基线InSAR距离向点数,为步骤3得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,|·|表示求绝对值。
采用公式
a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,计算得到方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,记为SUMm,m=1,2,…,MI。
步骤5、优化最大似然估计的多基线InSAR解缠相位,
比较SUM0和SUMm的大小:
如果SUMm<SUM0,m=1,2,…,MI,则令:a=1,2,…,Nx,r=1,2,…,Ny,优化方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数,Ny为步骤1初始化的多基线InSAR距离向点数,m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数,SUM0为步骤4得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,SUMm为步骤4得到的方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和。
如果SUMm≥SUM0,m=1,2,…,MI,的值不变。
步骤6、判断算法迭代条件
如果m<MI,迭代次数m加1,将代入
得到
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,m=1,2,…,MI),令:Am=0,重复步骤4至步骤6进行迭代计算,其中m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数,a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤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示意图
其中A0为多基线InSAR主天线,A1,A2,A3为多基线InSAR副天线,p为场景中心点目标,H为平台高度,θ为主天线雷达入射角,α为基线与水平方向夹角。
图3为本发明所提供方法的多基线InSAR相位解缠参数值
具体实施方式
本发明主要采用仿真实验的方法进行验证,所有步骤、结论都在MATLABR2013b软件上验证正确。具体实施步骤如下:
步骤1、初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数:
初始化迭代最大似然估计多基线InSAR相位解缠方法所需参数,包括:多基线InSAR方位向点数Nx=500;多基线InSAR距离向点数,记为Ny=500;多基线InSAR待解缠的缠绕相位,记为ω(a,r),a=1,2,…,Nx,r=1,2,…,Ny,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点;多基线InSAR非待解缠的缠绕相位组数t=6;多基线InSAR非待解缠的缠绕相位,记为β(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,其中i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位;多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ=0.9;迭代最大次数MI=6;相位加入的高斯噪声的标准差λ=0.8。
步骤2、将缠绕相位加入高斯噪声:
采用公式ψ(a,r)=λ·randn(Nx,Ny)+ω(a,r),a=1,2,…,Nx,r=1,2,…,Ny,Nx=500,Ny=500,计算得到加入高斯噪声的待解缠的缠绕相位,记为ψ(a,r),a=1,2,…,Nx,r=1,2,…,Ny,Nx=500,Ny=500,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数Nx=500,Ny为步骤1初始化的多基线InSAR距离向点数Ny=500,λ为步骤1初始化的相位加入的高斯噪声的标准差λ=0.8,ω(a,r)为步骤1初始化的多基线InSAR待解缠的缠绕相位,randn(Nx,Ny)表示MATLAB产生的均值为0、标准差为1的高斯噪声,并且为Nx×Ny的随机矩阵。
采用公式φ(i;a,r)=λ·randn(Nx,Ny)+β(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,Nx=500,Ny=500,t=6,计算得到加入高斯噪声的非待解缠的缠绕相位,记为φ(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,Nx=500,Ny=500,t=6,其中i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位,t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t=6,β(i;a,r)为步骤1初始化的多基线InSAR非待解缠的缠绕相位。
步骤3、利用最大似然函数估计多基线InSAR解缠相位:
多基线InSAR待解缠的缠绕相位的缠绕次数,记为k。采用公式
(a=1,2,…,Na,r=1,2,…,Nr,i=1,2,…,t,Nx=500,Ny=500,t=6),计算得到最大似然函数估计的缠绕次数,记为其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数Nx=500,Ny为步骤1初始化的多基线InSAR距离向点数Ny=500,i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位,t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t=6,εi,i=1,2,…,t,为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ=0.9,φ(i;a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位,ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位,∏(·)为累乘符号,·|表示求绝对值,cos(·)为余弦函数,cos-1(·)为反余弦函数,π为圆周率,argmax(·)表示求
取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。。
采用公式a=1,2,…,Nx,r=1,2,…,Ny,Nx=500,Ny=500,计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,记为a=1,2,…,Nx,r=1,2,…,Ny,Nx=500,Ny=500。
将代入公式
得到:
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,Nx=500,Ny=500,t=6),令:A0=0。
步骤4、计算迭代所需较大的似然函数值估计的解缠相位:
采用公式
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,m=1,2,…,MI,Nx=500,Ny=500,t=6,MI=6),计算得到第m次迭代的缠绕次数,记为m=1,2,…,MI,MI=6,其中m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数MI=6,a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数Nx=500,Ny为步骤1初始化的多基线InSAR距离向点数Ny=500,i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位,t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t=6,εi,i=1,2,…,t,为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数γ为步骤1初始化的多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数γ=0.9,φ(i;a,r)为步骤2得到的加入高斯噪声的非待解缠的缠绕相位,ψ(a,r)为步骤2得到的加入高斯噪声的待解缠的缠绕相位,∏(·)为累乘符号,·|表示求绝对值,cos(·)为余弦函数,cos-1(·)为反余弦函数,π为圆周率,argmax(·)表示求
取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值。
采用公式a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,Nx=500,Ny=500,MI=6,计算得到方位向第a个、距离向第r个点的第m次迭代的多基线InSAR解缠相位,记为a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,Nx=500,Ny=500,MI=6。
采用公式
a=1,2,…,Nx,r=1,2,…,Ny,Nx=500,Ny=500,计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,记为SUM0,其中a-1表示方位向第a-1个点,a+1表示方位向第a+1个点,r-1表示距离向第r-1个点,r+1表示距离向第r+1个点,Nx为步骤1初始化的多基线InSAR方位向点数Nx=500,Ny为步骤1初始化的多基线InSAR距离向点数Ny=500,为步骤3得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,|·|表示求绝对值。
采用公式
a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,Nx=500,Ny=500,MI=6,计算得到方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,记为SUMm,m=1,2,…,MI,MI=6。
步骤5、优化最大似然估计的多基线InSAR解缠相位,
比较SUM0和SUMm的大小:
如果SUMm<SUM0,m=1,2,…,MI,MI=6,则令:a=1,2,…,Nx,r=1,2,…,Ny,Nx=500,Ny=500,优化方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数Nx=500,Ny为步骤1初始化的多基线InSAR距离向点数Ny=500,m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数MI=6,SUM0为步骤4得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,SUMm为步骤4得到的方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和。
如果SUMm≥SUM0,m=1,2,…,MI,MI=6,则不进行任何计算,的值不变。
步骤6、判断算法迭代条件
如果m<MI,MI=6,迭代次数m加1,将代入
得到:
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,m=1,2,…,MI,Nx=500,Ny=500,t=6,MI=6),令:Am=0,重复步骤4至步骤6进行迭代计算,其中m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数MI=6,a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数Nx=500,Ny为步骤1初始化的多基线InSAR距离向点数Ny=500,i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位,t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数t=6,εi,i=1,2,…,t,为步骤1初始化的多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数γ为步骤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方位向点数,记为Nx;多基线InSAR距离向点数,记为Ny;多基线InSAR待解缠的缠绕相位,记为ω(a,r),a=1,2,…,Nx,r=1,2,…,Ny,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点;多基线InSAR非待解缠的缠绕相位组数,记为t;多基线InSAR非待解缠的缠绕相位,记为β(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,其中i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位;多基线InSAR待求的解缠相位和多基线InSAR非待求的解缠相位的相关系数,记为εi,i=1,2,…,t;多基线InSAR非待解缠的缠绕相位对应的副天线和多基线InSAR待求的解缠相位的相关系数,记为γ;迭代最大次数,记为MI;相位加入的高斯噪声的标准差,记为λ;
步骤2、将缠绕相位加入高斯噪声:
采用公式ψ(a,r)=λ·randn(Nx,Ny)+ω(a,r),a=1,2,…,Nx,r=1,2,…,Ny,计算得到加入高斯噪声的待解缠的缠绕相位,记为ψ(a,r),a=1,2,…,Nx,r=1,2,…,Ny,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数,Ny为步骤1初始化的多基线InSAR距离向点数,λ为步骤1初始化的相位加入的高斯噪声的标准差,ω(a,r)为步骤1初始化的多基线InSAR待解缠的缠绕相位,randn(Nx,Ny)表示MATLAB产生的均值为0、标准差为1的高斯噪声,并且为Nx×Ny的随机矩阵;
采用公式φ(i;a,r)=λ·randn(Nx,Ny)+β(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,计算得到加入高斯噪声的非待解缠的缠绕相位,记为φ(i;a,r),a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,其中i为正整数,表示多基线InSAR第i组非待解缠的缠绕相位,t为步骤1初始化的多基线InSAR非待解缠的缠绕相位组数,β(i;a,r)为步骤1初始化的多基线InSAR非待解缠的缠绕相位;
步骤3、利用最大似然函数估计多基线InSAR解缠相位;
多基线InSAR待解缠的缠绕相位的缠绕次数,记为k;采用公式
(a=1,2,…,Na,r=1,2,…,Nr,i=1,2,…,t),计算得到最大似然函数估计的缠绕次数,记为其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤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得到的加入高斯噪声的待解缠的缠绕相位,Π(·)为累乘符号,|·|表示求绝对值,cos(·)为余弦函数,cos-1(·)为反余弦函数,π为圆周率,argmax(·)表示求
取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值;
采用公式a=1,2,…,Nx,r=1,2,…,Ny,计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,记为a=1,2,…,Nx,r=1,2,…,Ny;
将取代公式中的k:
得到:
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t),令:A0=0;
步骤4、计算较大的似然函数值估计的解缠相位:
采用公式
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,m=1,2,…,MI),计算得到第m次迭代的缠绕次数,记为m=1,2,…,MI,其中m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数,a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤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得到的加入高斯噪声的待解缠的缠绕相位,Π(·)为累乘符号,|·|表示求绝对值,cos(·)为余弦函数,cos-1(·)为反余弦函数,π为圆周率,argmax(·)表示求
取最大值时的多基线InSAR待解缠的缠绕相位的缠绕次数k的值;
采用公式a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,计算得到方位向第a个、距离向第r个点的第m次迭代的多基线InSAR解缠相位,记为r=1,2,…,Ny,m=1,2,…,MI;
采用公式
a=1,2,…,Nx,r=1,2,…,Ny,计算得到方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,记为SUM0,其中a-1表示方位向第a-1个点,a+1表示方位向第a+1个点,r-1表示距离向第r-1个点,r+1表示距离向第r+1个点,Nx为步骤1初始化的多基线InSAR方位向点数,Ny为步骤1初始化的多基线InSAR距离向点数,为步骤3得到的方位向第a个、距离向第r个点的最大似然函数估计的 多基线InSAR解缠相位,|·|表示求绝对值;
采用公式
a=1,2,…,Nx,r=1,2,…,Ny,m=1,2,…,MI,计算得到方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,记为SUMm,m=1,2,…,MI;
步骤5、优化最大似然估计的多基线InSAR解缠相位:
比较SUM0和SUMm的大小:
如果SUMm<SUM0,m=1,2,…,MI,则令:a=1,2,…,Nx,r=1,2,…,Ny,优化方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位,其中a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤1初始化的多基线InSAR方位向点数,Ny为步骤1初始化的多基线InSAR距离向点数,m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数,SUM0为步骤4得到的方位向第a个、距离向第r个点的最大似然函数估计的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻的8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和,SUMm为步骤4得到的方位向第a个、距离向第r个点的第m次迭代优化的多基线InSAR解缠相位与方位向第a个、距离向第r个点相邻8个点的最大似然函数估计的多基线InSAR解缠相位之差的绝对值之和;
如果SUMm≥SUM0,m=1,2,…,MI,的值不变;
步骤6、判断算法迭代条件
如果m<MI,迭代次数m加1,将代入
得到
(a=1,2,…,Nx,r=1,2,…,Ny,i=1,2,…,t,m=1,2,…,MI),令:Am=0,重复步骤4至步骤6进行迭代计算,其中m为正整数,表示第m次迭代,MI为步骤1初始化的迭代最大次数,a和r为正整数,a表示方位向第a个点,r表示距离向第r个点,Nx为步骤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解缠相位即为最终的解缠相位。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510094461.4A CN104808203B (zh) | 2015-03-03 | 2015-03-03 | 一种迭代最大似然估计多基线InSAR相位解缠方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510094461.4A CN104808203B (zh) | 2015-03-03 | 2015-03-03 | 一种迭代最大似然估计多基线InSAR相位解缠方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104808203A CN104808203A (zh) | 2015-07-29 |
CN104808203B true CN104808203B (zh) | 2017-05-10 |
Family
ID=53693190
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510094461.4A Expired - Fee Related CN104808203B (zh) | 2015-03-03 | 2015-03-03 | 一种迭代最大似然估计多基线InSAR相位解缠方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104808203B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107193005B (zh) * | 2017-06-16 | 2020-03-17 | 桂林电子科技大学 | 一种无损卡尔曼滤波与粒子滤波相结合的相位展开算法 |
CN107102333B (zh) * | 2017-06-27 | 2020-01-10 | 北京航空航天大学 | 一种星载InSAR长短基线融合解缠方法 |
CN107390218B (zh) * | 2017-08-28 | 2020-06-05 | 西安电子科技大学 | 基于最小无穷范数的二维相位解缠绕方法 |
CN107783079B (zh) * | 2017-09-28 | 2021-10-26 | 淮海工学院 | 一种使用l0范数代价函数的二维相位快速解缠方法 |
CN108008383A (zh) * | 2017-11-09 | 2018-05-08 | 电子科技大学 | 一种迭代多基线高精度四次fft相位解缠方法 |
CN109633648B (zh) * | 2019-01-22 | 2022-04-29 | 北京航空航天大学 | 一种基于似然估计的多基线相位估计装置及方法 |
CN110161502B (zh) * | 2019-05-28 | 2020-10-27 | 北京邮电大学 | 一种星载多基线InSAR叠加数据的滤波方法及装置 |
CN114265062B (zh) * | 2021-11-11 | 2023-11-10 | 电子科技大学 | 一种基于相位梯度估计网络的InSAR相位解缠方法 |
CN113866768B (zh) * | 2021-12-02 | 2022-04-15 | 深圳大学 | 一种时序干涉雷达相位优化估计方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102508245A (zh) * | 2011-11-18 | 2012-06-20 | 北京航空航天大学 | 一种星载多频率与多基线InSAR高程估计精度等效性确定方法 |
CN103885059A (zh) * | 2014-01-26 | 2014-06-25 | 中国测绘科学研究院 | 一种多基线干涉合成孔径雷达三维重建方法 |
EP2784537A1 (en) * | 2013-05-15 | 2014-10-01 | Institute of Electronics, Chinese Academy of Sciences | Inversion method and apparatus based on polarimetric interferometric synthetic aperture radar |
-
2015
- 2015-03-03 CN CN201510094461.4A patent/CN104808203B/zh not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102508245A (zh) * | 2011-11-18 | 2012-06-20 | 北京航空航天大学 | 一种星载多频率与多基线InSAR高程估计精度等效性确定方法 |
EP2784537A1 (en) * | 2013-05-15 | 2014-10-01 | Institute of Electronics, Chinese Academy of Sciences | Inversion method and apparatus based on polarimetric interferometric synthetic aperture radar |
CN103885059A (zh) * | 2014-01-26 | 2014-06-25 | 中国测绘科学研究院 | 一种多基线干涉合成孔径雷达三维重建方法 |
Non-Patent Citations (1)
Title |
---|
利用L1范数的多基线InSAR相位解缠绕技术;于瀚雯 等;《西安电子科技大学学报(自然科学版)》;20130831;第40卷(第4期);第37-41页 * |
Also Published As
Publication number | Publication date |
---|---|
CN104808203A (zh) | 2015-07-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104808203B (zh) | 一种迭代最大似然估计多基线InSAR相位解缠方法 | |
CN102662171B (zh) | 一种sar层析三维成像方法 | |
CN103018730B (zh) | 分布式子阵波达方向估计方法 | |
CN106646564B (zh) | 一种基于低轨卫星增强导航方法 | |
CN104730519B (zh) | 一种采用误差迭代补偿的高精度相位解缠方法 | |
CN103885059B (zh) | 一种多基线干涉合成孔径雷达三维重建方法 | |
CN103675790A (zh) | 一种基于高精度DEM提高InSAR技术监测地表形变精度的方法 | |
JP2016507735A (ja) | 衛星レーダー干渉度のイオン歪み補正方法及びその装置 | |
CN107942331B (zh) | 基于谱分析的多通道sar系统通道偏差估计方法 | |
CN104007414A (zh) | 基于平面阵的二维波达方向估计方法和估计器 | |
CN102621549A (zh) | 多基线/多频段干涉相位解缠频域快速算法 | |
CN105158759B (zh) | 基于杂波相消的hrws sar通道相位偏差校正方法 | |
CN104020440A (zh) | 基于l型干涉式线性阵列的二维波达角估计方法 | |
CN103630898B (zh) | 对多基线干涉sar相位偏置进行估计的方法 | |
CN101887121A (zh) | 基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法 | |
CN104181513A (zh) | 一种雷达天线阵元位置的校正方法 | |
CN104251991A (zh) | 一种基于稀疏度估计的分维度阈值迭代稀疏微波成像方法 | |
CN103824287A (zh) | 一种基于稳健平面拟合的相位相关亚像素匹配方法 | |
CN103454636A (zh) | 基于多像素协方差矩阵的差分干涉相位估计方法 | |
CN104730521B (zh) | 一种基于非线性优化策略的SBAS‑DInSAR方法 | |
CN109696651B (zh) | 一种基于m估计的低快拍数下波达方向估计方法 | |
CN105572629B (zh) | 一种适用于任意阵列结构的低运算复杂度的二维测向方法 | |
CN104698448B (zh) | 运动平台下基于流形分离的共形阵列稳健估角方法 | |
CN105353374A (zh) | 一种用于自旋目标的单频雷达成像方法 | |
CN106932773A (zh) | 基于修正嵌入式容积卡尔曼滤波的相位展开算法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
EXSB | Decision made by sipo to initiate substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170510 Termination date: 20200303 |
|
CF01 | Termination of patent right due to non-payment of annual fee |