一种星载合成孔径雷达影像的区域网平面平差方法与流程

本发明涉及摄影测量技术领域，特别是涉及一种星载合成孔径雷达影像的区域网平面平差方法。

背景技术：

合成孔径雷达(syntheticapertureradar，以下简称为sar)是现代雷达技术的重要发展方向之一，星载sar作为其应用平台的一种，其应用领域正在快速的拓展和深入，包括环境监测、陆地测图、自然资源探测、灾害管理以及海洋测绘等。与光学、红外成像相比，星载sar具有多极化、全天时、全天候成像等独特优点。近来，国内外新一代的高分辨率星载sar卫星先后成功发射，其成像系统性能不断提升，星载sar影像摄影测量应用价值越来越高。当前，星载sar技术的研究与应用是我国对地观测领域重点发展的方向。我国的星载sar技术正处于一个追赶世界先进水平的阶段，因此研究如何利用区域网平差技术解决其影像之间的几何定位差异在实际的生产处理中具有重要的意义。

传统的星载sar影像的处理主要基于距离-多普勒模型(以下简称为r-d模型)，然而由于卫星平台，工作模式的不同，利用r-d模型在实际的生产工作中需要根据不同的卫星平台参数进行相应的设置，难以满足通用化统一处理流程的需求，制约其大规模的实际运用。而有理多项式模型(rationalpolynomialcoefficient，以下简称为rpc模型)，是一种通用的卫星遥感影像的几何模型，是在充分利用卫星遥感影像附带的辅助参数基础上，根据构建的严格成像几何模型进行拟合而得到的广义传感器模型。因此，利用rpc模型替代原有的r-d模型，能够实现不同卫星平台sar影像的通用处理，对于大规模运用sar影像解决实际生产需要，提升其利用价值具有重要意义。

技术实现要素：

为了克服上述现有技术的不足，本发明提供了一种星载合成孔径雷达影像的区域网平面平差方法。

本发明所采用的技术方案是：

步骤1，读取卫星ccd阵列影像，选择相邻卫星影像重叠区域的同名点作为影像间的连接点，根据控制点位信息在卫星影像上刺出控制点和检查点，并将连接点、控制点和检查点的坐标信息存成相应的点文件；

步骤2，读取卫星ccd阵列影像的rpc参数文件、连接点文件、控制点文件和检查点文件，得到卫星ccd阵列影像的rpc模型以及控制点和检查点的地面坐标、影像坐标和连接点的像点坐标；

步骤3，确定地面点坐标和像方变换模型的初值，在确定地面点坐标时采用dem约束的前方交会方法迭代得到；

步骤4，利用量测得到像点坐标和对应的地面点坐标，针对连接点和控制点分别对定向参数和地面点坐标求偏导并逐点构建误差方程，并将构建的误差方程进行法化；

步骤5，对法方程利用谱修正法迭代估计求解改正定向参数和地面点坐标；

步骤6，利用更新的影像连接点地面点坐标和影像定向参数，重复步骤3-步骤5，通过迭代过程不断更新地面点坐标和影像定向参数，当满足定向参数中的平移参数均小于阈值时，整个平差迭代结束；当不满足预设条件时，返回步骤3继续迭代计算，直到满足迭代收敛条件；如果迭代次数达到预设迭代次数，仍然不能收敛，那么迭代结束；

步骤7，每完成一次平差迭代时，计算出检查点对应的地面点坐标，并同时计算此次平差能达到的物方精度；计算连接点的像方误差，统计本次迭代能达到的像方精度；

步骤8，当整个平差迭代结束时，输出所计算得到的定向参数改正数，以及平差物方精度和像方精度的精度报告，根据精度报告对平差效果进行评价。

步骤2中利用影像的rpc文件里的模型参数，构建影像的rpc模型，

其中rpc模型定义如下：

式(1)中：

(p，l，h)为正则化地面点坐标，与非正则化的地面点坐标(latitude，longitude，height)，以下简写(lat，lon，h)的关系如式(2)；

(x，y)表示正则化影像坐标，与非正则化的影像坐标(s，1)的关系如式(3)；

numl(p，l，h)、denl(p，l，h)、nums(p，l，h)、dens(p，l，h)为三次多项式，不具有实际物理意义，下标l和s分别代表影像列和行：

numl(p，l，h)＝n1+n2l+n3p+n4h+n5lp+n6lh+n7ph+n8l²+n9p²+n10h²+n11plh+n12l³+n13lp²+n14lh²+n15l²p+n16p³+n17ph²+n18l²h+n19p²h+n20h³

denl(p，l，h)＝m1+m2l+m3p+m4h+m5lp+m6lh+m7ph+m8l²+m9p²+m10h²+m11plh+m12l³+m13lp²+m14lh²+m15l²p+m16p³+m17ph²+m18l²h+m19p²h+m20h³

nums(p，l，h)＝c1+c2l+c3p+c4h+c5lp+c6lh+c7ph+c8l²+c9p²+c10h²+c11plh+c12l³+c13lp²+c14lh²+c15l²p+c16p³+c17ph²+c18l²h+c19p²h+c20h³

dens(p，l，h)＝d1+d2l+d3p+d4h+d5lp+d6lh+d7ph+d8l²+d9p²+d10h²+d11plh+d12l³+d13lp²+d14lh²+d15l²p+d16p³+d17ph²+d18l²h+d19p²h+d20h³

三次多项式的系数n1，…，n20，m1，…，m20，c1，…，c20，d1，…，d20是rpc文件中提供的模型参数，b1和d1通常为1；

式(1)中的正则化地面点坐标定义为：

式(2)中：

lat_off、lat_scale、long_off、long_scale、height_off、height_scale为rpc文件中包含的地面点坐标正则化模型参数；

latitude表示经度、longitude表示纬度、height表示高程，此三项即可代表地面点的空间坐标；

式(1)中的正则化影像坐标定义为：

式(3)中：

samp_off、samp_scale、line_off、line_scale为rpc文件中包含的影像坐标正则化模型参数；

sample代表影像列坐标，其数值即为s；line代表影像行坐标，其数值即为l，即影像中的坐标由(s，l)表示。

进一步的，步骤3中地面点的三维空间坐标通过逐次迭代计算方法求得，具体步骤包括：

1)建立前方交会误差方程，求解地面点坐标改正数，联立式(1)、(2)、(3)，得到未正则化的地面坐标和影像坐标的关系式

将式(4)按照泰勒级数展开成线性形式：

其中s⁰和l⁰为地面点p⁰坐标(lat0，lon0，h0)带入式(1)中得到的投影点影像坐标；p⁰坐标(lat0，lon0，h0)为迭代的初值；

为各偏导项；

δlat，δlon，δh为地面点改正数，由于使用dem数据进行前方交会，因此这里δh＝0，在之后省去相关项；

改写式(5)得到误差方程如下：

式(6)这里(vs，vl)为影像像点坐标的误差项；

多个影像上的同名像点可以分别列出相应的误差方程，联立之后写成矩阵形式如下：

v＝ax-l，p(7)

这里的p矩阵为权矩阵，这里为单位矩阵；

由此可以得到地面点坐标(lat，lon，h)的改正数x＝(δlat，δlon)：

x＝(a^tpa)^-1a^tpl(8)

将式(8)解算得到的改正数x对地面点坐标进行改正，并将新的地面点坐标带入dem内插得到新的高程，即得到新的地面点p¹坐标(lat1，lon1，h1)；

2)利用新的地面点p¹带入1)中的进行解算，从而得到新的地面点改正数并更新地面点坐标，然后利用dem内插高程；

3)重复迭代过程，得到地面点p²，…，pⁿ的坐标，直到两次迭代的坐标更新值小于所设定的阈值后，结束迭代过程，即得到了前方交会后地面点的初值。

进一步的，步骤4的具体实现方式如下，

在上述rpc模型的基础上建立仿射变换模型

式中，δy和δx为控制点在影像坐标系中的量测坐标与真实坐标的差值，即改正数；a1，a2，a3和b1，b2，b3是影像的定向参数，(s，l)是控制点在影像坐标系中的坐标；

在此基础上，联立式(4)、式(9)并线性化展开建立误差方程：

式中，δdlat，δdlon，δdheight为地面点坐标改正数但是在平面平差中，由于高程是通过dem内插得到的，因此跟height相关的项都为0，之后省去；

δa1，δa2，δa3，δb1，δb2，δb3为影像定向参数改正数；

vx，vy为像点坐标改正数；

fx0，fy0为像点坐标近似值与像点坐标观测值之差；

为误差方程对定向参数所求的偏导数；

为误差方程对地面点坐标所求的偏导数；

误差方程写成矩阵形式记为：

v＝bt+ax-l(11)

同样可以对每个控制点建立如下线性方程：

误差方程记为矩阵形式如下：

v＝bt-l(13)

其中式(11)和式(13)中各参数为：

t＝(δa1δa2δa3δb1δb2δb3)^t

x＝(δdlatδdlon)^t

根据最小二乘平差原理，基于误差方程构建法方程如下：

进一步的，步骤5的具体实现方式如下，

将法方程两边同时加上得到：

式中，e为和阶数相同的单位矩阵，由于式子两边都含有未知参数所以只能采用迭代的方法求解，其迭代公式为：

将解算出来的带入式(16)解算出重复迭代过程解算出当两次迭代的改正数差值小于阈值时，停止迭代，即得到了改正参数t和地面点坐标改正数x。

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

(1)采用rpc模型替代原有的r-d模型，避免了针对不同卫星平台的分别参数设置，实现了各类影像类型的统一处理，有助于提升在实际生产中的生产效率；

(2)利用区域网平差技术对sar影像的几何定位不一致进行补偿，使得最终成果影像之间重叠区域的接边精度和与地理参考之间的绝对定位精度都能够满足精度要求，保障实际生产中产品数据的精度可靠性；

(3)由于sar影像在成像时与地面之间具有倾角，在转化成rpc模型之后，直接利用rpc模型进行地面点的前方交会会存在高程解算异常的问题，影响整个结果的最终精度。利用数字高程模型(digitalelevationmodel，以下简称dem)对前方交会过程进行高程方向的约束，能够解决高程解算异常的问题，保障整个结果可靠性。

附图说明

图1为本发明技术方案流程图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步说明。

如图1所示，本发明提供一种星载合成孔径雷达影像的区域网平面平差方法，包括以下步骤：

步骤1，读取sar影像，利用匹配技术或者人工采集等方式选择相邻星载sar影像重叠区域的同名点作为影像间的连接点，同时，根据控制点位信息在影像上刺出控制点和检查点的像方坐标，其中，控制点布设的原则应当尽量能够分布在测区的四角；整理rpc参数文件格式符合国际规范；同时整理连接点、控制点和检查点的坐标信息并存成相应点文件。

步骤2，读取rpc参数文件、连接点文件、控制点文件和检查点文件，得到sar影像的rpc模型以及控制点和检查点的地面坐标、影像坐标，以及连接点的像点坐标；利用影像的rpc文件里的模型参数，构建影像的rpc模型。

rpc模型定义如下：

式(1)中：

(p，l，h)为正则化地面点坐标，与非正则化的地面点坐标(latitude，longitude，height)，以下简写(lat，lon，h)的关系如式(2)；

(x，y)表示正则化影像坐标，与非正则化的影像坐标(s，l)的关系如式(3)；

numl(p，l，h)、denl(p，l，h)、nums(p，l，h)、dens(p，l，h)为三次多项式，不具有实际物理意义，下标l和s分别代表影像列和行：

numl(p，l，h)＝n1+n2l+n3p+n4h+n5lp+n6lh+n7ph+n8l²+n9p²+n10h²+n11plh+n12l³+n13lp²+n14lh²+n15l²p+n16p³+n17ph²+n18l²h+n19p²h+n20h³

denl(p，l，h)＝m1+m2l+m3p+m4h+m5lp+m6lh+m7ph+m8l²+m9p²+m10h²+m11plh+m12l³+m13lp²+m14lh²+m15l²p+m16p³+m17ph²+m18l²h+m19p²h+m20h³

nums(p，l，h)＝c1+c2l+c3p+c4h+c5lp+c6lh+c7ph+c8l²+c9p²+c10h²+c11plh+c12l³+c13lp²+c14lh²+c15l²p+c16p³+c17ph²+c18l²h+c19p²h+c20h³

dens(p，l，h)＝d1+d2l+d3p+d4h+d5lp+d6lh+d7ph+d8l²+d9p²+d10h²+d11plh+d12l³+d13lp²+d14lh²+d15l²p+d16p³+d17ph²+d18l²h+d19p²h+d20h³

三次多项式的系数n1，…，n20，m1，…，m20，c1，…，c20，d1，…，d20是rpc文件中提供的模型参数，b1和d1通常为1。

所谓正则化，是一项处理无限大、发散以及一些不合理表示式的方法，其方法透过引入一项辅助性的概念——正则化因子。在rpc模型中运用这种方法对地面点和影像点的坐标进行处理。

式(1)中的正则化地面点坐标定义为：

式(2)中：

lat_off、lat_scale、long_off、long_scale、height_off、height_scale为rpc文件中包含的地面点坐标正则化模型参数；

latitude表示经度、longitude表示纬度、height表示高程(某点沿铅垂线方向到大地水准面的距离，通常称为绝对高程或海拔，本技术领域简称高程)，此三项即可代表地面点的空间坐标。

式(1)中的正则化影像坐标定义为：

式(3)中：

samp_off、samp_scale、line_off、line_scale为rpc文件中包含的影像坐标正则化模型参数；

sample代表影像列坐标，其数值即为s；line代表影像行坐标，其数值即为l。即影像中的坐标由(s，l)表示。

基于rpc模型具有模拟精度高，通用性好，应用方便，计算量小等等优点，但是，该模型同时存在一大缺点，即参数没有严格的几何意义，在进行对地面目标几何定位处理时，无需建立起对应具有物理意义的严密模型。

步骤3，确定地面点坐标和仿射变换模型初值。其中所有地面点的坐标可以通过rpc模型对连接点进行直接前方交会得到地面点坐标(lat，lon，h)作为初值。在本技术中，由于sar影像直接进行前方交会高程解算容易收敛异常，因此在前方交会时使用dem数据进行高程约束。地面点的三维空间坐标可通过逐次迭代计算方法求得，具体步骤包括：

1)建立前方交会误差方程，求解地面点坐标改正数。联立式(1)、(2)、(3)，得到未正则化的地面坐标和影像坐标的关系式

将式(4)按照泰勒级数展开成线性形式：

其中s⁰和l⁰为地面点p⁰坐标(lat0，lon0，h0)带入式(1)中得到的投影点影像坐标；p⁰坐标(lat0，lon0，h0)为迭代的初值，一般可通过将左、右影像的地面点坐标正则化平移参数的平均值或者根据rpc参数的一次项部分进行前方交会所获得；

为各偏导项；

δlat，δlon，δh为地面点改正数，由于使用dem数据进行前方交会，因此这里δh＝0，在之后省去相关项；

改写式(5)得到误差方程如下：

式(6)这里(vs，vl)为影像像点坐标的误差项；

多个影像上的同名像点可以分别列出相应的误差方程，联立之后写成矩阵形式如下：

v＝ax-l，p(7)

这里的p矩阵为权矩阵，这里为单位矩阵；

由此可以得到地面点坐标(lat，lon，h)的改正数x＝(δlat，δlon)：

x＝(a^tpa)^-1a^tpl(8)

将式(8)解算得到的改正数x对地面点坐标进行改正，并将新的地面点坐标带入dem内插得到新的高程，即得到新的地面点p¹坐标(lat1，lon1，h1)；

2)利用新的地面点p¹带入1)中的进行解算，从而得到新的地面点改正数并更新地面点坐标，然后利用dem内插高程。

3)重复迭代过程，得到地面点p²，…，pⁿ的坐标，直到两次迭代的坐标更新值小于所设定的阈值后，(本文阈值设置为两次更新的向量的模差值小于0.001，但不限于此)，然后结束迭代过程，即得到了前方交会后地面点的初值。

像方变换模型初值一般设定为

步骤4，构建区域网平差模型，建立误差方程。

利用量测得到像点坐标和对应的地面点坐标，针对连接点和控制点分别对定向参数和地面点坐标求偏导并逐点构建误差方程并将构建的误差方程进行法化。具体过程为：

在上述rpc模型的基础上建立仿射变换模型

式中，δy和δx为控制点在影像坐标系中的量测坐标与真实坐标的差值，即改正数；a1，a2，a3和b1，b2，b3是影像的定向参数，(s，1)是控制点在影像坐标系中的坐标。

在此基础上，联立式(4)、式(9)并线性化展开建立误差方程：

式中，δdlat，δdlon，δdheight为地面点坐标改正数但是在平面平差中，由于高程是通过dem内插得到的，因此跟height相关的项都为0，之后省去；

δa1，δa2，δa3，δb1，δb2，δb3为影像定向参数改正数；

vx，vy为像点坐标改正数；

fx0，fy0为像点坐标近似值与像点坐标观测值之差；

为误差方程对定向参数所求的偏导数；

为误差方程对地面点坐标所求的偏导数。

误差方程写成矩阵形式记为：

v＝bt+ax-l(11)

同样可以对每个控制点建立如下线性方程：

误差方程记为矩阵形式如下：

v＝bt-l(13)

其中式(11)和式(13)中各参数为：

t＝(δa1δa2δa3δb1δb2δb3)^t

x＝(δdlatδdlon)^t

根据最小二乘平差原理，基于误差方程构建法方程如下：

步骤5，对误差方程法化，求解法方程，获取改正数，更新影像仿射变换系数。对法方程利用谱修正法迭代估计求解改正参数t和地面点坐标x，具体过程为：

将法方程两边同时加上得到：

式中，e为和阶数相同的单位矩阵，由于式子两边都含有未知参数所以只能采用迭代的方法求解，其迭代公式为：

将解算出来的带入式(16)解算出重复迭代过程解算出当两次迭代的改正数差值小于阈值时(本文阈值设置为两次更新的向量的模差值小于0.001，但不限于此)，停止迭代，即得到了改正参数t和地面点坐标改正数x。

步骤6，利用更新的影像连接点地面点坐标和影像定向参数，重复步骤3-步骤5，通过迭代过程不断更新地面点坐标和影像定向参数，直至影像的定向参数中的平移参数a1，b1小于阈值时(本实施例为0.1个像元pixel，但是不限于此)时，平差迭代结束；当不满足预设条件时，返回步骤3继续迭代计算，直到满足迭代收敛条件；如果迭代次数达到预设迭代次数(本实施例设定的迭代次数为20次，但是不限于此)，仍然不能收敛，那么平差失败退出，此时的平差精度可能有损失。

步骤7，平差迭代完成时，计算出检查点对应的地面点坐标，通过计算出的检查点的地面点坐标和已知的检查点地面坐标之差，即为检查点精度，也就是最终的平差后所能达到的物方精度。同时统计所有影像连接点的地面坐标通过rpc模型投影至影像上的投影坐标与其原始影像坐标的差值，即为像点精度，也就是平差后所能达到的像方精度。

步骤8，输出记录有计算得到的定向参数改正数，以及平差物方精度和像方精度的精度报告，根据精度报告可以对平差效果进行评价。如果通过精度报告得知平差精度较高则说明平差效果较好，平差参数可以用于后续的生产和测绘中，如果通过报告得知平差精度较低，可以从精度报告中查找原因，进而避免对后续的应用产生影响。

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