[发明专利]一种迭代曲面预测InSAR成像及高度估计方法

摘要

本发明公开了一种迭代曲面预测InSAR成像及高度估计方法，该方法利用后向投影成像算法，将InSAR原始回波数据投影到每次迭代估计得到的地形高度曲面上进行投影成像，从而得到该地形高程曲面下的干涉相位，再利用该干涉相位进行高程地形反演，然后对地形高度曲面进行高度补偿，通过多次迭代逼近处理，利用在得到的地形高度曲面上投影成像逐次迭代逼近估计真实的地形高度，从而实现了InSAR数据成像和高度估计一体化，提高了InSAR成像与高度估计精度，而且校正了地形高度起伏引起的几何畸变，即减小了地形几何失真，得到最终的InSAR成像结果和地形高程。

主权项

一种迭代曲面预测InSAR成像及高度估计方法，其特征是它包括以下几个步骤：步骤1、初始化InSAR成像系统参数：初始化InSAR成像系统参数，包括：雷达载波波长，记为λ，雷达平台主天线发射信号带宽，记为B，雷达发射脉冲时宽，记为Tr，雷达采样频率，记为Fs，主天线雷达入射角，记为θ，两组天线之间的基线长，记为Bl，基线与水平方向夹角，记为α，雷达脉冲重复频率，记为PRF，平台运动速度矢量，记为Vr，雷达系统距离向采样点数，记为Nr，雷达系统方位向采样点数，记为Na，雷达平台主天线初始位置，记为PM(0)，雷达平台副天线初始位置，记为PS(0)；上述参数中，雷达系统工作的信号波长λ、雷达平台主天线发射的信号带宽B、雷达平台主天线发射的脉冲时宽Tr、雷达平台接收系统的采样频率Fs、主天线雷达入射角θ、两组天线之间的基线长Bl、基线与水平方向夹角α，雷达脉冲重复频率PRF，在InSAR雷达系统设计过程中已经确定；平台运动速度矢量Vr，雷达系统距离向采样点数Nr、雷达系统方位向采样点数Na、雷达平台主天线初始位置PM(0)、雷达平台副天线初始位置PS(0)在InSAR雷达成像观测方案设计中已经确定；步骤2、获取InSAR原始回波数据：InSAR主天线在距离向第t个快时刻方位向第k个慢时刻的原始回波数据，记为Em(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，其中t和k为自然数，t表示距离向第t个快时刻，k表示方位向第k个慢时刻，Nr为步骤1初始化得到的雷达系统距离向采样点数，Na为步骤1初始化得到的雷达系统方位向采样点数；InSAR副天线在距离向第t个快时刻方位向第k个慢时刻的原始回波数据，记为Es(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na；在InSAR实际成像中，InSAR主天线在距离向第t个快时刻方位向第k个慢时刻的原始回波数据Em(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，和InSAR副天线在距离向第t个快时刻方位向第k个慢时刻的原始回波数据Es(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，由InSAR系统数据接收机提供；步骤3、对InSAR原始回波数据进行距离压缩：采用传统标准合成孔径雷达标准距离压缩方法对步骤2中得到的InSAR主天线在距离向第t个快时刻方位向第k个慢时刻的原始回波数据Em(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，进行距离压缩得到InSAR主天线在距离向第t个快时刻方位向第k个慢时刻的距离压缩后数据，记为Sm(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，其中t表示距离向第t个快时刻，k表示方位向第k个慢时刻；采用传统标准合成孔径雷达标准距离压缩方法对步骤2中得到的InSAR副天线在距离向第t个快时刻方位向第k个慢时刻的原始回波数据Es(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，进行距离压缩得到InSAR副天线在距离向第t个快时刻方位向第k个慢时刻的距离压缩后数据，记为Ss(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na；步骤4、初始化雷达投影成像曲面参数：InSAR雷达投影成像曲面选择为地表水平面坐标系，该坐标系水平横轴记为X轴，该坐标系水平纵轴记为Y轴，雷达投影成像曲面的中心坐标位于[0,0,0]，雷达投影成像曲面的X轴向分辨单元数，记为Nx，雷达投影成像曲面的Y轴向分辨单元数，记为Ny，雷达投影成像曲面的X轴向成像范围，记为Wx，雷达投影成像曲面的Y轴向成像范围，记为Wy，将投影成像曲面均匀等间隔进行划分，投影成像曲面中第a个X轴向和第r个Y轴向分辨单元的位置矢量，记为PT(a,r)＝[x(a,r),y(a,r),z(a,r)]，a＝1,…,Nx，r＝1,…,Ny，其中a和r均为自然数，a表示投影成像曲面的第a个X轴向分辨单元，r表示投影成像曲面的第r个Y轴向分辨单元，x(a,r)、y(a,r)和z(a,r)分别表示投影成像曲面中第a个X轴向第r个Y轴向分辨单元的X轴向位置、Y轴向位置和高度；步骤5、初始化迭代曲面预测InSAR成像及高度估计方法所需的参数：初始化迭代曲面预测InSAR成像及高度估计方法的参数，包括：成像处理迭代的最大次数，记为MI，成像处理迭代阈值，记为ε，第i次迭代过程中雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度，记为Zi(a,r)，a＝1,…,Nx，r＝1,…,Ny，i＝0,…,MI，其中i为自然数，表示成像算法的第i次迭代，Nx为步骤4初始化得到的雷达投影成像曲面的X轴向分辨单元数，Ny为步骤4中初始化得到的雷达投影成像曲面的Y轴向分辨单元数；Z0(a,r)为第0次迭代过程中雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度，即该步骤中初始化得到雷达投影成像曲面高度，初始化雷达投影成像曲面高度Z0(a,r)的值为Z0(a,r)＝0，a＝1,…,Nx，r＝1,…,Ny；步骤6、计算雷达投影成像曲面中分辨单元到平台的斜距：在成像处理第i次迭代中，若i＝1，则更新投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度为z(a,r)＝Z0(a,r)，更新得到投影成像曲面中第a个X轴向和第r个Y轴向分辨单元的位置矢量PT(a,r)＝[x(a,r),y(a,r),Z0(a,r)]；若i＞1，则更新投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度为z(a,r)＝Zi‑1(a,r)，得到PT(a,r)＝[x(a,r),y(a,r),Zi‑1(a,r)]，其中Z0(a,r)为步骤5中初始化得到的第0次迭代过程中雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度Z0(a,r)＝0，a＝1,…,Nx，r＝1,…,Ny，而Zi‑1(a,r)为成像处理第i‑1次迭代过程中雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度值；采用公式PM(k)＝PM(0)+Vr·k/PRF，k＝1,2,…,Na，计算得到雷达平台主天线第k个慢时刻的位置，记为PM(k)，k＝1,2,…,Na，其中PM(0)为步骤1初始化得到的雷达平台主天线初始位置，Vr为步骤1初始化得到的平台运动速度矢量，PRF为步骤1初始化得到的雷达脉冲重复频率，k表示方位向第k个慢时刻，Na为步骤1初始化得到的雷达系统方位向采样点数；采用公式PS(k)＝PS(0)+Vr·k/PRF，k＝1,2,…,Na，计算得到雷达平台副天线第k个慢时刻的位置，记为PS(k)，k＝1,2,…,Na，其中PS(0)为步骤1初始化得到的雷达平台副天线初始位置；采用公式Rm(k；a,r)＝||PM(k)‑PT(a,r)||2，a＝1,…,Nx，r＝1,…,Ny，k＝1,2,…,Na，计算得到雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元到InSAR主天线的斜距，记为Rm(k；a,r)，a＝1,…,Nx，r＝1,…,Ny，k＝1,2,…,Na，其中||·||2表示向量2范数运算符号，PT(a,r)为步骤4得到的投影成像曲面中第a个X轴向和第r个Y轴向分辨单元的位置矢量；采用公式Rs(k；a,r)＝||PS(k)‑PT(a,r)||2，a＝1,…,Nx，r＝1,…,Ny，k＝1,2,…,Na，计算得到雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元到InSAR副天线的斜距，记为Rs(k；a,r)，a＝1,…,Nx，r＝1,…,Ny，k＝1,2,…,Na；步骤7、采用标准后向投影成像算法进行成像：对步骤3得到的InSAR主天线在距离向第t个快时刻方位向第k个慢时刻的距离压缩后数据Sm(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，以及步骤6得到的雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元到InSAR主天线的斜距Rm(k；a,r)，a＝1,…,Nx，r＝1,…,Ny，采用传统标准后向投影成像算法进行成像，得到雷达平台主天线数据在投影成像曲面的成像结果，记为Bm(a,r)，a＝1,…,Nx，r＝1,…,Ny；对步骤3得到的InSAR副天线在距离向第t个快时刻方位向第k个慢时刻的距离压缩后数据Ss(t,k)，t＝1,2,…,Nr，k＝1,2,…,Na，以及步骤6得到的雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元到InSAR副天线的斜距Rs(k；a,r)，a＝1,…,Nx，r＝1,…,Ny，采用传统标准后向投影成像算法进行成像，得到雷达平台副天线数据在投影成像曲面的成像结果，记为Bs(a,r)，a＝1,…,Nx，r＝1,…,Ny；步骤8、InSAR主副天线干涉相位提取：采用公式φ(a,r)＝∠((Bm(a,r))*·Bs(a,r))，a＝1,…,Nx，r＝1,…,Ny，计算得到投影成像曲面下InSAR主副天线成像结果的干涉相位，记为φ(a,r)，a＝1,…,Nx，r＝1,…,Ny，其中Bm(a,r)和Bs(a,r)为步骤7中得到的雷达平台主天线与副天线数据在投影成像曲面的成像结果，∠(·)表示求解相位运算符号，上标*表示求解共轭运算符号；步骤9、InSAR干涉相位解缠：采用传统标准干涉相位解缠方法对步骤8中得到的投影成像曲面下InSAR主副天线成像结果的干涉相位φ(a,r)，a＝1,…,Nx，r＝1,…,Ny，进行干涉相位解缠处理，得到解缠后的InSAR干涉相位，记为a＝1,…,Nx，r＝1,…,Ny，其中a表示投影成像曲面第a个X轴向分辨单元，r表示投影成像曲面第r个Y轴向分辨单元；步骤10、计算雷达投影成像曲面条件下的地形高度：采用公式计算得到投影成像曲面第a个X轴向第r个Y轴向分辨单元P(a,r)的地形高度，记为H(a,r)，a＝1,…,Nx，r＝1,…,Ny，其中a表示投影成像曲面第a个X轴向分辨单元，r表示投影成像曲面第r个Y轴向分辨单元，λ为步骤1中初始化得到的雷达载波波长，为步骤6中时得到的雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元到InSAR主天线的斜距，ceil(·)表示四舍五入运算符号，sin(·)表示正弦函数运算符号，θ为步骤1中初始化得到的主天线雷达入射角，Bl为步骤1中初始化得到的两组天线之间的基线长，cos(·)表示余弦函数运算符号，α为步骤1中初始化得到的基线与水平方向夹角，为步骤9得到的解缠后的InSAR干涉相位，π为圆周率；步骤11、判定算法迭代条件：采用公式计算得到雷达投影成像曲面条件下地形高度绝对值的平均值，记为Hmean；采用公式Zi(a,r)＝H(a,r)+Zi‑1(a,r)，a＝1,…,Nx，r＝1,…,Ny，将地形高度H(a,r)，a＝1,…,Nx，r＝1,…,Ny，补偿到成像处理中第i次迭代过程中雷达投影成像曲面高度Zi(a,r)，a＝1,…,Nx，r＝1,…,Ny；如果平均值Hmean＞ε且i≤MI，将成像处理中迭代次数i加1得到i←i+1，再重复执行步骤6至步骤11；否则终止成像处理迭代，Zi(a,r)即为雷达投影成像曲面的最终地形高程，其中i表示成像处理第i次迭代，MI为步骤5中初始化得到算法迭代的最大次数，∑(·)表示向量求和运算符号，|·|表示求绝对值运算符号，←表示赋值运算符号，ε为步骤5中初始化得到的成像处理迭代阈值，H(a,r)，a＝1,…,Nx，r＝1,…,Ny，为步骤10中得到的投影成像曲面第a个X轴向第r个Y轴向分辨单元P(a,r)的地形高度，Zi‑1(a,r)为成像处理第i‑1次迭代过程中雷达投影成像曲面中第a个X轴向第r个Y轴向分辨单元的高度值。