CN113640797B - 一种用于参考条带模式InSAR的前斜视测高方法 - Google Patents

一种用于参考条带模式InSAR的前斜视测高方法 Download PDF

Info

Publication number
CN113640797B
CN113640797B CN202110909166.5A CN202110909166A CN113640797B CN 113640797 B CN113640797 B CN 113640797B CN 202110909166 A CN202110909166 A CN 202110909166A CN 113640797 B CN113640797 B CN 113640797B
Authority
CN
China
Prior art keywords
phase
antenna
grid
echo
calculating
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
CN202110909166.5A
Other languages
English (en)
Other versions
CN113640797A (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.)
Beihang University
Beijing Institute of Technology BIT
Original Assignee
Beihang University
Beijing Institute of Technology BIT
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 Beihang University, Beijing Institute of Technology BIT filed Critical Beihang University
Priority to CN202110909166.5A priority Critical patent/CN113640797B/zh
Publication of CN113640797A publication Critical patent/CN113640797A/zh
Application granted granted Critical
Publication of CN113640797B publication Critical patent/CN113640797B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/904—SAR modes
    • G01S13/9041—Squint 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

本发明公开了一种用于参考条带模式InSAR的前斜视测高方法，包括：建立三维直角坐标系，读取参考条带模式InSAR双通道回波数据、雷达系统参数及运动参数、地面控制点参数；在三维直角坐标系中的地平面对回波数据进行成像处理；将成像处理得到的两幅SAR图像进行干涉处理，并对干涉相位进行滤波；对干涉相位进行相位解缠绕，并基于地面控制点确定绝对相位；计算各像素点高程；对采样网格点上的数据进行插值重采样，得到几何校正后地面数字高程模型。本发明避免了传统InSAR测高处理方法在处理斜视数据时对主辅图像配准精度要求极高、难以保证干涉相位精度、在斜视模式下需要求解多普勒方程、参考条带模式InSAR的数字高程模型获取的问题。

Description

一种用于参考条带模式InSAR的前斜视测高方法
技术领域
本发明涉及干涉合成孔径雷达(InSAR)测高技术领域，具体涉及一种用于参考条带模式InSAR的前斜视测高方法。
背景技术
干涉合成孔径雷达(Interference Synthetic Aperture Radar,InSAR)测高是指对同一区域的两幅SAR图像进行干涉处理，从干涉相位中提取地面高程信息。目前，InSAR技术已广泛应用于机载、星载SAR系统，其应用领域涵盖地标形变探测、动目标检测、海洋测绘、森林制图、洪涝监测、交通检测和冰川研究等。用于测高的传统InSAR系统主要工作在正侧视模式下，用于干涉的两幅SAR图像中各像素点具有零多普勒中心频率或较小的多普勒中心频率；在此情况下，亚像素级的图像配准精度足以保证干涉相位精度和图像相干性。当InSAR系统工作在大斜视模式下时，受到多普勒中心频率影响，亚像素级的配准精度不再能够保证干涉相位精度。
参考条带模式SAR是指通过连续地沿俯仰向调整天线波束指向，以获取不与SAR平台运动轨迹平行的观测条带区域的一种新型SAR成像模式。该模式与传统的条带模式相比，可以更加灵活的设置观测区域，一定程度上摆脱平台运动轨迹对观测区域的限制。为减小数据量，便于设计时变的脉冲重复间隔以应对大距离徙动的影响，参考条带模式SAR需要保持波束指向垂直于条带走向，因而一般工作在大斜视模式下。
参考条带模式InSAR是指将干涉测高技术引入参考条带模式SAR中，通过单航过或重复航过获取同一区域的两幅SAR图像以完成干涉测高处理；为使参考条带模式InSAR具有干涉测高的能力，需要研究并设计一种可以用于大斜视模式的InSAR测高方法。
发明内容
为解决参考条带模式InSAR在大斜视模式下的测高问题，本发明提供一种用于参考条带模式InSAR的前斜视测高方法。
本发明公开了一种用于参考条带模式InSAR的前斜视测高方法，包括：
步骤1、建立三维直角坐标系，读取参考条带模式InSAR双通道回波数据、雷达系统参数及运动参数、地面控制点参数；
步骤2、基于所述雷达系统参数及运动参数，采用改进的BP算法在所述三维直角坐标系中的地平面对回波数据进行成像处理；
步骤3、将成像处理得到的两幅SAR图像进行干涉处理，并使用非线性相位滤波方法对干涉相位进行滤波；
步骤4、采用最小费用流方法对干涉相位进行相位解缠绕，并基于地面控制点确定绝对相位；
步骤5、采用前斜视模式SAR的高程反演公式计算各像素点高程；
步骤6、采用散点插值方法，对采样网格点上的数据进行插值重采样，得到几何校正后地面数字高程模型。
作为本发明的进一步改进，所述步骤1，具体包括：
步骤101、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系；其中，所述三维直角坐标系由满足右手定则的三个坐标轴X、Y、Z组成，X轴指向飞行方向，Y轴垂直于地平面指向上方，Z轴由右手定则确定；
步骤102、读取参考条带模式InSAR双通道回波数据、雷达系统参数及运动参数、地面控制点坐标；
所述参考条带模式InSAR双通道回波数据分别表示为EM和ES，分别对应InSAR系统主天线自发自收的回波信号矩阵和主天线发射辅天线接收的回波信号矩阵；EM和ES的都为Na×Nr的复矩阵，其中Na为方位向回波采样条数，Nr为距离向单条回波采样的点数；
所述雷达系统参数及运动参数包括雷达的工作波长λ、雷达的采样率fs、发射的线性调频信号的基带数字波形s(n)、光速c、各方位采样时刻的回波窗开启距离Rmin(na)、主天线相位中心在回波录取过程中的高度H及其空间位置坐标
Figure BDA0003203061240000021
辅天线相位中心在回波录取过程中的空间位置坐标
Figure BDA0003203061240000022
从主天线相位中心指向辅天线相位中心的基线矢量
Figure BDA0003203061240000023
雷达波束矢量在场景系下的水平方位角θy、雷达波束矢量在场景系下的下视角θp(na)、雷达天线方向图方位向3dB波束宽度θa3dB、雷达天线方向图俯仰向3dB波束宽度θr3dB、期望成像区域的零高度地面网格点坐标集合G(x,z)＝{(xj,0,zj)T}，其中n∈{1，2，…，Llfm}为基带数字波形的索引，数字波形长度设为Llfm，na∈{1，2，…，Na}为方位采样索引；
地面控制点坐标表示为
Figure BDA0003203061240000031
作为本发明的进一步改进，所述步骤2，具体包括：
步骤201、脉冲压缩；
通过在矩阵右端补零将回波信号矩阵EM扩充为Na×(Nr+Llfm)的矩阵EMp，通过基带数字波形s(n)的右端补零将其扩充为1×(Nr+Llfm)的向量sp(i)，其中i∈{1，2，…，Nr+Llfm}；对sp(i)进行快速傅里叶变换并取共轭得到匹配滤波器hlfm；依据如下公式逐行对EMp进行处理得到脉压后Na×(Nr+Llfm)的回波矩阵EMcp；
EMcp(i，：)＝IFFT[FFT[EMp(i，：)]*hlfm] (1)
其中FFT[·]与IFFT[·]分别表示快速傅里叶变换和快速傅里叶反变换，*表示矩阵的Hadamard积，(i，：)表示矩阵的第i行；
步骤202、初始化成像网格；
构造零初值的图像网格I(x，z)＝{p(xj，zj)}与地面网格点坐标集合G(x，z)一一对应；
步骤203、对第i行回波进行逆投影处理；
根据如下公式计算发射天线相位中心
Figure BDA0003203061240000032
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000033
Figure BDA0003203061240000034
根据如下公式计算天线系下发射天线相位中心
Figure BDA0003203061240000035
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000036
Figure BDA0003203061240000037
Figure BDA0003203061240000038
Figure BDA0003203061240000039
根据如下公式计算
Figure BDA00032030612400000310
在天线系下的方位角θa(i，j)和俯仰角θr(i，j)；
Figure BDA0003203061240000041
Figure BDA0003203061240000042
其中
Figure BDA0003203061240000043
表示
Figure BDA0003203061240000044
的第k个元素；
不同时满足以下条件的地面网格点的投影值v(i，j)＝0，对其余网格点进行后续步骤的处理；
Figure BDA0003203061240000045
Figure BDA0003203061240000046
若对主天线回波进行成像，则根据如下公式计算接收天线相位中心
Figure BDA0003203061240000047
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000048
Figure BDA0003203061240000049
若对辅天线回波进行成像，则根据如下公式计算接收天线相位中心
Figure BDA00032030612400000410
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA00032030612400000411
Figure BDA00032030612400000412
根据如下的公式计算地面网格点在回波中的索引idx(i，j)；
Figure BDA00032030612400000413
其中||·||表示向量的二范数；
不满足0≤idx(i，j)＜Nr的网格点的投影值v(i，j)＝0，对其余网格点进行后续步骤的处理；
根据如下公式计算第i个脉冲在成像网格点(xj，0，zj)T的投影值v(i，j)；
Figure BDA00032030612400000414
其中
Figure BDA00032030612400000415
表示向下取整，m为插值核点数，取正偶数；
步骤204、对全部方位采样时刻的回波进行步骤203；
步骤205、求和得到单视复图像；
根据如下公式求和计算图像网格I(x，z)＝{p(xj，zj)}；
p(xj，zj)＝∑iv(i，j) (14)
若对主天线成像则将成像结果表示为Im(x，z)＝{pm(xj，zj)}，若对辅天线成像，则将成像结果表示为Is(x，z)＝{ps(xj，zj)}；
步骤206、对辅天线回波矩阵EM重复步骤201到步骤204的操作。
作为本发明的进一步改进，所述步骤3，具体包括：
步骤301、干涉处理；
根据如下公式得到原始干涉相位；
φ(x，z)＝∠(Im(x，z)*conj(Is(x，z))) (15)
其中∠(·)表示取相位操作，conj(·)表示取共轭操作；
步骤302、非线性相位滤波；
根据王青松在IEEE Geoscience and Remote Sensing Letters期刊2011年11月第8卷第6期中的论文“An Efficient and Adaptive Approach for Noise Filtering ofSAR Interferometric Phase Images”中提出的非线性相位滤波方法完成相位滤波操作；
记滤波后相位为φfilt(x，z)。
作为本发明的进一步改进，所述步骤4，具体包括：
步骤401、最小费用流解相位缠绕；
根据Maori Costantini在Transactions on Geoscience and Remote Sensing期刊1998年5月第36卷第3期中的论文“A Novel Phase Unwrapping Method Based onNetwork Programming”中提出的最小费用最大流相位解缠绕方法完成相位解缠绕操作；
记解缠绕后相位为φunwrap(x，z)；
步骤402、计算地面控制点在主图像中聚焦位置；
根据如下公式计算地面控制点在主图像中的聚焦位置(xctrl，zctrl)；
xctrl＝pctrl_x (16)
Figure BDA0003203061240000051
其中，q∈{0，1}，当雷达右视时q＝0，当雷达左视时q＝1；
步骤403、计算地面控制点的绝对干涉相位；
根据如下公式计算控制点的合成孔径中心时刻主天线位置
Figure BDA0003203061240000052
Figure BDA0003203061240000053
Figure BDA0003203061240000054
根据如下公式计算控制点绝对干涉相位；
Figure BDA0003203061240000061
步骤404、计算解缠绕后的绝对干涉相位；
记
Figure BDA0003203061240000062
中离(xctrl，zctrl)最近的点为
Figure BDA0003203061240000063
根据如下公式计算解缠绕后的绝对干涉相位
Figure BDA0003203061240000064
Figure BDA0003203061240000065
Figure BDA0003203061240000066
其中round(·)表示四舍五入。
作为本发明的进一步改进，所述步骤5，具体包括：
步骤501、计算各地面像素网格对应的孔径中心时刻主天线相位中心位置；
根据如下公式计算各地面像素网格G(x，z)对应的孔径中心时刻主天线相位中心位置
Figure BDA0003203061240000067
Figure BDA0003203061240000068
步骤502、计算各像素网格高程值；
根据如下公式计算各像素网格坐标(xj，0，zj)T对应的高程值h(xj，zj)
Figure BDA0003203061240000069
作为本发明的进一步改进，所述步骤6，具体包括：
步骤601、计算各像素网格内散射点的空间坐标；
根据如下公式计算各像素网格点的空间坐标(x3j，y3j，z3j)T；
x3j＝xj (22)
y3j＝h(xj，zj) (23)
Figure BDA00032030612400000610
其中，q∈{0，1}，当雷达右视时q＝0，当雷达左视时q＝1；
步骤602、散点插值得到几何校正后成像区域数字高程模型；
根据Amidror在Journal of Electronic Imaging期刊2002年4月第11卷第2期中的论文“Scattered Data Interpolation Methods for Electronic Imaging Systems：aSurvey”中提出的方法，在上一步获取的空间点云{(x3j，y3j，z3j)T}上对成像网格坐标集合G(x，z)＝{(xj，0，zj)T}中空间坐标的第二维进行散点插值；得到期望成像区域的几何校正后的数字高程模型DEM＝{(xj,yj,zj)T}。
与现有技术相比，本发明的有益效果为：
本发明通过对两幅SAR图像进行干涉处理、相位滤波、相位解缠、高程反演、几何校正，得到目标区域的地面数字高程模型，实现大斜视条件下的InSAR测高处理；其避免了传统InSAR测高处理方法在处理斜视数据时对主辅图像配准精度要求极高、难以保证干涉相位精度的问题，避免了传统InSAR测高处理方法在斜视模式下需要求解多普勒方程的问题，能够较为快速、准确的解决参考条带模式InSAR的数字高程模型获取问题。
附图说明
图1为本发明一种实施例公开的用于参考条带模式InSAR的前斜视测高方法的流程图；
图2为本发明一种实施例公开的雷达发射信号数字波形图；
图3为本发明一种实施例公开的雷达回波窗开启距离变化图；
图4为本发明一种实施例公开的雷达波束指向下视角变化图；
图5为本发明一种实施例公开的雷达主天线BP成像结果幅度图；
图6为本发明一种实施例公开的雷达辅天线BP成像结果幅度图；
图7为本发明一种实施例公开的原始干涉相位图；
图8为本发明一种实施例公开的滤波后干涉相位图；
图9为本发明一种实施例公开的解缠绕后绝对干涉相位图；
图10为本发明一种实施例公开的几何校正前高程图；
图11为本发明一种实施例公开的数字高程模型图；
图12为本发明一种实施例公开的数字高程模型三维展示图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图对本发明做进一步的详细描述：
如图1所示，本发明提供一种用于参考条带模式InSAR的前斜视测高方法，包括：
步骤1、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系，读取参考条带模式InSAR的双通道回波数据，读取雷达系统参数及运动轨迹，读取地面控制点参数；
具体包括：
步骤101、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系；
三维直角坐标系由满足右手定则的三个坐标轴X、Y、Z组成，其中X轴指向飞行方向，Y轴垂直于地平面指向上方，Z轴由右手定则确定，涉及到三维空间坐标的数据全部定义在该坐标系下；
步骤102、读取参考条带模式InSAR双通道回波数据、读取雷达系统参数及运动参数、读取地面控制点三维坐标；
数据、参数、坐标在本实施例中具体为：参考条带模式InSAR双通道回波数据分别表示为EM和ES，分别对应InSAR系统主天线自发自收的回波信号矩阵和主天线发射辅天线接收的回波信号矩阵。EM和ES的都为Na×Nr的复矩阵，其中Na＝15334为方位向回波采样条数，Nr＝3201为距离向单条回波采样的点数。雷达系统参数包括雷达的工作波长λ＝0.02m、雷达的采样率fs＝160MHz、发射的线性调频信号的基带数字波形s(n)如图2所示、光速c＝3e8m/s、各方位采样时刻的回波窗开启距离Rmin(na)如图3所示、主天线相位中心在回波录取过程中的高度H＝10km及其空间位置坐标
Figure BDA0003203061240000081
辅天线相位中心在回波录取过程中的空间位置坐标
Figure BDA0003203061240000082
从主天线相位中心指向辅天线相位中心的基线矢量
Figure BDA0003203061240000083
雷达波束矢量在场景系下的水平方位角θy＝45deg、雷达波束矢量在场景系下的下视角θp(na)如图4所示、雷达天线方向图方位向3dB波束宽度θa3dB＝1deg、雷达天线方向图俯仰向3dB波束宽度θr3dB＝2deg、期望成像区域的零高度地面网格点坐标集合G(x，z)＝{(xj，0，zj)T}为以坐标(10606.601718，0，10606.601718)T为中心的1000×1000的0.1m间隔网格，其中n∈{1，2，…，Llfm}为基带数字波形的索引，数字波形长度设为Llfm＝1601，na∈{1，2，…，Na}为方位采样索引；地面控制点坐标表示为
Figure BDA0003203061240000091
步骤2、基于雷达系统参数和运动轨迹，采用改进的BP算法在所建三维直角坐标系中的地平面完成成像处理；
具体包括：
步骤201、脉冲压缩；
通过在矩阵右端补零将回波信号矩阵EM扩充为Na×(Nr+Llfm)的矩阵EMp，通过基带数字波形s(n)的右端补零将其扩充为1×(Nr+Llfm)的向量sp(i)，其中i∈{1，2，…，Nr+Llfm}；对sp(i)进行快速傅里叶变换并取共轭得到匹配滤波器hlfm；依据如下公式逐行对EMp进行处理得到脉压后Na×(Nr+Llfm)的回波矩阵EMcp；
EMcp(i，：)＝IFFT[FFT[EMp(i，：)]*hlfm] (25)
其中FFT[·]与IFFT[·]分别表示快速傅里叶变换和快速傅里叶反变换，*表示矩阵的Hadamard积，(i，：)表示矩阵的第i行。
步骤202、初始化成像网格；
构造零初值的图像网格I(x，z)＝{p(xj，zj)}与地面网格点坐标集合G(x，z)一一对应；
步骤203、对第i行回波进行逆投影处理；
根据如下公式计算发射天线相位中心
Figure BDA0003203061240000092
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000093
Figure BDA0003203061240000094
根据如下公式计算天线系下发射天线相位中心
Figure BDA0003203061240000095
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000096
Figure BDA0003203061240000097
Figure BDA0003203061240000098
Figure BDA0003203061240000099
根据如下公式计算
Figure BDA00032030612400000910
在天线系下的方位角θa(i，j)和俯仰角θr(i，j)；
Figure BDA0003203061240000101
Figure BDA0003203061240000102
其中
Figure BDA0003203061240000103
表示
Figure BDA0003203061240000104
的第k个元素；
不同时满足以下条件的地面网格点的投影值v(i，j)＝0，对其余网格点进行后续步骤的处理；
|θa(i，j)|≤0.5deg (32)
|θr(i，j)|≤1deg (33)
若对主天线回波进行成像，则根据如下公式计算接收天线相位中心
Figure BDA0003203061240000105
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000106
Figure BDA0003203061240000107
若对辅天线回波进行成像，则根据如下公式计算接收天线相位中心
Figure BDA0003203061240000108
到成像网格点(xj，0，zj)T的矢量表示
Figure BDA0003203061240000109
Figure BDA00032030612400001010
根据如下的公式计算地面网格点在回波中的索引idx(i，j)；
Figure BDA00032030612400001011
其中||·||表示向量的二范数；
不满足0≤idx(i，j)＜Nr的网格点的投影值v(i，j)＝0，对其余网格点进行后续步骤的处理；
根据如下公式计算第i个脉冲在成像网格点(xj，0，zj)T的投影值v(i，j)；
Figure BDA00032030612400001012
其中
Figure BDA00032030612400001013
表示向下取整，m＝16为插值核点数；
步骤204、对全部方位采样时刻的回波进行步骤203；
步骤205、求和得到单视复图像；
根据如下公式求和计算图像网格I(x，z)＝{p(xj，zj)}；
p(xj，zj)＝∑iv(i，j) (38)
若对主天线成像则将成像结果表示为Im(x，z)＝{pm(xj，zj)}，若对辅天线成像，则将成像结果表示为Is(x，z)＝{ps(xj，zj)}；
步骤206、对辅天线回波矩阵EM重复步骤201到步骤204的操作；
经过步骤201到步骤206处理的主辅天线复图像幅值如图5、图6所示。
步骤3、将成像处理得到的两幅SAR图像进行干涉处理，并使用非线性相位滤波方法对干涉相位进行滤波；
具体包括：
步骤301、干涉处理；
根据如下公式得到原始干涉相位如图7所示；
φ(x，z)＝∠(Im(x，z)*conj(Is(x，z))) (39)
其中∠(·)表示取相位操作，conj(·)表示取共轭操作；
步骤302、非线性相位滤波
根据王青松在IEEE Geoscience and Remote Sensing Letters期刊2011年11月第8卷第6期中的论文“An Efficient and Adaptive Approach for Noise Filtering ofSAR Interferometric Phase Images”中提出的非线性相位滤波方法完成相位滤波操作；
记滤波后相位为φfilt(x，z)，如图8所示；
步骤4、采用最小费用流方法对干涉相位进行相位解缠绕，并基于地面控制点确定绝对相位；
具体包括：
步骤401、最小费用流解相位缠绕；
根据Maori Costantini在Transactions on Geoscience and Remote Sensing期刊1998年5月第36卷第3期中的论文“A Novel Phase Unwrapping Method Based onNetwork Programming”中提出的最小费用最大流相位解缠绕方法完成相位解缠绕操作；
记解缠绕后相位为φunwrap(x，z)，如图9所示；
步骤402、计算地面控制点在主图像中聚焦位置；
根据如下公式计算地面控制点在主图像中的聚焦位置(xctrl，zctrl)；
xctrl＝pctrl_x＝10606.5518m (40)
Figure BDA0003203061240000111
其中，q∈{0，1}，当雷达右视时q＝0，当雷达左视时q＝1；
步骤403、计算地面控制点的绝对干涉相位；
根据如下公式计算控制点的合成孔径中心时刻主天线位置
Figure BDA0003203061240000121
Figure BDA0003203061240000122
Figure BDA0003203061240000123
根据如下公式计算控制点绝对干涉相位；
Figure BDA0003203061240000124
步骤404、计算解缠绕后的绝对干涉相位；
记
Figure BDA0003203061240000125
中离(xctrl，zctrl)最近的点为
Figure BDA0003203061240000126
根据如下公式计算解缠绕后的绝对干涉相位
Figure BDA0003203061240000127
Figure BDA0003203061240000128
Figure BDA0003203061240000129
其中round(·)表示四舍五入；
步骤5、采用前斜视模式SAR的高程反演公式计算各像素点高程；
具体包括：
步骤501、计算各地面像素网格对应的孔径中心时刻主天线相位中心位置；
根据如下公式计算各地面像素网格G(x，z)对应的孔径中心时刻主天线相位中心位置
Figure BDA00032030612400001210
Figure BDA00032030612400001211
步骤502、计算各像素网格高程值；
根据如下公式计算各像素网格坐标(xj，0，zj)T对应的高程值h(xj，zj)如图10所示；
Figure BDA00032030612400001212
步骤6、采用散点插值方法，对采样网格点上的数据进行插值重采样，得到几何校正后地面数字高程模型；
具体包括：
步骤601、计算各像素网格内散射点的空间坐标；
根据如下公式计算各像素网格点的空间坐标(x3j，y3j，z3j)T；
x3j＝xj (48)
y3j＝h(xj，zj) (49)
Figure BDA0003203061240000131
其中，q∈{0，1}，当雷达右视时q＝0，当雷达左视时q＝1；
步骤602、散点插值得到几何校正后成像区域数字高程模型；
根据Amidror在Journal of Electronic Imaging期刊2002年4月第11卷第2期中的论文“Scattered Data Interpolation Methods for Electronic Imaging Systems：aSurvey”中提出的方法，在上一步获取的空间点云{(x3j，y3j，z3j)T}上对成像网格坐标集合G(x，z)＝{(xj，0，zj)T}中空间坐标的第二维进行散点插值；得到期望成像区域的几何校正后的数字高程模型DEM＝{(xj，yj，zj)T}，如图11与图12所示。
本发明的优点为：
本发明通过对两幅SAR图像进行干涉处理、相位滤波、相位解缠、高程反演、几何校正，得到目标区域的地面数字高程模型，实现大斜视条件下的InSAR测高处理；其避免了传统InSAR测高处理方法在处理斜视数据时对主辅图像配准精度要求极高、难以保证干涉相位精度的问题，避免了传统InSAR测高处理方法在斜视模式下需要求解多普勒方程的问题，能够较为快速、准确的解决参考条带模式InSAR的数字高程模型获取问题。
以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种用于参考条带模式InSAR的前斜视测高方法，其特征在于，包括：
步骤1、建立三维直角坐标系，读取参考条带模式InSAR双通道回波数据、雷达系统参数及运动参数、地面控制点参数；具体包括：
步骤101、沿平行和垂直于飞行速度的方向及天向建立三维直角坐标系；其中，所述三维直角坐标系由满足右手定则的三个坐标轴X、Y、Z组成，X轴指向飞行方向，Y轴垂直于地平面指向上方，Z轴由右手定则确定；
步骤102、读取参考条带模式InSAR双通道回波数据、雷达系统参数及运动参数、地面控制点坐标；
所述参考条带模式InSAR双通道回波数据分别表示为EM和ES，分别对应InSAR系统主天线自发自收的回波信号矩阵和主天线发射辅天线接收的回波信号矩阵；EM和ES的都为Na×Nr的复矩阵，其中Na为方位向回波采样条数，Nr为距离向单条回波采样的点数；
所述雷达系统参数及运动参数包括雷达的工作波长λ、雷达的采样率fs、发射的线性调频信号的基带数字波形s(n)、光速c、各方位采样时刻的回波窗开启距离Rmin(na)、主天线相位中心在回波录取过程中的高度H及其空间位置坐标
Figure FDA0003515627870000011
辅天线相位中心在回波录取过程中的空间位置坐标
Figure FDA0003515627870000012
从主天线相位中心指向辅天线相位中心的基线矢量
Figure FDA0003515627870000013
雷达波束矢量在场景系下的水平方位角θy、雷达波束矢量在场景系下的下视角θp(na)、雷达天线方向图方位向3dB波束宽度θa3dB、雷达天线方向图俯仰向3dB波束宽度θr3dB、期望成像区域的零高度地面网格点坐标集合G(x,z)＝{(xj,0,zj)T}，其中n∈{1,2,…,Llfm}为基带数字波形的索引，数字波形长度设为Llfm，na∈{1,2,…,Na}为方位采样索引；
地面控制点坐标表示为
Figure FDA0003515627870000014
步骤2、基于所述雷达系统参数及运动参数，在所述三维直角坐标系中的地平面对回波数据进行成像处理；具体包括：
步骤201、脉冲压缩；
通过在矩阵右端补零将回波信号矩阵EM扩充为Na×(Nr+Llfm)的矩阵EMp，通过基带数字波形s(n)的右端补零将其扩充为1×(Nr+Llfm)的向量sp(i)，其中i∈{1,2,…,Nr+Llfm}；对sp(i)进行快速傅里叶变换并取共轭得到匹配滤波器hlfm；依据如下公式逐行对EMp进行处理得到脉压后Na×(Nr+Llfm)的回波矩阵EMcp；
EMcp(i,:)＝IFFT[FFT[EMp(i,:)]*hlfm] (1)
其中FFT[·]与IFFT[·]分别表示快速傅里叶变换和快速傅里叶反变换，*表示矩阵的Hadamard积，(i,:)表示矩阵的第i行；
步骤202、初始化成像网格；
构造零初值的图像网格I(x,z)＝{p(xj,zj)}与地面网格点坐标集合G(x,z)一一对应；
步骤203、对第i行回波进行逆投影处理；
根据如下公式计算发射天线相位中心
Figure FDA0003515627870000021
到成像网格点(xj,0,zj)T的矢量表示
Figure FDA0003515627870000022
Figure FDA0003515627870000023
根据如下公式计算天线系下发射天线相位中心
Figure FDA0003515627870000024
到成像网格点(xj,0,zj)T的矢量表示
Figure FDA0003515627870000025
Figure FDA0003515627870000026
Figure FDA0003515627870000027
Figure FDA0003515627870000028
根据如下公式计算
Figure FDA0003515627870000029
在天线系下的方位角θa(i,j)和俯仰角θr(i,j)；
Figure FDA00035156278700000210
Figure FDA00035156278700000211
其中
Figure FDA00035156278700000212
表示
Figure FDA00035156278700000213
的第k个元素；
不同时满足以下条件的地面网格点的投影值v(i,j)＝0，对其余网格点进行后续步骤的处理；
Figure FDA00035156278700000214
Figure FDA00035156278700000215
若对主天线回波进行成像，则根据如下公式计算接收天线相位中心
Figure FDA00035156278700000216
到成像网格点(xj,0,zj)T的矢量表示
Figure FDA0003515627870000031
Figure FDA0003515627870000032
若对辅天线回波进行成像，则根据如下公式计算接收天线相位中心
Figure FDA0003515627870000033
到成像网格点(xj,0,zj)T的矢量表示
Figure FDA0003515627870000034
Figure FDA0003515627870000035
根据如下的公式计算地面网格点在回波中的索引idx(i,j)；
Figure FDA0003515627870000036
其中‖·‖表示向量的二范数；
不满足0≤idx(i,j)<Nr的网格点的投影值v(i,j)＝0，对其余网格点进行后续步骤的处理；
根据如下公式计算第i个脉冲在成像网格点(xj,0,zj)T的投影值v(i,j)；
Figure FDA0003515627870000037
其中
Figure FDA0003515627870000038
表示向下取整，m为插值核点数，取正偶数；
步骤204、对全部方位采样时刻的回波进行步骤203；
步骤205、求和得到单视复图像；
根据如下公式求和计算图像网格I(x,z)＝{p(xj,zj)}；
p(xj,zj)＝∑iv(i,j) (14)
步骤206、将主天线成像则将成像结果表示为Im(x,z)＝{pm(xj,zj)}；
步骤207、对辅天线回波矩阵EM重复步骤201到步骤205的操作，将辅天线成像结果表示为Is(x,z)＝{ps(xj,zj)}；
步骤3、将成像处理得到的两幅SAR图像进行干涉处理，并使用非线性相位滤波方法对干涉相位进行滤波；
步骤4、采用最小费用流方法对干涉相位进行相位解缠绕，并基于地面控制点确定绝对相位；
步骤5、采用前斜视模式SAR的高程反演公式计算各像素点高程；
步骤6、采用散点插值方法，对采样网格点上的数据进行插值重采样，得到几何校正后地面数字高程模型。
2.如权利要求1所述的前斜视测高方法，其特征在于，所述步骤3，具体包括：
步骤301、干涉处理；
根据如下公式得到原始干涉相位；
φ(x,z)＝∠(Im(x,z)*conj(Is(x,z))) (15)
其中∠(·)表示取相位操作，conj(·)表示取共轭操作；
步骤302、非线性相位滤波，记滤波后相位为φfilt(x,z)。
3.如权利要求2所述的前斜视测高方法，其特征在于，所述步骤4，具体包括：
步骤401、最小费用流解相位缠绕，记解缠绕后相位为φunwrap(x,z)；
步骤402、计算地面控制点在主图像中聚焦位置；
根据如下公式计算地面控制点在主图像中的聚焦位置(xctrl,zctrl)；
xctrl＝pctrl_x (16)
Figure FDA0003515627870000041
其中，q∈{0,1}，当雷达右视时q＝0，当雷达左视时q＝1；
步骤403、计算地面控制点的绝对干涉相位；
根据如下公式计算控制点的合成孔径中心时刻主天线位置
Figure FDA0003515627870000049
Figure FDA00035156278700000410
Figure FDA0003515627870000042
根据如下公式计算控制点绝对干涉相位；
Figure FDA0003515627870000043
步骤404、计算解缠绕后的绝对干涉相位；
记
Figure FDA0003515627870000044
中离(xctrl,zctrl)最近的点为
Figure FDA0003515627870000045
根据如下公式计算解缠绕后的绝对干涉相位
Figure FDA0003515627870000046
Figure FDA0003515627870000047
Figure FDA0003515627870000048
其中round(·)表示四舍五入。
4.如权利要求3所述的前斜视测高方法，其特征在于，所述步骤5，具体包括：
步骤501、计算各地面像素网格对应的孔径中心时刻主天线相位中心位置；
根据如下公式计算各地面像素网格G(x,z)对应的孔径中心时刻主天线相位中心位置
Figure FDA0003515627870000051
Figure FDA0003515627870000052
步骤502、计算各像素网格高程值；
根据如下公式计算各像素网格坐标(xj,0,zj)T对应的高程值h(xj,zj)
Figure FDA0003515627870000053
5.如权利要求4所述的前斜视测高方法，其特征在于，所述步骤6，具体包括：
步骤601、计算各像素网格内散射点的空间坐标；
根据如下公式计算各像素网格点的空间坐标(x3j,y3j,z3j)T；
x3j＝xj (22)
y3j＝h(xj,zj) (23)
Figure FDA0003515627870000054
其中，q∈{0,1}，当雷达右视时q＝0，当雷达左视时q＝1；
步骤602、散点插值得到几何校正后成像区域数字高程模型；
在获取的空间点云{(x3j,y3j,z3j)T}上对成像网格坐标集合G(x,z)＝{(xj,0,zj)T}中空间坐标的第二维进行散点插值；得到期望成像区域的几何校正后的数字高程模型DEM＝{(xj,yj,zj)T}。
CN202110909166.5A 2021-08-09 2021-08-09 一种用于参考条带模式InSAR的前斜视测高方法 Active CN113640797B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110909166.5A CN113640797B (zh) 2021-08-09 2021-08-09 一种用于参考条带模式InSAR的前斜视测高方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110909166.5A CN113640797B (zh) 2021-08-09 2021-08-09 一种用于参考条带模式InSAR的前斜视测高方法

Publications (2)

Publication Number Publication Date
CN113640797A CN113640797A (zh) 2021-11-12
CN113640797B true CN113640797B (zh) 2022-04-12

Family

ID=78420335

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110909166.5A Active CN113640797B (zh) 2021-08-09 2021-08-09 一种用于参考条带模式InSAR的前斜视测高方法

Country Status (1)

Country Link
CN (1) CN113640797B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114609633B (zh) * 2022-03-17 2023-09-01 电子科技大学 一种圆周聚束模式干涉sar测高方法

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102508245A (zh) * 2011-11-18 2012-06-20 北京航空航天大学 一种星载多频率与多基线InSAR高程估计精度等效性确定方法
CN102788979A (zh) * 2012-07-20 2012-11-21 电子科技大学 一种基于后向投影InSAR成像配准的GPU实现方法
CN102854506A (zh) * 2012-09-10 2013-01-02 电子科技大学 一种基于后向投影算法的动基线干涉sar基线补偿方法
CN103018741A (zh) * 2012-12-11 2013-04-03 电子科技大学 一种基于后向投影的InSAR成像去平地一体化方法
CN103487809A (zh) * 2013-09-23 2014-01-01 中国科学院电子学研究所 一种基于BP算法和时变基线的机载InSAR数据处理方法
CN103616682A (zh) * 2013-09-27 2014-03-05 电子科技大学 一种基于曲面投影的多基线InSAR处理方法
WO2014074631A1 (en) * 2012-11-07 2014-05-15 Neva Ridge Technologies Sar point cloud generation system
CN104698457A (zh) * 2014-09-02 2015-06-10 电子科技大学 一种迭代曲面预测InSAR成像及高度估计方法
CN105182337A (zh) * 2015-09-09 2015-12-23 北京航空航天大学 一种基于曲面后向投影算法的形变反演方法
WO2016086699A1 (zh) * 2014-12-01 2016-06-09 中国科学院电子学研究所 一种结合局部频率估计的小波域InSAR干涉相位滤波方法
CN105929399A (zh) * 2016-04-25 2016-09-07 电子科技大学 一种干涉sar数据成像及高程估计方法
CN106707281A (zh) * 2017-01-05 2017-05-24 北京航空航天大学 一种基于多频数据处理的机载D‑InSAR形变检测方法
CN112882030A (zh) * 2021-01-12 2021-06-01 中国科学院空天信息创新研究院 InSAR成像干涉一体化处理方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8154435B2 (en) * 2008-08-22 2012-04-10 Microsoft Corporation Stability monitoring using synthetic aperture radar
TWI486556B (zh) * 2013-01-04 2015-06-01 Univ Nat Central Integration of Radar and Optical Satellite Image for Three - dimensional Location

Patent Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102508245A (zh) * 2011-11-18 2012-06-20 北京航空航天大学 一种星载多频率与多基线InSAR高程估计精度等效性确定方法
CN102788979A (zh) * 2012-07-20 2012-11-21 电子科技大学 一种基于后向投影InSAR成像配准的GPU实现方法
CN102854506A (zh) * 2012-09-10 2013-01-02 电子科技大学 一种基于后向投影算法的动基线干涉sar基线补偿方法
WO2014074631A1 (en) * 2012-11-07 2014-05-15 Neva Ridge Technologies Sar point cloud generation system
CN103018741A (zh) * 2012-12-11 2013-04-03 电子科技大学 一种基于后向投影的InSAR成像去平地一体化方法
CN103487809A (zh) * 2013-09-23 2014-01-01 中国科学院电子学研究所 一种基于BP算法和时变基线的机载InSAR数据处理方法
CN103616682A (zh) * 2013-09-27 2014-03-05 电子科技大学 一种基于曲面投影的多基线InSAR处理方法
CN104698457A (zh) * 2014-09-02 2015-06-10 电子科技大学 一种迭代曲面预测InSAR成像及高度估计方法
WO2016086699A1 (zh) * 2014-12-01 2016-06-09 中国科学院电子学研究所 一种结合局部频率估计的小波域InSAR干涉相位滤波方法
EP3229038A1 (en) * 2014-12-01 2017-10-11 Institute of Electronics, Chinese Academy of Sciences Wavelet domain insar interferometric phase filtering method in combination with local frequency estimation
CN105182337A (zh) * 2015-09-09 2015-12-23 北京航空航天大学 一种基于曲面后向投影算法的形变反演方法
CN105929399A (zh) * 2016-04-25 2016-09-07 电子科技大学 一种干涉sar数据成像及高程估计方法
CN106707281A (zh) * 2017-01-05 2017-05-24 北京航空航天大学 一种基于多频数据处理的机载D‑InSAR形变检测方法
CN112882030A (zh) * 2021-01-12 2021-06-01 中国科学院空天信息创新研究院 InSAR成像干涉一体化处理方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"InSAR DEM reconstruction based on backprojection algorithm in two converse flights";Xiaoning Hu等;《2019 6th Asia-pacific conference on synthetic aperture radar(APSAR)》;20191231;第1-4页 *
"斜视情况下的分布式卫星insar处理方法";郭交;《系统工程与电子技术》;20110630;第33卷(第6期);第1243-1246页 *

Also Published As

Publication number Publication date
CN113640797A (zh) 2021-11-12

Similar Documents

Publication Publication Date Title
CN109031219B (zh) 基于相位测距的宽带雷达弹道目标微动几何参数估计方法
CN107229048B (zh) 一种高分宽幅sar动目标速度估计与成像方法
CN109856635B (zh) 一种csar地面动目标重聚焦成像方法
CN105929399B (zh) 一种干涉sar数据成像及高程估计方法
CN103487809B (zh) 一种基于BP算法和时变基线的机载InSAR数据处理方法
CN105137430B (zh) 一种前视阵列sar的回波稀疏获取及其三维成像方法
CN104597447B (zh) 一种子孔径SAR大斜视改进Omega‑K成像方法
CN105223572B (zh) 一种基于pfa算法的正前视双基sar成像处理方法
CN110673143A (zh) 一种子孔径大斜视sar俯冲成像的两步处理方法
CN113589285B (zh) 一种飞行器sar实时成像方法
CN114545411B (zh) 一种基于工程实现的极坐标格式多模高分辨sar成像方法
CN103018739A (zh) 一种校正多通道幅相误差的微波三维成像方法
CN105607055A (zh) 一种基于天线方向图的机载雷达单脉冲前视成像方法
CN113640797B (zh) 一种用于参考条带模式InSAR的前斜视测高方法
CN107589421B (zh) 一种阵列前视sar成像方法
CN110879391B (zh) 基于电磁仿真和弹载回波仿真的雷达图像数据集制作方法
CN109633639A (zh) Topsar干涉数据的高精度快速配准方法
CN110554377B (zh) 基于多普勒中心偏移的单通道sar二维流场反演方法及系统
CN115685200A (zh) 一种高精度大前斜视sar成像运动补偿与几何校正方法
CN108132466A (zh) 一种机载阵列天线下视三维成像方法和系统
CN112558070B (zh) 圆周扫描地基sar的频域成像方法及装置
CN108107427A (zh) 基于超分辨技术的机载/弹载阵列雷达前视成像方法
CN109799502B (zh) 一种适用于滤波反投影算法的两维自聚焦方法
CN107765240A (zh) 一种运动状态的判断方法、装置和电子设备
CN105572648B (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