基于Hermite的SAR卫星轨道插值干涉测量应用研究

马 飞,隋立春,姚顽强

(1.长治学院 计算机系,山西 长治 046011; 2.长安大学 地质与测绘工程学院,西安 710054;3.西安科技大学 测绘科学与技术学院,西安 710054)

0 引言

在合成孔径雷达(SAR)卫星技术的应用中，卫星轨道参数是合成孔径雷达干涉测量(InSAR, interferometric synthetic aperture radar)数据处理过程中重要参数之一，可直接影响主辅SAR图像配准、DEM测高、基线估算等精度。在干涉测量数据处理过程中，最初主影像和辅影像粗配准过程中会引入卫星轨道参数作为影像初始偏移量，去除平地相位过程中会利用卫星轨道参数模拟计算平地相位，生成DEM过程中会利用卫星轨道参数估算基线长度，因此卫星轨道参数是干涉合成孔径雷达的基础，它直接决定着合成孔径雷达干涉测量系统的特性优劣[1]。

大部分SAR卫星运营商都能提供卫星轨道精密参数，例如欧空局发射的ERS卫星、ENVISAT卫星由荷兰代夫特对地空间研究中心(DEOS， delft institute earth—oriented space research)提供的精密轨道数据精度优于6 cm[2]；哨兵卫星POD精密定轨星历数据定位精度优于5厘米[3]；TerraSAR X卫星搭载全球定位系统(GPS)接收机和激光精密定轨反射器，轨道确定精度达到厘米级[4]。因此关于SAR卫星轨道插值的文献研究较少，很多学者专注于研究InSAR数学模型、轨道精度对InSAR测量结果的影响等，分析InSAR干涉测量的数学模型，包括轨道精度、距离向频谱、干涉临界基线距、模糊高度、差分相位对形变的敏感度等问题，但未就如何改进轨道精度提出解决方案。针对卫星轨道插值问题更多的文献专注于研究全球定位系统(GPS, global positioning system)、 北斗导航卫星定位系统(BDS, beiDou navigation satellite system)和伽利略(Galileo)卫星定位系统的轨道插值问题，插值方法包括埃尔米特插值、切比雪夫插值、三次样条插值、拉格朗日插值算法等。

针对未提供精密轨道文件的SAR卫星数据，严重制约着该卫星影像数据处理精度。例如加拿大的RADARSAT卫星，日本的JERS卫星、ALOS卫星和ALOS-2卫星等等，其中ALOS卫星搭载的PALSAR传感器扫描地面影像的持续时间大概为13秒，在影像的参数文件中提供了卫星从生成第一个像素到最后一个像素期间的所有卫星状态矢量，包括卫星三维位置、三维移动速率等数据，其状态矢量采样间隔为60秒；因此，主辅SAR影像在进行干涉测量时只有一个卫星状态矢量参与了计算。这样在利用轨道参数进行去除椭球参考相位、像素点匹配、去除平地相位时都会遇到困难，干涉相位图中会有残余条纹，造成最终获取的DEM或差分干涉图中混有系统误差。因此本文拟采用编程工具对卫星轨道数据进行模拟，利用埃尔米特插值法对其轨道参数进行等距插值，从而估计其轨道参数，缩小卫星轨道状态矢量采样间隔，提高轨道精度。

1 埃尔米特插值和SAR干涉测量基本原理

1.1 插值算法的简介

插值算法是一种关于函数逼近的简单又实用的数学方法，主要用于某函数在有限个点处取值状况，估计其在其他点处的值。插值算法是数值分析中一种基本算法。根据被插值的函数自变量个数，插值法可分为一维插值、二维插值、和多维插值。根据选用的分段直线、多项式或样条函数插值函数，插值又可分线性、多项式和样条插值[5]。

目前较为常用的插值算法有Lagrange插值、牛顿插值、埃尔米特插值、三次样条插值、高维插值等。Lagrange插值、牛顿插值主要对函数进行直接插值，对函数的一阶或多阶导数没有限制，没有考虑插值点的导数值，即没有考虑卫星轨道的瞬时速度矢量，不适用于SAR卫星轨道插值；埃尔米特插值及三次样条插值不仅对原函数进行限制，根据限制条件的不同，可以对插值函数的倒数进行限制，如埃尔米特插值可以对函数的一阶导数进行限制，而三次样条插值法可以根据不同的条件，可以对一阶或二阶导数进行限制，从而可以保证插值函数的连续性，这两种插值方法既可以拟合SAR卫星轨道的位置矢量，也可以计算卫星在某一位置的瞬时速度矢量，因此适用于SAR卫星轨道插值。根据本次对卫星轨道参数插值的条件，即在已知点位和速度矢量在X、Y、Z方向的分量条件下，在时间上进行插值，故选择的插值方法有埃尔米特插值和三次样条插值法。武汉大学吴宏安教授[6]等人对InSAR数学模型研究发现，利用原始轨道参数文件进行埃尔米特插值，可以有效地加密轨道文件参数，利用新的加密后的轨道参数文件进行像素点配准、基线估算、椭球参考相位估算等，可以提高干涉测量精度，去除平地残余相位。本文利用埃尔米特插值法估算SAR卫星轨道参数，并进行相关的实验研究。

1.2 埃尔米特轨道参数的插值算法

卫星描述轨道模型的星历数据包括卫星的位置矢量(X,Y,Z)和瞬时的速度矢量(Vx,Vy,Vz)，用下式来表示：

(1)

(2)

其中:(X,Y,Z)分别为固定笛卡尔参考坐标系的三个坐标轴，固定笛卡尔参考坐标系统是以地心为原点，X、Y轴垂直位于赤道平面，X轴由地心指向格林威治子午线，Z轴垂直于赤道平面且符合右手定律，与地球椭球极轴重合指向北极。

卫星轨道参数文件中的状态矢量可以进行多项式拟合，将(X,Y,Z)和(Vx,Vy,Vz)作为时间T的函数，构建矩阵如下式：

(3)

式中,a,b,c分别为X,Y,Z维的多项式系数。

当求一个插值多项式H(x)，使其满足条件为：

H(xi)=yi,i=1,2,…,n

(4)

(5)

根据函数的插值条件，共有m+n+2个。所构造函数H(x)次数一般小于或等于m+n+1次，可由牛顿插值的思想来构造H(x)为：

H(x)=Nn(x)+Pm(x){(x-x0)(x-x1)…(x-xn)}

(6)

其中:Nn(X)为牛顿基本的差值多项式，Pm(X)为函数对应的m次多项式。根据式(6)和已知数据可以得到:

H(xi)=Nn(xi)=yi,i=1,2,…,n

(7)

在确定Pm(X)，需先对H(x)求导：

(8)

其中:ωn+1(x)=(x-x0)(x-x1)…(x-xn)，令x=xi,i=0，1，2，3，…,m，且把H′(xi)=y′,i=1,2,3,…,m代入式上式可得:

(9)

进一步可得Pm(X):

(10)

对于求Pm(X)转化成为已知Pm(x)的函数表，如表1所示。

表1 Pm(x)的数据表

当Pm(X)不超过m的插值多项式。Pm(x)满足式:

Pm(x)=Pm(x0)+Pm[x0,x1](x-x1)+…+

Pm[x0,…,xm](x-x0)…(x-xm-1)=

(11)

其中:Pm[x0,…,xm]表示为第m阶差商，形式为:

(12)

再将Nn(X)、Pm(x)带入式(6)，可得到插值函数H(x)。

本论文采用埃尔米特等距插值法，即要求在给定节点处函数值和一阶导数均相等[7]。满足下式：

yi=f(xi)

(13)

(14)

以上为埃尔米特插值算法基本原理，本文基于软件编程实现对SAR卫星轨道参数插值计算。

1.3 SAR干涉测量技术简介

InSAR技术测量是通过对地表目标点的两个视向的观测实现的，其工作方式包括距离向干涉测量(Across-track Interferometry)、方位向干涉测量(Along-track Interferometry)和重轨干涉测量(repeat-pass Interferometry)。

距离向干涉测量通过在飞行平台装载两个天线，一个天线用于发射并接收雷达波，另外一个仅用于接收雷达波。通过这种干涉测量方式，只要确定平台位置，就较容易地获得高质量的干涉测量数据及较高精度的DEM。当前这种干涉测量方式常用于航空平台上，即为前面所提到了机载InSAR。在采用这种测量方式，要求飞行平台较为稳定，而且需要飞机姿态的影响。

方位向干涉测量，该方式将两个天线安装在平台的前端和后端，主要适用于机载，已应用在运动目标的探测、海上波浪谱线等方面的测量工作。

本次用于研究的影像为基于卫星的重轨干涉测量获取的影像。重轨干涉测量，其作业方式仅需要一个天线，利用卫星在短时间回到大致相同的轨道上，获取同一范围内的主辅影像数据。该测量适合于几百公里上的卫星作业方式。一般而言，星载SAR要获得较好的影像数据，需满足以下3个条件：首先是在两次观测期间，地面没有发生形变或仅发生微小的形变；其次平台与地表目标间具有较稳定的观测几何关系，这需要搭载平台的稳定和平台的准确定位参数，在获取相同的影像，当卫星有精密的轨道参数时，处理后的影像精度更高；最后为SAR处理器能够对做运动的补偿信号保留住内在的相位信息。

数据预处理，该阶段主要从主辅影像轨道参数文件等生成轨道参数文件和单视复影像。当研究区域为全影像的一部分时，可以对原始影像进行裁剪，将得到裁剪区域用于后面的处理研究。

数据配准，该过程可分为去噪、控制点自动搜索、重采样与配准质量评价等，主要目的是将辅影像与主影像进行匹配，使得两景影像坐标系统一致。在处理影像过程中，需采用多级配准方案来提高影像的配准精度。一般分为粗配准、精配准及子像元配准。

干涉图生成，需进行相位干涉，即对两景影像的匹配点的数值做共轭相乘，得到地面任意点的干涉后影像图的数值。通常计算得到的相位差由干涉图的灰度表示，并以条纹形状表现。

重采样与配准质量评价。重采样指对对原始影像数据的实部和虚部做内插，插值方法有双线性和双三次样条方法。为了获得高精度的相位值，对重采样提出较高的要求。配准结果的评价是对配准结果的综合评价，对干涉图的生成有重要影响。质量评价的量化标准为相干系数，其取值范围[-1,1]，为配置质量评价量化标准。当相干系数越接近0，代表着两景数据相关程度越小；当相干系数的绝对值越接近1时，表示两景数据的相关程度越大。

去平地效应，当地面某两点高程不同，其相位差可能有两部分造成的，其一为平地两点的相位差；其二为两点间存在的高差而造成相位差。通过去平地效应，使得到干涉图的相位差只由地形引起，从而得到更为准确的地形数据。

相位解缠，即从干涉图得到的相位差，取值范围为(-π,π]，代表了各个相位差的主值，为了得到真实相位差需在这个主值的基础上减去2π或加上2π的整数倍。相位解缠的算法可分为两大类，一类为基于路径积分相位解缠方法(主要为枝切法)；另一类为基于最小范数的解缠算法。

地理编码，由干涉图得到高程值后，只是得到各点在SAR坐标系下的高程值。为了得到对应坐标系下的数字高程模型，需将SAR坐标系下的每一点高程值进行坐标转换，得到最终的地面DEM。

2 实验分析与研究

2.1 巴姆地区Envisat卫星实验结果验证

巴姆地区位于伊朗的东南部，中心位置为北纬29.01°，东经58.26°。本文利用覆盖该地区的Envisat卫星ASAR传感器影像和DEOS提供的精密轨道文件对埃尔米特轨道插值实验进行分析验证。根据雷达干涉测量原理[15]，对两幅Envisat卫星主辅影像进行干涉处理，先进行主辅影像配准，计算主辅影像上所有的像素点在参考托球面上的投影点，对所有配准的像素点进行相位干涉，获取各像素的相位值，该相位值中包含地形起伏引起的相位、地表形变引起的相位变化、椭球参考面相位、噪声相位等等。本次数据处理过程中像素点配准和椭球参考面相位都受到轨道文件精度的影响，因此干涉相位图的参考面相位的条纹规律性和相干性值是轨道参数文件质量定性和定量评判的重要标准。

如图1所示，图(a)为基于原始轨道文件进行干涉测量得到的干涉图，图(b)为基于埃尔米特插值轨道文件进行干涉测量得到的干涉图，图(c)为基于DEOS提供的精密轨道文件干涉测量得到的干涉图。

图1 基于不同轨道参数得到的干涉图

对比图(a),(b)两幅干涉图可知，两幅图中均出现了清晰且规律性很强的干涉条纹，其中心位置出现的蝴蝶状条纹是地震形变引起的相位变化，属于形变条纹；其他规律性条纹主要是参考面相位条纹；在图中左下角的蓝线区域的干涉条纹对比具有较大的差距，该区域地形起伏较大，属于地形相位条纹和参考面相位条纹叠加的结果。图(b)中经轨道插值后得到的干涉图干涉条纹更清晰，在地形复杂的地区仍具有较高的相干性，在形变区的干涉条纹具有较强的连贯性。分析图(b),(c)可知，由插值得到的干涉图和利用卫星精密轨道得到的干涉图有较强的一致性,尤其是在左下角地形起伏较大的区域，干涉条纹基本无差异，说明利用DEOS精密轨道文件和埃尔米特插值轨道文件进行干涉测量处理后相位无差异，埃尔米特插值轨道文件的精度与DEOS精密轨道文件精度相当。

通过对以上结果分析，可以定性地判断基于埃尔米特卫星轨道参数的插值算法，可以提高SAR影像的干涉测量效果，有效地提高轨道文件精度。

由于欧空局会提供Envisat卫星的DEOS精密轨道文件，所以使用埃尔米特插值对原始文件插值计算精密轨道文件的研究价值并不在于为Envisat卫星轨道插值，而在于为无精密轨道文件的SAR卫星插值轨道数据，例如日本ALOS卫星。因此本文再次使用埃尔米特插值算法对ALOS卫星轨道进行插值实验。

2.2 陕西地区ALOS卫星轨道插值结果

本文选取的ALOS/PALSAR数据为覆盖陕西咸阳彬县和长武县的彬长矿区的影像。该矿储煤丰富，质量优良，矿区的西北、北部以甘肃省为界，东西距离为46 km，南北长约为36 km，矿区地表覆盖面积近1 000 km2。该矿区位于黄土高原的塬梁沟壑地带，地表覆有黄土层。气候干燥，降水量少，年平均气温在10°，高差多达200 m。该地区的属于泾河流域，其中流经该矿区的最大支流为黑河。彬长矿区的地表植被相对较少，位于开采沉陷区域内的村落较少，有利于获取高质量的SAR影像。ALOS/PALSAR数据的原始参数文件中包含有11个轨道状态矢量，但由于其轨道矢量时间间隔为60秒，每景影像像素扫描时间为13秒，因此干涉测量过程中有可能无轨道状态矢量参与。本文采用埃尔米特插值法进行等距插值，设置时间间隔为10秒、5秒、2秒分别计算插值后的轨道状态矢量，在插值过程中以时间T为自变量，卫星的位置矢量作为函数值，速度矢量作为时间的一阶导数值。

具体实验数据影像参数如表2所示，利用编程工具对实验数据计算加密后的轨道参数，再将新的轨道参数文件替换原始文件进行干涉测量数据处理。

表2 轨道插值运算实验数据

在ALOS/PALSAR数据的原始参数文件中，包含有影像起始时间(start_time)、结束时间(end_time)、中点位置时间(center_time)，11个采样时间间隔为60秒的轨道矢量(X,Y,Z,Vx,Vy,Vz,T)，本文以这11个轨道状态矢量为已知点构建埃尔米特基函数，以第一个轨道状态矢量采样时间为起始时间，以最后一个轨道状态矢量采样时间为结束，分别进行10秒、5秒、2秒等间距插值，利用轨道矢量插值后干涉测量结果如图2所示。

图2 不同轨道参数采样间隔的地理编码影像干涉相位测量图

通过观察图2中的(a),(b),(c),(d)四幅图发现图中的干涉条纹基本无差异，通过分析InSAR数学模型[7-8]可知，日本ALOS卫星的PALSAR传感器使用L波段对地扫描，L波段波长相对C波段较长，干涉测量后的干涉条纹对地形和形变的敏感度变低，地形高程变化178米才会引起一个周期的相位条纹变化。因此在(a),(b),(c),(d)四幅图中的干涉条纹变化周期较图1中的ENVISAT卫星的C波段变化周期较少，四幅图中的干涉条纹肉眼基本观察不到差异。需要引入定量分析方法分析不同轨道采样间隔下干涉测量的效果又和不同。

图3 不同卫星轨道参数采样间隔的影像的相干图

为了进一步比较四种轨道参数采样间隔对影像相干系数的影响，埃尔米特插值法得到的卫星轨道参数采样间隔为60秒、10秒、5秒、2秒的相干图，选用编程工具软件分别对该相干图做灰度均值运算，求得每幅相干图的平均系数，如表3所示。

表3 不同采样间隔下影像相干图的平均相干系数表

通过表3可知当卫星轨道参数采样间隔缩小时，两幅SAR影像干涉得到平均相关系数比原始采样间隔60秒的影像相干图的平均相干系数均高，其代表对原始卫星轨道参数进行插值的方法对提高SAR影像的相干性具有可行性，并取得一定成效。再进一步分析，随着轨道参数采样间隔的不断减少，提高影像的相干性并非呈正相关，本次实验中，当采样间隔为5秒时，得到相干图的平均相干系数最高，而其采样间隔为2秒时，其相干图的平均相关系数略低。这说明随着轨道参数间隔的减小，其相干性呈先提高后降低的趋势，以此估计得到矿区地表点也呈该沉降趋势。本次研究中，原始轨道参数间隔为60秒，经插值算法采样间隔为10秒、5秒、2秒，初步判断当采样间隔为5秒时，得到相干性最好的差分干涉影像图。

3 结束语

1)本文基于埃尔米特插值多项式拟合SAR卫星轨道，对SAR卫星轨道模型的中的位置矢量和瞬时的速度矢量进行插值计算，可以有效地加密卫星轨道星历数据，提高卫星轨道参数文件精度，减少干涉测量过程中的系统误差。

2)通过对Envisat卫星原始轨道文件进行埃尔米特插值计算，得到加密后的卫星轨道星历数据。利用原始轨道文件、加密轨道文件、DEOS精密轨道文件分别进行干涉测量，定性研究发现，利用加密轨道文件和精轨文件得到的干涉图中干涉条纹清晰且一致，说明埃尔米特插值算法可以有效加密SAR卫星轨道星历数据，提高干涉测量精度。

3)通过对ALOS卫星原始粗轨进行埃尔米特插值计算，轨道星历数据采样间隔分别设置为10秒、5秒、2秒，结果表明进行轨道插值后得到的干涉图相干性比基于粗轨得到的干涉图相干性显著提升；随着轨道参数间隔的缩小，其相干性呈先提高后降低的趋势，采样间隔为5秒时相干性最高。