CN108983239A - 星载干涉sar数字高程模型重建方法 - Google Patents

Abstract

一种星载干涉SAR数字高程模型重建方法，通过选取符合SAR立体测量条件的两幅SAR图像，利用SAR图像的目标散射特性选取SAR图像的稳定散射点，并利用SAR立体测量技术实现稳定散射点的高效同名点匹配和高精度三维重建，从而可以将稳定散射点作为地面控制点，用于干涉参数的定标，进而实现数字高程模型重建。本发明克服了星载干涉SAR数字高程模型重建时地面控制点获取困难的问题，实现了高精度的数字高程模型重建。

Description

星载干涉SAR数字高程模型重建方法

技术领域

本发明涉及干涉合成孔径雷达信号处理技术领域，尤其涉及一种星载干涉SAR数字高程模型重建方法。

背景技术

干涉合成孔径雷达(Interferometric Synthetic Aperture Radar，InSAR)利用同一区域多次观测的SAR复数据获取干涉相位信息，从而提取地表的高程信息或变化信息的一项技术，它将SAR的测量拓展到三维空间，具有全天时、全天候、高精度的特点，因此在诸多领域都有广泛的应用，地形测绘是其主要应用之一。

星载干涉SAR由于其覆盖范围广、平台稳定的优点，有利于进行大面积、高精度的地形测绘。然而，在数字高程模型(Digital Elevation Model，DEM)重建过程中，由于系统测量误差、信号处理误差等因素的存在，使得重建模型中的相关参数的取值并不准确，需要利用地面控制点(Ground Control Points，GCP)已知的三维位置信息，来修正DEM重建模型中的参数偏差，从而提高DEM重建的精度，这一过程称为干涉定标。对于星载干涉SAR系统，具体包括两方面：一是，解缠后的干涉相位与绝对干涉相位仍然相差一个常数，即2π的整数倍，不能直接用于DEM的反演，需要通过GCP解算解缠相位和真实相位之间的绝对模糊数，生成与DEM对应的绝对相位；二是，星上GPS测量的基线值存在一定误差，需要利用GCP进行定标，以提高最终DEM产品的精度。由此可见，高精度的GCP信息的获取是实现星载干涉SAR高精度DEM的一个必要条件。

目前，干涉SAR定标最常规的方法是建立定标场，在其中人工布设一定数量的角反射器等定标设备，并测量出其准确的三维位置信息，作为GCP进行干涉参数的定标。然而，实际工程中，大量定标器的布设非常耗费人力物力，尤其是在复杂的地形条件下，这一过程难以实现。因而，通常利用定标场进行两次定标之间的间隔时间长，难以实现常态化定标；而且，定标场的设置范围有限，对于星载干涉SAR全球测绘的目标而言，难以实现所有区域的高精度定标。

发明内容

有鉴于此，本发明的主要目的在于提供一种星载干涉SAR数字高程模型重建方法，以期至少部分地解决上述提及的技术问题中的至少之一。

为实现上述目的，本发明的技术方案如下：

一种星载干涉SAR数字高程模型重建方法，包括以下步骤：

步骤A：对星载干涉SAR的主、辅图像对进行干涉处理，得到初始绝对干涉相位、初始基线长度和初始基线倾角；

步骤B：获取待测绘数字高程模型区域的具有相同数据源的两幅SAR复图像，其中一幅为所述星载干涉SAR的主图像，在所述两幅SAR复图像中选取N对散射特征稳定的同名点，其中N≥4；

步骤C：对每一对选出的同名点对，分别在所述两幅SAR复图像中进行点目标分析，通过插值得到各自亚像素级的点目标峰值位置；

步骤D：根据每一对同名点对在所述两幅SAR复图像中的点目标峰值位置以及所述两幅SAR复图像的成像参数文件，获取以下参数：点目标在两幅SAR复图像中的斜距和方位成像时刻，相应成像时刻的SAR传感器速度矢量和位置矢量，以及两幅SAR复图像的成像参数文件中读取的多普勒中心频率；

步骤E：根据步骤D获取的参数，建立距离-多普勒方程组，并对所述距离-多普勒方程组进行最小二乘迭代求解，得到每一对同名点的三维位置矢量；

步骤F：在所述干涉SAR的主图像中，将获得三维位置矢量的N对同名点作为地面控制点，根据星载干涉SAR数字高程模型重建几何关系式，以初始绝对干涉相位、初始基线长度和初始基线倾角为初始参数建立敏感度方程，通过对敏感度方程迭代求解来对星载干涉SAR的基线长度、基线倾角和绝对干涉相位进行定标；

步骤G：将定标后的基线长度、基线倾角和绝对干涉相位代入星载干涉SAR数字高程模型重建几何关系式，反演出星载干涉SAR图像中每一个像素对应的高程值。

从上述技术方案可以看出，本发明的星载干涉SAR数字高程模型重建方法具有以下有益效果：

(1)利用干涉SAR图像中散射特性稳定的点作为GCP，一方面避免了在定标场中布设定标器获取GCP，从而大大降低了DEM测绘中定标过程的人力物力成本，有利于实现常态化定标；另一方面，克服了定标场设置范围有限，不利于高精度全球测绘的缺点；

(2)充分利用丰富的星载SAR数据，构成满足SAR立体测量的数据对，利用SAR立体测量技术，能够获取稳定散射点的高精度三维位置信息，有利于实现高精度干涉参数的定标，进而实现高精度的DEM重建。

附图说明

图1是本发明实施例星载干涉SAR数字高程模型重建方法流程图；

图2是本发明实施例星载于涉SAR成像几何示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

为解决现有技术的问题，本发明提供了一种星载干涉SAR数字高程模型重建方法，选择干涉SAR图像中散射特性稳定的点作为GCP，并引入SAR立体测量技术，充分利用目前丰富的星载SAR数据，构成满足SAR立体测量的数据对，从而获取GCP的三维位置信息，实现干涉参数的定标和DEM的高精度重建，克服了星载干涉SAR DEM重建时获取GCP困难的问题。

本发明提供了一种星载干涉SAR数字高程模型重建方法，图1为本发明实施例星载干涉SAR数字高程模型重建方法流程图，参照图1，该方法包括：

步骤1：对星载干涉SAR主、辅图像对进行干涉处理，得到初始绝对干涉相位、初始基线长度和初始基线倾角。具体包括：

(11)对主、辅图像进行复图像配准；

(12)对配准后的主、辅图像进行共轭相乘，对共轭相乘的结果取相位，生成干涉相位图；

(13)对干涉相位图进行去平地效应，得到去平地后的干涉相位图和平地相位图，这一步骤可根据(12)得到的干涉相位图条纹的疏密情况进行取舍，若条纹比较稀疏，则可以不进行去平地；

(14)若包括步骤(13)，则对去平地后的干涉相位图进行滤波去噪，否则对步骤(2)生成的干涉相位图进行滤波去噪；

(15)对去噪后的干涉相位图进行相位解缠，得到解缠相位。若包括步骤(13)，则将步骤(13)得到的平地相位与解缠相位相加，作为初始绝对干涉相位，否则直接将解缠相位作为初始绝对干涉相位；

(16)根据主、辅图像对应的轨道参数计算初始基线长度、初始基线倾角。

步骤2：获取待测绘DEM区域的具有相同数据源的两幅SAR复图像，作为SAR立体测量的数据对，其中一幅为星载干涉SAR的主图像，在两幅SAR复图像中选取N对散射特征稳定的同名点，其中N≥4。具体方法为：

(21)将星载干涉SAR的主图像作为SAR立体测量的第一幅图像；

(22)选择该测绘区域相同数据源的另一幅SAR图像作为SAR立体测量的第二幅图像。具体而言，若测绘区域为地形平坦的场景，则选择观测方向与第一幅图像相反的数据，即若第一幅图像为升轨观测，则第二幅图像选择降轨观测，反之亦然；若测绘区域为地形起伏的场景，则选择观测方向与第一幅图像相同的数据，即第一幅图像和第二幅图像同为升轨观测或同为降轨观测；

(23)在第一幅图像中通过人工解译识别的方法选择N个分布均匀的稳定散射点，根据第一幅图像的距离-多普勒-椭球方程，计算该像素P₁(x₁，y₁)在地面上的粗略三维位置P_粗，N应大于需要定标的参数个数，这里需要对基线长度、基线倾角、绝对干涉相位3个参数进行定标，因此取N≥4；

根据SAR图像的目标散射特性，稳定散射点主要从以下目标中选择：桥梁、道路、铁路、建筑物、角反射器、裸露岩石等；

(24)然后再利用第二幅图像距离-多普勒-椭球方程计算地面上的粗略三维位置P_粗在第二幅图像中的像素坐标P₂(x₂，y₂)，实现了稳定散射点在第二幅图像中的粗匹配。

其中，根据SAR复图像中的像素坐标结合距离-多普勒-椭球方程来解算对应的地面点三维位置、以及根据该地面点三维位置结合距离-多普勒-椭球方程来解算SAR复图像中的像素坐标，均为本领域的公知常识，在此提供距离-多普勒-椭球方程如下：

R＝|S-P_粗|

其中，R为目标P的斜距，λ为波长，f_dc为多普勒中心频率，P_粗为目标P的三维位置，S和V分别为目标P成像时刻SAR传感器的位置和速度矢量，R_eq和R_p分别为地球的赤道半径和极地半径，h_ave为成像区域相对于大地水准面的平均高度。

(25)在第二幅图像中粗匹配的位置附近区域，再次通过人工解译识别的方法找出与第一幅图像中的稳定散射点对应的同名点。

步骤3：对每一对选出的同名点对，分别在两幅SAR复图像中进行点目标分析，通过插值得到各自亚像素级的点目标峰值位置(i₁，j₁)，(i₂，J₂)，其中l₁、J₁为第一幅图像中点目标点方位向位置和距离向位置，i₂、J₂为第二幅图像中点目标点方位向位置和距离向位置。

步骤4：根据每一对同名点对在两幅SAR复图像中的点目标峰值位置以及两幅SAR复图像的成像参数文件，获取以下参数：点目标在两幅SAR复图像中的斜距和方位成像时刻，相应成像时刻的SAR传感器速度矢量和位置矢量，以及两幅SAR复图像的成像参数文件中读取的多普勒中心频率；

具体地，通过以下公式计算出点目标在两幅SAR复图像中的斜距R₁、R₂和方位成像时刻t₁、t₂：

其中，R_{1_near}，R_{2_near}分别为两幅SAR复图像的近距，c为光速，f_s1，f_s2分别为两幅SAR复图像的采样率，t_{1_start}，t_{2_start}分别为两幅SAR复图像的起始成像时刻，PRF₁，PRF₂分别为两幅SAR复图像的脉冲重复频率。

从两幅SAR复图像的成像参数文件中读取各自的多普勒中心频率，。

根据两幅SAR复图像的方位向成像时刻t₁、t₂和成像参数文件中的轨道状态矢量拟合出相应成像时刻的SAR传感器速度矢量V₁(t₁)，V₂(t₂)和位置矢量S₁(t1)，S₂(t₂)，具体的拟合方法为本领域的公知常识，在此不作赘述。

步骤5：对每一对选出的同名点对，根据步骤4获取的参数，建立如下的距离-多普勒方程组，对该方程组进行最小二乘迭代求解，得到每一对同名点的三维位置矢量P：

R₁＝|S₁(t₁)-P|

R₂＝|S₂(t₂)-P|

其中，λ为波长。

步骤6：在干涉SAR的主图像中，将获得三维位置矢量的点作为GCP，根据星载干涉SAR数字高程模型重建几何关系式，以初始绝对干涉相位、初始基线长度和初始基线倾角为初始参数建立敏感度方程，通过对敏感度方程迭代求解来对星载干涉SAR的基线长度、基线倾角和绝对干涉相位进行定标。

图2为星载干涉SAR几何关系的示意图。结合图2，星载干涉SAR每一像素DEM重建几何关系式如下：

R_m ²+(H+R_e)²-(h+R_t)²＝2R_m(H+R_e)cosθ

其中，φ为绝对干涉相位，B和α分别为基线长度和基线倾角，R_m，R_s分别为主、辅图像中目标T的斜距，H为平台高度，θ为主图像的观测视角，Q为系数，对于星载重轨干涉SAR模式，Q＝2，对于星载双站干涉SAR模式，Q＝1，B和α分别为基线长度和基线倾角，R_e为星下点的地球半径，R_t为目标区域的地球半径，h为反演出的目标高程值。

本步骤的具体方法为：

(61)根据上述DEM重建几何关系式式，以初始绝对干涉相位、初始基线长度和初始基线倾角为初始参数，获取N个地面控制点重建的高程值与真实高程值之间的偏差，并建立敏感度方程：

A·ΔX＝Δh；

其中，为敏感度矩阵，ΔX＝[ΔB，Δα，Δφ]^T为各参数的偏差，Δh＝[Δh₁，...，Δh_N]^T为各个地面控制点重建的高程值与真实高程值之间的偏差。

具体地，第i个地面控制点重建的高程值与真实高程值之间的偏差的计算公式为：

Δh_i＝h_{i_ideal}-h_i(φ₀，B₀，α₀)；

其中，i＝1，2，...N；n_{i_ideal}为第i个控制点的三维位置矢量通过坐标转换得到的真实高程值，h_i(φ₀，B₀，α₀)为将初始绝对干涉相位φ₀、初始基线长度B₀和初始基线倾角α₀代入星载干涉SAR DEM重建几何关系式，而反演出的第i个控制点重建的高程值。

(62)求解上述超定线性方程组的最小二乘解即可获得各定标参数的修正量ΔX，各定标参数包括基线长度、基线倾角和绝对干涉相位。

(63)对定标参数修正后，重新计算各GCP的高程偏差，并更新敏感度矩阵，再次求解各定标参数的修正量，如此迭代，直至定标参数收敛，即完成了参数定标过程。

步骤7：将定标后的基线长度B、基线倾角α和绝对干涉相位φ代入步骤6中的DEM重建几何关系式，反演出星载干涉SAR图像中每一个像素对应的高程值，从而完成星载干涉SAR数字高程模型的重建。

本发明的关键点在于，根据SAR图像的目标散射特性，选出图像中的稳定散射点，并通过选取符合SAR立体测量条件的数据对，利用SAR立体测量技术，实现稳定散射点的高效同名点匹配和高精度三维重建，从而可以将稳定散射点作为GCP，用于干涉参数的定标，进而实现高精度的DEM重建，克服了星载干涉SAR DEM重建时GCP获取困难的问题。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种星载干涉SAR数字高程模型重建方法，其特征在于，包括以下步骤：

步骤A：对星载干涉SAR的主、辅图像对进行干涉处理，得到初始绝对干涉相位、初始基线长度和初始基线倾角；

步骤B：获取待测绘数字高程模型区域的具有相同数据源的两幅SAR复图像，其中一幅为所述星载干涉SAR的主图像，在所述两幅SAR复图像中选取N对散射特征稳定的同名点，其中N≥4；

步骤C：对每一对选出的同名点对，分别在所述两幅SAR复图像中进行点目标分析，通过插值得到各自亚像素级的点目标峰值位置；

步骤D：根据每一对同名点对在所述两幅SAR复图像中的点目标峰值位置以及所述两幅SAR复图像的成像参数文件，获取以下参数：点目标在两幅SAR复图像中的斜距和方位成像时刻，相应成像时刻的SAR传感器速度矢量和位置矢量，以及两幅SAR复图像的成像参数文件中读取的多普勒中心频率；

步骤E：根据步骤D获取的参数，建立距离-多普勒方程组，并对所述距离-多普勒方程组进行最小二乘迭代求解，得到每一对同名点的三维位置矢量；

步骤F：在所述干涉SAR的主图像中，将获得三维位置矢量的N对同名点作为地面控制点，根据星载干涉SAR数字高程模型重建几何关系式，以初始绝对干涉相位、初始基线长度和初始基线倾角为初始参数建立敏感度方程，通过对敏感度方程迭代求解来对星载干涉SAR的基线长度、基线倾角和绝对干涉相位进行定标；

步骤G：将定标后的基线长度、基线倾角和绝对干涉相位代入星载干涉SAR数字高程模型重建几何关系式，反演出星载干涉SAR图像中每一个像素对应的高程值，从而完成星载干涉SAR数字高程模型的重建。

2.根据权利要求1所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤A包括：

子步骤A1：对所述主、辅图像进行复图像配准；

子步骤A2：对配准后的主、辅图像进行共轭相乘，对共轭相乘的结果取相位，生成干涉相位图；

子步骤A3：判断子步骤A2得到的干涉相位图条纹的疏密情况，若条纹密集，则对干涉相位图进行去平地效应，得到去平地后的干涉相位图和平地相位图，并对去平地后的干涉相位图进行滤波去噪，反之，则直接对子步骤A2得到的干涉相位图进行滤波去噪；

子步骤A4：对去噪后的干涉相位图进行相位解缠，得到解缠相位，若子步骤A3包括去平地的步骤，则将子步骤A3得到的平地相位与解缠相位相加，作为初始绝对干涉相位，否则直接将解缠相位作为初始绝对干涉相位；

子步骤A5：根据所述主、辅图像对应的轨道参数计算初始基线长度、初始基线倾角。

3.根据权利要求1所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤B中获取具有相同数据源的两幅SAR复图像的步骤包括：

子步骤B1：将星载干涉SAR的主图像作为SAR立体测量的第一幅图像；

子步骤B2：选择该测绘区域相同数据源的另一幅SAR图像作为SAR立体测量的第二幅图像，其中，若测绘区域为地形平坦的场景，则选择的第二幅图像的观测方向与第一幅图像相反；若测绘区域为地形起伏的场景，则选择的第二幅图像的观测方向与第一幅图像相同。

4.根据权利要求3所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤B中选取N对散射特征稳定的同名点的步骤包括：

子步骤B3：在第一幅图像中通过人工解译识别的方法选择N个分布均匀的稳定散射点，根据第一幅图像的距离-多普勒-椭球方程和每个稳定散射点的像素坐标，计算每个稳定散射点在地面上的粗略三维位置；

子步骤B4：然后利用第二幅图像的距离-多普勒-椭球方程和所述粗略三维位置计算每个稳定散射点在第二幅图像中的像素坐标P₂(x₂，y₂)，实现了稳定散射点在第二幅图像中的粗匹配；

子步骤B5：在第二幅图像中粗匹配的位置附近区域，再次通过人工解译识别的方法找出与第一幅图像中的稳定散射点对应的同名点。

5.根据权利要求4所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述距离-多普勒-椭球方程如下：

R＝|S-P_粗|

其中，R为目标P的斜距，λ为波长，f_dc为多普勒中心频率，P_粗为目标P的三维位置，S和V分别为目标P成像时刻SAR传感器的位置和速度矢量，R_eq和R_p分别为地球的赤道半径和极地半径，h_ave为成像区域相对于大地水准面的平均高度。

6.根据权利要求5所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤B中稳定散射点选自桥梁、道路、铁路、建筑物、角反射器或裸露岩石。

7.根据权利要求5所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤D中，所述点目标在两幅SAR复图像中的斜距R₁、R₂和方位成像时刻t₁、t₂的计算方法为：

其中，R_{1_near}、R_{2_near}分别为两幅SAR复图像的近距，c为光速，f_s1、f_s2分别为两幅SAR复图像的采样率，t_{1_start}、t_{2_start}分别为两幅SAR复图像的起始成像时刻，PRF₁、PRF₂分别为两幅SAR复图像的脉冲重复频率，i₁、j₁分别为第一幅图像中同名点方位向位置和距离向位置，i₂、j₂分别为第二幅图像中同名点方位向位置和距离向位置。

8.根据权利要求1所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤F和步骤G中星载干涉SAR数字高程模型重建几何关系式为：

R_m ²+(H+R_e)²-(h+R_t)²＝2R_m(H+R_e)cosθ

其中，φ为绝对干涉相位，B和α分别为基线长度和基线倾角，R_m，R_s分别为主、辅图像中以同名点为目标T的斜距，H为平台高度，θ为主图像的观测视角，Q为系数，对于星载重轨干涉SAR模式，Q＝2，对于星载双站干涉SAR模式，Q＝1，B和α分别为基线长度和基线倾角，R_e为星下点的地球半径，R_t为目标区域的地球半径，h为反演出的目标高程值。

9.根据权利要求8所述的星载干涉SAR数字高程模型重建方法，其特征在于，所述步骤F中所述对星载干涉SAR的基线长度、基线倾角和绝对干涉相位进行定标的步骤具体包括：

子步骤F1：根据星载干涉SAR数字高程模型重建几何关系式，以初始绝对干涉相位、初始基线长度和初始基线倾角为初始参数，获取N个地面控制点重建的高程值与真实高程值之间的偏差，并建立敏感度方程：

A·ΔX＝Δh：

其中，为敏感度矩阵，ΔX＝[ΔB，Δα，Δφ]^T为各参数的偏差，Δh＝[Δh₁，…，Δh_N]^T为各个地面控制点重建的高程值与真实高程值之间的偏差；

子步骤F2：通过求解所述敏感度方程的最小二乘解来获得定标参数的修正量ΔX，所述定标参数包括基线长度、基线倾角和绝对干涉相位；

子步骤F3：对所述定标参数修正后，重新计算N个地面控制点的高程偏差，并更新敏感度矩阵，再次求解各定标参数的修正量，如此迭代，直至定标参数收敛，即完成了参数定标过程。

10.根据权利要求9所述的星载干涉SAR数字高程模型重建方法，其特征在于，子步骤F1中，第i个地面控制点重建的高程值与真实高程值之间的偏差的计算公式为：

Δh_i＝h_{i_ideal}-h_i(φ₀，B₀，α₀)；

其中，i＝1，2，...N；n_{i_ideal}为第i个控制点的三维位置矢量通过坐标转换得到的真实高程值，h_i(φ₀，B₀，α₀)为将初始绝对干涉相位φ₀、初始基线长度B₀和初始基线倾角α₀代入星载干涉SAR数字高程模型重建几何关系式，而反演出的第i个控制点重建的高程值。