CN105929398A - 结合外控点的InSAR高精度高分辨率DEM获取方法 - Google Patents

Abstract

本发明公开了一种结合外控点的InSAR高精度高分辨率DEM获取方法，充分融合了InSAR监测手段分辨率高、相对精度高、覆盖范围大和地面控制点绝对精度高的特点，既解决InSAR获取DEM绝对精度不高，又克服了传统DEM获取方法例如星载和机载光学及激光测高手段获取高分辨测图效费比低和受云层影响，地面测图费时费力且无法得到整个野外山区高分辨地貌等缺点。整个流程结构清晰，具有实现简单、费用低、监测精度高、监测范围大、自动化程度高等优点。

Description

结合外控点的InSAR高精度高分辨率DEM获取方法

技术领域

本发明属于基于遥感影像的大地测量领域，是一种利用InSAR(合成孔径干涉雷达)技术和外部控制点获取高分辨率高精度DEM的方法。

背景技术

由于具备全天候，大范围，高精度的观测能力，合成孔径干涉雷达测量(Interferometric SyntheticAperture radar，InSAR)已经成为国内外热点关注技术。随着各种SAR传感器的升空，例如JERS、ERS-1/2、ENVISAT、ALOS/PALSAR、ALOS-2/PALSAR-2、RadarSat-1/2、TerraSAR-X、COSMO-SkyMed、Sentinel-1等，可用SAR数据的时间和空间分辨率不断提高，该技术在监测由于地下水抽取、煤矿开采、石油开采、地下铁路修建、填海、地震、冰川等引起的地表形变方面得到了空前的发展和应用。InSAR监测地表形变主要是借助于覆盖同一个地区的两幅或两幅以上SAR图像，利用包含在SAR图像中的差分相位信息获取来实现。令Δφ_def为地表形变引起的相位变化，λ为雷达波长，则形变值ΔD可通过下面的模型得到：

$\mathrm{\Δ}D=\frac{\mathrm{\λ}}{4\mathrm{\π}}{\mathrm{\Δ\φ}}_{def}---\left(1\right)$

当Δφ_def＝2π时，对应的形变大小，即形变模糊度为：

$\mathrm{\Δ}D=\frac{\mathrm{\λ}}{2}---\left(2\right)$

除监测形变以外，InSAR技术另一个基本应用是通过在同一地区不同角度拍摄的SAR图像的相位差别来提取地形起伏。令Δφ_top为地表起伏Δh引起的相位变化，R为雷达视距，θ⁰为雷达入射角，B_⊥为垂直基线，则相位差到地形起伏的转换模型为:

$\mathrm{\Δ}h=-\frac{\mathrm{\λ}}{4\mathrm{\π}}\·\frac{{\mathrm{Rsin\θ}}^0}{B_{\⊥}}{\mathrm{\Δ\φ}}_{top}---\left(3\right)$

干涉相位图上2π相位变化引起的高程变化，即模糊高为

$h_{2\mathrm{\π}}=\left|\frac{{\mathrm{\λRsin\θ}}^0}{2B_{\⊥}}\right|---\left(4\right)$

λ一般为几厘米级到几分米，B_⊥一般为几十米到上千米，而R为几百千米，显然，地形起伏对相位变化的敏感度要远大于形变对相位变化的敏感度。常规的重复轨道观测相位变化往往包含了两次观测间大气变化引起的相关部分和去相干引起的噪声。这些相位误差会极大地影响InSAR高程提取结果。因此，除美国航空局NASA的SRTM项目(以获取DEM为目标的双天线单次观测)以外，其他单天线重复观测SAR项目长期以来以形变监测应用为主。DEM即为数字高程模型，鉴于全球多个领域对高精度高分辨率DEM的需求，德国宇航局DLR开发了具有创新性的一发双收观测模式，即双基(Bistatic)模式。由串联飞行的TerraSAR-X或TanDEM-X卫星发射雷达信号，而TerraSAR-X与TanDEM-X卫星同时接收回波信号。两个卫星接收回波信号的时间差极短(一分钟内)，这样保证了两次成像的大气状态和地表散射特性的高度相似性，进而保证的相位观测的精度。双基模式下相位差到地形起伏的转换模型为:

$\mathrm{\Δ}h=-\frac{\mathrm{\λ}}{2\mathrm{\π}}\·\frac{{\mathrm{Rsin\θ}}^0}{B_{\⊥}}{\mathrm{\Δ\φ}}_{top}---\left(5\right)$

目前TerraSAR-X和TanDEM-X获取的影像对(即TanDEM-XCoSSC影像对)像元大小可以达到3m以内，通过迭代差分干涉的方式能满足高分辨DEM生产像元大小这一块的要求。但是通过差分干涉方法获取的DEM是相对于外部DEM的高程值，其绝对精度与解缠参考点的选取有很大关系，需要外部控制点来进行绝对高程校正；而且，通过InSAR手段获取的初始DEM，其平面位置精度主要由地理编码决定，具体来说是通过真实SAR影像与外部DEM模拟的SAR影像之间的强度匹配来实现(配准法则是窗口内强度相关系数达到最大)。匹配的结果一般是通过拟合得到的表征系统偏移量的二次多项式来检核，也就是基于内符合精度。尽管目前SAR影像星历参数和轨道数据精度均较高，但山区坡度起伏大的特点决定即使细小的平面偏差也会引起较大的高程误差。因此需要一种匹配法则不同于强度相关最大的配准方法来检核和校正初始DEM的平面和垂直精度。此外，目前通过轨道来估计基线一般都包含残差。基线残差影响传递到DEM产品表现为平面上的系统误差趋势；SAR的斜距成像使得影像不可避免地出现透视收缩。透视收缩的强度与局部坡度和坡向密切相关。因此还要考虑初始DEM产品与坡度和坡向相关的系统误差。上述系统误差都要借助外部控制点来进行最小二乘平差来削弱或者消除。

综上所述，目前采用TanDEM-X CoSSC数据和InSAR差分干涉手段可以得到高分辨的DEM产品，但其高程是相对的，而且受各种现有的InSAR干涉处理过程中无法完全去除的误差影响，需要外部控制点来进行校正以提高绝对的平面和垂直精度。

发明内容：

为了克服上述背景技术的缺陷，本发明提供一种结合外控点的InSAR高精度高分辨率DEM获取方法，具有实现简单、监测精度高范围大、自动化程度高的优点。

为了解决上述技术问题本发明的所采用的技术方案为：

一种结合外控点的InSAR高精度高分辨率DEM获取方法，包括：

步骤1，通过双基InSAR差分干涉处理单轨影像对，获取单轨DEM；

步骤2，采取不同轨道的影像对分别获取DEM，并进行逐点融合；

步骤3，将经融合的DEM与地面控制点进行配准；

步骤4，校正DEM与坡度、平面位置和坡向相关的系统误差。

较佳地，步骤1包括：

步骤1.1，根据DEM需求区域坐标范围定制影像对，影像对包括主影像和辅影像；

步骤1.2，将主影像和辅影像共轭相乘生成初始干涉图；

步骤1.3通过轨道数据估计基线，根据基线由外部DEM模拟包含平地相位的地形相位图；

步骤1.4将地形相位图与初始干涉相位图相减得到差分干涉图；

步骤1.5对差分干涉图进行相位滤波和相位解缠，将解缠相位转化为高程差；

步骤1.6将高程差与外部DEM叠加得到新DEM；

步骤1.7通过地理编码将新DEM从SAR坐标系转换到外部DEM的地理坐标系得到单轨DEM。

较佳地，还包括：将步骤1.7所得单轨DEM作为外部DEM，对步骤1.1至步骤1.7进行迭代，直至获取原始影像像元大小的DEM作为单轨DEM。

较佳地，步骤1.5对差分干涉图进行相位滤波采用的滤波方式包括Nonlocal均值滤波方式，滤波模型包括：

$u\left(s\right)=\frac{\underset{t}{\Σ}w(s,t)\mathrm{\φ}\left(t\right)}{\underset{t}{\Σ}w(s,t)}$

s为滤波窗口内的中心像元，t为滤波窗口内的其他任意像元，真实的干涉相位为u，含噪声的原始干涉相位为φ，相位估计值u(s)为滤波结果，w(s,t)为由像元t的邻域像元与像元s的邻域像元之间的SAR影像强度和相干性的综合相似度而决定的权重。

较佳地，步骤2的具体步骤包括：

步骤21，将不同轨道的影像对统一配准到相同坐标系之下，不同轨道的影像对包括升轨和降轨的影像对；

步骤22，采用融合模型进行融合，融合模型包括其中，H_fusion为融合后高程，B_a、W_a和H_a分别为升轨的垂直基线值，叠影或阴影掩膜值及高程值，B_d、W_d和H_d分别为降轨的垂直基线值，叠影或阴影掩膜值及高程值。

较佳地，步骤3包括：

步骤3.1，将外部控制点的空间直角坐标转为大地经纬度坐标；

步骤3.2，采用双线性插值提取外部控制点所在DEM高程值；

步骤3.3，将经融合的DEM与数个外部控制点进行统计配准，配准模型包括

其中，β为控制点位置的局部坡向，α为控制点位置的局部坡度，为所有控制点位置坡度的平均值，a为DEM平面偏移矢量的模，b为DEM偏移矢量的方向角，ΔZ为DEM垂直偏移，Δh为DEM任一点的高程值与控制点高程值的差。

较佳地，对配准过程进行迭代，直至高程差标准差改善率低于2％。

较佳地，步骤4对DEM结果与平面位置相关的误差进行校正，是通过二次曲面拟合来完成的，二次曲面拟合模型包括

Δh_(x,y)＝a₁x²+a₂y²+a₃xy+a₄x+a₅y+a₆

其中x和y表示某个像元的横坐标和纵坐标，Δh_(x,y)为该点的DEM值与控制点高程差，a₁……a₆为待估模型参数。采用解算控制点位的位置与高程差作为观测值，通过最小二乘平差方法来求解待估参数a₁……a₆，然后将模型求取的误差趋势值与经步骤3改正后DEM进行相减,得到进一步改正后的DEM。

较佳地，步骤4对DEM结果与坡度相关系统误差进行校正是通过二次或三次多项式模型来完成的，

二次或三次多项式包括Δh_(slope)＝b₁α²+b₂α+b₃或Δh_(slope)＝b₁α³+b₂α²+b₃α+b₄，其中α表示某个像元点的坡度值，Δh_(slope)为该点的DEM值与控制点高程差b₁……b₄为待估模型参数，

采用解算控制点位的坡度值与高程差作为观测值，通过最小二乘平差方法来求解待估参数b₁……b₄，然后将模型求取的误差趋势与经过平面位置误差改正后的DEM进行相减，得到进一步改正后的DEM。

较佳地，步骤4对DEM结果与坡向相关系统误差进行校正是通过正弦模型来完成的，正弦模型为Δh_aspect＝A₁sin(ωβ+φ)+C₁

其中β为某个像元点的坡向值，Δh_aspect为该点的DEM值与控制点高程差，A₁、ω、φ和C₁为待估模型参数。

采用解算控制点位的坡向值与高程差作为观测值，通过最小二乘平差方法来求解待估参数A₁、ω、φ和C₁，然后将模型求取的误差趋势与经过坡度误差改正后的DEM进行相减，得到进一步改正后的DEM。

本发明的有益效果在于：本发明所述方法在差分干涉InSAR获取大范围高分辨率DEM的基础上，针对差分干涉获取的高程是相对值，SAR基线估计和地理编码包含误差，SAR成像受透视收缩影响等主要问题，采用外部控制点来进行约束和改善，提高InSAR差分干涉技术获取的高分辨率DEM的精度。充分融合了InSAR监测手段分辨率高、相对精度高、覆盖范围大和地面控制点绝对精度高的特点，既解决InSAR获取DEM绝对精度不高，又克服了传统DEM获取方法例如星载和机载光学及激光测高手段获取高分辨测图效费比低和受云层影响，地面测图费时费力且无法得到整个野外山区高分辨地貌等缺点。整个流程结构清晰，具有实现简单、费用低、监测精度高、监测范围大、自动化程度高等优点。

附图说明

图1是本发明实施例的方法流程图。

图2是本发明实施例统计配准时平面偏移矢量大小在坡向为β的山坡处投影S_xy的示意图；

图3是本发明实施例统计配准时DEM某点位的高程与真实高程之间之差Δh的示意图；

图4是本发明实施例获取的DEM与地面水准检查点的对比结果图。

具体实施方式

下面结合附图和实施例对本发明做进一步的说明。

一种结合外控点的InSAR高精度高分辨率DEM获取方法，所述外控点即为外部控制点，如图1所示，包括：

步骤1，通过双基InSAR差分干涉处理单轨影像对，获取单轨DEM；

步骤1.1，根据DEM需求区域坐标范围定制影像对，影像对包括主影像和辅影像；

步骤1.2，将主影像和辅影像共轭相乘生成初始干涉图；

步骤1.3，通过轨道数据估计基线，根据基线由外部DEM模拟包含平地相位的地形相位图；

步骤1.4，将地形相位图与初始干涉相位图相减得到差分干涉图；

步骤1.5，对差分干涉图进行相位滤波和相位解缠，将解缠相位转化为高程差；

对差分干涉图进行相位滤波采用的滤波方式包括Nonlocal均值滤波方式，所述的Nonlocal均值滤波方式即为非局部均值滤波方式，滤波模型包括：

$u\left(s\right)=\frac{\underset{t}{\Σ}w(s,t)\mathrm{\φ}\left(t\right)}{\underset{t}{\Σ}w(s,t)}---\left(6\right)$

s为滤波窗口内的中心像元，t为滤波窗口内的其他任意像元，真实的干涉相位为u，含噪声的原始干涉相位为φ，相位估计值u(s)为滤波结果，w(s,t)为由像元t的邻域像元与像元s的邻域像元之间的SAR影像强度和相干性的综合相似度而决定的权重。定义|z|、|z′|分别为两幅SAR影像的强度。s,k、t,k分别表示两个以像元s和t为中心的小窗口中处于相同位置的第k个像元。令P(k)为像元s,k与t,k的相似度。而P(k)的计算公式为：

$P\left(k\right)={\sqrt{\frac{C}{B}}}^3(\frac{A+B}{A}\sqrt{\frac{B}{A-B}}-\arcsin \sqrt{\frac{B}{A}})---\left(7\right)$

A＝(|z_s,k|²+|z_t,k|²+|z'_s,k|²+|z'_t,k|²) (8)

B＝4(|z_s,k|²|z'_s,k|²+|z_t,k|²|z'_t,k|²+2|z_s,k||z'_s,k||z_t,k||z'_t,k|cos(φ_s,k-φ_t,k)) (9)

C＝|z_s,k||z'_s,k||z_t,k||z'_t,k| (10)

$w(s,t)=\underset{k}{\Π}p{\left(k\right)}^{1/h}---\left(11\right)$

其中h为平滑参数，h越大权值的差异性越小，滤波的强度则越大。

步骤1.6，将高程差与外部DEM叠加得到新DEM；

步骤1.7，通过地理编码将新DEM从SAR坐标系转换到外部DEM的地理坐标系得到单轨DEM。

步骤1.8，将步骤1.7所得单轨DEM作为外部DEM，对步骤1.1至步骤1.7进行迭代，直至获取原始影像像元大小的DEM作为单轨DEM。

一次标准的差分干涉InSAR提取DEM的过程即完成。因为TanDEM-X影像初始分辨率(3m左右)与外部DEM(例如30m分辨率的SRTM数据)相差较大，在利用外部DEM进行地理编码和相位模拟时，如果采用较大过采样倍数会引入较大误差，所以初次生成DEM时需对TanDEM-X CoSSC影像多视以缩小外部DEM的过采样倍数。通过迭代的方式来逐级提高DEM产品的分辨率，即将上一次生成的新DEM作为外部DEM重复上述InSAR差分干涉过程，直到获取TanDEM-X CoSSC单视下的DEM产品。

在上述过程中，相位滤波可以极大地影响相位图的质量，进而影响生成DEM的精度，因此至关重要。本套方法采用的是考虑滤波窗口内像元异质性的Non-Local均值滤波，该方法可以有效抑制相位噪声而不损害相位的原本真实信息。Non-Local均值滤波方法的核心是对滤波窗口内的像素值进行加权平均后赋予中心像元点。给予与中心像元相似度高的像元更高的权值，反之与中心像元差异性越大更低的权值。因此权值的计算是否合理直接关系到该估计方法的可靠性和精度。

步骤2，采取不同轨道的影像对分别获取DEM，并进行逐点融合以降低SAR影像透视收缩和阴影的影响。

步骤21，将不同轨道的影像对统一配准到相同坐标系之下，不同轨道的影像对主要是时间间隔较短的升轨和降轨影像对；

步骤22，采用融合模型进行融合，融合模型包括

$H_{fusion}=\frac{B_a\·W_a\·H_a+B_d\·W_d\·H_d}{B_a\·W_a+B_d\·W_d}---\left(12\right)$

其中，H_fusion为融合后高程，B_a、W_a和H_a分别为升轨的垂直基线值，叠影或阴影掩膜值(叠影或阴影为真则该值为0，否则为1)及高程值，B_d、W_d和H_d分别为降轨的垂直基线值，叠影或阴影掩膜值及高程值。

步骤3，将经融合的DEM与地面控制点进行配准；

步骤3.1，将外部控制点的空间直角坐标转为大地经纬度坐标；

国内常见的控制点平面基准为北京54平面高斯投影坐标系和1985黄海高程基准，而InSAR生成的DEM为WGS84大地经纬度坐标系和大地高(WGS84椭球高)。

首先将外部控制点的北京54平面高斯投影坐标转换到北京54大地经纬度坐标，即(x,y,h)_beijing54→(B,L,H)_beijing54。转换过程中需要已知高斯投影参数如中央子午线。

其次将外部控制点的北京54大地经纬度坐标转换到北京54空间直角坐标，即(B,L,H)_beijing54→(X,Y,Z)_beijing54。

然后再将外部控制点的北京54空间直角坐标转换到WGS84空间直角坐标，即(X,Y,Z)_beijing54→(X,Y,Z)_WGS84具体转换过程如下：

采用布尔莎七参数模型将外部控制点的北京54空间直角坐标转换到WGS84空间直角坐标，公式如下：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS84}=(1+m){\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{beijing54}+\left[\begin{array}{ccc}0& +{\mathrm{\ϵ}}_z& -{\mathrm{\ϵ}}_y\\ -{\mathrm{\ϵ}}_z& 0& +{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y& -{\mathrm{\ϵ}}_x& 0\end{array}\right]\·{\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{beijing54}+\left[\begin{array}{c}\mathrm{\Δ}{X}_0\\ \mathrm{\Δ}{Y}_0\\ {\mathrm{\ΔZ}}_0\end{array}\right]---\left(13\right)$

其中ΔX₀，ΔY₀，ΔZ₀是平移参数；m是尺度比参数，ε_x,ε_y,ε_z是旋转参数。

对其进行变换得到如下式子：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS84}={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{beijing54}+\left[\begin{array}{ccccccc}1& 0& 0& X_{beijing54}& 0& +{\mathrm{\ϵ}}_z& -{\mathrm{\ϵ}}_y\\ 0& 1& 0& Y_{beijing54}& -{\mathrm{\ϵ}}_z& 0& +{\mathrm{\ϵ}}_x\\ 0& 0& 1& Z_{beijing54}& {\mathrm{\ϵ}}_y& -{\mathrm{\ϵ}}_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δ}{X}_0\\ \mathrm{\Δ}{Y}_0\\ {\mathrm{\ΔZ}}_0\\ m\\ {\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_z\end{array}\right]---\left(14\right)$

解算这七参数，至少需要三个已知的公共控制点的北京54空间直角坐标(X,Y,Z)_beijing54和WGS-84空间直角坐标(X,Y,Z)_WGS84，采用间接平差模型最小二乘准则进行解算，用矩阵形式表示为：

$V=B\hat{x}-l---\left(15\right)$

$B=\left[\begin{array}{ccccccc}1& 0& 0& X_{beijing54}& 0& +{\mathrm{\ϵ}}_z& -{\mathrm{\ϵ}}_y\\ 0& 1& 0& Y_{beijing54}& -{\mathrm{\ϵ}}_z& 0& +{\mathrm{\ϵ}}_x\\ 0& 0& 1& Z_{beijing54}& {\mathrm{\ϵ}}_y& -{\mathrm{\ϵ}}_x& 0\end{array}\right]---\left(16\right)$

$l={\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS84}-{\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{beijing54}---\left(17\right)$

其中V为观测值改正数向量，为待估参数残差值，采用最小二乘迭代方法进行参数残差求解，

$\hat{x}={\left(B^{T}B\right)}^{-1}{B}^{T}l---\left(18\right)$

解得七参数，除了本身解算时的精度评定，还可以采用多余的已知控制点对七参数准确性进行检核。根据解得的七参数，就可以求出每个外部控制点的WGS84空间直角坐标。

最后外部控制点的WGS84空间直角坐标转为WGS84大地经纬度坐标，即(X,Y,Z)_WGS84→(B,L,H)_WGS84。

步骤3.2，采用双线性插值提取外部控制点所在DEM高程值；

外部控制点在DEM格网内对应P点位置，Q12，Q22，Q11和Q21为P点四周的DEM格网点。欲得到P点的高程值，采用P点四周的四个格网点的高程值进行双线性插值。首先在x轴方向上，对R1和R2两个点进行插值。然后根据R1和R2对P点进行插值。

一般来讲，欲求函数f在点P＝(x,y)的值，假设我们已知函数f在Q₁₁＝(x₁,y₁)，Q₁₂＝(x₁,y₂)，Q₂₁＝(x₂,y₁)及Q₂₂＝(x₂,y₂)四个点的值，则首先在x方向进行线性插值，得到

$\begin{array}{c}f\left(R_1\right)\≈\frac{x_2-x}{x_2-x_1}f\left(Q_{11}\right)+\frac{x-x_1}{x_2-x_1}f\left(Q_{21}\right),R_1=(x,y_1)\\ f\left(R_2\right)\≈\frac{x_2-x}{x_2-x_1}f\left(Q_{12}\right)+\frac{x-x_1}{x_2-x_1}f\left(Q_{22}\right),R_2=(x,y_2)\end{array}---\left(19\right)$

然后在y方向进行线性插值，得到

$f\left(P\right)\≈\frac{y_2-y}{y_2-y_1}f\left(R_1\right)+\frac{y-y_1}{y_2-y_1}f\left(R_2\right)---\left(20\right)$

这样就得到所要的结果f(x,y)，

$\begin{array}{c}f(x,y)\≈\frac{f\left(Q_{11}\right)}{(x_2-x_1)(y_2-y_1)}(x_2-x)(y-y_1)+\frac{f\left(Q_{21}\right)}{(x_2-x_1)(y_2-y_1)}(x-x_1)(y_2-y)\\ +\frac{f\left(Q_{12}\right)}{(x_2-x_1)(y_2-y_1)}(x_2-x)(y-y_1)+\frac{f\left(Q_{22}\right)}{(x_2-x_1)(y_2-y_1)}(x-x_1)(y-y_1)\end{array}---\left(21\right)$

步骤3.3，将经融合的DEM与数个外部控制点进行统计配准，所述配准模型包括

其中，β为控制点位置的局部坡向，α为控制点位置的局部坡度，为所有控制点位置坡度的平均值，a为DEM平面偏移矢量的模，b为DEM偏移矢量的方向角，ΔZ为DEM垂直偏移，Δh为DEM任一点的高程值与控制点高程值的差。

上述配准模型的推导如下：

假设InSAR生成的DEM的平面偏移矢量，AZ_S为平面偏移矢量的方位向，则平面偏移矢量大小在坡向为β的山坡处的投影S_xy的大小为(如图2所示)：

S_xy＝|S|cos(AZ_S-β) (22)

DEM某点位的高程与真实高程之间的差Δh为(如图3所示)：

Δh＝|S_xy|tan_α+|S_z| (23)

α为局部坡度，S_z为DEM在垂直向上的偏移矢量，上式可以进一步转为

$\frac{\mathrm{\Δ}h}{tan\mathrm{\α}}=\left|S\right|cos({\mathrm{AZ}}_S-\mathrm{\β})+\frac{\left|S_z\right|}{tan\mathrm{\α}}---\left(24\right)$

将上式简化为

$\frac{\mathrm{\Δ}h}{tan\mathrm{\α}}=acos(b-\mathrm{\β})+c---\left(25\right)$

c等于则公式(24)成为一个干净的余弦模型，Δh，α，β为观测值，a、b、c为待估模型参数。

在所有控制点中选取一部分在坡度、坡向和空间上较为均匀分布的控制点作为参与模型求解的解算点，剩余的点作为误差改进效果的检查点。

采用解算控制点高程差、坡度和坡向作为观测值来进行最小二乘平差求解配准模型参数，然后根据配准模型参数对DEM进行平移达到配准目的，完成一次统计配准后求取在检查控制点位高程差值的标准方差改善幅度。

步骤3.4，对步骤3.1至3.3的配准过程进行迭代，直至高程差标准差改善率低于2％。

因为步骤3.3的统计配准中用替代tanα会引入跟局部坡度相关的误差。坡度分布在0～90之内，所以坡度小的区域配准改善度偏小，而坡度大的区域配准改善度过大。需用一个与局部坡度相关的正切模型来消除这个引入的误差：

Δh＝a₇tan(α)+b₇ (26)

采用解算控制点位的坡度值与高程差作为观测值，通过最小二乘平差方法来求解待估参数a₇和b₇，然后将(26)所述的模型求取的误差趋势值与经过步骤3.3改正后的DEM进行相减，得到进一步改正后的DEM。

步骤4，校正DEM与坡度、平面位置和坡向相关的系统误差。

步骤4.1，对DEM结果与平面位置相关的误差进行校正，通过二次曲面拟合来完成。

在InSAR差分干涉处理过程中的基线估计往往会有误差，基线残差在DEM产品的体现为平面趋势。因此可以通过对DEM结果进行二次曲面拟合来消除这部分与平面位置相关的误差，二次曲面拟合模型包括

Δh_(x,y)＝a₁x²+a₂y²+a₃xy+a₄x+a₅y+a₆ (27)

其中x和y表示某个像元的横坐标和纵坐标，Δh_(x,y)为该点的DEM值与控制点高程差，a₁……a₆为待估模型参数。采用解算控制点位的位置与高程差作为观测值，通过最小二乘平差方法来求解待估参数a₁……a₆，然后将公式(27)所述的模型求取的误差趋势值与经过步骤3改正后DEM进行相减，得到进一步改正后的DEM。

步骤4.2，对DEM结果与坡度相关系统误差进行校正，通过二次或三次多项式模型来完成。

由于SAR斜距成像特点，使得SAR影像不可避免的出现透视收缩。不同平台下得到的DEM融合可以消除一部分透视收缩影响，但仍有可能具有不可忽视的残留误差。因为透视收缩与局部坡度大小有直接关系，所以设定一个二次或三次多项式模型来消除，具体选择哪个模型视高程差的标准差改善效果而定，

二次多项式为

Δh_(slope)＝b₁α²+b₂α+b₃ (28)

三次多项式为

Δh_(slope)＝b₁α³+b₂α²+b₃α+b₄ (29)

其中α表示某个像元点的坡度值，Δh_(slope)为该点的DEM值与控制点高程差b₁……b₄为待估模型参数，

采用解算控制点位的坡度值与高程差作为观测值，通过最小二乘平差方法来求解待估参数b₁……b₄，然后将公式(28)或(29)所述的模型求取的误差趋势值与经过平面位置误差改正后的DEM进行相减，得到进一步改正后的DEM。

步骤4.3，对DEM结果与坡向相关系统误差进行校正，通过正弦模型来完成的。

除局部坡度以外，SAR影像透视收缩还与局部坡向有直接关系，所以设定一个正弦模型来消除与坡向相关的系统误差，所述正弦模型为

Δh_aspect＝A₁sin(ωβ+φ)+C₁ (30)

其中β为某个像元点的坡向值，Δh_aspect为该点的DEM值与控制点高程差，A₁、ω、φ和C₁为待估模型参数。为方便最小二乘平差解算，该正弦模型可进一步变换为一阶傅立叶模型

Δh_aspect＝A₁sinωβ+B₁cosωβ+C₁ (31)

其中β为某个像元点的坡向值，Δh_aspect为该点的DEM值与控制点高程差，A₁、B₁、ω和C₁为待估模型参数，

采用解算控制点位的坡向值与高程差作为观测值，通过最小二乘平差方法来求解待估参数A₁、B₁、ω和C₁，然后将公式(31)所述的模型求取的误差趋势与经过坡度误差改正后的DEM进行相减，得到进一步改正后的DEM。

本实施例所述方法充分融合了InSAR监测手段分辨率高、相对精度高、覆盖范围大和地面控制点绝对精度高的特点，既解决InSAR获取DEM绝对精度不高，又克服了传统DEM获取方法例如星载和机载光学及激光测高手段获取高分辨测图效费比低和受云层影响，地面测图费时费力且无法得到整个野外山区高分辨地貌等缺点。整个流程结构清晰，具有实现简单、费用低、监测精度高、监测范围大、自动化程度高等优点。通过实地测试，经本流程得到的地形复杂区域DEM绝对精度可达到2m左右，其中72％的点精度在1.5m以内(如图4所示)。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (10)

1.一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，包括：

步骤1，通过双基InSAR差分干涉处理单轨影像对，获取单轨DEM；

步骤2，采取不同轨道的影像对分别获取DEM，并进行逐点融合；

步骤3，将经融合的DEM与地面控制点进行配准；

步骤4，校正DEM与坡度、平面位置和坡向相关的系统误差。

2.根据权利要求1所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤1包括：

步骤1.1，根据DEM需求区域坐标范围定制影像对，所述影像对包括主影像和辅影像；

步骤1.2，将所述主影像和所述辅影像共轭相乘生成初始干涉图；

步骤1.3通过轨道数据估计基线，根据基线由外部DEM模拟包含平地相位的地形相位图；

步骤1.4将所述地形相位图与所述初始干涉相位图相减得到差分干涉图；

步骤1.5对差分干涉图进行相位滤波和相位解缠，将解缠相位转化为高程差；

步骤1.6将高程差与外部DEM叠加得到新DEM；

步骤1.7通过地理编码将所述新DEM从SAR坐标系转换到外部DEM的地理坐标系得到所述单轨DEM。

3.根据权利要求2所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，还包括：将所述步骤1.7所得所述单轨DEM作为所述外部DEM，对步骤1.1至步骤1.7进行迭代，直至获取原始影像像元大小的DEM作为所述单轨DEM。

4.根据权利要求2所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤1.5对差分干涉图进行相位滤波采用的滤波方式包括Nonlocal均值滤波方式，滤波模型包括：

$u\left(s\right)=\frac{\underset{t}{\Σ}w(s,t)\mathrm{\φ}\left(t\right)}{\underset{t}{\Σ}w(s,t)}$

s为滤波窗口内的中心像元，t为滤波窗口内的其他任意像元，真实的干涉相位为u，含噪声的原始干涉相位为φ，相位估计值u(s)为滤波结果，w(s,t)为由像元t的邻域像元与像元s的邻域像元之间的SAR影像强度和相干性的综合相似度而决定的权重。

5.根据权利要求1所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤2的具体步骤包括：

步骤21，将所述不同轨道的影像对统一配准到相同坐标系之下，所述不同轨道的影像对包括升轨和降轨的影像对；

步骤22，采用融合模型进行融合，所述融合模型包括其中，H_fusion为融合后高程，B_a、W_a和H_a分别为升轨的垂直基线值，叠影或阴影掩膜值及高程值，B_d、W_d和H_d分别为降轨的垂直基线值，叠影或阴影掩膜值及高程值。

6.根据权利要求1所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤3包括：

步骤3.1，将外部控制点的空间直角坐标转为大地经纬度坐标；

步骤3.2，采用双线性插值提取外部控制点所在DEM高程值；

步骤3.3，将经融合的DEM与数个外部控制点进行统计配准，所述配准模型包括

其中，β为控制点位置的局部坡向，α为控制点位置的局部坡度，为所有控制点位置坡度的平均值，a为DEM平面偏移矢量的模，b为DEM偏移矢量的方向角，ΔZ为DEM垂直偏移，Δh为DEM任一点点的高程值与控制点高程值的差。

7.根据权利要求1或6所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于：对所述配准过程进行迭代，直至高程差标准差改善率低于2％。

8.根据权利要求1所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤4对经步骤三改正后DEM与平面位置相关的误差进行校正，是通过二次曲面拟合来完成的，二次曲面拟合模型包括

Δh_(x,y)＝a₁x²+a₂y²+a₃xy+a₄x+a₅y+a₆

其中x和y表示任一个像元的横坐标和纵坐标，Δh_(x,y)为所述像元点的DEM值与控制点高程差，a₁……a₆为待估模型参数，采用解算控制点位的位置与高程差作为观测值，通过最小二乘平差方法来求解待估参数a₁……a₆，然后将模型求取的误差趋势值与经步骤3改正后DEM进行相减,得到进一步改正后的DEM。

9.根据权利要求1所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤4对DEM结果与坡度相关系统误差进行校正是通过二次或三次多项式模型来完成的，

二次或三次多项式包括Δh_(slope)＝b₁α²+b₂α+b₃或Δh_(slope)＝b₁α³+b₂α²+b₃α+b₄，其中α表示某个像元点的坡度值，Δh_(slope)为该点的DEM值与控制点高程差b₁……b₄为待估模型参数，

采用解算控制点位的坡度值与高程差作为观测值，通过最小二乘平差方法来求解待估参数b₁……b₄，然后将模型求取的误差趋势与经过平面位置误差改正后的DEM进行相减，得到进一步改正后的DEM。

10.根据权利要求1所述的一种结合外控点的InSAR高精度高分辨率DEM获取方法，其特征在于，所述步骤4对DEM结果与坡向相关系统误差进行校正是通过正弦模型来完成的，所述正弦模型为Δh_aspect＝A₁sin(ωβ+φ)+C₁

其中β为某个像元点的坡向值，Δh_aspect为该点的DEM值与控制点高程差，A₁、ω、φ和C₁为待估模型参数。

采用解算控制点位的坡向值与高程差作为观测值，通过最小二乘平差方法来求解待估参数A₁、ω、φ和C₁，然后将模型求取的误差趋势与经过坡度误差改正后的DEM进行相减，得到进一步改正后的DEM。