CN104112078A - 一种自适应几何关系强弱的有理函数模型区域网平差方法 - Google Patents

Abstract

本发明涉及一种自适应几何关系强弱的有理函数模型区域网平差方法，它包括以下步骤：1)读取卫星遥感影像的有理函数模型参数、归一化因子，以及连接点、控制点和检查点的影像坐标；2)根据读取的有理函数模型参数、归一化因子和连接点坐标，计算每一个连接点的所有交会角，并取每一个连接点的最大的交会角作为立体几何关系强弱的量化指标；3)根据每一个连接点最大交会角大小，将连接点分为第一类连接点和第二类连接点；4)根据连接点的最大交会角大小，算法自适应确定连接点的物方大地坐标初始值；5)建立误差方程，进行区域网平差；6)输出有理函数模型精化参数，完成基于RFM的区域网平差。本方法特别适用于宽视角、大范围的遥感影像平差以及平差数据复杂、数据来源广泛的多源数据间的区域网平差。

Description

一种自适应几何关系强弱的有理函数模型区域网平差方法

技术领域

本发明涉及摄影成像技术领域，涉及一种有理函数模型的区域网平差方法，特别是关于一种适用于航天及深空探测卫星遥感影像的自适应立体成像几何关系强弱的有理函数模型的区域网平差技术。

背景技术

随着遥感卫星的不断发展，各种新的轨道器和传感器不断涌现，传感器成像形式的多样化使得严密成像几何模型的形式日益复杂，构建难度也不断增加，这样加大了用户的使用难度，提高了用户使用的专业门槛，从而限制卫星影像的使用和发展。另外，构建严密成像几何模型所需的卫星平台和有效载荷等关键技术参数是航天大国的技术机密，因此遥感卫星的几何处理往往采用有理函数模型(RFM)这一通用的几何模型。

RFM具体定义关系式如下：

$X=\frac{{\mathrm{Num}}_L(P,L,H)}{{\mathrm{Den}}_L(P,L,H)}$    (1)

$Y=\frac{{\mathrm{Num}}_s(P,L,H)}{{\mathrm{Den}}_s(P,L,H)}$

式中，(X,Y)为归一化的影像坐标，(P,L,H)为归一化的地面点坐标，Num_L(P,L,H)、Den_L(P,L,H)、Num_s(P,L,H)、Den_s(P,L,H)是关于P、L和H的三次多项式：

Num_L(P,L,H)＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²

            +a₁₀H²+a₁₁PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³+a₁₇PH²

            +a₁₈L²H+a₁₉P²H+a₂₀H³

Den_L(P,L,H)＝b₁+b₂L+b₃P+b₄H+b₅LP+b₆LH+b₇PH+b₈L²+b₉P²

            +b₁₀H²+b₁₁PLH+b₁₂L³+b₁₃LP²+b₁₄LH²+b₁₅L²P+b₁₆P³+b₁₇PH²

            +b₁₈L²H+b₁₉P²H+b₂₀H³

Num_s(P,L,H)＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²

            +c₁₀H²+c₁₁PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³+c₁₇PH²

            +c₁₈L²H+c₁₉P²H+c₂₀H³

Den_s(P,L,H)＝d₁+d₂L+d₃P+d₄H+d₅LP+d₆LH+d₇PH+d₈L²+d₉P²

            +d₁₀H²+d₁₁PLH+d₁₂L³+d₁₃LP²+d₁₄LH²+d₁₅L²P+d₁₆P³+d₁₇PH²

            +d₁₈L²H+d₁₉P²H+d₂₀H³

式中，三次多项式的系数a₁,…,a₂₀，b₁,…,b₂₀，c₁,…,c₂₀，d₁,…,d₂₀是有理函数模型参数(RPC)，由RPC文件提供。b₁和d₁通常为1。

所谓归一化，是一项数据范围相差很大、避免计算中舍入误差的方法，其方法引入一项辅助性的概念—“归一化因子”，在RFM中运用这种方法对地面点和影像点的坐标进行处理。归一化地面点坐标定义公式为

$P=\frac{\mathrm{Latitude}-\mathrm{LAT}\_\mathrm{OFF}}{\mathrm{LAT}\_\mathrm{SCALE}}$

$L=\frac{\mathrm{Longitude}-\mathrm{LONG}\_\mathrm{OFF}}{\mathrm{LONG}\_\mathrm{SCALE}}---\left(2\right)$

$H=\frac{\mathrm{Height}-\mathrm{HEIGHT}\_\mathrm{OFF}}{\mathrm{HEIGHT}\_\mathrm{SCALE}}$

式中，归一化因子LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF、HEIGHT_SCALE是RPC文件中包含的地面点坐标归一化模型参数。Latitude表示经度、Longitude表示纬度、Height表示高程(某点沿铅垂线方向到大地水准面的距离，通常称为绝对高程或海拔)此三项代表地面点的空间坐标。

归一化影像坐标定义公式为

$X=\frac{\mathrm{Sample}-\mathrm{SAMP}\_\mathrm{OFF}}{\mathrm{SAMP}\_\mathrm{SCALE}}$    (3)

$Y=\frac{\mathrm{Line}-\mathrm{LINE}\_\mathrm{OFF}}{\mathrm{LINE}\_\mathrm{SCALE}}$

式中，归一化因子SAMP_OFF、SAMP_SCALE、LINE_OFF、LINE_SCALE也是RPC文件中包含的影像坐标归一化模型参数，Sample表示影像列坐标，其数值记为s，line表示影像行坐标，其数值记为l。

基于RFM具有模拟精度高，通用性好，应用方便，计算量小等等优点，是目前航天摄影测量中应用较为广泛的一种通用成像几何模型。

RFM不仅可以用于单幅遥感影像的几何精校正与正射校正，还可以用于多景遥感影像的立体测图。在利用RFM进行多像立体测图时，不同平台、不同传感器、不同分辨率的影像均可由统一的模型求解，并利用区域网平差的方法在稀少控制点的情况下提高定位精度，进行模型精化。

然而，在利用有理函数模型进行区域网平差时，由于数据的来源广泛、数据特征各异，利用传统的区域网平差，以RFM的精化参数和各连接点的三维坐标为未知参数进行求解时往往出现方程病态、求解不稳定、求解精度降低等情况。如在进行同轨影像平差时，影像间的连接点即为成像时的同一成像点，连接点无立体信息；相邻轨影像间同样也会产生交会条件弱的问题；在多类型多源数据进行联合平差时，由于数据的来源各异，连接点的交会条件也有所不同，在利用传统的方式进行区域网平差时，往往会使得求解不稳定，解算精度低。例如，我国高分一号卫星的16米分辨率影像由四个在一条线上排列的相机同时推扫式成像，相机间有较小的重叠，从而形成大幅宽的影像。由于相机无法按理想的方式排列成像，且成像时间同步性、姿态轨道精度等因素的影响，使得四个相机影像利用RFM正射纠正后无法无缝镶嵌。此时利用传统的RFM区域网平差方法求解出现病态，求解无法收敛，而利用只求解大地坐标中的经纬度进行平差时对多时相多角度数据精度不够，现有处理方法的不足成了这类新型传感器、多源卫星遥感数据联合平差使用时的瓶颈问题。

发明内容

针对上述问题，本发明的目的在于解决现有技术的不足，提出一种自适应成像立体几何关系强弱的有理函数模型区域网平差方法，该方法能够通过卫星遥感影像间几何交会关系，自适应确定求解连接点的未知数，并根据卫星覆盖区域的地形特性，自动选择连接点高程初始值来源，最后进行整体平差求解，对RFM模型进行精化。该方法能够有效解决因区域网平差中部分卫星影像间几何关系弱而引起的求解不稳定、精度差的弱点，并兼顾几何关系较强的卫星遥感影像，使得平差过程中不损伤原有的几何精度，处理后影像一致性较好。

为实现上述目的，本发明采取以下技术方案：一种自适应几何关系强弱的有理函数模型区域网平差方法，它包括以下步骤：

1)读取卫星遥感影像的有理函数模型参数、归一化因子，以及连接点、控制点和检查点的影像坐标；

2)根据读取的有理函数模型参数、归一化因子和连接点坐标，计算每一个连接点的所有交会角，并取每一个连接点的最大的交会角作为立体几何关系强弱的量化指标；

3)根据每一个连接点最大交会角大小，将连接点分为第一类连接点和第二类连接点；

4)根据连接点的最大交会角大小，算法自适应确定连接点的物方大地坐标初始值，具体过程为：针对第一类连接点：首先获取影像范围内的高差dh＝Height_max-Height_min，式中的Height_max表示最大高程；Height_min表示最小高程；然后计算因高差引起的最大投影差ds＝dh×tanθ_max，当ds≤k×r时，高程初值取值为参数HEIGHT_OFF，经纬度初值用有理函数模型和高程初值直接求解，式中的k表示精度要求；r表示影像分辨率；当ds≥k×r时，利用初始高程HEIGHT_OFF和有理函数模型参数，直接求解出初始的大地坐标，然后利用该大地坐标经纬度在辅助数字高程模型上内插出新的高程值，比较新的高程值与原始高程的差值，如果差值较小，则该经纬度与高程作为解算初值，如果差值较大，则利用新的高程替代HEIGHT_OFF为初值迭代求解；针对第二类连接点：利用有理函数模型空间前方交会的方式直接得到初始值；

5)建立误差方程，进行区域网平差，其中：

对于控制点，其线性化误差方程为

$F_x=F_{x0}+\frac{{\∂F}_x}{{\∂e}_0}\·{\mathrm{\Δe}}_0+\frac{{\∂F}_x}{{\∂e}_1}\·{\mathrm{\Δe}}_1+\frac{{\∂F}_x}{{\∂e}_2}\·{\mathrm{\Δe}}_2=0$

$F_y=F_{y0}+\frac{{\∂F}_y}{{\∂f}_0}\·{\mathrm{\Δf}}_0+\frac{{\∂F}_y}{{\∂f}_1}\·{\mathrm{\Δf}}_1+\frac{{\∂F}_y}{{\∂f}_2}\·{\mathrm{\Δf}}_2=0$

对于第一类连接点，其线性化误差方程为

$F_x=F_{x0}+\frac{{\∂F}_x}{{\∂e}_0}\·{\mathrm{\Δe}}_0+\frac{{\∂F}_x}{{\∂e}_1}\·{\mathrm{\Δe}}_1+\frac{{\∂F}_x}{{\∂e}_2}\·{\mathrm{\Δe}}_2+\frac{{\∂F}_x}{\∂\mathrm{lat}}\·\mathrm{\Δlat}+\frac{{\∂F}_x}{\∂\mathrm{lon}}\·\mathrm{\Δlon}=0$

$F_y=F_{y0}+\frac{{\∂F}_y}{{\∂f}_0}\·{\mathrm{\Δf}}_0+\frac{{\∂F}_y}{{\∂f}_1}\·{\mathrm{\Δf}}_1+\frac{{\∂F}_y}{{\∂f}_2}\·{\mathrm{\Δf}}_2+\frac{{\∂F}_y}{\∂\mathrm{lat}}\·\mathrm{\Δlat}+\frac{{\∂F}_y}{\∂\mathrm{lon}}\·\mathrm{\Δlon}=0$

对于第二类连接点，其线性化误差方程为

$F_x=F_{x0}+\frac{{\∂F}_x}{{\∂e}_0}\·{\mathrm{\Δe}}_0+\frac{{\∂F}_x}{{\∂e}_1}\·{\mathrm{\Δe}}_1+\frac{{\∂F}_x}{{\∂e}_2}\·{\mathrm{\Δe}}_2+\frac{{\∂F}_x}{\∂\mathrm{lat}}\·\mathrm{\Δlat}+\frac{{\∂F}_x}{\∂\mathrm{lon}}\·\mathrm{\Δlon}+\frac{{\∂F}_x}{\∂h}\·\mathrm{\Δh}=0$

$F_y=F_{y0}+\frac{{\∂F}_y}{{\∂f}_0}\·{\mathrm{\Δf}}_0+\frac{{\∂F}_y}{{\∂f}_1}\·{\mathrm{\Δf}}_1+\frac{{\∂F}_y}{{\∂f}_2}\·{\mathrm{\Δf}}_2+\frac{{\∂F}_y}{\∂\mathrm{lat}}\·\mathrm{\Δlat}+\frac{{\∂F}_y}{\∂\mathrm{lon}}\·\mathrm{\Δlon}+\frac{{\∂F}_y}{\∂h}\·\mathrm{\Δh}=0$

上式中，F_x表示方程在列方向的误差，F_y表示方程在行方向的误差；F_x0、F_y0表示泰勒级数0次项，e₀,e₁,e₂和f₀,f₁,f₂为仿射变换参数；line、sample为影像行列坐标，x、y为精化后的影像行列坐标；表示纬度、经度和高程的一阶偏导数；

6)输出有理函数模型精化参数，完成基于RFM的区域网平差。

所述步骤2)中任一连接点的最大交会角的计算方法如下：

①选取该连接点在第一张影像上的坐标(line₁,sample₁)；

②利用第一张影像的有理函数模型中的归一化参数求解该影像范围内的最大高程Height_max和最小高程Height_min，求解公式如下：

Height_max＝HEIGHT_OFF+HEIGHT_SCALE

Height_min＝HEIGHT_OFF-HEIGHT_SCALE

③根据有理函数模型的正变换，由连接点在第一张影像上影像坐标(line₁,sample₁)、影像的最大高程Height_max、最小高程Height_min，求解该像点在最大高程面上的大地经纬度坐标(Latitude_maxh,Longitude_maxh,Height_max)和该像点在最小高程面上的大地经纬度坐标(Latitude_minh,Longitude_minh,Height_min)；

④将像点在最大高程面上的大地经纬度坐标(Latitude_maxh,Longitude_maxh,Height_max)和在最小高程面上的大地经纬度坐标(Latitude_minh,Longitude_minh,Height_min)转换为星固直角坐标系下的直角坐标(X_maxh,Y_maxh,Z_max)，(X_minh,Y_minh,Z_min)；

⑤计算影像上最大高程与最小高程对应星固直角坐标的三维矢量[X_maxh-X_minh,Y_maxh-Y_minh,Z_max-Z_min]；

⑥选取该连接点在下一张影像上的坐标(line_i,sample_i)，并重复步骤2)～5)，直至获得连接点在所有的影像上的最大高程与最小高程对应星固直角坐标的三维矢量；

⑦利用矢量夹角公式求取每两个矢量的夹角，选取最大矢量夹角θ_max即最大交会角。

所述步骤3)中，第一类连接点和第二类连接点的分类方法如下：

①将各连接点按照最大交会角的大小由小到大排列；

②从第一个连接点开始，进行空间前方交会解算，如果第一个连接点解算过程中迭代收敛，能解算出稳定的解，则取第一个连接点的最大交会角的值作为阈值θ₁，结束计算；如果不能收效，则对下一个连接点进行空间前方交会解算，直至找出能迭代收敛，解算正确稳定的连接点为止，并取这个点的最大交会角的值为阈值θ₁；如果所有点都不能迭代收敛，则取最后一个连接点的最大交会角值加上1度作为阈值θ₁；

③当连接点的最大交会角值小于阈值θ₁时，该连接点分类为第一类连接点，否则分类为第二类连接点。

本发明由于采取以上技术方案，其具有以下优点：1、本发明提出了一种稳健的基于有理函数模型的连接点交会角计算方法，通过引入连接点的最大交会角，使得连接点几何关系强弱可量化。2、本发明通过对不同几何关系强弱连接点自适应的未知数确定，使得平差求解方案更为稳健，提高了平差的精度及适用范围。3、本发明根据不同条件下精度要求，自适应确定连接点平差时大地坐标及高程的初值，使得平差求解稳定性和精度提高。4、本发明通过以上三条技术的引入，解决了因几何关系强弱不一致而导致的现有区域网平差技术求解不稳定、精度不高的问题，提高了平差解算精度和稳定性，扩大了平差解算的适用范围。本方法特别适用于宽视角、大范围的遥感影像平差以及平差数据复杂、数据来源广泛的多源数据间的区域网平差。

附图说明

图1本发明方法的流程示意图；

图2是本发明立体成像交会角示意图；其中，图2(a)表示弱的交会关系；图2(b)表示强的交会关系。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

如图1所示，本发明的自适应几何关系强弱的有理函数模型区域网平差方法，包括以下步骤：

1、读取卫星遥感影像的有理函数模型参数(RPC)、归一化因子，以及连接点、控制点和检查点的影像坐标。

卫星影像产品中附带RPC文件(以文档的形式给出)，而RPC文件中包含多项式参数和归一化因子(背景技术中已详细给出各参数含义)。

连接点含义：在参与平差的影像中，选取各影像上对同一物体成像的点，即为连接点，连接点的影像坐标是连接点在每一幅其出现的影像上的像点坐标。用(line_i,sample_i)表示连接点在第i幅影像上的影像坐标。连接点可由影像匹配或人工选点的方式获取，两种方式均为本领域现有的成熟方法，故不予赘述。

控制点含义：已知影像坐标和其对应地物的地面点的三维坐标(由大地经纬度坐标(Latitude,Longitude)和高程Height构成)，在平差中作为真值提高平差精度的点。

大地经纬度坐标的含义：在大地坐标系中某点的大地经度是通过该点的子午面与本初子午面所构成的二面角。由本初子午面起算，向东为正，称为东经，向西为负，称为西经，该点的法线与赤道面的夹角，叫做大地纬度，由赤道面起算，向北为正，称为北纬，向南为负，称为南纬。由大地经度和大地纬度所表示的地面点坐标为大地经纬度坐标，根据大地经纬度坐标和高程可以标注地面点的三维坐标。

检查点含义：已知影像坐标和其对应地物的地面点的三维坐标，不直接参与平差，只作为真值检查平差精度的点。

2、根据读取的有理函数模型参数(RPC)、归一化因子和连接点坐标，计算每一个连接点的所有交会角(交会角为摄影光线夹角)，并取每一个连接点的最大的交会角作为立体几何关系强弱的量化指标。

立体几何关系的强弱含义：利用多影像立体量测，求解摄影物体的三维坐标时，如果影像空间前方交会深度方向误差椭圆小，求解稳定，则表示成像立体几何关系强，如果影像成像光线立体交会求解深度方向误差椭圆大，则求解出现病态，则表示成像立体几何关系弱。例如，如图2所示，P_l，P_r为两张影像，分别对P点进行摄影，两根光线的夹角为θ即为交会角，交会点处的椭圆为求解三维坐标时坐标的不确定性示意椭圆。(a)表示弱的交会关系，(b)表示强的交会关系。在航空摄影测量中，可以用航片的基高比(摄影基线与摄影高程的比值)或重叠度(两张航片的重叠区域与航片大小的比值，一般分航向重叠度和旁向重叠度)来量化立体成像几何关系的强弱。根据多源、复杂卫星遥感影像有理函数模型平差的特点，本发明利用连接点的最大交会角θ_max的值表示立体成像几何关系的强弱。

下面给出一种基于RFM求解某一连接点的最大交会角的方法：

1)针对某一连接点，选取该连接点在第一张影像上的坐标(line₁,sample₁)。

2)利用这张影像的有理函数模型中的归一化参数求解该影像范围内的最大高程Height_max和最小高程Height_min，求解公式如下：

Height_max＝HEIGHT_OFF+HEIGHT_SCALE

Height_min＝HEIGHT_OFF-HEIGHT_SCALE

3)根据有理函数模型的正变换(已知影像坐标line、sample和高程Height，求解大地经纬度Latitude、Longitude的过程为正变换，该求解过程为公开的成熟算法)，由连接点在第一张影像上影像坐标(line₁,sample₁)、影像的最大高程Height_max、最小高程Height_min，求解该像点在最大高程面上的大地经纬度坐标(Latitude_maxh,Longitude_maxh,Height_max)和该像点在最小高程面上的大地经纬度坐标(Latitude_minh,Longitude_minh,Height_min)。

4)将像点在最大高程面上的大地经纬度坐标(Latitude_maxh,Longitude_maxh,Height_max)和在最小高程面上的大地经纬度坐标(Latitude_minh,Longitude_minh,Height_min)转换为星固直角坐标系下的直角坐标(X_maxh,Y_maxh,Z_max)，(X_minh,Y_minh,Z_min)。

星固直角坐标系：采用笛卡尔直角坐标系，以星球质心为原点，星球赤道面为基本平面，采用右手螺旋法则，Z轴垂直于基本平面指向北极，X轴指向经度原点方向，Y轴垂直于X轴。

大地经纬度坐标系与星固直角坐标系之间的转换公式为：

$\left\{\begin{array}{c}X=(R+\mathrm{Height})\cos \left(\mathrm{Latitude}\right)\cos \left(\mathrm{Longitude}\right)\\ Y=(R+\mathrm{Height})\cos \left(\mathrm{Latitude}\right)\sin \left(\mathrm{Longitude}\right)\\ Z=(R+\mathrm{Height})\sin \left(\mathrm{Latitude}\right)\end{array}\right.---\left(4\right)$

上式中，R表示星球半径。

5)计算影像上最大高程与最小高程对应星固直角坐标的三维矢量[X_maxh-X_minh,Y_maxh-Y_minh,Z_max-Z_min]。

6)选取该连接点在下一张(即第i张,i＝2,…,n,n为该求解连接点所在摄影的影像张数)影像上的坐标(line_i,sample_i)，并重复步骤2)～5)，直至获得连接点在所有的影像上的最大高程与最小高程对应星固直角坐标的三维矢量。

7)利用矢量夹角公式求取每两个矢量的夹角，选取最大矢量夹角θ_max(即最大交会角)作为成像几何关系强弱的量化指标。

针对每一连接点重复运用上述方法，即可获得所有连接点的最大交会角。

3、根据每一个连接点最大交会角大小，将连接点分为两类，具体方法为：

1)将各连接点按照最大交会角的大小由小到大排列。

2)从第一个连接点开始，进行空间前方交会解算(空间前方交会是指利用同名光线的交会确定模型点空间位置的方法，该方法为公开的成熟算法)，如果第一个连接点解算过程中迭代收敛，能解算出正确稳定的解。则取第一个连接点的最大交会角的值作为阈值θ₁(在理论上存在一阈值θ₁：当最大交会角值小于阈值θ₁时，空间前方交会求解不稳定)，结束计算；如果不能收效，则对下一个连接点进行空间前方交会解算，直至找出能迭代收敛，解算正确稳定的连接点为止，并取这个点的最大交会角的值为阈值θ₁；如果所有点都不能迭代收敛，则取最后一个连接点的最大交会角值加上1度作为阈值θ₁。

3)当连接点的最大交会角值小于阈值θ₁时(θ_max<θ₁)，该连接点分类为第一类连接点，否则分类为第二类连接点。

4、根据连接点的最大交会角θ_max大小，算法自适应确定连接点的物方大地坐标初始值，具体方法为：

针对第一类连接点：

首先获取影像范围内的高差dh＝Height_max-Height_min，其中，Height_max、Height_min的求解方法见步骤2。

然后计算因高差引起的最大投影差ds＝dh×tanθ_max。

当ds≤k×r(k表示精度要求，r表示影像分辨率)时，高程初值取值为参数HEIGHT_OFF，经纬度初值用有理函数模型和高程初值直接求解。

当ds≥k×r时，利用初始高程HEIGHT_OFF和有理函数模型参数，直接求解出初始的大地坐标，然后利用该大地坐标经纬度在辅助数字高程模型(DEM)上内插出新的高程值，比较新的高程值与原始高程的差值，如果差值较小，则该经纬度与高程作为解算初值，如果差值较大，则利用新的高程替代HEIGHT_OFF为初值迭代求解。

针对第二类连接点：利用有理函数模型空间前方交会的方式直接得到初始值。

5、建立误差方程，进行区域网平差。

以在像方添加仿射变换参数为例，进行区域网平差。误差方程建立及求解过程如下：

1)根据像方模型，进行线性化。

以像方防射变换为例，像方模型为：

F_x＝e₀+e₁·sample+e₂·line-x   (5)

F_y＝f₀+f₁·sample+f₂·line-y

上式中，F_x表示方程在列方向的误差，F_y表示方程在行方向的误差；e₀,e₁,e₂和f₀,f₁,f₂为仿射变换参数；line、sample为影像行列坐标，x、y为精化后的影像行列坐标。

对上式进行一阶泰勒级数展开，可以得到线性化的误差方程：

$F_x=F_{x0}+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_0}\&CenterDot;{\mathrm{\&Delta;e}}_0+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_1}\&CenterDot;{\mathrm{\&Delta;e}}_1+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_2}\&CenterDot;{\mathrm{\&Delta;e}}_2=0$    (6)

$F_y=F_{y0}+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_0}\&CenterDot;{\mathrm{\&Delta;f}}_0+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_1}\&CenterDot;{\mathrm{\&Delta;f}}_1+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_2}\&CenterDot;{\mathrm{\&Delta;f}}_2=0$

$F_x=F_{x0}+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_0}\&CenterDot;{\mathrm{\&Delta;e}}_0+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_1}\&CenterDot;{\mathrm{\&Delta;e}}_1+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_2}\&CenterDot;{\mathrm{\&Delta;e}}_2+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}=0$    (7)

$F_y=F_{y0}+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_0}\&CenterDot;{\mathrm{\&Delta;f}}_0+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_1}\&CenterDot;{\mathrm{\&Delta;f}}_1+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_2}\&CenterDot;{\mathrm{\&Delta;f}}_2+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}=0$

$F_x=F_{x0}+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_0}\&CenterDot;{\mathrm{\&Delta;e}}_0+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_1}\&CenterDot;{\mathrm{\&Delta;e}}_1+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_2}\&CenterDot;{\mathrm{\&Delta;e}}_2+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}+\frac{{\&PartialD;F}_x}{\&PartialD;h}\&CenterDot;\mathrm{\&Delta;h}=0$    (8)

$F_y=F_{y0}+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_0}\&CenterDot;{\mathrm{\&Delta;f}}_0+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_1}\&CenterDot;{\mathrm{\&Delta;f}}_1+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_2}\&CenterDot;{\mathrm{\&Delta;f}}_2+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}+\frac{{\&PartialD;F}_y}{\&PartialD;h}\&CenterDot;\mathrm{\&Delta;h}=0$

上式中，F_x0、F_y0表示泰勒级数0次项，分别表示仿射变换参数的一阶偏导数；表示纬度、经度和高程的一阶偏导数；

对于控制点来说，未知数为仿射变换参数(e₀,e₁,e₂,f₀,f₁,f₂)，因此线性化误差方程为式(6)；对第一类连接点，未知数包括仿射变换参数和相应的物方经纬度坐标，因此线性化误差方程为式(7)。而第二类连接点，未知数包括了和经纬度坐标、高程，因此线性化误差方程为式(8)。

2)对每个控制点和连接点分别按1)中的误差方程，逐点计算法方程。

3)对法方程求解，并对未知数进行改正。并判断平差精度，如平差精度满足要求(如当未知数改正数小于0.00001时)，则输出平差结果，如平差精度不满足要求，则重复步骤1)，2)，3)，迭代求解。

6、输出有理函数模型(RFM)精化参数(如像方仿射变换参数(e₀,e₁,e₂,f₀,f₁,f₂))，完成基于RFM的区域网平差。

上述各实施例仅用于说明本发明，其中方法的实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (3)

1.一种自适应几何关系强弱的有理函数模型区域网平差方法，它包括以下步骤：

1)读取卫星遥感影像的有理函数模型参数、归一化因子，以及连接点、控制点和检查点的影像坐标；

2)根据读取的有理函数模型参数、归一化因子和连接点坐标，计算每一个连接点的所有交会角，并取每一个连接点的最大的交会角作为立体几何关系强弱的量化指标；

3)根据每一个连接点最大交会角大小，将连接点分为第一类连接点和第二类连接点；

4)根据连接点的最大交会角大小，算法自适应确定连接点的物方大地坐标初始值，具体过程为：

针对第一类连接点：首先获取影像范围内的高差dh＝Height_max-Height_min，式中的Height_max表示最大高程；Height_min表示最小高程；然后计算因高差引起的最大投影差ds＝dh×tanθ_max，当ds≤k×r时，高程初值取值为参数HEIGHT_OFF，经纬度初值用有理函数模型和高程初值直接求解，式中的k表示精度要求；r表示影像分辨率；当ds≥k×r时，利用初始高程HEIGHT_OFF和有理函数模型参数，直接求解出初始的大地坐标，然后利用该大地坐标经纬度在辅助数字高程模型上内插出新的高程值，比较新的高程值与原始高程的差值，如果差值较小，则该经纬度与高程作为解算初值，如果差值较大，则利用新的高程替代HEIGHT_OFF为初值迭代求解；

针对第二类连接点：利用有理函数模型空间前方交会的方式直接得到初始值；

5)建立误差方程，进行区域网平差，其中：

对于控制点，其线性化误差方程为

$F_x=F_{x0}+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_0}\&CenterDot;{\mathrm{\&Delta;e}}_0+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_1}\&CenterDot;{\mathrm{\&Delta;e}}_1+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_2}\&CenterDot;{\mathrm{\&Delta;e}}_2=0$

$F_y=F_{y0}+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_0}\&CenterDot;{\mathrm{\&Delta;f}}_0+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_1}\&CenterDot;{\mathrm{\&Delta;f}}_1+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_2}\&CenterDot;{\mathrm{\&Delta;f}}_2=0$

对于第一类连接点，其线性化误差方程为

$F_x=F_{x0}+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_0}\&CenterDot;{\mathrm{\&Delta;e}}_0+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_1}\&CenterDot;{\mathrm{\&Delta;e}}_1+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_2}\&CenterDot;{\mathrm{\&Delta;e}}_2+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}=0$

$F_y=F_{y0}+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_0}\&CenterDot;{\mathrm{\&Delta;f}}_0+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_1}\&CenterDot;{\mathrm{\&Delta;f}}_1+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_2}\&CenterDot;{\mathrm{\&Delta;f}}_2+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}=0$

对于第二类连接点，其线性化误差方程为

$F_x=F_{x0}+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_0}\&CenterDot;{\mathrm{\&Delta;e}}_0+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_1}\&CenterDot;{\mathrm{\&Delta;e}}_1+\frac{{\&PartialD;F}_x}{{\&PartialD;e}_2}\&CenterDot;{\mathrm{\&Delta;e}}_2+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_x}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}+\frac{{\&PartialD;F}_x}{\&PartialD;h}\&CenterDot;\mathrm{\&Delta;h}=0$

$F_y=F_{y0}+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_0}\&CenterDot;{\mathrm{\&Delta;f}}_0+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_1}\&CenterDot;{\mathrm{\&Delta;f}}_1+\frac{{\&PartialD;F}_y}{{\&PartialD;f}_2}\&CenterDot;{\mathrm{\&Delta;f}}_2+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lat}}\&CenterDot;\mathrm{\&Delta;lat}+\frac{{\&PartialD;F}_y}{\&PartialD;\mathrm{lon}}\&CenterDot;\mathrm{\&Delta;lon}+\frac{{\&PartialD;F}_y}{\&PartialD;h}\&CenterDot;\mathrm{\&Delta;h}=0$

上式中，F_x表示方程在列方向的误差，F_y表示方程在行方向的误差；F_x0、F_y0表示泰勒级数0次项，e₀,e₁,e₂和f₀,f₁,f₂为仿射变换参数；line、sample为影像行列坐标，x、y为精化后的影像行列坐标；表示纬度、经度和高程的一阶偏导数；

6)输出有理函数模型精化参数，完成基于RFM的区域网平差。

2.如权利要求1所述的一种自适应几何关系强弱的有理函数模型区域网平差方法，其特征在于：所述步骤2)中任一连接点的最大交会角的计算方法如下：

①选取该连接点在第一张影像上的坐标(line₁,sample₁)；

②利用第一张影像的有理函数模型中的归一化参数求解该影像范围内的最大高程Height_max和最小高程Height_min，求解公式如下：

Height_max＝HEIGHT_OFF+HEIGHT_SCALE

Height_min＝HEIGHT_OFF-HEIGHT_SCALE

③根据有理函数模型的正变换，由连接点在第一张影像上影像坐标(line₁,sample₁)、影像的最大高程Height_max、最小高程Height_min，求解该像点在最大高程面上的大地经纬度坐标(Latitude_maxh,Longitude_maxh,Height_max)和该像点在最小高程面上的大地经纬度坐标(Latitude_minh,Longitude_minh,Height_min)；

④将像点在最大高程面上的大地经纬度坐标(Latitude_maxh,Longitude_maxh,Height_max)和在最小高程面上的大地经纬度坐标(Latitude_minh,Longitude_minh,Height_min)转换为星固直角坐标系下的直角坐标(X_maxh,Y_maxh,Z_max)，(X_minh,Y_minh,Z_min)；

⑤计算影像上最大高程与最小高程对应星固直角坐标的三维矢量[X_maxh-X_minh,Y_maxh-Y_minh,Z_max-Z_min]；

⑥选取该连接点在下一张影像上的坐标(line_i,sample_i)，并重复步骤2)～5)，直至获得连接点在所有的影像上的最大高程与最小高程对应星固直角坐标的三维矢量；

⑦利用矢量夹角公式求取每两个矢量的夹角，选取最大矢量夹角θ_max即最大交会角。

3.如权利要求1或2所述的一种自适应几何关系强弱的有理函数模型区域网平差方法，其特征在于：所述步骤3)中，第一类连接点和第二类连接点的分类方法如下：

①将各连接点按照最大交会角的大小由小到大排列；

②从第一个连接点开始，进行空间前方交会解算，如果第一个连接点解算过程中迭代收敛，能解算出稳定的解，则取第一个连接点的最大交会角的值作为阈值θ₁，结束计算；如果不能收效，则对下一个连接点进行空间前方交会解算，直至找出能迭代收敛，解算正确稳定的连接点为止，并取这个点的最大交会角的值为阈值θ₁；如果所有点都不能迭代收敛，则取最后一个连接点的最大交会角值加上1度作为阈值θ₁；

③当连接点的最大交会角值小于阈值θ₁时，该连接点分类为第一类连接点，否则分类为第二类连接点。