一种基于多频数据处理的机载D-InSAR形变检测方法
技术领域
[0001] 本
发明涉及雷达技术领域，尤其涉及一种基于机载平台，利用低频回
波数据反演地形DEM，将其作为高频回波数据的BP成像参考平面得到相应的SAR图像，通过差分干涉处理获取高
精度地形形变的检测方法。
背景技术
[0002] 机载D-InSAR系统相较于星载D-InSAR技术具有轨道灵活、重访周期短、检测精度高等优势，但同时机载平台存在飞行控制能
力低、航迹误差大等问题，因此如何提高系统检测精度一直是它的研究热点与难点。传统单频D-InSAR系统通过SAR
图像配准、干涉处理、去平地操作、滤波处理、
相位解缠处理、
相位差分处理等步骤得到地形的形变信息。
[0003] 对于星载D-InSAR系统可近似主辅天线到目标的瞬时斜距夹
角为0，因此在形变检测公式推导时可近似主辅天线的斜距差为基线在LOS方向的
水平分量。同时由于星载平台轨道固定，且精度较高，可以采用基于轨道参数的去平地方法，保证平地相位计算精度的同时减小干涉条纹
密度，降低相位解缠难度。但是机载平台的飞行高度较低，此时主辅天线的瞬时斜距夹角不可近似为0，如果仍将斜距差近似为基线的水平分量，会引入较大的系统检测误差。同时机载平台的航迹误差较大，这给去平地操作的具体方法选择带来了一定困难。
[0004] 此外，D-InSAR系统还存在检测精度与数据处理难度相矛盾的问题，突出表现为相位解缠难度较大，因此解决上述问题是提高机载D-InSAR系统检测精度与稳定度的关键。
发明内容
[0005] 针对上述问题中存在的不足之处，本发明提供一种基于多频数据处理的机载D-InSAR形变检测方法。
[0006] 为实现上述目的，本发明提供一种基于多频数据处理的机载D-InSAR形变检测方法，包括：
[0007] 步骤1、低频回波数据反演高程DEM信息：将两个低频回波数据L1、L2经BP成像处理得到两幅SAR图像SLC1、SLC2，SLC1与SLC2通过图像配准、共轭相乘、滤波处理、相位解缠处理、DEM重建与地理编码操作得到目标地形的DEM信息；
[0008] 步骤2、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2干涉处理：将所述DEM信息作为H1与H2的BP成像参考网格，通过相位补偿与相干
叠加得到两幅SAR图像SLC3、SLC4；SLC3与SLC4通过图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含地形信息的解缠相位；
[0009] 步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理：将所述DEM信息作为H3的BP成像参考网格，通过相位补偿与相干叠加得到一幅SAR图像SLC5；SLC5与SLC3通过图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含地形信息、地形形变信息的解缠相位；
[0010] 步骤4、两个解缠相位反演地形形变信息：将步骤2、步骤3中得到的两个解缠相位通过差分处理得到包含地形形变信息的相位，再通过反演处理得到相应的地形形变量。
[0011] 作为本发明的进一步改进，所述步骤1包括：
[0012] 步骤11、BP成像：设置一个1×Nr矩阵的频域
滤波器F1，其中Nr与回波矩阵距离向长度相同，利用距离向调
频率Kr、快时间τ、瞬时斜距R，如公式(1)设置滤波器：
[0013]
[0014] 调用函数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中找到j4πR(η)/λ相应的位置点，补偿相位因子e ，相干叠加处理后得到SAR图像矩阵SLC1和SLC2；
[0015] 步骤12、图像配准：图像矩阵SLC1和SLC2粗配准，在SLC1中心位置取A×B的匹配窗，在SLC2中心位置取C×D搜索窗，其中A＜C，B＜D；在搜索窗内逐行逐列以A×B对应窗口遍历计算其与匹配窗的相干系数，相干系数最大时的偏移量为SLC1、SLC2两图像的相对偏移量，通过矩阵切割处理获取公共部分矩阵Re1、Re2；然后进行亚
像素处理：
[0016] a、对矩阵Re1、Re2进行频域补0处理，插值间隔为0.01个像素，在信号频域将原始信号以中心点分割为四等份，依次位列新插值矩阵的四角边缘区域，其余值填充为0，得到插值后的矩阵In1、In2；
[0017] b、主图像In1匹配窗设置为32×32，辅图像In2搜索窗设置为64×64，通过矩阵遍历计算偏移量，此时主图像的匹配窗数组为100个，偏移量拟合；
[0018] c、辅图像Re2经插值处理、校正偏移量、
采样处理得到配准后的辅图像矩阵Re3；
[0019] 步骤13、共轭相乘：矩阵Re1与Re3的共轭矩阵相乘，得到干涉相位矩阵Minf；
[0020] 步骤14、滤波处理：设置一个1×31的窗口在矩阵Minf中遍历移动，每个信号点的值用邻域各点值的中值代替，得到滤波后的矩阵Mfilter；
[0021] 步骤15、相位解缠处理：调用函数unwrap(Mfilter,[],2)，得到解缠相位矩阵Munwrap；
[0022] 步骤16、DEM重建与地理编码：采用公式(2)反演地形高程：
[0023]
[0024] 其中H为平台高度，R为多普勒斜距信息，B为基线长度，α为基线角，利用轨道参数，基线信息以及系统控制点的瞬时斜距信息，将斜距
坐标系的高程矩阵Z转换到地理坐标系，得到目标地形真实的DEM信息。
[0025] 作为本发明的进一步改进，所述步骤2包括：
[0026] 步骤21、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2基于DEM信息进行BP成像：设置距离向频域滤波器，然后对滤波处理后的矩阵进行距离向插值处理，与步骤1中的操作相同；BP成像的参考网格为DEM信息，以每个网格为相位中心两侧取1/2天线合成孔径长度为雷达成像的方位区域，计算网格与每个雷达方位位置的瞬时斜距，在回波插值矩阵中找到对应的矩阵单元，补偿相位因子并相干叠加处理后得到两个SAR图像矩阵SLC3和SLC4；
[0027] 步骤22、图像矩阵SLC3和SLC4执行步骤12～步骤15，经图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含
检测区域的地形信息的解缠相位矩阵φ12。
[0028] 作为本发明的进一步改进，所述步骤3包括：
[0029] 步骤31、包含地形信息、地形形变信息的高频辅天线回波矩阵H3基于DEM信息进行BP成像：具体操作与步骤21相同，处理得到图像矩阵SLC5；
[0030] 步骤32、图像矩阵SLC3和SLC5执行步骤12～步骤15，经图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含检测区域的地形信息、地形形变信息的解缠相位矩阵φ13。
[0031] 作为本发明的进一步改进，所述步骤4包括：
[0032] 步骤41、将矩阵φ12和φ13代入公式(3)，得到仅包含地形形变信息的相位矩阵φd，其中 为形变基线与垂直基线水平分量的比值；
[0033]
[0034] 步骤42、将矩阵φd带入公式(4)，得到形变场矩阵ΔD；
[0035]
[0036] 与
现有技术相比，本发明的有益效果为：
[0037] 1、本发明提供的方法可以控制系统近似计算引入的误差，提高系统检测精度；
[0038] 2、本发明提供的方法能克服传统D-InSAR系统检测精度与数据处理复杂度之间的矛盾，有效获取高频数据的信息，降低相位解缠难度，提高系统检测精度与稳定度；
[0039] 3、本发明提供的方法省略了去平地操作，简化了干涉数据处理流程。
附图说明
[0040] 图1为本发明一种
实施例公开的基于多频数据处理的机载D-InSAR形变检测方法的
流程图；
[0041] 图2为本发明利用的将DEM信息作为BP成像处理参考网格的示意图；
[0042] 图3(a)为检测区域的地形信息；
[0043] 图3(b)为检测区域的形变场信息；
[0044] 图4(a)为低频地形信息干涉相位；
[0045] 图4(b)为低频地形信息滤波后相位；
[0046] 图4(c)为低频地形信息解缠相位；
[0047] 图5为低频回波数据反演的DEM信息；
[0048] 图6(a)为基于多频处理的地形信息干涉相位；
[0049] 图6(b)为基于多频处理的地形信息滤波后相位；
[0050] 图6(c)为基于多频处理的地形信息解缠相位；
[0051] 图7(a)为基于多频处理的形变后地形干涉相位；
[0052] 图7(b)为基于多频处理的形变后地形滤波相位；
[0053] 图7(c)为基于多频处理的形变后地形解缠相位；
[0054] 图8(a)为沿LOS方向的形变检测结果；
[0055] 图8(b)为形变检测误差场。
具体实施方式
[0056] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0057] 为了解决机载D-InSAR系统近似误差大、去平地操作误差大以及相位解缠难度大的问题，本发明根据系统近似误差产生的原因以及干涉条纹密度对相位解缠操作的影响，采用高低两种频率的回波数据融合处理的方式提高系统的检测精度。首先利用低频天线干涉处理获得地形的粗精度DEM，然后将该DEM信息作为高频天线BP成像处理的参考平面，以控制系统的近似误差与高频天线的干涉条纹密度，最后利用得到的高频SAR图像进行差分干涉处理获取地形的形变信息，并在处理流程中省略去平地操作，简化数据处理流程的同时进一步控制数据处理误差，实现了一种能获取高精度地表形变信息的机载D-InSAR检测方法。
[0058] 下面结合附图对本发明做进一步的详细描述：
[0059] 如图1所示，本发明提供一种基于多频数据处理的机载D-InSAR形变检测方法，包括：
[0060] 步骤1、低频回波数据反演DEM信息：将两个低频回波矩阵L1和L2经BP成像处理得到两个SAR图像矩阵SLC1和SLC2，将矩阵SLC1与SLC2通过图像配准、共轭相乘、滤波处理、相位解缠处理得到一个解缠相位矩阵Munwrap，最后通过DEM重建与地理编码操作得到目标地形真实的DEM信息。
[0061] 步骤11、BP成像：设置一个1×Nr矩阵的频域滤波器F1，其中Nr与回波矩阵距离向长度相同，利用距离向调频率Kr、快时间τ、瞬时斜距R，如公式(1)设置滤波器：
[0062]
[0063] 调用函数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(η)，在j4πR(η)/λ矩阵M2中找到相应的位置点，补偿相位因子e ，相干叠加处理后得到SAR图像矩阵SLC1和SLC2；
[0064] 步骤12、图像配准：图像矩阵SLC1和SLC2粗配准，在SLC1中心位置取A×B的匹配窗，在SLC2中心位置取C×D搜索窗，其中A＜C，B＜D；在搜索窗内逐行逐列以A×B对应窗口遍历计算其与匹配窗的相干系数，相干系数最大时的偏移量为SLC1、SLC2两图像的相对偏移量，通过矩阵切割处理获取公共部分矩阵Re1、Re2；然后进行亚像素处理，主要有三步：
[0065] a、对矩阵Re1、Re2进行频域补0处理，插值间隔为0.01个像素，在信号频域将原始信号以中心点分割为四等份，依次位列新插值矩阵的四角边缘区域，其余值填充为0，得到插值后的矩阵In1、In2；
[0066] b、主图像In1匹配窗设置为32×32，辅图像In2搜索窗设置为64×64，通过矩阵遍历计算偏移量，此时主图像的匹配窗数组为100个，偏移量拟合；
[0067] c、辅图像Re2经插值处理、校正偏移量、采样处理得到配准后的辅图像矩阵Re3；
[0068] 步骤13、共轭相乘：矩阵Re1与Re3的共轭矩阵相乘，得到干涉相位矩阵Minf；
[0069] 步骤14、滤波处理：设置一个1×31的窗口在矩阵Minf中遍历移动，每个信号点的值用邻域各点值的中值代替，得到滤波后的矩阵Mfilter；
[0070] 步骤15、相位解缠处理：调用函数unwrap(Mfilter,[],2)，得到解缠相位矩阵Munwrap；
[0071] 步骤16、DEM重建与地理编码：采用公式(2)反演地形高程Z：
[0072]
[0073] 其中H为平台高度，R为多普勒斜距信息，B为基线长度，α为基线角，利用轨道参数，基线信息以及系统控制点的瞬时斜距信息，将斜距坐标系的高程矩阵Z转换到地理坐标系，得到目标地形真实的DEM信息。
[0074] 步骤2、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2干涉处理：对矩阵H1和H2做BP成像，将DEM信息作为参考网格，通过计算雷达与网格之间的瞬时斜距得到其在回波插值矩阵中的具体位置，进而通过相位补偿与相干叠加得到SAR图像矩阵SLC3和SLC4，两矩阵SLC3和SLC4通过图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含地形信息的解缠相位矩阵φ12。
[0075] 步骤21、高频主天线回波数据H1与包含地形信息的高频辅天线回波数据H2基于DEM信息进行BP成像：设置距离向频域滤波器，然后对滤波处理后的矩阵进行距离向插值处理，与步骤1中的操作相同；BP成像的参考网格为DEM信息，以每个网格为相位中心两侧取1/2天线合成孔径长度为雷达成像的方位区域，计算网格与每个雷达方位位置的瞬时斜距，在回波插值矩阵中找到对应的矩阵单元，补偿相位因子并相干叠加处理后得到两个SAR图像矩阵SLC3和SLC4，具体实现方法如图2所示；
[0076] 步骤22、图像矩阵SLC3和SLC4执行步骤12～步骤15，经图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含检测区域的地形信息的解缠相位矩阵φ12。
[0077] 步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理：对H3进行BP成像，参考网格为DEM信息，计算雷达与网格之间的瞬时斜距得到其在回波插值矩阵中的对应矩阵单元，进行相位补偿与相干叠加处理得到SAR图像矩阵SLC5，矩阵SLC3与SLC5通过图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含地形信息、地形形变信息的解缠相位矩阵φ13。
[0078] 步骤31、包含地形信息、地形形变信息的高频辅天线回波矩阵H3基于DEM信息进行BP成像：具体操作与步骤21相同，处理得到图像矩阵SLC5；
[0079] 步骤32、图像矩阵SLC3和SLC5执行步骤12～步骤15，经图像配准、共轭相乘、滤波处理、相位解缠处理，得到包含检测区域的地形信息、地形形变信息的解缠相位矩阵φ13。
[0080] 步骤4、两个解缠相位反演地形形变信息：将步骤2、步骤3中得到的两个解缠相位通过差分处理得到包含地形形变信息的相位，再通过反演处理得到相应的地形形变量。
[0081] 步骤41、将矩阵φ12和φ13代入公式(3)，得到仅包含地形形变信息的相位矩阵φd，其中 为形变基线与垂直基线水平分量的比值；
[0082]
[0083] 步骤42、将矩阵φd带入公式(4)，得到形变场矩阵ΔD；
[0084]
[0085] 实施例1：
[0086] 针对一个前后坡场景、凹陷的圆锥形变场进行地形形变检测，雷达成像方式为正侧视。其中形变场信息如图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。
[0087] 步骤1、低频回波数据反演DEM信息。
[0088] (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；
[0089] (b)图像配准：输入矩阵SLC1和SLC2，在SLC1中心位置取50×100的匹配窗，在SLC2中心位置取200×400搜索窗，在搜索窗内逐行逐列以50×100大小的对应窗口遍历计算相干系数，相干系数最大时偏移量为0。亚像素级配准时匹配窗的大小设置为20×20，搜索窗的大小设置为50×50，匹配窗的个数为100个，通过偏移量拟合后，对辅图像进行插值、校正、采样处理，得到配准后的辅图像矩阵Re3，大小为625×4000；
[0090] (c)共轭相乘：矩阵SLC1与Re3的共轭矩阵相乘，得到干涉矩阵Minf如图4(a)所示；
[0091] (d)滤波处理：设置一个1×31的窗口在矩阵Minf中遍历移动，每个信号点的值用邻域各点值的中值代替，滤波后的矩阵Mfilter如图4(b)所示；
[0092] (e)相位解缠处理：调用函数unwrap(Mfilter,[],2)，得到解缠相位矩阵Munwrap，如图4(c)所示；
[0093] (f)DEM重建与地理编码，利用公式(2)得到高程信息矩阵Z，地理编码后得到DEM信息，如图5所示。
[0094] 步骤2、高频主天线回波矩阵H1与包含地形信息的高频辅天线回波矩阵H2做干涉处理。
[0095] (a)矩阵H1和H2将DEM信息作为BP成像处理的参考网格。设置距离向匹配滤波器，滤波后的矩阵距离向频域补0插值，与步骤1中的操作相同。合成孔径长度为177m，以DEM信息中每个单元为方位中心取两侧88.5m范围内对应的各个雷达位置，计算该矩阵单元到每个雷达位置的瞬时斜距，确定该矩阵单元在距离向插值信号中对应的具体位置，通过相位补偿处理和相干叠加后得到SAR图像矩阵SLC3和SLC4。
[0096] (b)图像矩阵SLC3和SLC4执行步骤1中的操作(b)～(e)，依次得到包含地形信息的干涉相位矩阵如图6(a)、滤波矩阵如图6(b)；和包含地形信息的解缠相位矩阵φ12如图6(c)。
[0097] 步骤3、高频主天线回波数据H1与包含地形信息、地形形变信息的高频辅天线回波数据H3干涉处理。
[0098] (a)矩阵H3将DEM信息作为BP成像处理的参考网格。设置距离向匹配滤波器，滤波后的矩阵距离向频域补0处理，与步骤2中(a)的操作相同。合成孔径长度为177m，以DEM信息中每个单元为方位中心取两侧88.5m范围内对应的各个雷达位置，计算该矩阵单元到每个雷达位置的瞬时斜距，确定该矩阵单元在距离向插值信号中对应的具体位置，通过相位补偿处理和相干叠加后得到SAR图像矩阵SLC5。
[0099] (b)图像矩阵SLC3和SLC5执行步骤1中的操作(b)～(e)，依次得到包含地形信息、地形形变信息的干涉相位矩阵如图7(a)、滤波矩阵如图7(b)和包含地形信息、地形形变信息的解缠相位矩阵φ13如图7(c)。
[0100] 步骤4、矩阵φ12和φ13差分处理得到检测区域的地形形变量。
[0101] (a)将矩阵φ12和φ13代入公式(3)，得到仅包含地形形变信息的相位矩阵φd；
[0102] (b)将矩阵φd带入公式(4)，得到形变场矩阵ΔD，如图8(a)所示，根据理论形变场对检测结果进行分析，得到的误差场如图8(b)所示，误差均值为4.866×10-4m，最大误差为0.43mm，均方差为4.65×10-2mm，达到了mm级的检测精度且稳定度较高，检测结果理想。
[0103] 本发明主要针对机载D-InSAR系统近似误差大、去平地操作不理想、相位解缠难度大的
缺陷，利用高低两种频率的回波数据融合处理，通过低频回波数据干涉处理反演得到高程DEM信息，将该DEM作为高频回波BP成像处理时的参考网格，得到相应的SAR图像，在差分干涉处理中控制系统近似误差，省略去平地操作，同时降低了干涉条纹密度，减小数据数理难度，提高系统检测精度与稳定度，实现了一种高精度的基于多频数据融合处理的机载D-InSAR形变检测方法。其具有以下优点：
[0104] 1、本发明提供的方法可以控制系统近似计算引入的误差，提高系统检测精度；
[0105] 2、本发明提供的方法能克服传统D-InSAR系统检测精度与数据处理复杂度之间的矛盾，有效获取高频数据的信息，降低相位解缠难度，提高系统检测精度与稳定度；
[0106] 3、本发明提供的方法省略了去平地操作，简化了干涉数据处理流程。
[0107] 以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
