一种高精度星载SAR几何定标方法与流程

本发明属于微波遥感数据处理技术领域，尤其涉及一种高精度星载sar几何定标方法。

背景技术：

随着日本7m分辨率的alos‐palsar、加拿大lm分辨率的radarsat‐2、意大利lm分辨率的cosmo‐skymed、德国lm分辨率的军民两用卫星terrasar‐x和tandem‐x、日本1m分辨率的alos‐2、中国1m分辨率的高分三号等sar卫星的相继成功发射，标志着高分辨率星载sar时代的到来。随着sar遥感影像的分辨率提高，对其几何定位能力的需求也相应提高。由于sar卫星入轨后，环境变化、设备老化等因素严重影响其定位性能，因此，需要通过高精度的几何定标方法，提升sar系统的无控定位能力。

星载sar几何定标方法的核心是，利用控制点在sar影像上的像素坐标与利用严密成像几何模型反算的像素坐标存在偏移量，根据偏移量再反演星载sar系统的几何定标参数，进而实现提升星载sar系统的无控定位精度。sar信号的发射和接收都是在时间尺度上完成的，主要包括快时间(距离向)和慢时间(方位向)。这二维的时间误差主要影响sar影像在距离向和方位向的几何定位误差，是星载sar几何定位的主要误差源。经典的星载sar几何定标主要是针对方位向起始时间和距离向起始斜距两个参数进行标定。由于系统设备时间控制单元的误差引起，系统设备时间控制单元用于记录所有事件的系统日期，它有一定的时间精度及一定的记录频率，它们会对方位向的开始时间精度带来一定影响。由于雷达信号经过信号通道的各个器件时产生的系统时延，不同类型(时宽和带宽组合)雷达信号的系统时延也不同；另外，雷达信号穿越大气到达地面再从地面返回大气时的大气路径双向延迟。因此，雷达信号的时宽和带宽以及大气传播延迟均影响着几何定标精度，而经典的几何定标方法均为充分考虑这两个影响因素。

技术实现要素：

针对现有技术存在的不足，本发明顾及雷达信号的时宽和带宽特点和大气传播延迟影响，提出了一种高精度星载sar几何定标方法，包括如下步骤：

步骤1，根据sar系统设计的波位表，查找出与某一种时宽和带宽组合模式相同的波位号；

步骤2，通过对sar卫星系统的成像任务规划，以查找出的任意波位对几何定标场区域进行成像；

步骤3，远程调整几何定标场区域的所有自动角反射器指向，使角反射器的法线方向与入射波的方向保持一致，根据自动角反射器调整的方位角和俯仰角，精确计算成像时刻各个角反射器顶点的地理坐标，同时，对各个角反射器顶点的地理坐标进行固体潮改正补偿；

步骤4，对几何定标场区域的回波数据进行成像处理，获得sar影像和辅助参数文件；

步骤5，利用基于质心法的亚像元点提取方法，对每景sar影像中的角反射器点进行精确提取，获得角反射器点的像方坐标；

步骤6，针对每景sar影像，计算雷达信号在天线相位中心到各个角反射器点之间的大气传播延迟改正值，其中大气传播延迟包括中性大气延迟和电离层延迟；

步骤7，分别计算每景sar影像的几何定标参数，所述几何定标参数包括距离向起始时间改正值、方位向起始时间改正值、采样频率改正值、脉冲重复频率改正值；

步骤8，根据几何定标场区域多景sar影像的几何定标参数，取均值，获取该成像模式的几何定标参数。

进一步，还包括以下步骤，对几何定标后的几何精度验证场区域sar影像进行精度评价，

步骤9，针对步骤1中的时宽和带宽组合模式，对几何精度验证场区域进行成像任务规划、拍摄；

步骤10，根据步骤8计算出的几何定标参数，对几何精度验证场区域的sar回波数据进行补偿，并进行成像处理，获取sar影像和辅助参数文件；

步骤11，获取几何精度验证场的地面控制点的像方坐标，并重复进行步骤6的操作，获得每个地面控制点的大气传播延迟改正值；

步骤12，根据sar影像中地面控制点的像方坐标和大气传播延迟改正值，利用rd模型的直接定位方法，计算sar影像中每个控制点的物方坐标，并与每个地面控制点的真实坐标进行对比，所得的差值即为该点的定位误差；统计所有地面控制点定位误差的标准差，即得到sar系统在该成像模式下经几何定标后无控制点情况的几何定位精度。

进一步的，所述步骤3中根据自动角反射器调整的方位角和俯仰角，精确计算成像时刻各个角反射器顶点的地理坐标，同时，对各个角反射器顶点的地理坐标进行固体潮改正补偿，其实现方式如下，

假设以云台垂直方向的转动轴中心为坐标系原点，以真北方向为坐标系的x轴，建立左手直角坐标系o-xyz，设p为角反射器顶点，p'为p在xoy平面的投影，p'a和p'b分别为p'到y轴和x轴的距离，ω为点p的方位角，为点p的俯仰角；

根据高斯投影的方法，将地面控制点wgs‐84系下的大地坐标(lat0,lon0,h0)转换成高斯平面坐标(x0,y0,h0)，由此可知，坐标系原点坐标o为(x0,y0,h0+h)，点p的高斯平面坐标为其中，h为云台垂直方向的转动中心相对于地面控制点的高度，po是角反射器顶点与云台垂直方向的转动中心之间的距离；根据高斯投影关系，反算出点p在wgs‐84系下的大地坐标(latp,lonp,hp)，即获得高精度角反射器顶点的地理坐标；

采用以下公式对角反射器点的重力固体潮进行校正，

δg(t)＝δthg(t)-δfc(1)

式中，δth为潮汐因子，δfc是地心纬度的函数，即

式中，为地理纬度的函数，即cm为地心至月心的平均距离，rm为地心至月心的距离，cs为地心至日心的平均距离，rs为地心至日心的距离，zm为月球对地面点的地心天顶距，zs为太阳对地面点的地心天顶距。

进一步，所述步骤5的实现方式如下，

1)通过sar图像间接定位模型，由地面大地坐标计算角反射器对应的sar图像像点坐标，再由影像地面点坐标和轨道参数获得该地面点对应的方位向时间和角反射器点的星地距离后，解算该点对应的行列号；其中地面点坐标，由外业gps测量获得，轨道参数由星上gps获得，下传给地面使用；

2)在步骤1)中计算的行列号附近区域，寻找在此区域像素灰度值最大的亮点，再通过人工目视判读，比较周围地物，确定该点是否为角反射器在sar图像上的粗略位置，初步确定角反射器的粗位置；

3)以步骤2)中的粗略位置所对应像素点为中心，选取其邻域为目标质心寻找区，将该领域内所有像点对应的坐标和灰度强度信息按顺序提取出来，并把该灰度赋予每个像素的几何中心；

4)对该目标区内每一像素进行的细分，细分尺寸根据不同的精度要求划分，然后采用双线性插值方法或sinc函数插值方法求得每一细分点的灰度；

5)计算图像回波强度质心。

在计算质心时，采用逐行分析计算，取平均值法计算质心，计算方式如下：

计算每一行质心的坐标

其中，n为被插值点的个数，m为细分后点的个数，xk为被插值点的坐标，xu为细分后点的行坐标，ik为被插值点的灰度值，f(u)为细分后点的灰度值；

然后，以每一行质心坐标和对应灰度值为参数，求整个目标区的质心，

其中，为每一行质心坐标，yu为细分后点的列坐标，为细分后的点对应平均灰度值；计算得到的(xc,yc)即为目标区域的质心，也即角反射器点对应的像方坐标。

进一步的，步骤6中所述的中性大气延迟改正值的计算方式如下，

上式中地表压强psruf＝pd+e，pw为大气可降水量；md、mw分别为干、湿大气分子量，md＝28.9644kg/kmol，mw＝18.0152kg/kmol；rm为摩尔气体常量，rm＝8.31451j/(mol·k)；

k1、k2值与卫星载荷发射电磁波信号的波长λ相关，求解公式：

其中，地表压强和大气可降水量从全球大气模型中获得。

进一步的，步骤6中所述的电离层延迟改正值的计算方式如下，

设目标点t垂直方向上的全部电子都压缩于一点t'处，雷达信号在传播路径上的电离层延迟为，

式中k＝40.28m³/s²，θ'是压缩点t'处的入射角，tec的单位是1016个电子/米²，计算公式如下：

第一，在ionex文件中取得全球的tec含量图，然后以穿刺点的经纬度(lat,lon)为参数，从它附近的网格点做内插，具体的公式为，

式中，i0,0，i1,0，i1,1，i0,1为穿刺点相邻四角网格点的tec含量，lat0为左下角网格点的纬度，lon0为左下角网格点的精度，δlat为电离层网格点的纬度间隔，δlon为电离层网格点的经度间隔；

第二，得到空间内插后进行相应时间点t的时间双线性内插，具体的公式为，

式中，ti和ti+1为时间点t两个临近时间格网的时间；

第三，取得穿刺点对应时刻垂直方向上的tec后，根据公式(7)获得雷达信号在传播路径上的电离层延迟值。

进一步的，步骤7中每景sar影像的几何定标参数计算方式如下，

构建星载sar几何定标模型，

式中，tf、ts分别为距离向的快时间和方位向的慢时间，tf0、ts0分别为距离向起始时间的测量值和方位向起始时间测量值，tdelay为大气传播延迟时间，为大气传播延迟值除以光速，δtf、δts为系统时延误差，x、y为像素坐标，fs为采样频率，fp为脉冲重复频率；

对ts0进行补偿，

式中，ts'0为星上记录的方位向起始时间，n为雷达信号从发射到接收经历整周期的个数，tsample_delay为星上记录的采样时间延迟；

星载sar几何定标模型(10)可以表示成如下形式：

式(12)的误差方程为

v＝bx-l(13)

其中，

设初始值δtf＝0、δts＝0，fs和fp由星上下传的辅助参数文件获得，采用迭代运算的方法获取几何定标参数x＝[δtfδtsδfsδfp]^t，实现方式如下，

a)根据sar影像间接定位算法，通过角反射器点的地面坐标和轨道参数，可以计算出该地面角反射器点对应的sar影像方位向时间ts和该点的星地距离rst，星地距离rst即距离向时间tf；

b)根据公式(10)，利用几何定标场区域地面布设的n个角反射器点，组建n个方程组；通过基于谱修正法迭代法的最小二乘平差，精确计算方位向起始时间延迟δts、距离向起始时间延迟δtf、脉冲重复频率改正值δfp＝fp'-fp和采样频率改正值δfs＝fs'-fs；

c)更新四个几何定标参数，重新执行步骤a)，再次计算四个几何定标参数，判断两次计算的几何定标参数是否小于一个预先设定的极小值，若小于则迭代终止，否则转向a)继续迭代运算；

d)迭代运算结束后的结果为该景sar影像的几何定标参数。

与现有技术相比，本发明的有益效果：

现有的星载sar几何定标技术发展仍不成熟，尚无系统的定标方法和解决方案。本发明首先将影响斜距测量精度的主要影响因素进行了分类，主要分为sar系统时延误差和大气传播延迟误差。sar系统时延误差主要是雷达信号在sar载荷内部经过各个部件时产生的固有的时延误差，即从振荡器处产生到sar天线发射端、再从sar天线接收端到回波数据的接收采样，主要受雷达信号的带宽和脉冲宽度影响；大气传播延迟误差主要是雷达信号在传播路径上受到大气环境的折射等影响产生的时变的时延误差。由此，根据雷达信号的带宽和脉冲宽度组合，分组进行星载sar的几何定标，无需考虑成像模式、波位、左右侧视等因素的影响，大大降低了星载sar几何定标的工作量；考虑了时变的大气环境对斜距测量精度的影响，提高了定标精度。

附图说明

图1为高精度星载sar几何定标流程图；

图2为自动角反射器顶点坐标补偿示意图；

图3为单层电离层模型内插示意图。

具体实施方式

本发明提供了一种高精度星载sar几何定标方法，其流程图如图1所示。以一种时宽和带宽组合的成像模式为例作为一种实施方式，具体实现步骤如下：

步骤1，根据sar系统设计的波位表，查找出与该时宽和带宽组合模式相同的波位号。

步骤2，通过对sar卫星系统的成像任务规划，以查找出的任意波位对几何定标场区域进行成像，至少成像3次；其中，几何定标场区域的sar影像用于计算sar卫星系统的几何定标参数。

步骤3，远程调整几何定标场的所有自动角反射器指向，使角反射器的法线方向与入射波的方向保持一致。根据自动角反射器调整的方位角和俯仰角，精确计算成像时刻各个角反射器顶点的地理坐标。同时，对各个角反射器顶点的地理坐标进行固体潮改正补偿。

由于自动角反射器设备在水平和垂直方向上转动，其顶点位置会随着方位角和俯仰角的变化而变化。为了获取角反射器顶点的精确地面坐标，需要对角反射器顶点的坐标进行补偿修正。假设以云台垂直方向的转动轴中心为坐标系原点，以真北方向为坐标系的x轴，建立左手直角坐标系o-xyz，如图2所示。其中，p为角反射器顶点，p'为p在xoy平面的投影，p'a和p'b分别为p'到y轴和x轴的距离，ω为点p的方位角，为点p的俯仰角。

根据高斯投影的方法，将地面控制点wgs‐84系下的大地坐标(lat0,lon0,h0)转换成高斯平面坐标(x0,y0,h0)。由此可知，坐标系原点坐标o为(x0,y0,h0+h)，点p的高斯平面坐标为其中，h为云台垂直方向的转动中心相对于地面控制点的高度，po是角反射器顶点与云台垂直方向的转动中心之间的距离。再根据高斯投影关系，反算出点p在wgs‐84系下的大地坐标(latp,lonp,hp)，即获得高精度角反射器顶点的地理坐标。

由于地球自传，地球上任意点相对于月球和太阳的位置发生周期地改变，故重力固体潮随时间而周期变化。在进行角反射器点的重力固体潮校正时，一般采用下面公式计算：

δg(t)＝δthg(t)-δfc(1)

式中，δth为潮汐因子，一般采用1.15或1.16；δfc是地心纬度的函数，即

式中，为地理纬度的函数，即cm为地心至月心的平均距离，rm为地心至月心的距离，cs为地心至日心的平均距离，rs为地心至日心的距离，zm为月球对地面点的地心天顶距，zs为太阳对地面点的地心天顶距。

步骤4，对几何定标场区域的回波数据进行成像处理，获得sar影像和辅助参数文件。

步骤5，利用基于质心法的亚像元点提取方法，对每景sar影像中的角反射器点进行精确提取，获得角反射器点的像方坐标。

亚像元方法是利用目标像点的灰度分布特性通过内插细分算法确定出目标像点位置，从而使测量精度可以达到亚像元级。它是通过构造一个尽可能精确反映目标区域像素灰度值位置和目标像点的质心位置之间关系的数学模型，从而实现对像点位置的精确估计。为了提高求得的质心位置的精度，需要将读入的图像进行细分，通过对图像各点的灰度值的插值，得到插值曲线，求得细分后各点的灰度值。根据不同的精度要求，图像细分比例可取0.1、0.01、0.001等。

其算法具体步骤如下：

1)通过sar图像间接定位模型，由地面大地坐标计算角反射器对应的sar图像像点坐标，再由影像地面点坐标和轨道参数获得该地面点对应的方位向时间和角反射器点的星地距离后，解算该点对应的行列号；其中地面点坐标，由外业gps测量获得；轨道参数由星上gps获得，下传给地面使用。

2)在步骤1)中计算的行列号附近区域(如10*10像素)，通过计算机寻找在此区域像素灰度值最大的亮点，再通过人工目视判读，比较周围地物，确定该点是否为角反射器在sar图像上的粗略位置，初步确定角反射器的粗位置。

3)以步骤2)中的粗略位置所对应像素点为中心，选取其邻域(例如选取5×5领域)为目标质心寻找区。将该领域内所有像点对应的坐标和灰度强度信息按顺序提取出来，并把该灰度赋予每个像素的几何中心，例如对像素坐标点(8000，5000)，其中心(7999.5，4999.5)的灰度为该像素的灰度大小。

4)对该目标区内每一像素进行的细分，细分尺寸根据不同的精度要求可以是10×10或100×100或1000×1000，当细分尺寸为10×10时，插值点步长为0.1，它的精度即可达到0.1个像素。采用一定的数学函数方法(如双线性插值方法、sinc函数插值方法等)求得每一细分点的灰度。

5)计算图像回波强度质心。

在计算质心时，采用逐行分析计算，取平均值法计算质心，其实现原理为：

对目标区域，逐行对像素图像(行列都以细分后的单位为基准，如细分尺寸为100×100，则行列基本单位为原图像的0.01，即0.01行/列)进行处理(包括插值、图像分割、质心计算等)，并求出每一行的质心坐标和平均灰度值，质心坐标也就是x轴方向坐标，最后将每一行的(取决于图像的高度)质心坐标求平均值，即得原图像质心位移。

计算每一行质心的坐标x_u公式为：

其中，n为被插值点的个数，m为细分后点的个数，xk为被插值点的坐标，xu为细分后点的行坐标，ik为被插值点的灰度值，f(u)为细分后点的灰度值。

然后，以每一行质心坐标和对应灰度值为参数，求整个目标区的质心公式为：

其中，x_u为每一行质心坐标，yu为细分后点的列坐标，f__(_u_)为细分后的点对应平均灰度值。

计算得到的(xc,yc)即为目标区域的质心，也即角反射器点对应的像方坐标。计算质心时，根据插值点步长(如有效步长为0.001)的选取和像素的计算精度要求，精度取到0.0001即可。

步骤6，针对每景sar影像，计算雷达信号在天线相位中心到各个角反射器点之间的大气传播延迟改正值，大气传播延迟主要分为中性大气延迟和电离层延迟。

jamesc.owens结合流体静力学方程和非理想气体公式推导出中性大气中干分量和湿分量的距离改正模型，如下式：

上式中地表压强psruf＝pd+e，pw为大气可降水量；md、mw分别为干、湿大气分子量，md＝28.9644kg/kmol，mw＝18.0152kg/kmol；rm为摩尔气体常量，rm＝8.31451j/(mol·k)。

k1、k2值与卫星载荷发射电磁波信号的波长λ相关，求解公式：

由式(5)可知，干大气和湿大气的延迟量求解还需要获取大气的地表压强和大气可降水量，所述的地表压强和大气可降水量通常从全球大气模型中获得。全球大气模型(globalatmosphericmodels,gam)可以提供由卫星和气象基站记录的地表气象数据，包括温度、气压、水汽含量、风速风向等。本发明以美国国家环境预报中心(ncep)的大气分析模型作为外部数据，模型提供了每隔6小时1°×1°的经纬网格点存储的等压面数据。

图3为单层电离层模型示意图，可看作把目标点t垂直方向上的全部电子都压缩于一点t'处，电离层延迟是压缩点t'处的tec含量和视向与压缩点t'处单层电离层法线方向的夹角有很大关系。由此，可以得到雷达信号在传播路径上的电离层延迟：

式中k＝40.28m³/s²，tec的单位是1016个电子/米²，θ'是压缩点t'处的入射角。由上式可知，电离层延迟与tec和电磁波频率fc有关(已知)。

欧洲定轨中心(code)以ionex(ionosphereexchangeformat)文件格式每天给出，一天中从utc零时到24时，每过2个小时会产生一幅全球电离层tec图，每天给出13幅(从2014年10月19日开始，以1小时为间隔每天给出25幅)，igs计算电离层的格网划分按照全球经度方向间隔5°，纬度方向间隔2.5°，经度范围为‐180°w至180°e，纬度范围为87.5°s至‐87.5°n，共有5183个格网点。如果要获得某一特定时间穿刺点的延迟，要经过以下步骤：首先对历元时刻的tec含量采用双线性内插原理进行四站网格空间内插；接着按实际的时间对历元之间做内插，具体过程如下：

第一，在ionex文件中取得全球的tec含量图，然后以穿刺点的经纬度(lat,lon)为参数，从它附近的网格点做内插，具体的公式为：

式中，i0,0，i1,0，i1,1，i0,1为穿刺点相邻四角网格点的tec含量，lat0为左下角网格点的纬度，lon0为左下角网格点的精度，δlat为电离层网格点的纬度间隔，δlon为电离层网格点的经度间隔。

第二，得到空间内插后进行相应时间点t的时间双线性内插，具体的公式是：

式中，ti和ti+1为时间点t两个临近时间格网的时间。

第三，取得穿刺点对应时刻垂直方向上的tec后，根据公式(7)就能获得雷达信号在传播路径上的电离层延迟值。

步骤7，分别计算每景sar影像的几何定标参数，即距离向起始时间改正值δtf、方位向起始时间改正值δts、采样频率改正值δf、脉冲重复频率改正值δfp。

sar信号的发射和接收都是在时间尺度上完成的，主要包括快时间(距离向)和慢时间(方位向)。这二维的时间误差主要影响sar影像在距离向和方位向的几何定位误差，是星载sar几何定位的主要误差源。由此，构建星载sar几何定标模型：

式中，tf、ts分别为距离向的快时间和方位向的慢时间，tf0、ts0分别为距离向起始时间的测量值和方位向起始时间测量值，tdelay为大气传播延迟时间，为大气传播延迟值除以光速，δtf、δts为系统时延误差，x、y为像素坐标，fs为采样频率，fp为脉冲重复频率(prf)。

星上记录的方位向起始时间ts'0是雷达信号的接收时刻，然而sar卫星是在持续运动中发射和接收雷达信号，等效的sar成像时刻ts0为雷达信号的发射时刻与接收时刻的中间时刻，故需进行补偿：

式中，n为雷达信号从发射到接收经历整周期的个数，tsample_delay为星上记录的采样时间延迟。

星载sar几何定标模型(10)可以表示成如下形式：

式(12)的误差方程为

v＝bx-l(13)

其中，

设初始值δtf＝0、δts＝0，fs和fp由星上下传的辅助参数文件获得，采用迭代运算的方法获取几何定标参数x＝[δtfδtsδfsδfp]^t。几何定标参数解算的具体步骤如下：

1)根据sar影像间接定位算法，通过角反射器点的地面坐标和轨道参数，可以计算出该地面角反射器点对应的sar影像方位向时间ts和该点的星地距离rst(即距离向时间tf)。

2)根据公式(10)，利用定标场地面布设的n个角反射器点，可以组建n个方程组。通过基于谱修正法迭代法的最小二乘平差，可精确计算方位向起始时间改正值δts、距离向起始时间改正值δtf、prf改正值δfp＝fp'-fp和采样频率改正值δfs＝fs'-fs。

3)更新四个几何定标参数，重新执行步骤1)，再次计算四个几何定标参数，判断两次计算的几何定标参数是否小于一个预先设定的极小值，若小于则迭代终止，否则转向1)继续迭代运算。

4)迭代运算结束后的结果，就是该景sar影像的几何定标参数。

步骤8，根据几何定标场区域多景sar影像的几何定标参数，取均值，获取该成像模式的几何定标参数。

步骤9，针对该时宽和带宽组合模式，对几何精度验证场区域进行成像任务规划、拍摄，至少2个精度验证区域；其中，几何精度验证场区域的sar影像用于对上述步骤获取的几何定标参数的正确性进行验证。

步骤10，根据步骤8计算出的几何定标参数，对几何精度验证场区域的sar回波数据进行补偿，并进行成像处理，获取sar影像和辅助参数文件。

步骤11，获取几何精度验证场的地面控制点(角反射器点、gnss测量点等)的像方坐标，并重复进行步骤6的操作，获得每个地面控制点的大气传播延迟改正值。

步骤12，对几何定标后的几何精度验证场区域sar影像进行精度评价。根据sar影像中地面控制点的像方坐标和大气传播延迟改正值，利用rd模型的直接定位方法，计算sar影像中每个控制点的物方坐标，并与每个地面控制点的真实坐标进行对比，所得的差值即为该点的定位误差。统计所有地面控制点定位误差的标准差，即得到sar系统在该成像模式下经几何定标后无控制点情况的几何定位精度。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。