CN102183761A - 星载干涉合成孔径雷达数字高程模型重建方法 - Google Patents

Abstract

本发明一种提供干涉合成孔径雷达数字高程模型重建方法。技术方案是：首先，对原始绝对干涉相位图进行粗采样。然后，拟合出粗采样绝对干涉相位图每一像素处数字高程模型重建三维坐标与绝对干涉相位的多项式。最后，通过对粗采样绝对干涉相位图每一像素处的多项式进行双线性内插，得到原始绝对干涉相位图每一像素处的多项式，将该像素处的绝对干涉相位值代入其对应的多项式即可得到相应的三维坐标。本发明计算量小，能够在现有计算能力下在较短的时间内完成；同时本方法的计算精度能够满足干涉合成孔径雷达系统数字高程模型重建的要求。

Description

星载干涉合成孔径雷达数字高程模型重建方法

技术领域

本发明属于微波遥感和信号处理的交叉技术领域，特别涉及一种干涉合成孔径雷达数字高程模型重建方法。

背景技术

干涉合成孔径雷达能够全天候、全天时地获取具有较高空间分辨率和高程精度的数字高程模型。利用数字高程模型，可监测地球陆地表面和冰雪表面的变化情况，进行地震、火山爆发、滑坡和大洪水等灾害预测，为农业、林业、渔业等生产提供帮助，为军事活动提供信息支持。总之，利用干涉合成孔径雷达获得的数字高程模型具有重要的应用价值和广泛的应用领域。

早期的干涉合成孔径雷达主要用于机载系统，在应用范围和应用空间上有较大的局限性。20世纪90年代随着技术的进步，干涉合成孔径雷达成功应用于星载系统，为全球提供了丰富的星载干涉合成孔径雷达数据。2000年2月11日美国国家航空航天局和国家影像与测绘局联合进行了航天飞机雷达地形测绘任务，这次基于干涉合成孔径雷达技术的试验首次实现了真正意义上的全球地形三维测绘，试验的部分数据和结果已经公布以供各科研机构处理和评估。2007年，德国发射了分辨率达1米的X波段合成孔径雷达卫星TerraSAR-X，这是迄今为止分辨率最高的民用合成孔径雷达卫星。德国于2010年6月21日又发射了第二颗几乎一样的卫星，两颗合成孔径雷达卫星将组成双星编队的TanDEM-X系统进行全球地形测绘，获取全球高精度数字高程模型数据，公布的预期指标将高于美国的航天飞机雷达地形测绘任务。目前德国已经将部分的TerraSAR-X的重复轨道干涉合成孔径雷达数据在互联网公布供全球研究机构使用。

数字高程模型重建是干涉合成孔径雷达数据处理过程中的关键步骤。干涉合成孔径雷达系统获取的原始回波数据经成像、配准、干涉、滤波、解缠等处理后得到绝对干涉相位，数字高程模型重建即利用该绝对干涉相位及雷达系统的轨道参数来反演雷达图像每一分辨单元所对应的地球表面三维坐标。

数字高程模型重建的耗时与雷达图像像素的数目成正比。星载干涉合成孔径雷达系统的一个重要任务是进行全球测绘，获取全球高精度数字高程模型产品。德国空间局计划利用TanDEM-X系统在三年内完成全球数字高程模型测量，预计获取的数据量将大于1300TB(2000年美国航天飞机雷达地形测绘任务的全球数据量约9.8TB)。面对如此海量实测数据，在保证精度的情况下，如何提高数据处理速度是一个亟待解决的重要问题。因此研究基于星载干涉合成孔径雷达系统的高精度、快速数字高程模型重建方法具有重要意义。

在干涉合成孔径雷达数字高程模型重建方法的研究中，Goblirsch于1997年提出了一种正侧视、平地假设下的数字高程模型重建方法，该方法重建得到的三维坐标误差较大；Giovanni Nico于2002年基于距离球、多普勒锥和相位双曲面推导了数字高程模型重建的解析解，该解析解的推导是基于两幅干涉合成孔径雷达图像的发射斜距重合或者重复轨道干涉测量模式的假设；Marcus于1998年研究了如何基于产品坐标格网快速获取高程的方法，该方法不能获得雷达图像每一分辨单元的三维坐标；在2005年，Eineder基于最大似然估计方法来同时进行相位解缠和定位处理，并且融合不同观测几何情况下的干涉相位图得到数字高程模型。上述数字高程模型重建方法在处理海量实测数据时往往难以克服处理速度与重建精度的矛盾。

发明内容

本发明要解决的技术问题是提供干涉合成孔径雷达数字高程模型重建方法，本方法计算量小，能够在现有计算能力下在较短的时间内完成；同时本方法的计算精度能够满足干涉合成孔径雷达系统数字高程模型重建的要求。

本发明技术方案的思路是：首先，对原始绝对干涉相位图进行粗采样。然后，拟合出粗采样绝对干涉相位图每一像素处数字高程模型重建三维坐标与绝对干涉相位的多项式。最后，通过对粗采样绝对干涉相位图每一像素处的多项式进行双线性内插，得到原始绝对干涉相位图每一像素处的多项式，将该像素处的绝对干涉相位值代入其对应的多项式即可得到相应的三维坐标。

本发明技术方案是：

利用已知的原始绝对干涉相位图及雷达系统的轨道参数，完成以下步骤：

第一步：衍生绝对干涉相位偏差图

对原始绝对干涉相位图按照方位向间隔k个像素、距离向间隔l个像素进行粗采样(k、l的取值由最终希望得到的数字高程模型精度要求决定，精度越高取值越小，可取[3，20]区间某一整数值)，得到粗采样绝对干涉相位图，将n个相位偏差分别叠加至粗采样绝对干涉相位图，衍生出n幅绝对干涉相位偏差图，其中n的取值也由最终希望得到的数字高程模型精度要求决定，精度越高取值越大，可取[4，10]区间内某一整数值。当n为奇数时，相位偏差Vφ_i，i＝1，2，Ln的取值如下：

$\left\{\begin{array}{cc}{\mathrm{V\φ}}_i=-i\·\frac{C\·2\mathrm{\π}\sqrt{{(f_a\·k)}^2+{(f_r\·l)}^2}}{n-1}& i=1,...\frac{n-1}{2},\\ V{\mathrm{\φ}}_i=0& i=\frac{n+1}{2},\\ V{\mathrm{\φ}}_i=-V{\mathrm{\φ}}_{i-\frac{n+1}{2}}& i=\frac{n+1}{2}+1,...n.\end{array}\right.$

当n为偶数时，相位偏差Vφ_i，i＝1，2，Ln的取值如下：

$\left\{\begin{array}{cc}{\mathrm{V\φ}}_i=-(2i-1)\·\frac{C\·2\mathrm{\π}\sqrt{{(f_a\·k)}^2+{(f_r\·l)}^2}}{2(n-1)}& i=1,...\frac{n}{2},\\ V{\mathrm{\φ}}_i=-V{\mathrm{\φ}}_{i-\frac{n}{2}}& i=\frac{n}{2}+1,...n.\end{array}\right.$

其中，C＝1.5，f_a、f_r分别为原始绝对干涉相位图的方位向频率和距离向频率。

第二步：绝对干涉相位偏差图的数字高程模型重建

基于每一幅绝对干涉相位偏差图，结合干涉合成孔径雷达系统的斜距方程、多普勒方程和干涉相位方程进行数字高程模型重建，得到n幅数字高程模型重建结果。数字高程模型重建具体方法参见硕士论文《天基InSAR理想干涉量的仿真与应用研究》(国防科技大学，2009年，王青松，第21页)2.3.2节。

第三步：多项式拟合

基于前述n幅绝对干涉相位偏差图以及相应的n幅数字高程模型重建结果，在粗采样绝对干涉相位图每一像素位置(a_s，r_s)处拟合出三个以绝对干涉相位值φ为自变量的m阶多项式其中m∈[1，n-1]，根据需要的计算精度选取，多项式

表示在像素位置(a_s，r_s)处重建的数字高程模型中的x轴坐标取值与绝对干涉相位值φ的关系，

表示在像素位置(a_s，r_s)处重建的数字高程模型中的y轴坐标取值与绝对干涉相位值φ的关系，

表示在像素位置(a_s，r_s)处重建的数字高程模型中的z轴坐标取值与绝对干涉相位值φ的关系。

第四步：多项式内插

对多项式

进行双线性内插可以求得原始绝对干涉相位图每一像素位置(a，r)处的三个多项式p′_x(a，r)(φ)，p′_y(a，r)(φ)，p′_z(a，r)(φ)：

$p_{x(a,r)}^{\′}\left(\mathrm{\φ}\right)=\mathrm{bilinear}\left(\begin{array}{c}\left[\begin{array}{cc}{p}_{x(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{x(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\\ p_{x(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{x(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\end{array}\right]\end{array}\right)$

$p_{y(a,r)}^{\′}\left(\mathrm{\φ}\right)=\mathrm{bilinear}\left(\begin{array}{c}\left[\begin{array}{cc}{p}_{y(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{y(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\\ p_{y(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{y(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\end{array}\right]\end{array}\right)$

$p_{z(a,r)}^{\′}\left(\mathrm{\φ}\right)=\mathrm{bilinear}\left(\begin{array}{c}\left[\begin{array}{cc}{p}_{z(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{z(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\\ p_{z(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{z(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\end{array}\right]\end{array}\right)$

其中，bilinear()表示双线性内插操作，floor()表示向小取整操作。

第五步：数字高程模型计算

将原始绝对干涉相位图每一像素位置(a，r)处的原始绝对干涉相位值φ(a，r)代入多项式p′_x(a，r)(φ)，p′_y(a，r)(φ)，p′_z(a，r)(φ)，即可得到该像素位置所对应的地面目标点三维坐标(x，y，z)，即最终得到的数字高程模型重建结果。

$\left\{\begin{array}{c}x=p_{x(a,r)}^{\′}\left(\mathrm{\φ}(a,r)\right)\\ y=p_{y(a,r)}^{\′}\left(\mathrm{\φ}(a,r)\right)\\ z=p_{z(a,r)}^{\′}\left(\mathrm{\φ}(a,r)\right)\end{array}\right.$

采用本发明可取得以下技术效果：

本发明的数字高程模型重建方法中，衍生出的多幅绝对干涉相位偏差图的数字高程模型重建过程是独立的，可以分别完成，因此本发明的数字高程模型重建方法十分有利于进行并行计算，可以充分发挥集群计算机的优势，设计并行程序大幅提高程序运行效率。

本发明的数字高程模型重建方法利用多项式拟合极大地减少了重复迭代求解三维非线性方程组的次数，有效降低了计算量。因此，采用本发明可以高效地完成数字高程模型重建，实现在现有计算能力下，快速完成数字高程模型重建，同时精度能够满足干涉合成孔径雷达系统测绘任务的要求。

附图(表)说明

图1本发明提供的干涉合成孔径雷达数字高程模型重建方法流程示意图；

图2为实际数据处理中使用的原始绝对干涉相位图；

图3为传统逐点重建方法得到的数字高程模型图；

图4为利用本发明得到的数字高程模型图。

具体实施方式

图1为本发明提供的干涉合成孔径雷达数字高程模型重建方法流程示意图。整个流程分为五步。第一步，衍生绝对干涉相位偏差图；第二步，对衍生的绝对干涉相位偏差图进行数字高程模型重建；第三步，多项式拟合，拟合出粗采样绝对干涉相位图每一像素处数字高程模型重建三维坐标与绝对干涉相位值的多项式；第四步，多项式内插，通过对拟合得到的多项式进行双线性内插可以求得原始绝对干涉相位图每一像素处的多项式；第五步，数字高程模型计算，将原始绝对干涉相位图每一像素处的原始绝对干涉相位值代入内插得到的多项式中即可得到最终的数字高程模型重建结果。

图2为数据处理实验中输入的原始绝对干涉相位图。该数据的点数为4381×3386，本幅图描述的场景范围约9000米×7000米。

图3为传统逐点重建方法得到的数字高程模型图，处理时间为2149.016秒。

图4为利用本发明的具体实施方式在k＝8，l＝8，n＝4，m＝3情况下得到的数字高程模型图，处理时间为231.012秒。相对于图3的重建结果，图4在x轴的坐标误差均方根值是0.03745米，在y轴的坐标误差均方根值是0.18804米，在z轴的坐标误差均方根值是0.17574米，相对于原始场景，重建精度能够满足干涉合成孔径雷达系统测绘任务的要求。

Claims (1)

1.一种星载干涉合成孔径雷达数字高程模型重建方法，其特征在于，利用已知的原始绝对干涉相位图及雷达系统的轨道参数，完成以下步骤：

第一步：衍生绝对干涉相位偏差图；

对原始绝对干涉相位图按照方位向间隔k个像素、距离向间隔l个像素进行粗采样，k、l的取值由最终希望得到的数字高程模型精度要求决定，精度越高取值越小，得到粗采样绝对干涉相位图；将n个相位偏差分别叠加至粗采样绝对干涉相位图，衍生出n幅绝对干涉相位偏差图，其中n的取值也由最终希望得到的数字高程模型精度要求决定，精度越高取值越大，可取[4，10]区间内某一整数值。当n为奇数时，相位偏差Vφ_i，i＝1，2，Ln的取值如下：

当n为偶数时，相位偏差Vφ_i，i＝1，2，Ln的取值如下：

其中，C＝1.5，f_a、f_r分别为原始绝对干涉相位图的方位向频率和距离向频率。

第二步：绝对干涉相位偏差图的数字高程模型重建

基于每一幅绝对干涉相位偏差图，结合干涉合成孔径雷达系统的斜距方程、多普勒方程和干涉相位方程进行数字高程模型重建，得到n幅数字高程模型重建结果。 

第三步：多项式拟合

基于前述n幅绝对干涉相位偏差图以及相应的n幅数字高程模型重建结果，在粗采样绝对干涉相位图每一像素位置(a_s，r_s)处拟合出三个以绝对干涉相位值φ为自变量的m阶多项式 

其中m∈[1，n-1]，根据需要的计算精度选取，多项式 

表示在像素位置(a_s，r_s)处重建的数字高程模型中的x轴坐标取值与绝对干涉相位值φ的关系， 

表示在像素位置(a_s，r_s)处重建的数字高程模型中的y轴坐标取值与绝对干涉相位值φ的关系， 

表示在像素位置(a_s，r_s)处重建的数字高程模型中的z轴坐标取值与绝对干涉相位值φ的关系。

第四步：多项式内插

对多项式 

进行双线性内插可以求得原始绝对干涉相位图每一像素位置(a，r)处的三个多项式p′_x(a，r)(φ)，p′_y(a，r)(φ)，p′_z(a，r)(φ)：

其中，bilinear()表示双线性内插操作，floor()表示向小取整操作。 

第五步：数字高程模型计算

将原始绝对干涉相位图每一像素位置(a，r)处的原始绝对干涉相位值φ(a，r)代入多项式p′_x(a，r)(φ)，p′_y(a，r)(φ)，p′_z(a，r)(φ)，即可得到该像素位置所对应的地面目标点三维坐标(x，y，z)，即最终得到的数字高程模型重建结果。

k、l的取值由最终希望得到的数字高程模型精度要求决定，精度越高取值越小，可取[3，20]区间某一整数值

n的取值也由最终希望得到的数字高程模型精度要求决定，精度越高取值越大，可取[4，10]区间内某一整数值。 

Images