CN106707281B - 一种基于多频数据处理的机载D-InSAR形变检测方法 - Google Patents
一种基于多频数据处理的机载D-InSAR形变检测方法 Download PDFInfo
- Publication number
- CN106707281B CN106707281B CN201710006318.4A CN201710006318A CN106707281B CN 106707281 B CN106707281 B CN 106707281B CN 201710006318 A CN201710006318 A CN 201710006318A CN 106707281 B CN106707281 B CN 106707281B
- Authority
- CN
- China
- Prior art keywords
- matrix
- slc
- processing
- phase
- data
- 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
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
-
- 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/885—Radar or analogous systems specially adapted for specific applications for ground probing
-
- 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/9004—SAR image acquisition techniques
- G01S13/9017—SAR image acquisition techniques with time domain processing of the SAR signals in azimuth
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)
- Signal Processing (AREA)
- Image Processing (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提供一种基于多频数据处理的机载D‑InSAR形变检测方法,通过低频回波数据干涉处理反演得到高程DEM信息,将该DEM信息作为高频回波BP成像处理时的参考网格,得到相应的SAR图像,对SAR图像进行干涉处理分别得到包含地形信息的干涉相位和包含地形、地形形变信息的干涉相位;将两个干涉相位通过差分处理、反演处理得到相应的地形形变量。本发明提供的方法可以控制***近似计算引入的误差,提高***检测精度;能克服传统D‑InSAR***检测精度与数据处理复杂度之间的矛盾,有效获取高频数据的信息,降低相位解缠难度,提高***检测精度与稳定度;省略了去平地操作,简化了干涉数据处理流程。
Description
技术领域
本发明涉及雷达技术领域,尤其涉及一种基于机载平台,利用低频回波数据反演地形DEM,将其作为高频回波数据的BP成像参考平面得到相应的SAR图像,通过差分干涉处理获取高精度地形形变的检测方法。
背景技术
机载D-InSAR***相较于星载D-InSAR技术具有轨道灵活、重访周期短、检测精度高等优势,但同时机载平台存在飞行控制能力低、航迹误差大等问题,因此如何提高***检测精度一直是它的研究热点与难点。传统单频D-InSAR***通过SAR图像配准、干涉处理、去平地操作、滤波处理、相位解缠处理、相位差分处理等步骤得到地形的形变信息。
对于星载D-InSAR***可近似主辅天线到目标的瞬时斜距夹角为0,因此在形变检测公式推导时可近似主辅天线的斜距差为基线在LOS方向的水平分量。同时由于星载平台轨道固定,且精度较高,可以采用基于轨道参数的去平地方法,保证平地相位计算精度的同时减小干涉条纹密度,降低相位解缠难度。但是机载平台的飞行高度较低,此时主辅天线的瞬时斜距夹角不可近似为0,如果仍将斜距差近似为基线的水平分量,会引入较大的***检测误差。同时机载平台的航迹误差较大,这给去平地操作的具体方法选择带来了一定困难。
此外,D-InSAR***还存在检测精度与数据处理难度相矛盾的问题,突出表现为相位解缠难度较大,因此解决上述问题是提高机载D-InSAR***检测精度与稳定度的关键。
发明内容
针对上述问题中存在的不足之处,本发明提供一种基于多频数据处理的机载D-InSAR形变检测方法。
为实现上述目的,本发明提供一种基于多频数据处理的机载D-InSAR形变检测方法,包括:
步骤1、低频回波数据反演高程DEM信息:将两个低频回波数据L1、L2经BP成像处理得到两幅SAR图像SLC1、SLC2,SLC1与SLC2通过图像配准、共轭相乘、滤波处理、相位解缠处理、DEM重建与地理编码操作得到目标地形的DEM信息;
步骤2、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2干涉处理:将所述DEM信息作为H1与H2的BP成像参考网格,通过相位补偿与相干叠加得到两幅SAR图像SLC3、SLC4;SLC3与SLC4通过图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含地形信息的解缠相位;
步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理:将所述DEM信息作为H3的BP成像参考网格,通过相位补偿与相干叠加得到一幅SAR图像SLC5;SLC5与SLC3通过图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含地形信息、地形形变信息的解缠相位;
步骤4、两个解缠相位反演地形形变信息:将步骤2、步骤3中得到的两个解缠相位通过差分处理得到包含地形形变信息的相位,再通过反演处理得到相应的地形形变量。
作为本发明的进一步改进,所述步骤1包括:
步骤11、BP成像:设置一个1×Nr矩阵的频域滤波器F1,其中Nr与回波矩阵距离向长度相同,利用距离向调频率Kr、快时间τ、瞬时斜距R,如公式(1)设置滤波器:
调用函数fft(L1,[],2)得到回波矩阵的距离向频域信号,调用repmat(F1,Na,1)得到Na×Nr的滤波矩阵,频域相乘后得到矩阵M1,对M1频域补零得到矩阵M2;BP成像参考网格距离向间距设置为m,取距离向分辨率的1/3,方位向参考网格间距设置为n,取方位向分辨率的1/3,其中c为光速,Br为调频信号带宽,La为方位向天线长度,网格中心为天线波束中心;以每一个网格为方位向中心,左右各取1/2天线合成孔径长度作为雷达成像的方位向区域,计算网格与每个雷达位置的瞬时斜距R(η),在矩阵M2中找到相应的位置点,补偿相位因子ej4πR(η)/λ,相干叠加处理后得到SAR图像矩阵SLC1和SLC2;
步骤12、图像配准:图像矩阵SLC1和SLC2粗配准,在SLC1中心位置取A×B的匹配窗,在SLC2中心位置取C×D搜索窗,其中A<C,B<D;在搜索窗内逐行逐列以A×B对应窗口遍历计算其与匹配窗的相干系数,相干系数最大时的偏移量为SLC1、SLC2两图像的相对偏移量,通过矩阵切割处理获取公共部分矩阵Re1、Re2;然后进行亚像素处理:
a、对矩阵Re1、Re2进行频域补0处理,插值间隔为0.01个像素,在信号频域将原始信号以中心点分割为四等份,依次位列新插值矩阵的四角边缘区域,其余值填充为0,得到插值后的矩阵In1、In2;
b、主图像In1匹配窗设置为32×32,辅图像In2搜索窗设置为64×64,通过矩阵遍历计算偏移量,此时主图像的匹配窗数组为100个,偏移量拟合;
c、辅图像Re2经插值处理、校正偏移量、采样处理得到配准后的辅图像矩阵Re3;
步骤13、共轭相乘:矩阵Re1与Re3的共轭矩阵相乘,得到干涉相位矩阵Minf;
步骤14、滤波处理:设置一个1×31的窗口在矩阵Minf中遍历移动,每个信号点的值用邻域各点值的中值代替,得到滤波后的矩阵Mfilter;
步骤15、相位解缠处理:调用函数unwrap(Mfilter,[],2),得到解缠相位矩阵Munwrap;
步骤16、DEM重建与地理编码:采用公式(2)反演地形高程:
其中H为平台高度,R为多普勒斜距信息,B为基线长度,α为基线角,利用轨道参数,基线信息以及***控制点的瞬时斜距信息,将斜距坐标系的高程矩阵Z转换到地理坐标系,得到目标地形真实的DEM信息。
作为本发明的进一步改进,所述步骤2包括:
步骤21、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2基于DEM信息进行BP成像:设置距离向频域滤波器,然后对滤波处理后的矩阵进行距离向插值处理,与步骤1中的操作相同;BP成像的参考网格为DEM信息,以每个网格为相位中心两侧取1/2天线合成孔径长度为雷达成像的方位区域,计算网格与每个雷达方位位置的瞬时斜距,在回波插值矩阵中找到对应的矩阵单元,补偿相位因子并相干叠加处理后得到两个SAR图像矩阵SLC3和SLC4;
步骤22、图像矩阵SLC3和SLC4执行步骤12~步骤15,经图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含检测区域的地形信息的解缠相位矩阵φ12。
作为本发明的进一步改进,所述步骤3包括:
步骤31、包含地形信息、地形形变信息的高频辅天线回波矩阵H3基于DEM信息进行BP成像:具体操作与步骤21相同,处理得到图像矩阵SLC5;
步骤32、图像矩阵SLC3和SLC5执行步骤12~步骤15,经图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含检测区域的地形信息、地形形变信息的解缠相位矩阵φ13。
作为本发明的进一步改进,所述步骤4包括:
步骤41、将矩阵φ12和φ13代入公式(3),得到仅包含地形形变信息的相位矩阵φd,其中为形变基线与垂直基线水平分量的比值;
步骤42、将矩阵φd带入公式(4),得到形变场矩阵ΔD;
与现有技术相比,本发明的有益效果为:
1、本发明提供的方法可以控制***近似计算引入的误差,提高***检测精度;
2、本发明提供的方法能克服传统D-InSAR***检测精度与数据处理复杂度之间的矛盾,有效获取高频数据的信息,降低相位解缠难度,提高***检测精度与稳定度;
3、本发明提供的方法省略了去平地操作,简化了干涉数据处理流程。
附图说明
图1为本发明一种实施例公开的基于多频数据处理的机载D-InSAR形变检测方法的流程图;
图2为本发明利用的将DEM信息作为BP成像处理参考网格的示意图;
图3(a)为检测区域的地形信息;
图3(b)为检测区域的形变场信息;
图4(a)为低频地形信息干涉相位;
图4(b)为低频地形信息滤波后相位;
图4(c)为低频地形信息解缠相位;
图5为低频回波数据反演的DEM信息;
图6(a)为基于多频处理的地形信息干涉相位;
图6(b)为基于多频处理的地形信息滤波后相位;
图6(c)为基于多频处理的地形信息解缠相位;
图7(a)为基于多频处理的形变后地形干涉相位;
图7(b)为基于多频处理的形变后地形滤波相位;
图7(c)为基于多频处理的形变后地形解缠相位;
图8(a)为沿LOS方向的形变检测结果;
图8(b)为形变检测误差场。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
为了解决机载D-InSAR***近似误差大、去平地操作误差大以及相位解缠难度大的问题,本发明根据***近似误差产生的原因以及干涉条纹密度对相位解缠操作的影响,采用高低两种频率的回波数据融合处理的方式提高***的检测精度。首先利用低频天线干涉处理获得地形的粗精度DEM,然后将该DEM信息作为高频天线BP成像处理的参考平面,以控制***的近似误差与高频天线的干涉条纹密度,最后利用得到的高频SAR图像进行差分干涉处理获取地形的形变信息,并在处理流程中省略去平地操作,简化数据处理流程的同时进一步控制数据处理误差,实现了一种能获取高精度地表形变信息的机载D-InSAR检测方法。
下面结合附图对本发明做进一步的详细描述:
如图1所示,本发明提供一种基于多频数据处理的机载D-InSAR形变检测方法,包括:
步骤1、低频回波数据反演DEM信息:将两个低频回波矩阵L1和L2经BP成像处理得到两个SAR图像矩阵SLC1和SLC2,将矩阵SLC1与SLC2通过图像配准、共轭相乘、滤波处理、相位解缠处理得到一个解缠相位矩阵Munwrap,最后通过DEM重建与地理编码操作得到目标地形真实的DEM信息。
步骤11、BP成像:设置一个1×Nr矩阵的频域滤波器F1,其中Nr与回波矩阵距离向长度相同,利用距离向调频率Kr、快时间τ、瞬时斜距R,如公式(1)设置滤波器:
调用函数fft(L1,[],2)得到回波矩阵的距离向频域信号,调用repmat(F1,Na,1)得到Na×Nr的滤波矩阵,频域相乘后得到矩阵M1,对M1频域补零得到矩阵M2;BP成像参考网格距离向间距设置为m,取距离向分辨率的1/3,方位向参考网格间距设置为n,取方位向分辨率的1/3,其中c为光速,Br为调频信号带宽,La为方位向天线长度,网格中心为天线波束中心;以每一个网格为方位向中心,左右各取1/2天线合成孔径长度作为雷达成像的方位向区域,计算网格与每个雷达位置的瞬时斜距R(η),在矩阵M2中找到相应的位置点,补偿相位因子ej4πR(η)/λ,相干叠加处理后得到SAR图像矩阵SLC1和SLC2;
步骤12、图像配准:图像矩阵SLC1和SLC2粗配准,在SLC1中心位置取A×B的匹配窗,在SLC2中心位置取C×D搜索窗,其中A<C,B<D;在搜索窗内逐行逐列以A×B对应窗口遍历计算其与匹配窗的相干系数,相干系数最大时的偏移量为SLC1、SLC2两图像的相对偏移量,通过矩阵切割处理获取公共部分矩阵Re1、Re2;然后进行亚像素处理,主要有三步:
a、对矩阵Re1、Re2进行频域补0处理,插值间隔为0.01个像素,在信号频域将原始信号以中心点分割为四等份,依次位列新插值矩阵的四角边缘区域,其余值填充为0,得到插值后的矩阵In1、In2;
b、主图像In1匹配窗设置为32×32,辅图像In2搜索窗设置为64×64,通过矩阵遍历计算偏移量,此时主图像的匹配窗数组为100个,偏移量拟合;
c、辅图像Re2经插值处理、校正偏移量、采样处理得到配准后的辅图像矩阵Re3;
步骤13、共轭相乘:矩阵Re1与Re3的共轭矩阵相乘,得到干涉相位矩阵Minf;
步骤14、滤波处理:设置一个1×31的窗口在矩阵Minf中遍历移动,每个信号点的值用邻域各点值的中值代替,得到滤波后的矩阵Mfilter;
步骤15、相位解缠处理:调用函数unwrap(Mfilter,[],2),得到解缠相位矩阵Munwrap;
步骤16、DEM重建与地理编码:采用公式(2)反演地形高程Z:
其中H为平台高度,R为多普勒斜距信息,B为基线长度,α为基线角,利用轨道参数,基线信息以及***控制点的瞬时斜距信息,将斜距坐标系的高程矩阵Z转换到地理坐标系,得到目标地形真实的DEM信息。
步骤2、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2干涉处理:对矩阵H1和H2做BP成像,将DEM信息作为参考网格,通过计算雷达与网格之间的瞬时斜距得到其在回波插值矩阵中的具***置,进而通过相位补偿与相干叠加得到SAR图像矩阵SLC3和SLC4,两矩阵SLC3和SLC4通过图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含地形信息的解缠相位矩阵φ12。
步骤21、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2基于DEM信息进行BP成像:设置距离向频域滤波器,然后对滤波处理后的矩阵进行距离向插值处理,与步骤1中的操作相同;BP成像的参考网格为DEM信息,以每个网格为相位中心两侧取1/2天线合成孔径长度为雷达成像的方位区域,计算网格与每个雷达方位位置的瞬时斜距,在回波插值矩阵中找到对应的矩阵单元,补偿相位因子并相干叠加处理后得到两个SAR图像矩阵SLC3和SLC4,具体实现方法如图2所示;
步骤22、图像矩阵SLC3和SLC4执行步骤12~步骤15,经图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含检测区域的地形信息的解缠相位矩阵φ12。
步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理:对H3进行BP成像,参考网格为DEM信息,计算雷达与网格之间的瞬时斜距得到其在回波插值矩阵中的对应矩阵单元,进行相位补偿与相干叠加处理得到SAR图像矩阵SLC5,矩阵SLC3与SLC5通过图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含地形信息、地形形变信息的解缠相位矩阵φ13。
步骤31、包含地形信息、地形形变信息的高频辅天线回波矩阵H3基于DEM信息进行BP成像:具体操作与步骤21相同,处理得到图像矩阵SLC5;
步骤32、图像矩阵SLC3和SLC5执行步骤12~步骤15,经图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含检测区域的地形信息、地形形变信息的解缠相位矩阵φ13。
步骤4、两个解缠相位反演地形形变信息:将步骤2、步骤3中得到的两个解缠相位通过差分处理得到包含地形形变信息的相位,再通过反演处理得到相应的地形形变量。
步骤41、将矩阵φ12和φ13代入公式(3),得到仅包含地形形变信息的相位矩阵φd,其中为形变基线与垂直基线水平分量的比值;
步骤42、将矩阵φd带入公式(4),得到形变场矩阵ΔD;
实施例1:
针对一个前后坡场景、凹陷的圆锥形变场进行地形形变检测,雷达成像方式为正侧视。其中形变场信息如图3所示,其中(a)为地形信息,(b)为形变场信息。前后坡地形两侧的平地区域距离向长度均为150m,方位向范围1000m,距离向范围1200m,斜坡高度80m,形变锥半径为173m,高线为10mm。***参数设置为:机载平台高度12km,平台速度150m/s,下视角60°,调频信号脉宽500us,调频信号带宽160MHz,方位向天线长度为6m,方位向采样率95,基线B1为10m,基线角α1为0,基线B2为5m,基线角为α20,低频频率4GHz,高频频率15GHz。
步骤1、低频回波数据反演DEM信息。
(a)BP成像处理:低频回波矩阵的距离向频域匹配滤波器为1024×2048矩阵,矩阵行向量为滤波处理后得到矩阵M1,M1进行距离向频域补0处理,插值倍数为16倍,得到矩阵M2,距离向分辨率为0.94m,方位向分辨率为3m,网格间隔m设置为0.3m,n设置为1m,网格范围设置为距离向1200m,方位向1000m,网格中心为天线的波束中心。合成孔径长度为177m,以每个网格为方位中心两侧取88.5m范围内对应的各个雷达成像的方位向位置,计算该网格到每个雷达位置的瞬时斜距,确定网格点在距离向插值信号中具体矩阵单元,补偿相位ej4πR(η)/λ与相干叠加后得到SAR图像矩阵SLC1和SLC2,大小为625×4000;
(b)图像配准:输入矩阵SLC1和SLC2,在SLC1中心位置取50×100的匹配窗,在SLC2中心位置取200×400搜索窗,在搜索窗内逐行逐列以50×100大小的对应窗口遍历计算相干系数,相干系数最大时偏移量为0。亚像素级配准时匹配窗的大小设置为20×20,搜索窗的大小设置为50×50,匹配窗的个数为100个,通过偏移量拟合后,对辅图像进行插值、校正、采样处理,得到配准后的辅图像矩阵Re3,大小为625×4000;
(c)共轭相乘:矩阵SLC1与Re3的共轭矩阵相乘,得到干涉矩阵Minf如图4(a)所示;
(d)滤波处理:设置一个1×31的窗口在矩阵Minf中遍历移动,每个信号点的值用邻域各点值的中值代替,滤波后的矩阵Mfilter如图4(b)所示;
(e)相位解缠处理:调用函数unwrap(Mfilter,[],2),得到解缠相位矩阵Munwrap,如图4(c)所示;
(f)DEM重建与地理编码,利用公式(2)得到高程信息矩阵Z,地理编码后得到DEM信息,如图5所示。
步骤2、高频主天线回波矩阵H1与包含地形信息的高频辅天线回波矩阵H2做干涉处理。
(a)矩阵H1和H2将DEM信息作为BP成像处理的参考网格。设置距离向匹配滤波器,滤波后的矩阵距离向频域补0插值,与步骤1中的操作相同。合成孔径长度为177m,以DEM信息中每个单元为方位中心取两侧88.5m范围内对应的各个雷达位置,计算该矩阵单元到每个雷达位置的瞬时斜距,确定该矩阵单元在距离向插值信号中对应的具***置,通过相位补偿处理和相干叠加后得到SAR图像矩阵SLC3和SLC4。
(b)图像矩阵SLC3和SLC4执行步骤1中的操作(b)~(e),依次得到包含地形信息的干涉相位矩阵如图6(a)、滤波矩阵如图6(b);和包含地形信息的解缠相位矩阵φ12如图6(c)。
步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理。
(a)矩阵H3将DEM信息作为BP成像处理的参考网格。设置距离向匹配滤波器,滤波后的矩阵距离向频域补0处理,与步骤2中(a)的操作相同。合成孔径长度为177m,以DEM信息中每个单元为方位中心取两侧88.5m范围内对应的各个雷达位置,计算该矩阵单元到每个雷达位置的瞬时斜距,确定该矩阵单元在距离向插值信号中对应的具***置,通过相位补偿处理和相干叠加后得到SAR图像矩阵SLC5。
(b)图像矩阵SLC3和SLC5执行步骤1中的操作(b)~(e),依次得到包含地形信息、地形形变信息的干涉相位矩阵如图7(a)、滤波矩阵如图7(b)和包含地形信息、地形形变信息的解缠相位矩阵φ13如图7(c)。
步骤4、矩阵φ12和φ13差分处理得到检测区域的地形形变量。
(a)将矩阵φ12和φ13代入公式(3),得到仅包含地形形变信息的相位矩阵φd;
(b)将矩阵φd带入公式(4),得到形变场矩阵ΔD,如图8(a)所示,根据理论形变场对检测结果进行分析,得到的误差场如图8(b)所示,误差均值为4.866×10-4m,最大误差为0.43mm,均方差为4.65×10-2mm,达到了mm级的检测精度且稳定度较高,检测结果理想。
本发明主要针对机载D-InSAR***近似误差大、去平地操作不理想、相位解缠难度大的缺陷,利用高低两种频率的回波数据融合处理,通过低频回波数据干涉处理反演得到高程DEM信息,将该DEM作为高频回波BP成像处理时的参考网格,得到相应的SAR图像,在差分干涉处理中控制***近似误差,省略去平地操作,同时降低了干涉条纹密度,减小数据数理难度,提高***检测精度与稳定度,实现了一种高精度的基于多频数据融合处理的机载D-InSAR形变检测方法。其具有以下优点:
1、本发明提供的方法可以控制***近似计算引入的误差,提高***检测精度;
2、本发明提供的方法能克服传统D-InSAR***检测精度与数据处理复杂度之间的矛盾,有效获取高频数据的信息,降低相位解缠难度,提高***检测精度与稳定度;
3、本发明提供的方法省略了去平地操作,简化了干涉数据处理流程。
以上仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (5)
1.一种基于多频数据处理的机载D-InSAR形变检测方法,其特征在于,包括:
步骤1、低频回波数据反演高程DEM信息:将两个低频回波数据L1、L2经BP成像处理得到两幅SAR图像SLC1、SLC2,SLC1与SLC2通过图像配准、共轭相乘、滤波处理、相位解缠处理、DEM重建与地理编码操作得到目标地形的DEM信息;
步骤2、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2干涉处理:将所述DEM信息作为H1与H2的BP成像参考网格,通过相位补偿与相干叠加得到两幅SAR图像SLC3、SLC4;SLC3与SLC4通过图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含地形信息的解缠相位;
步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理:将所述DEM信息作为H3的BP成像参考网格,通过相位补偿与相干叠加得到一幅SAR图像SLC5;SLC5与SLC3通过图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含地形信息、地形形变信息的解缠相位;
步骤4、两个解缠相位反演地形形变信息:将步骤2、步骤3中得到的两个解缠相位通过差分处理得到包含地形形变信息的相位,再通过反演处理得到相应的地形形变量。
2.如权利要求1所述的基于多频数据处理的机载D-InSAR形变检测方法,其特征在于,所述步骤1包括:
步骤11、BP成像:设置一个1×Nr矩阵的频域滤波器F1,其中Nr与回波矩阵距离向长度相同,利用距离向调频率Kr、快时间τ、瞬时斜距R,如公式(1)设置滤波器:
调用函数fft(L1,[],2)得到回波矩阵的距离向频域信号,其中[]表示快速傅里叶变换的数据长度取的列向量个数;调用repmat(F1,Na,1)得到Na×Nr的滤波矩阵,其中Na与回波矩阵方位向长度相同;频域相乘后得到矩阵M1,对M1频域补零得到矩阵M2;BP成像参考网格距离向间距设置为m,取距离向分辨率的1/3,方位向参考网格间距设置为n,取方位向分辨率的1/3,其中c为光速,Br为调频信号带宽,La为方位向天线长度,网格中心为天线波束中心;以每一个网格为方位向中心,左右各取1/2天线合成孔径长度作为雷达成像的方位向区域,其中Rc为天线到目标点的距离;计算网格与每个雷达位置的瞬时斜距R(η),在矩阵M2中找到相应的位置点,补偿相位因子ej4πR(η)/λ,相干叠加处理后得到SAR图像矩阵SLC1和SLC2;
步骤12、图像配准:图像矩阵SLC1和SLC2粗配准,在SLC1中心位置取A×B的匹配窗,在SLC2中心位置取C×D搜索窗,其中A<C,B<D;在搜索窗内逐行逐列以A×B对应窗口遍历计算其与匹配窗的相干系数,相干系数最大时的偏移量为SLC1、SLC2两图像的相对偏移量,通过矩阵切割处理获取公共部分矩阵Re1、Re2;然后进行亚像素处理:
a、对矩阵Re1、Re2进行频域补0处理,插值间隔为0.01个像素,在信号频域将原始信号以中心点分割为四等份,依次位列新插值矩阵的四角边缘区域,其余值填充为0,得到插值后的矩阵In1、In2;
b、主图像In1匹配窗设置为32×32,辅图像In2搜索窗设置为64×64,通过矩阵遍历计算偏移量,此时主图像的匹配窗数组为100个,偏移量拟合;
c、辅图像Re2经插值处理、校正偏移量、采样处理得到配准后的辅图像矩阵Re3;
步骤13、共轭相乘:矩阵Re1与Re3的共轭矩阵相乘,得到干涉相位矩阵Minf;
步骤14、滤波处理:设置一个1×31的窗口在矩阵Minf中遍历移动,每个信号点的值用邻域各点值的中值代替,得到滤波后的矩阵Mfilter;
步骤15、相位解缠处理:调用函数unwrap(Mfilter,[],2),得到解缠相位矩阵Munwrap,其中[]表示相位解缠的跃变容差取默认值π;
步骤16、DEM重建与地理编码:采用公式(2)反演地形高程:
其中H为平台高度,R为多普勒斜距信息,B为基线长度,α为基线角,利用轨道参数,基线信息以及***控制点的瞬时斜距信息,将斜距坐标系的高程矩阵Z转换到地理坐标系,得到目标地形真实的DEM信息。
3.如权利要求2所述的基于多频数据处理的机载D-InSAR形变检测方法,其特征在于,所述步骤2包括:
步骤21、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2基于DEM信息进行BP成像:设置距离向频域滤波器,然后对滤波处理后的矩阵进行距离向插值处理,与步骤1中的操作相同;BP成像的参考网格为DEM信息,以每个网格为相位中心两侧取1/2天线合成孔径长度为雷达成像的方位区域,计算网格与每个雷达方位位置的瞬时斜距,在回波插值矩阵中找到对应的矩阵单元,补偿相位因子并相干叠加处理后得到两个SAR图像矩阵SLC3和SLC4;
步骤22、图像矩阵SLC3和SLC4执行步骤12~步骤15,经图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含检测区域的地形信息的解缠相位矩阵φ12。
4.如权利要求3所述的基于多频数据处理的机载D-InSAR形变检测方法,其特征在于,所述步骤3包括:
步骤31、包含地形信息、地形形变信息的高频辅天线回波矩阵H3基于DEM信息进行BP成像:具体操作与步骤21相同,处理得到图像矩阵SLC5;
步骤32、图像矩阵SLC3和SLC5执行步骤12~步骤15,经图像配准、共轭相乘、滤波处理、相位解缠处理,得到包含检测区域的地形信息、地形形变信息的解缠相位矩阵φ13。
5.如权利要求4所述的基于多频数据处理的机载D-InSAR形变检测方法,其特征在于,所述步骤4包括:
步骤41、将矩阵φ12和φ13代入公式(3),得到仅包含地形形变信息的相位矩阵φd,其中为形变基线与垂直基线水平分量的比值;
步骤42、将矩阵φd带入公式(4),得到形变场矩阵ΔD;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710006318.4A CN106707281B (zh) | 2017-01-05 | 2017-01-05 | 一种基于多频数据处理的机载D-InSAR形变检测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710006318.4A CN106707281B (zh) | 2017-01-05 | 2017-01-05 | 一种基于多频数据处理的机载D-InSAR形变检测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106707281A CN106707281A (zh) | 2017-05-24 |
CN106707281B true CN106707281B (zh) | 2019-01-11 |
Family
ID=58906738
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710006318.4A Active CN106707281B (zh) | 2017-01-05 | 2017-01-05 | 一种基于多频数据处理的机载D-InSAR形变检测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106707281B (zh) |
Families Citing this family (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108088358B (zh) * | 2017-12-18 | 2019-08-20 | 电子科技大学 | 一种基于多基线雷达轨道形变检测方法 |
CN108363056B (zh) * | 2018-01-25 | 2020-08-18 | 北京无线电测量研究所 | 一种多频点干涉成像方法及*** |
CN108594222B (zh) * | 2018-03-21 | 2020-05-08 | 中国科学院电子学研究所 | 一种双频干涉合成孔径雷达的高程重建方法和装置 |
CN108983231B (zh) * | 2018-06-06 | 2021-12-31 | 电子科技大学 | 一种基于视频合成孔径雷达的干涉视频测量方法 |
CN108765353A (zh) * | 2018-06-06 | 2018-11-06 | 中国电子科技集团公司第二十九研究所 | 一种面向sar影像的时序滤波方法 |
CN109839635B (zh) * | 2019-03-13 | 2022-12-27 | 武汉大学 | 一种通过CryoSat-2 SARIn模式L1b级波形数据提取测高脚点高程的方法 |
CN110490827B (zh) * | 2019-08-30 | 2022-04-12 | 三亚中科遥感研究所 | 机载InSAR数据的快速实时处理方法及*** |
CN111562574B (zh) * | 2020-05-22 | 2022-08-16 | 中国科学院空天信息创新研究院 | 基于后向投影的mimo探地雷达三维成像方法 |
CN111983608A (zh) * | 2020-07-07 | 2020-11-24 | 西安瑞得空间信息技术有限公司 | 一种快速dem生成方法 |
CN113640797B (zh) * | 2021-08-09 | 2022-04-12 | 北京航空航天大学 | 一种用于参考条带模式InSAR的前斜视测高方法 |
CN113835089A (zh) * | 2021-08-27 | 2021-12-24 | 中国空间技术研究院 | 一种基于差频InSAR的缓变高程反演方法 |
CN114353709B (zh) * | 2021-12-16 | 2023-08-04 | 武汉滨湖电子有限责任公司 | 一种多阵面天线的面精度调整方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
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处理方法 |
CN105182337A (zh) * | 2015-09-09 | 2015-12-23 | 北京航空航天大学 | 一种基于曲面后向投影算法的形变反演方法 |
-
2017
- 2017-01-05 CN CN201710006318.4A patent/CN106707281B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
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处理方法 |
CN105182337A (zh) * | 2015-09-09 | 2015-12-23 | 北京航空航天大学 | 一种基于曲面后向投影算法的形变反演方法 |
Also Published As
Publication number | Publication date |
---|---|
CN106707281A (zh) | 2017-05-24 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106707281B (zh) | 一种基于多频数据处理的机载D-InSAR形变检测方法 | |
CN103487809B (zh) | 一种基于BP算法和时变基线的机载InSAR数据处理方法 | |
CN111208512B (zh) | 一种基于视频合成孔径雷达的干涉测量方法 | |
CN105487065B (zh) | 一种时序星载雷达数据处理方法和装置 | |
ES2671924T3 (es) | Método para el procesamiento de datos brutos Spotlight SAR de alta resolución obtenidos espacialmente | |
CN109375222B (zh) | 一种合成孔径雷达干涉测量电离层相位估计和补偿方法 | |
CN107741586A (zh) | 基于DBF‑TOPS加权的星载Ka InSAR信号处理方法 | |
CN106249236B (zh) | 一种星载InSAR长短基线图像联合配准方法 | |
CN110488288B (zh) | 一种机载sar高分辨层析成像方法 | |
CN103713287B (zh) | 一种基于互质多基线的高程重建方法及装置 | |
CN112711021B (zh) | 一种多分辨率InSAR交互干涉时序分析方法 | |
CN110007302A (zh) | 一种星载双天线斜视顺轨干涉sar洋流速度测量方法 | |
CN103323840B (zh) | 干涉sar回波数据与平台运动及姿态数据的时间对准方法 | |
CN102866393B (zh) | 一种基于pos与dem数据的sar多普勒参数估计方法 | |
CN102608576B (zh) | 一种用于大前斜合成孔径雷达回波图像的几何校正方法 | |
CN105372657A (zh) | 基于回波数据的视频合成孔径雷达运动补偿成像方法 | |
CN109188432A (zh) | 一种平行双基聚束sar快速bp成像方法 | |
CN115685200A (zh) | 一种高精度大前斜视sar成像运动补偿与几何校正方法 | |
CN107515397A (zh) | 基于InSAR技术高频率地震地区换流站边坡沉降监测方法 | |
CN109239710A (zh) | 雷达高程信息的获取方法及装置、计算机可读存储介质 | |
CN109633639A (zh) | Topsar干涉数据的高精度快速配准方法 | |
KR20170014167A (ko) | 멀티간섭기법 기반의 이온왜곡 보정 방법 및 그 장치 | |
KR20120009186A (ko) | Sar데이터를 이용하여 수치고도모형을 제작하기 위한 방법 | |
CN104330796B (zh) | 一种基于子图像相干合成的地基合成孔径雷达快速成像方法 | |
CN103630898A (zh) | 对多基线干涉sar相位偏置进行估计的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20170524 Assignee: SPACETY Co.,Ltd. (CHANGSHA) Assignor: BEIHANG University Contract record no.: X2023990000783 Denomination of invention: An Airborne D-InSAR Deformation Detection Method Based on Multifrequency Data Processing Granted publication date: 20190111 License type: Common License Record date: 20230829 |
|
EE01 | Entry into force of recordation of patent licensing contract |