CN113835090B - 一种基于多通道sar系统的高精度干涉相位获取方法 - Google Patents

一种基于多通道sar系统的高精度干涉相位获取方法 Download PDF

Info

Publication number
CN113835090B
CN113835090B CN202111009890.9A CN202111009890A CN113835090B CN 113835090 B CN113835090 B CN 113835090B CN 202111009890 A CN202111009890 A CN 202111009890A CN 113835090 B CN113835090 B CN 113835090B
Authority
CN
China
Prior art keywords
imaging
azimuth
scene
iteration
distance
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.)
Active
Application number
CN202111009890.9A
Other languages
English (en)
Other versions
CN113835090A (zh
Inventor
张晓玲
王楠
陈益飞
詹旭
韦顺军
师君
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202111009890.9A priority Critical patent/CN113835090B/zh
Publication of CN113835090A publication Critical patent/CN113835090A/zh
Application granted granted Critical
Publication of CN113835090B publication Critical patent/CN113835090B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • 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/904—SAR modes
    • G01S13/9052—Spotlight mode

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)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于多通道SAR系统的高精度干涉相位获取方法，它是通过多通道SAR系统可获取一幅高分辨率和一幅低分辨率图像，本发明通过减少副天线通道为单通道，得到两幅分辨率不同的SAR图像，再由两幅不同分辨率图像获得高精度干涉相位，利用二维频域干涉相位的高稀疏特性，建立了从低分辨图像到二维频域干涉相位的矩阵重构方程，为求解该方程提出了一种矩阵形式稀疏重构算法，最终获得了观测区域的干涉相位。本发明与传统的干涉测量方法比较，有效提高了InSAR复图像干涉相位的精度，且极大地降低了内存和计算量的消耗，大幅度提高计算效率。

Description

一种基于多通道SAR系统的高精度干涉相位获取方法
技术领域
本发明属于雷达技术领域，它特别涉及了一种多通道合成孔径雷达(SAR)系统及其干涉相位获取方法。
背景技术
作为一种工作在微波波段的有源雷达，合成孔径雷达(Synthetic ApertureRadar,SAR)具有全天时、全天候的成像能力，即无论是白天或黑夜、晴天还是雷雨风雪天气，都可以随时随地成像，克服了光学和红外系统不能在晚上和复杂天气条件进行成像的缺点。
干涉合成孔径雷达(InSAR)测量技术是在合成孔径雷达(SAR)的基础上，利用两个或者多个位置不同的天线观测同一个目标场景，根据目标到不同天线的斜距差获得测量数据的干涉相位，再通过平台与地面观测场景的几何关系反演出地面场景的数字高程信息的技术。由于具有全天时、全天候的特点，InSAR已经成为当前提取大面积地表三维图像和地形高程变化信息的一项重要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。
高分辨率SAR图像不仅能提高SAR系统对观测区域的检测能力，而且对于获取高精度的干涉相位和地形高程反演提供了有力的保障，其重要程度不言而喻。扫描式SAR模式是通过牺牲方位向分辨率来换取距离向的宽测绘带，从而会降低SAR成像质量，进而影响InSAR的干涉相位和高程反演精度。多通道SAR模式(沿方位向布设)是解决宽测绘带和高分辨率的一种有效工作模式，它的核心是通过沿平台运动方向布设多通道，以空间采样等效时间采样，在保障图像分辨率的同时，也能兼顾成像测绘的幅宽。该模式要具有InSAR功能，可以在垂直于方位向增设一组多通道，但是会极大提高SAR系统的硬件成本和实现难度，也会大大提高SAR系统对数据采集、传输和存储等方面的压力，因此为了降低系统复杂度和硬件成本等，可在垂直于方位向增设一个通道来实现干涉测量，但会形成不同分辨率的主副SAR图像，而传统的干涉相位测量是使用两幅分辨率相同的SAR图像共轭相乘得到，如果直接对高分辨率进行抽样或对低分辨率进行插值在采用传统的共轭相乘算法都会影响干涉相位精度，无法保证干涉相位的准确提取，进而影响干涉相位的精度。因此本发明在二维频域干涉相位稀疏重构的技术基础上提出了一种基于多通道SAR系统的高精度干涉相位获取方法。
发明内容
本发明基于主天线多通道副天线单通道的SAR系统及二维频域稀疏重构技术，提出了一种基于多通道SAR系统的高精度干涉相位获取方法，通过该多通道SAR系统可获取一幅高分辨率和一幅低分辨率图像，再利用二维频域稀疏重构方法可获取目标场景的高精度干涉相位，该系统在降低了双天线多通道的硬件成本的基础上保证了高精度的干涉相位图，不仅降低了系统成本而且同时降低了SAR系统对数据采集、传输和存储等方面的压力
针对如何从不同分辨率SAR图像获得高精度干涉相位问题，建立了从低分辨图像到频域干涉相位的矢量重构方程。由于方程求解存在克罗内克积计算，使得内存和计算量消耗巨大，针对该问题，利用二维频域干涉相位的高稀疏特性，建立了从低分辨图像到二维频域干涉相位的矩阵重构方程，为求解该方程提出了一种矩阵形式稀疏重构算法，该算法通过建立矩阵式稀疏重构方程，和传统的干涉测量方法比较，有效提高了InSAR复图像干涉相位的精度，且极大地降低了内存和计算量的消耗。
为了方便描述本发明的内容，首先作以下术语定义：
定义1、合成孔径雷达(SAR)
合成孔径雷达是将雷达固定于载荷运动平台上，结合平台的运动以合成等效阵列以实现阵列向的分辨率，再利用雷达波束向回波延时实现距离一维成像，从而实现对观测目标二维成像的一种合成孔径雷达技术，详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。
定义2、标准合成孔径雷达回波数据距离向脉冲压缩
标准合成孔径雷达回波数据距离向脉冲压缩是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行信号聚焦成像的过程。详见文献“雷达成像技术”，保铮，邢孟道，王彤，电子工业出版社，2005。
定义3、范数
设X是复数域上线性空间，其中/>表示复数域，若它满足如下性质：||X||≥0，且当||X||＝0仅有X＝0；||aX||＝|a|||X||，其中a为任意常数；||X1+X2||≤||X1||+||X2||，则称||X||为X空间上的范数(norm)，其中X1和X2为X空间上的任意两个值。对于定义1中的N×1维离散信号向量X＝[x1,x2,…,xN]T，向量X的LP范数表达式为/>其中xi为向量X的第i个元素，∑|·|表示绝对值求和运算符号，向量X的L1范数表达式为向量X的L2范数表达式为/>向量X的L0范数表达式为且xi≠0。详见文献“矩阵理论”，黄廷祝等编著，高等教育出版社出版。
定义4、方位向、距离向
将雷达平台运动的方向叫做方位向，将垂直于方位向的方向叫做距离向。详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。
定义5、散射特性相位
散射特性相位由目标点的散射特性决定，记为详见文献“InSAR技术在地表形变监测中的应用”赵佳曼，2019.
定义6、二维离散傅里叶变换
若f(x,y)∈R×C为一个二维图像，二维离散傅里叶变换是利用公式将二维图像从空间域转换到频域的变换方法，其中 为距离向离散傅里叶变换矩阵， 为方位向离散傅里叶变换矩阵,R为二维图像f(x,y)的行数，C为二维图像f(x,y)的列数。
定义7、哈达玛积
哈达玛积是矩阵的一类运算，若A＝(aij)和B＝(bij)是两个同阶矩阵，若cij＝aij×bij则称矩阵C＝(cij)为A和B的哈达玛积，即
定义8、压缩感知稀疏重构理论
如果一个信号是稀疏的或可压缩的，那么该信号就可以用远低于奈奎斯特采样定理所要求的采样率来无失真的重构出该信号。如果信号稀疏，并且测量矩阵满足不相干和RIP属性，使用压缩感知恢复的信号稀疏重建可以通过解决以下最优化问题来实现：
其中，α是估计信号，y是测量信号，Θ是测量矩阵，ε是噪声门限。详见文献“阵列三维合成孔径雷达稀疏成像技术研究”韦顺军，2013。
定义9、合成孔径雷达原始回波仿真方法
合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定系统参数条件下具有合成孔径雷达回波信号特性的原始信号的方法，详见文献“张朋，合成孔径雷达回波信号仿真研究，西北工业大学博士论文，2004”。
定义10、线阵SAR的快时刻和慢时刻
线阵SAR运动平台飞过一个方位向合成孔径长度所需要的时间称为慢时间，雷达系统以一定时间长度的重复周期发射接收脉冲，因此慢时间可以表示为一个以脉冲重复周期为步长的离散化时间变量，其中每一个脉冲重复周期离散时间变量值为一个慢时刻。快时刻是指在一个脉冲重复周期内，距离向采样回波信号的时间间隔变量。详见文献“合成孔径雷达成像原理”，皮一鸣等编著，电子科技大学出版社出版。
定义11、传统的迭代软阈值迭代算法((Iterative Soft ThresholdingAlgorithm,ISTA)
迭代软阈值算法(ISTA)是求解范数稀疏性正则化的反问题目标函数，其要求所要求解的原始反问题是信号的稀疏表示，并且测量矩阵满足不相干和RIP属性。用来解决如下优化问题：
其中,x是估计信号，y是测量信号，Θ是测量矩阵，λ为正则化参数，其标准软阈值函数是
传统迭代软阈值算法(ISTA)详见文献“K.Bredies,D.A.Lorenz.Linearconvergence of iterative soft-thresholding,2008”。
定义12、传统理论成像分辨率
线阵SAR成像传统理论分辨率是指利用经典匹配滤波理论成像算法得到线阵SAR系统在距离向、方位向和切航迹向的成像分辨率。对于收发共用天线，线阵SAR距离向的分辨率记为ρr，近似表达式为其中C为电磁波在空气中传播的速度，Br为线阵SAR发射信号带宽；方位向分辨率记为ρa，近似表达为/>其中Da为天线在方位向的真实孔径。详见文献“合成孔径雷达成像原理，杨建宇等编著，电子科技大学出版社”。
定义13、传统的BP算法
BP算法首先将原始回波数据进行距离压缩，然后通过精确补偿成像空间中每个采样点到合成孔径长度内每个天线相位中心的回波时延相位并进行相干累加，从而恢复场景散射点的目标函数。传统的BP算法详见文献“阵列三维合成孔径雷达稀疏成像技术研究，韦顺军，2013”。
定义14、方位向抽样矩阵
取一个M阶单位矩阵IM的第行到第/>行构成一个抽样率为β的方位向抽样矩阵，记为Ωβ，其中round(·)为取整运算符号。
定义15、降序排列
descend(·)表示对·进行从大到小的一个排序。
本发明提出的一种基于多通道SAR系统的高精度干涉相位获取方法，它包括如下步骤：
步骤1、多通道SAR系统构建：
构建一个具有两幅天线的多通道SAR系统，包括雷达发射机、接收机、控制和处理计算机；其中主天线是单发多收模式，即一个发射机多个接收机，而副天线采用单发单收模式，即一个接收机一个发射机，采用控制和处理计算机来控制SAR系统信号的发射、接收、采集功能，并完成后续采集数据的处理及成像；
步骤2、初始化SAR系统参数：
初始化SAR系统参数包括：平台速度矢量记为vr为方位向速度；两天线的基线长度，记为H；雷达信号载频，记为fc；脉冲重复时间记为PRI；雷达系统的脉冲重复频率为PRF；雷达发射信号带宽记为Br；电磁波在空气中的传播速度记为c；方位向慢时刻记为τ，τ＝1,…,α，α为方位向慢时刻总数；距离向快时间记为t，t＝1,…,T，T为距离向快时刻总数；
上述参数均为SAR系统标准参数，其中雷达信号载频fc，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽Br；平台速度矢量在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；
步骤3、划分SAR的成像场景空间：
以雷达波束照射场区域地平面的单位向量所构成的平面直角坐标系作为InSAR的成像场景目标空间；初始化距离向成像场空间长度为Lx，方位向成像场空间长度为Ly；将成像场景目标空间进行划分为大小相等的平面单元格，主天线成像场景平面的距离向、方位向单元格数分别为N，L，副天线成像场景平面的距离向、方位向单元格数分别为I，J；
步骤4、主副天线成像：
在实际InSAR成像中，原始回波数据由数据接收机提供，主天线在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为zm(τ,t)，副天线在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为zs(τ,t)；采用定义13中的传统的BP算法分别对zm(τ,t)和zs(τ,t)进行成像处理，得到成像处理后的图像为ym(mxm,mym)和ys(mxs,mys),其中mxm＝1,…,N，N为步骤3中主天线成像场景平面的距离向单元格数，mym＝1,…,L，L为步骤3中主天线成像场景平面的方位向单元格数，mxs＝1,…,I，I为步骤3中副天线的成像场景平面的距离向单元格数，mys＝1,…,J，J为步骤3中副天线成像场景平面的方位向单元格数；
采用公式和/>计算得到ym(mxm,mym)和ys(mxs,mys)的矩阵形式，记为Ys和Ym，其中(·)T为转置运算符号；
步骤5、建立矩阵重构方程：
步骤5.1、采用公式计算得到Ys的二维离散傅里叶变换，记为其中/>为定义6中定义的二维傅里叶变换，Ys为步骤4中ys(mxs,mys)矩阵形式，FI∈I×I为定义6中定义的距离向离散傅里叶变换，FJ∈J×J为定义6中定义的方位向离散傅里叶变换，I为步骤3中副天线成像场景平面的距离向单元格数，J为步骤3中副天线成像场景平面的方位向单元格数；
步骤5.2、采用公式Θm＝angle(Ym)，计算得到矩阵形式的主天线相位，记为Θm，其中Ym为步骤4中ym(mxm,mym)矩阵形式，angle(·)为相位运算符号；
采用公式计算得到干涉图像相位，记为U，其中/>为步骤5.1中Ys的二维离散傅里叶变换，Ys为步骤4中ys(τ,t)矩阵形式，Ωβ∈J×L为定义14中定义的抽样率为β的方位向抽样矩阵，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换，(·)H为矩阵的共轭转置运算符号，N为步骤3中主天线成像场景平面的距离向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
步骤6、重构干涉相位：
步骤6.1、初始化参数：
初始化基于定义11中传统迭代软阈值算法的重构干涉相位算法的参数初始化：初始化误差迭代终止门限为ε，初始化残差η(0)，初始化范数项系数为p，初始化最大迭代次数为Tmax，初始化场景的非稀疏解稀疏解X(0)，图像差值W(0)，场景目标估计数K，背景相对幅值参数ξ，“噪声”矩阵的标准差δ0；
步骤6.2、估计非稀疏解：
采用公式计算得到第t+1次迭代的场景非稀疏解，记为/>其中t＝0,…,Tmax，Tmax为步骤6.1中的最大迭代次数，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换，Ωβ∈J×L为定义14中定义的抽样率为β的方位向抽样矩阵，(·)H为步骤5.2中矩阵的共轭转置变换运算符号，Θm∈N×L为步骤5.2中矩阵形式的主天线相位，X(t)为第t次迭代的场景稀疏解，N为步骤3中主天线成像场景平面的距离向单元格数，J为步骤3中副天线成像场景平面的方位向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
步骤6.3、估计“噪声”标准差：
采用公式计算得到第t+1次迭代的“噪声”矩阵的标准差，记为δt+1，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，K为步骤6.1中的场景目标估计数，descend(·)为定义15中定义的降序排列，|·|为取模运算符号；
步骤6.4、计算图像差值：
采用公式计算得到第t+1次迭代的图像差值，记为W(t +1)，其中/>为步骤5.2中的副天线SAR图像的二维离散傅里叶变换，Ωβ∈J×L为步骤5.2中的方位向抽样矩阵，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换，/>为定义7中定义的哈达玛积运算，(·)T为矩阵的转置变换运算符号,N为步骤3中主天线成像场景平面的距离向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
步骤6.5、估计稀疏解：
采用公式计算得到第t+1次迭代的场景稀疏解，记为X(t+1)，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，δt+1为步骤6.3中的第t+1次迭代的“噪声”矩阵的标准差，ξ为步骤6.1中的背景相对幅值参数，soft(·；·)为定义11中定义的ISTA算法中的标准软阈值函数；
步骤6.6、迭代终止判断：
采用公式计算得到第t+1次迭代的残差，记为η(t+1)，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，X(t)为步骤6.5中的第t次迭代的场景稀疏解，||·||p为定义3中定义的范数，p为步骤6.1中的范数项系数；
如果t<Tmax且η(t+1)>ε，继续执行步骤6.2～步骤6.4，且t＝t+1，X(t)＝X(t+1)，W(t)＝W(t+1),其中，ε为步骤6.1中的误差迭代终止门限，Tmax为步骤6.1中的最大迭代次数；
若不满足t<Tmax或η(t+1)>ε，算法终止迭代，则输出非稀疏解和稀疏解
步骤7、获取时域干涉相位：
采用公式计算得到非稀疏解的时域干涉图，记为Xnon，其中/>步骤6.6中的非稀疏解，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换,N为步骤3中主天线成像场景平面的距离向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
采用公式计算得到稀疏解的时域干涉图，记为X，其中/>步骤6.6中的稀疏解；
至此，我们全场景InSAR的干涉相位结果，整个重构方法结束。
本发明的创新点是：为获取高分辨率图像，主副天线均为多通道的InSAR系统带来了硬件成本大和实现难度的问题，本发明通过减少副天线通道为单通道，得到两幅分辨率不同的SAR图像，再由两幅不同分辨率图像获得高精度干涉相位，由于方程求解存在克罗内克积计算，使得内存和计算量消耗巨大，针对该问题，利用二维频域干涉相位的高稀疏特性，建立了从低分辨图像到二维频域干涉相位的矩阵重构方程，为求解该方程提出了一种矩阵形式稀疏重构算法，最终获得了观测区域的干涉相位。
本发明优点在于利用了多通道SAR系统，在极大降低了系统复杂度和硬件成本的优势上，得到了高精度的干涉相位，并且针对矢量方程内存和计算量消耗大的缺点，建立了从两幅不同分辨率图像到频域干涉相位的矩阵重构方程，该算法极大地提高算法的运算效率。
附图说明
图1为本发明流程图；
图2数据传输框图；
图3SAR系统模型图
其中，(a)几何示意图，(b)天线相位中心图，其中A1表示多通道主天线，A2表示单通道副天线；图4为系统参数表；
具体实施方式
本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-2017b上验证正确。具体实施步骤如下：
步骤1、多通道SAR系统构建：
构建一个具有两幅天线的多通道SAR系统，包括雷达发射机、接收机、控制和处理计算机；其中主天线是单发多收模式，即一个发射机多个接收机，而副天线采用单发单收模式，即一个接收机一个发射机，采用控制和处理计算机来控制SAR系统信号的发射、接收、采集等功能，并完成后续采集数据的处理及成像；
步骤2、初始化SAR系统参数：
初始化SAR系统参数包括：平台速度矢量记为vr＝50m/s为方位向速度；两天线的基线长度，记为H＝1m；雷达发射信号载频，记为fc＝35GHz；脉冲重复时间记为PRI＝5ms；雷达系统的脉冲重复频率为PRF＝200Hz；雷达发射信号带宽记为Br＝0.4GHz；电磁波在空气中的传播速度记为c＝299792458m/s；方位向慢时刻记为τ，τ＝1,…,α，α＝1503为方位向慢时刻总数；距离向快时间记为t，t＝1,…,T，T＝800为距离向快时刻总数
上述参数均为SAR系统标准参数，其中雷达信号载频fc，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽Br；平台速度矢量在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；
步骤3、划分SAR的成像场景空间：
以雷达波束照射场区域地平面的单位向量所构成的平面直角坐标系作为InSAR的成像场景目标空间；初始化距离向成像场空间长度为Lx＝300，方位向成像场空间长度为Ly＝300；将成像场景目标空间进行划分为大小相等的平面单元格，主天线成像场景平面的距离向、方位向单元格数分别为N＝1601，L＝801，副天线成像场景平面的距离向、方位向单元格数分别为I＝1601，J＝401；
步骤4、主副天线成像：
在实际InSAR成像中，原始回波数据由数据接收机提供，主天线在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为zm(τ,t)，副天线在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为zs(τ,t)；采用定义13中的传统的BP算法分别对zm(τ,t)和zs(τ,t)进行成像处理，得到成像处理后的图像为ym(mxm,mym)和ys(mxs,mys),其中mxm＝1,…,N，N＝1601为步骤3中主天线成像场景平面的距离向单元格数，mym＝1,…,L，L＝801为步骤3中主天线成像场景平面的方位向单元格数，mxs＝1,…,I，I＝1601为步骤3中副天线的成像场景平面的距离向单元格数，mys＝1,…,J，J＝401为步骤3中副天线成像场景平面的方位向单元格数；
采用公式和/>计算得到ym(mxm,mym)和ys(mxs,mys)的矩阵形式，记为Ys和Ym，其中(·)T为转置运算符号；
步骤5、建立矩阵重构方程：
步骤5.1、采用公式计算得到Ys的二维离散傅里叶变换，记为其中/>为定义6中定义的二维傅里叶变换，Ys为步骤4中ys(mxs,mys)矩阵形式，FI∈I×I为定义6中定义的距离向离散傅里叶变换，FJ∈J×J为定义6中定义的方位向离散傅里叶变换，I＝1601为步骤3中副天线成像场景平面的距离向单元格数，J＝401为步骤3中副天线成像场景平面的方位向单元格数；
步骤5.2、采用公式Θm＝angle(Ym)，计算得到矩阵形式的主天线相位，记为Θm，其中Ym为步骤4中ym(mxm,mym)矩阵形式，angle(·)为相位运算符号；
采用公式计算得到干涉图像相位，记为U，其中/>为步骤5.1中Ys的二维离散傅里叶变换，Ys为步骤4中ys(τ,t)矩阵形式，Ωβ∈J×L为定义14中定义的抽样率为β的方位向抽样矩阵，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换，(·)H为矩阵的共轭转置运算符号，N＝1601为步骤3中主天线成像场景平面的距离向单元格数，L＝801为步骤3中主天线成像场景平面的方位向单元格数；
步骤6、重构干涉相位：
步骤6.1、初始化参数：
初始化基于定义11中传统迭代软阈值算法(ISTA)的重构干涉相位算法的参数初始化：初始化误差迭代终止门限为ε＝1×10-3，初始化残差η(0)＝1×1014，初始化范数项系数为p＝2，初始化最大迭代次数为Tmax＝50000，初始化场景的非稀疏解稀疏解X(0)＝0，图像差值W(0)＝0，场景目标估计数K＝127×2500，背景相对幅值参数ξ＝1，“噪声”矩阵的标准差δ0＝0；
步骤6.2、估计非稀疏解：
采用公式计算得到第t+1次迭代的场景非稀疏解，记为/>其中t＝0,…,Tmax，Tmax＝50000为步骤6.1中的最大迭代次数，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换，Ωβ∈J×L为定义14中定义的抽样率为β＝2的方位向抽样矩阵，(·)H为步骤5.2中矩阵的共轭转置变换运算符号，Θm∈N×L为步骤5.2中矩阵形式的主天线相位，X(t)为第t次迭代的场景稀疏解，N＝1601为步骤3中主天线成像场景平面的距离向单元格数，J＝401为步骤3中副天线成像场景平面的方位向单元格数，L＝801为步骤3中主天线成像场景平面的方位向单元格数；
步骤6.3、估计“噪声”标准差：
采用公式计算得到第t+1次迭代的“噪声”矩阵的标准差，记为δt+1，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，K＝127×2500为步骤6.1中的场景目标估计数，descend(·)为定义15中定义的降序排列，|·|为取模运算符号；
步骤6.4、计算图像差值：
采用公式计算得到第t+1次迭代的图像差值，记为W(t +1)，其中/>为步骤5.3中的副天线SAR图像的二维离散傅里叶变换，Ωβ∈J×L为步骤5.3中的方位向抽样矩阵，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换，/>为定义7中定义的哈达玛积运算，(·)T为矩阵的转置变换运算符号,N＝1601为步骤3中主天线成像场景平面的距离向单元格数，L＝801为步骤3中主天线成像场景平面的方位向单元格数；
步骤6.5、估计稀疏解：
采用公式计算得到第t+1次迭代的场景稀疏解，记为X(t+1)，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，δt+1为步骤6.3中的第t+1次迭代的“噪声”矩阵的标准差，ξ＝1为步骤6.1中的背景相对幅值参数，soft(·；·)为定义11中定义的ISTA算法中的标准软阈值函数；
步骤6.6、迭代终止判断：
采用公式计算得到第t+1次迭代的残差，记为η(t+1)，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，X(t)为步骤6.5中的第t次迭代的场景稀疏解，||·||p为定义3中定义的范数，p＝2为步骤6.1中的范数项系数；
如果t<Tmax且η(t+1)>ε，继续执行步骤6.2～步骤6.4，且t＝t+1，X(t)＝X(t+1)，W(t)＝W(t+1),其中，ε＝1×10-3为步骤6.1中的误差迭代终止门限，Tmax＝50000为步骤6.1中的最大迭代次数；
若不满足t<Tmax或η(t+1)>ε，算法终止迭代，则输出非稀疏解和稀疏解
步骤7、获取时域干涉相位：
采用公式计算得到非稀疏解的时域干涉图，记为Xnon，其中/>步骤6.6中的非稀疏解，FN∈N×N为定义6中定义的距离向离散傅里叶变换，FL∈L×L为定义6中定义的方位向离散傅里叶变换,N＝1601为步骤3中主天线成像场景平面的距离向单元格数，L＝801为步骤3中主天线成像场景平面的方位向单元格数；
采用公式计算得到稀疏解的时域干涉图，记为X，其中/>步骤6.6中的稀疏解；
至此，我们全场景InSAR的干涉相位结果，整个重构方法结束。
经过计算机仿真及实测数据结果证明，本发明首先，使用主天线多通道、副天线单通道技术和InSAR技术相结合实现了宽测绘InSAR，在高分宽幅多通道SAR模型基础上，建立了一种宽测绘InSAR模型。通过该InSAR系统得到了两幅分辨率不同的SAR图像，然后建立了从两幅不同分辨率的图像到频域干涉相位的矩阵重构方程，有效提高了InSAR成像干涉图的精度并降低了运算量，极大地提升了InSAR成像的运算效率，并且极大程度上降低了系统的硬件成本。

Claims (1)

1.一种基于多通道SAR系统的高精度干涉相位获取方法，其特征是它包括以下几个步骤：
步骤1、多通道SAR系统构建：
构建一个具有两幅天线的多通道SAR系统，包括雷达发射机、接收机、控制和处理计算机；其中主天线是单发多收模式，即一个发射机多个接收机，而副天线采用单发单收模式，即一个接收机一个发射机，采用控制和处理计算机来控制SAR系统信号的发射、接收、采集功能，并完成后续采集数据的处理及成像；
步骤2、初始化SAR系统参数：
初始化SAR系统参数包括：平台速度矢量记为vr为方位向速度；两天线的基线长度，记为H；雷达信号载频，记为fc；脉冲重复时间记为PRI；雷达系统的脉冲重复频率为PRF；雷达发射信号带宽记为Br；电磁波在空气中的传播速度记为c；方位向慢时刻记为τ，τ＝1,…,α，α为方位向慢时刻总数；距离向快时间记为t，t＝1,…,T，T为距离向快时刻总数；
上述参数均为SAR系统标准参数，其中雷达信号载频fc，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽Br；平台速度矢量在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；
步骤3、划分SAR的成像场景空间：
以雷达波束照射场区域地平面的单位向量所构成的平面直角坐标系作为InSAR的成像场景目标空间；初始化距离向成像场空间长度为Lx，方位向成像场空间长度为Ly；将成像场景目标空间进行划分为大小相等的平面单元格，主天线成像场景平面的距离向、方位向单元格数分别为N，L，副天线成像场景平面的距离向、方位向单元格数分别为I，J；
步骤4、主副天线成像：
在实际InSAR成像中，原始回波数据由数据接收机提供，主天线在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为zm(τ,t)，副天线在第τ个方位向慢时刻和第t个距离向快时刻中InSAR的原始回波数据记为zs(τ,t)；采用传统的BP算法分别对zm(τ,t)和zs(τ,t)进行成像处理，得到成像处理后的图像为ym(mxm,mym)和ys(mxs,mys),其中mxm＝1,…,N，N为步骤3中主天线成像场景平面的距离向单元格数，mym＝1,…,L，L为步骤3中主天线成像场景平面的方位向单元格数，mxs＝1,…,I，I为步骤3中副天线的成像场景平面的距离向单元格数，mys＝1,…,J，J为步骤3中副天线成像场景平面的方位向单元格数；
采用公式和/>计算得到ym(mxm,mym)和ys(mxs,mys)的矩阵形式，记为Ys和Ym，其中(·)T为转置运算符号；
步骤5、建立矩阵重构方程：
步骤5.1、采用公式计算得到Ys的二维离散傅里叶变换，记为/>其中/>为二维傅里叶变换，Ys为步骤4中ys(mxs,mys)矩阵形式，FI∈I×I为距离向离散傅里叶变换，FJ∈J×J为方位向离散傅里叶变换，I为步骤3中副天线成像场景平面的距离向单元格数，J为步骤3中副天线成像场景平面的方位向单元格数；
步骤5.2、采用公式Θm＝angle(Ym)，计算得到矩阵形式的主天线相位，记为Θm，其中Ym为步骤4中ym(mxm,mym)矩阵形式，angle(·)为相位运算符号；
采用公式计算得到干涉图像相位，记为U，其中/>为步骤5.1中Ys的二维离散傅里叶变换，Ys为步骤4中ys(τ,t)矩阵形式，Ωβ∈J×L为抽样率为β的方位向抽样矩阵，FN∈N×N为距离向离散傅里叶变换，FL∈L×L为方位向离散傅里叶变换，(·)H为矩阵的共轭转置运算符号，N为步骤3中主天线成像场景平面的距离向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
步骤6、重构干涉相位：
步骤6.1、初始化参数：
初始化基于传统迭代软阈值算法的重构干涉相位算法的参数初始化：初始化误差迭代终止门限为ε，初始化残差η(0)，初始化范数项系数为p，初始化最大迭代次数为Tmax，初始化场景的非稀疏解稀疏解X(0)，图像差值W(0)，场景目标估计数K，背景相对幅值参数ξ，“噪声”矩阵的标准差δ0；
步骤6.2、估计非稀疏解：
采用公式计算得到第t+1次迭代的场景非稀疏解，记为其中t＝0,…,Tmax，Tmax为步骤6.1中的最大迭代次数，FN∈N×N为距离向离散傅里叶变换，FL∈L×L为方位向离散傅里叶变换，Ωβ∈J×L为抽样率为β的方位向抽样矩阵，(·)H为步骤5.2中矩阵的共轭转置变换运算符号，Θm∈N×L为步骤5.2中矩阵形式的主天线相位，X(t)为第t次迭代的场景稀疏解，N为步骤3中主天线成像场景平面的距离向单元格数，J为步骤3中副天线成像场景平面的方位向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
步骤6.3、估计“噪声”标准差：
采用公式计算得到第t+1次迭代的“噪声”矩阵的标准差，记为δt+1，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，K为步骤6.1中的场景目标估计数，descend(·)为降序排列，|·|为取模运算符号；
步骤6.4、计算图像差值：
采用公式计算得到第t+1次迭代的图像差值，记为W(t+1)，其中/>为步骤5.2中的副天线SAR图像的二维离散傅里叶变换，Ωβ∈J×L为步骤5.2中的方位向抽样矩阵，FN∈N×N为距离向离散傅里叶变换，FL∈L×L为方位向离散傅里叶变换，/>为哈达玛积运算，(·)T为矩阵的转置变换运算符号,N为步骤3中主天线成像场景平面的距离向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
步骤6.5、估计稀疏解：
采用公式计算得到第t+1次迭代的场景稀疏解，记为X(t+1)，其中为步骤6.2中的第t+1次迭代的场景非稀疏解，δt+1为步骤6.3中的第t+1次迭代的“噪声”矩阵的标准差，ξ为步骤6.1中的背景相对幅值参数，soft(·；·)为ISTA算法中的标准软阈值函数；
步骤6.6、迭代终止判断：
采用公式计算得到第t+1次迭代的残差，记为η(t+1)，其中/>为步骤6.2中的第t+1次迭代的场景非稀疏解，X(t)为步骤6.5中的第t次迭代的场景稀疏解，||·||p为范数，p为步骤6.1中的范数项系数；
如果t<Tmax且η(t+1)>ε，继续执行步骤6.2～步骤6.4，且t＝t+1，X(t)＝X(t+1)，W(t)＝W(t+1),其中，ε为步骤6.1中的误差迭代终止门限，Tmax为步骤6.1中的最大迭代次数；
若不满足t<Tmax或η(t+1)>ε，算法终止迭代，则输出非稀疏解和稀疏解
步骤7、获取时域干涉相位：
采用公式计算得到非稀疏解的时域干涉图，记为Xnon，其中/>步骤6.6中的非稀疏解，FN∈N×N为距离向离散傅里叶变换，FL∈L×L为方位向离散傅里叶变换,N为步骤3中主天线成像场景平面的距离向单元格数，L为步骤3中主天线成像场景平面的方位向单元格数；
采用公式计算得到稀疏解的时域干涉图，记为X，其中/>步骤6.6中的稀疏解；
至此，我们全场景InSAR的干涉相位结果，整个重构方法结束。
CN202111009890.9A 2021-08-31 2021-08-31 一种基于多通道sar系统的高精度干涉相位获取方法 Active CN113835090B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111009890.9A CN113835090B (zh) 2021-08-31 2021-08-31 一种基于多通道sar系统的高精度干涉相位获取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111009890.9A CN113835090B (zh) 2021-08-31 2021-08-31 一种基于多通道sar系统的高精度干涉相位获取方法

Publications (2)

Publication Number Publication Date
CN113835090A CN113835090A (zh) 2021-12-24
CN113835090B true CN113835090B (zh) 2024-04-12

Family

ID=78961623

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111009890.9A Active CN113835090B (zh) 2021-08-31 2021-08-31 一种基于多通道sar系统的高精度干涉相位获取方法

Country Status (1)

Country Link
CN (1) CN113835090B (zh)

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103439693A (zh) * 2013-08-16 2013-12-11 电子科技大学 一种线阵sar稀疏重构成像与相位误差校正方法
CN103616682A (zh) * 2013-09-27 2014-03-05 电子科技大学 一种基于曲面投影的多基线InSAR处理方法
CN103698763A (zh) * 2013-12-12 2014-04-02 电子科技大学 基于硬阈值omp的线阵sar稀疏成像方法
CN103698764A (zh) * 2013-12-27 2014-04-02 中国科学院电子学研究所 一种稀疏采样条件下的干涉合成孔径雷达成像方法
CN103869316A (zh) * 2014-03-27 2014-06-18 西安电子科技大学 基于稀疏表征的前视阵列sar超分辨成像方法
CN107957574A (zh) * 2017-12-28 2018-04-24 桂林电子科技大学 基于ifft和混合匹配追踪的时分地基mimo滑坡雷达成像方法
CN108226927A (zh) * 2017-12-14 2018-06-29 电子科技大学 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法
CN109061642A (zh) * 2018-07-13 2018-12-21 电子科技大学 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法
CN110133656A (zh) * 2019-06-06 2019-08-16 电子科技大学 一种基于互质阵列的分解与融合的三维sar稀疏成像方法
CN111145337A (zh) * 2019-12-13 2020-05-12 电子科技大学 基于分辨率逼近的快速稀疏重构的线阵sar三维成像方法
CN111679277A (zh) * 2020-05-28 2020-09-18 电子科技大学 一种基于sbrim算法的多基线层析sar三维成像方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103439693A (zh) * 2013-08-16 2013-12-11 电子科技大学 一种线阵sar稀疏重构成像与相位误差校正方法
CN103616682A (zh) * 2013-09-27 2014-03-05 电子科技大学 一种基于曲面投影的多基线InSAR处理方法
CN103698763A (zh) * 2013-12-12 2014-04-02 电子科技大学 基于硬阈值omp的线阵sar稀疏成像方法
CN103698764A (zh) * 2013-12-27 2014-04-02 中国科学院电子学研究所 一种稀疏采样条件下的干涉合成孔径雷达成像方法
CN103869316A (zh) * 2014-03-27 2014-06-18 西安电子科技大学 基于稀疏表征的前视阵列sar超分辨成像方法
CN108226927A (zh) * 2017-12-14 2018-06-29 电子科技大学 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法
CN107957574A (zh) * 2017-12-28 2018-04-24 桂林电子科技大学 基于ifft和混合匹配追踪的时分地基mimo滑坡雷达成像方法
CN109061642A (zh) * 2018-07-13 2018-12-21 电子科技大学 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法
CN110133656A (zh) * 2019-06-06 2019-08-16 电子科技大学 一种基于互质阵列的分解与融合的三维sar稀疏成像方法
CN111145337A (zh) * 2019-12-13 2020-05-12 电子科技大学 基于分辨率逼近的快速稀疏重构的线阵sar三维成像方法
CN111679277A (zh) * 2020-05-28 2020-09-18 电子科技大学 一种基于sbrim算法的多基线层析sar三维成像方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"Compressive sensing for MIMO radar";Y. Yao 等;《2009 IEEE International Conference on Acoustics, Speech and Signal Processing》;第3017-3020页 *
"Parameter Selection in Sparsity-Driven SAR Imaging";O. Batu 等;《IEEE Transactions on Aerospace and Electronic Systems》;第47卷(第4期);第3040-3050页 *

Also Published As

Publication number Publication date
CN113835090A (zh) 2021-12-24

Similar Documents

Publication Publication Date Title
CN109061642B (zh) 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法
CN104111458B (zh) 基于双重稀疏约束的压缩感知合成孔径雷达成像方法
CN107037429B (zh) 基于门限梯度追踪算法的线阵sar三维成像方法
CN111679277B (zh) 一种基于sbrim算法的多基线层析sar三维成像方法
CN108226927B (zh) 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法
CN112099008B (zh) 基于cv-admmn的sa-isar成像与自聚焦方法
CN110244303B (zh) 基于sbl-admm的稀疏孔径isar成像方法
CN103399315B (zh) 实孔径相控阵雷达高分辨探测成像方法
CN112099007B (zh) 适用于非理想天线方向图的方位向多通道sar模糊抑制方法
CN105842693B (zh) 一种基于压缩感知的双通道sar动目标检测的方法
CN104950305A (zh) 一种基于稀疏约束的实波束扫描雷达角超分辨成像方法
CN107493106B (zh) 一种基于压缩感知的频率和角度联合估计的方法
CN111505639A (zh) 基于变重频采样模式的合成孔径雷达宽幅稀疏成像方法
CN107576961A (zh) 一种互质降采样间歇合成孔径雷达稀疏成像方法
CN111722227B (zh) 基于近似观测矩阵的聚束sar压缩感知成像方法
CN109031299B (zh) 低信噪比条件下基于相位差分的isar平动补偿方法
CN113608218B (zh) 一种基于后向投影原理的频域干涉相位稀疏重构方法
CN107656271B (zh) 基于压缩感知重构的太赫兹雷达成像算法
Liu et al. Ambiguities Suppression for Azimuth Multichannel SAR Based on ${L_ {2, q}} $ Regularization With Application to Gaofen-3 Ultra-Fine Stripmap Mode
CN113484859B (zh) 一种基于融合技术的二维超分辨雷达成像方法
CN107561534B (zh) 一种基于全极化高轨sar的电离层时变tec测量方法
CN113835090B (zh) 一种基于多通道sar系统的高精度干涉相位获取方法
CN108931770B (zh) 基于多维贝塔过程线性回归的isar成像方法
CN110133656B (zh) 一种基于互质阵列的分解与融合的三维sar稀疏成像方法
CN115453523A (zh) 一种扫描雷达稀疏目标批处理超分辨方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant