CN113324527B - 一种同轨激光测高点与三线阵立体影像联合测绘处理方法 - Google Patents

Abstract

本发明公开了一种同轨激光测高点与三线阵立体影像联合测绘处理方法，该方法首先将激光测高点作为高程控制点，利用激光测高点和同轨获取的正视影像相对平面误差较小的特性，获取分布于各种地形类型(尤其是丘陵和山地地形)区域的激光测高点在三线阵立体影像上的准确像方坐标；然后以轨为单位将整轨三线阵立体影像构建平差区域网，开展激光测高点与立体影像的联合区域网平差；最后根据平差结果更新立体影像的成像几何模型。本发明充分利用激光测高点与光学影像同轨获取且高程精度极高的优势，在无地面控制点条件下，有效提升了三线阵立体影像的高程精度，为利用卫星影像开展全球范围高精度地理信息资源建设提供了技术支撑。

Description

一种同轨激光测高点与三线阵立体影像联合测绘处理方法

技术领域

本发明属于遥感与摄影测量技术领域，涉及一种激光测高点与立体影像联合几何处理方法，尤其涉及一种激光测高点与同轨获取的三线阵立体影像联合测绘处理方法。

背景技术

测图卫星立体影像目前已经成为世界各国开展多尺度测绘地理信息产品生成的主要数据源。利用卫星立体影像开展测绘应用的最大难点是保障影像的几何精度，尤其是高程精度。虽然通过使用足够数量的高精度地面控制资料可以有效提升卫星影像几何精度，然而开展全球范围测图应用时，受政治和经济等因素限制，很多区域根本无法获取控制资料。此外，完全不使用或者减少使用地面控制点，也是开展测绘应用时降低测图成本、提高测图效率的迫切需求。

2020年7月25日发射的资源三号03星(ZY3-03)，主要服务于全球范围1:5万比例尺的立体测图。卫星搭载了三台全色线阵推扫式光学相机和一台包含近红外、红、绿和蓝4个谱段的多光谱相机，三台全色相机分别为沿轨向前倾22°的前视相机、近乎垂直对地的正视相机和沿轨向后倾22°的后视相机，可以在一次过境摄影成像时，获取三线阵立体影像和4谱段多光谱影像。其中正视相机的地面分辨率均为2.1米，前视和后视相机的地面分辨率2.5米，多光谱相机的地面分辨率5.8米。幅宽均大于51千米。根据卫星影像应用实践情况来看，ZY3-03星三线阵立体影像的原始高程中误差约为6-10米、原始平面中误差8-20米，其高程精度无法满足全球范围1:5万比例尺测图应用精度要求。为此，ZY3-03星专门搭载了一台单波束激光测高仪，用于在获取立体影像的同时，以2Hz的对地观测频率同步获取立体影像范围内的高精度激光测高点，激光测高点没有足印影像，激光光斑在地面的直径约50米。激光测高仪作为一种高精度测距仪器，因具有方向性好、测距精度高等特点，获取的地表高程信息精度极高。根据卫星的设计指标，以及卫星应用实践情况来看，ZY3-03星激光测高点的高程中误差在平原地区优于0.5米、丘陵地区优于0.8米、山地地区优于1米。ZY3-03星是世界首颗同时搭载业务化激光测高仪和三线阵立体相机的测图卫星。

将激光测高点作为立体影像的高程控制点并开展联合测绘处理，来提升立体影像高程精度时，其核心步骤是获取激光测高点在立体影像上的准确像方坐标，若激光测高点在立体影像上的像点位置与实际像点位置存在较大偏差，在丘陵和山地地形区域(即地表起伏较大时)，等效于为激光测高点引入了额外的高程误差，将会严重影响联合测绘处理效果。而如何克服ZY3-03星激光测高点没有足印影像的影像，实现激光侧高点在立体影像上的准确定位，进而设计一种将激光测高点和三线阵立体影像进行联合几何处理的有效方法，实现卫星立体影像的高程精度大幅度提升的目标，已经成为了ZY3-03星(及类似卫星)影像实现全球范围(即稀少控制或无地面控制条件下)1:5万及更大比例尺立体测图应用的迫切需求，也将填补三线阵立体影像与激光测高点联合处理方面的技术空白。

发明内容

本发明主要解决的技术问题是：为了解决无地面控制(如境外区域无法获取地面控制点)条件下卫星影像的几何精度(尤其是高程精度)无法满足高精度测图应用等不足，充分利用激光测高点高程精度较高的特点，以及同轨获取激光测高点和正视影像之间平面相对误差较小的特性，提出了一种同轨激光测高点与三线阵立体影像联合测绘处理方法，能够大幅提升各类地形区域的三线阵立体影像的高程精度。

为解决上述技术问题，本发明的目的通过以下的技术方案来实现：

一种同轨激光测高点与三线阵立体影像联合测绘处理方法，包括如下步骤：

步骤1，以卫星成像摄影轨道为单位，将同一轨多景三线阵立体影像构建一个平差区域网，并采用基于像方仿射变换补偿有理函数模型RFM(Rational Function Model)的区域网平差模型；

步骤2，在各三线阵立体影像内部布设密集的连接点，在相邻的三线阵立体影像之间布设足量的公共连接点；

步骤3，利用同轨获取的激光测高点和正视影像相对平面误差较小等特性，根据激光测高点的经纬度坐标和高程值，通过正视影像的RFM，计算获取激光测高点在正视影像上的像点坐标；

步骤4，分别将前视和后视影像与正视影像开展高精度的同名点匹配，获取正视影像上的激光测高点像点在前视和后视影像上的同名点，亦即获取了激光测高点在前视影像和后视影像上的像点坐标；

步骤5，针对每个激光测高点像点，将影像RFM的像方仿射变换参数和激光测高点像点对应的地面经纬度坐标作为未知数，而将激光测高点的高程值作为已知量，逐点构建误差方程；针对每个连接点，将影像RFM的像方仿射变换参数和连接点对应的地面三维坐标作为未知数，逐点构建误差方程；

步骤6，对误差方程进行法化，得到法方程；采用最小二乘原理，进行整体平差解算，解算每张影像RFM的仿射变换参数、连接点对应的地面三维坐标改正数和激光测高像点对应的地面经纬度坐标改正数；

步骤7，每完成一次平差解算时，判断是否满足迭代收敛条件(即所有影像RFM的仿射变换参数中平移参数小于阈值)；如果满足，则平差结束；反之，则利用平差解算获取的改正数更新影像RFM、连接点对应地面三维坐标和激光测高像点对应的地面经纬度坐标，重复步骤5—步骤6，再次进行平差迭代解算，不断更新影像RFM的仿射变换参数、连接点对应的地面三维坐标和激光测高像点对应的地面经纬度坐标，直到满足迭代收敛条件为止。最后输出每张影像RFM的仿射变换参数；

步骤8，利用最后一次平差迭代计算得到的每张影像RFM的仿射变换参数,对RFM模型进行像方仿射变换补偿，为每张影像生成得到新的RFM。

与现有技术相比，本发明的一个或多个实施例可以具有如下优点：

(1)通过将激光测高点作为高程控制，针对立体影像开展区域网平差，大幅提升了三线阵立体影像的高程精度，使其满足1:5万乃至更大比例尺测图精度要求，使得全球范围无地面控制条件下高精度测图成为了可能。

(2)本发明充分利用同轨获取的激光测高点与正视影像之间平面相对误差较小的特性，获取的激光测高点在立体影像上像点位置的精度极高，有效避免了地表起伏较大区域(主要是丘陵和山地地形区域)，因激光测高点在立体影像上的量测位置与实际位置偏差较大，进而影响最终区域网平差精度。保证了丘陵和山地地形区域激光测高点的有效利用。

(3)利用激光测高点高程精度极高的优势，在测图应用中实现了对地面高精度高程控制点的替代，避免外业实地测量，大大降低工作成本。

附图说明

图1是同轨激光测高点与三线阵立体影像联合测绘处理方法流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合实施例及附图对本发明作进一步详细的描述。

如图1所示，为同轨激光测高点与三线阵立体影像联合测绘处理方法流程，包括以下步骤：

步骤1，以卫星成像摄影轨道为单位，将同一轨多景三线阵立体影像构建一个平差区域网，并采用基于像方仿射变换补偿有理函数模型(RFM)的区域网平差模型；

RFM是指利用有理多项式建立影像的像素坐标与其对应的地面点的大地坐标之间的数学映射关系，其基本方程定义如下：

式中，(X_n,Y_n,Z_n)和(r_n,c_n)分别为地面点坐标(X,Y,Z)和影像坐标(r,c)经过平移和缩放后的归一化坐标，取值在[-1,1]之间，归一化的目的是为避免计算过程中参数数值量级差别过大而引入舍入误差，以增强参数求解的稳定性。R₀,C₀为影像坐标正则化的平移参数、X₀,Y₀,Z₀为地面点坐标正则化的平移参数。R_s,C_s为影像坐标正则化的比例系数、X_s,Y_s,Z_s为地面点坐标正则化的比例系数。P_i(i＝1,2,3,4)分别表示一般多项式，式中各变量X_n,Y_n,Z_n的幂均不超过3次，所有变量的幂之和也不超过3次，其形式如下：

式中，a_ij(i＝1,2,3,4；j＝0,1,…,19)为有理多项式系数；

影像RFM的误差补偿采用像方仿射变换补充模型，则公式(1)中RFM描述的像点坐标(r,c)和地面点坐标(X,Y,Z)之间的关系被修正为：

式中，(Δr,Δc)为像点坐标(r,c)的系统误差像方补偿值，其值为：

式中，(a₀,a₁,a₂,b₀,b₁,b₂)表示补偿RFM系统误差的仿射变换参数。

步骤2：在各三线阵立体影像内部布设密集的连接点，在相邻的三线阵立体影像之间布设足量的公共连接点。

本步骤进一步包括：

步骤2.1，首先采用影像自动匹配连接点，常见的影像自动匹配算法有基于像方灰度的影像匹配算法(如最小二乘法)、基于物方的影像匹配算法、基于像方特征的跨接法影像匹配、金字塔多级影像匹配、SIFT等。要求区域网中每个三线阵立体模型内的连接点应分布均匀，且连接点数量应大于25个；相邻三线阵立体模型之间的连接点数量应大于5个；区域内90％的连接点维度(即该连接点连接影像的数量)应等于重叠区影像的总数量。

步骤2.2，当影像自动匹配的连接点分布不均匀、数量不足或维度不够时，应采用人工判读/量测方法补测连接点，人工判读/量测的连接点应位于影像清晰、特征明显、反差较大、易于转刺和量测的固定目标上。

步骤3：利用同轨获取的激光测高点和正视影像相对平面误差较小等特性，根据激光测高点的经纬度坐标和高程值，通过正视影像的RFM，计算获取激光测高点在同轨获取的正视影像上的像点坐标。

由于激光测高点与正视影像的获取时间几乎相同，因此两者所受的外部误差的类型和量级也基本相同，且激光测高仪的指向角度和正视相机主光轴指向角均垂直对地，因此外部误差对两者平面精度的影响效果也非常相似。这些因素使得激光测高点和同轨的正视影像的平面相对误差较小，必然小于激光测高点的地面足印半径，因此可以直接将激光测高点的地面经纬度坐标和高程值，代入正视影像RFM公式(公式1)，计算激光测高点在正视影像上的像点坐标。

步骤4，分别将前视和后视影像与正视影像开展高精度的同名点匹配，获取正视影像上的激光测高点像点在前视和后视影像上的同名点，亦即获取了激光测高点在前视影像和后视影像上的像点坐标。

本步骤进一步包括：

步骤4.1，计算正视影像上激光测高点像点坐标位置(整像素位置)的SIFT(Scale-invariant feature transform，尺度不变特征变换)特征，SIFT作为图像处理领域的一种局部特征描述子，在一定程度上这种描述具有尺度不变性，它对物体的尺度变化，刚体变换，光照强度和遮挡都具有较好的稳定性。保存获取的SIFT特征描述符；

步骤4.2，根据激光测高点的地面经纬度坐标和高程值，利用前视影像的RFM，计算激光测高点在前视影像上的大致像素位置，并以此为中心，向四周外扩一定像素范围作为匹配的搜索范围；计算前视影像搜索范围内每个像素的SIFT特征，并保存其特征描述符；

步骤4.3，使用快速最近邻逼近搜索函数库(Fast Approximate NearestNeighbor Search Library)，搜索步骤4.2获取的多个特征描述符中，与正视影像上激光测高点像点的SIFT特征描述符最近邻值，获取该特征描述符对应的前视影像上的像点坐标，亦即前视影像上激光测高点像点的整像素坐标值。

步骤4.4，以前视影像上激光测高点像点的整像素坐标位置为中心，与正视影像上激光测高点像点开展最小二乘匹配，获取前视影像上激光测高点像点的精确坐标值。

步骤4.5，重复步骤4.2-步骤4.4的方法，同样获取后视影像上激光测高点像点的精确坐标值。

步骤5，针对每个激光测高点像点，将影像RFM的像方仿射变换参数和激光测高点像点对应的地面经纬度坐标作为未知数，而将激光测高点的高程值作为已知量，逐点构建误差方程；针对每个连接点，将影像RFM的像方仿射变换参数和连接点对应的地面三维坐标作为未知数，逐点构建误差方程。

本步骤进一步包括：

步骤5.1，确定连接点对应的地面三维坐标初值、激光测高点像点对应的地面经纬度坐标初值和影像RFM仿射变换参数初值；

其中，通过RFM模型对连接点进行直接前方交会得到地面点坐标作为连接点对应的地面三维坐标初值。通过RFM模型对激光测高点像点进行直接前方交会得到地面点经纬度坐标作为激光测高点像点对应的地面经纬度坐标初值。

本实施例中，仿射变换模型初值一般设定为

步骤5.2，设置连接点和激光测高点的误差方程如下：

式中，第一行误差方程为立体影像上连接点列出的误差方程，将影像RFM的仿射变换参数和连接点对应的地面三维坐标作为未知数；第二行误差方程为激光高点在立体影像上的像点坐标列出的误差方程，将影像RFM的仿射变换参数和激光测高点在立体影像上像点对应的地面经纬度坐标作为未知数。V₁＝[v_R1 v_C1]^T和V₂＝[v_R2 v_C2]^T分别为连接点和激光测高点在影像上行、列坐标观测值残差向量；t＝[Δa₀ Δa₁ Δa₂ Δb₀ Δb₁ Δb₂]^T为影像RFM的像方坐标系统误差的仿射变化补偿参数的改正数向量；x₁＝[ΔX_tie ΔY_tie ΔZ_tie]^T为连接点对应地面大地坐标的改正数向量；x₂＝[ΔX_laser ΔY_laser]^T为激光测高点像点对应的地面经纬度坐标的改正数向量；

是未知数t对应的偏导数系数矩阵；

分别是未知数x₁和x₂对应的偏导数系数矩阵；

为利用初值代入观测方程计算得到的常数项，(r⁰,c⁰)为利用未知数的近似值通过RFM计算出的连接点像点行、列坐标；P_i(i＝1,2)为各观测值的权矩阵，通常各类观测值的权值由他们的先验信息确定，一般可取观测值10倍先验标准差，并在每次平差迭代计算后重新计算各类观测值的权。

步骤6，对误差方程进行法化，得到法方程；采用最小二乘原理，进行整体平差解算，解算每张影像RFM的仿射变换参数、连接点对应的地面三维坐标改正数和激光测高点像点对应的地面经纬度坐标改正数。

本步骤进一步包括：

步骤6.1：根据最小二乘平差原理，对误差方程进行法化，得到的法方程形式如下：

记为：

步骤6.2：由于影像上连接点数量较多，法方程中待求解的未知数个数过多，通过对误差方程进行变换消去X，只求解其中的仿射变换参数的改正数，然后通过再次前方交会的方式更新连接点对应的地面三维坐标，以提升解算的效率，将上式(8)变化为如下形式：

Nt＝G (9)

其中，

N中

为2阶对角矩阵，在解算时通过对每个小矩阵求逆的方式得到；这一步可以多线程并行处理，能够大大节省解算时间，提升计算效率。

步骤6.3：求解出

之后便能够分别得到N、G，针对公(9)的方程，利用数学中的共轭梯度下降法进行迭代求解，在两次求解得到的t的差值小于设定的阈值，或者求解次数超过设定的次数之后结束迭代，输出得到最终的t，也就是仿射变换参数改正数。

步骤7，每完成一次平差解算时，判断所有影像RFM的仿射变换参数中平移参数a₀,b₀是否小于阈值(本实施例为0.1个像元pixel，但是不限于此)；如果满足，则平差结束；反之，则利用平差解算获取的改正数更新影像RFM、连接点对应地面三维坐标和激光测高点像点对应的地面经纬度坐标，重复步骤5-步骤6，再次进行平差迭代解算，不断更新影像RFM的仿射变换参数、连接点对应的地面三维坐标和激光测高点像点对应的地面经纬度坐标，直到满足迭代收敛条件为止。最后输出每张影像RFM的仿射变换参数。

步骤8，利用最后一次平差迭代计算得到的每张影像的RFM仿射变换参数,对影像的RFM进行像方仿射变换补偿，为每张影像生成得到新的RFM。

本步骤进一步包括：

步骤8.1：利用影像原有RFM及仿射变换补偿参数，计算影像四个角点对应的地面范围；根据美国地质调查局提供的全球1km分辨率DEM(Global 30-arc-second DigitalElevation Model)，计算该地区的最大最小椭球高。然后，在高程方向以一定的间隔分层(为了防止设计矩阵状态恶化，一般高程方向分层的层数超过2，本实施例平面分为10层，但是不限于此)，在平面上，以一定的格网大小建立地面规则格网(本实施例平面分为15×15格网，就是将该影像对应影像范围分成15×15的格子，共有16×16个格网点，但是不限于此)，生成控制点地面坐标，最后利用影像原有RFM的反变换及仿射变换补偿参数，计算控制点的影像坐标。加密控制格网和层，建立独立检查点。

步骤8.2：采用最小二乘平差原理求解RFM参数，首先要将RFM的方程(公式1)变形为：

则得线性化误差方程为：

V＝Bx-l,W (11)

式中，

x＝[a_1ia_2j a_3i a_4j]^T，W为权矩阵。

根据最小二乘平差原理，可以求解:

x＝(B^TB)^-1B^Tl (12)

经过变形的RFM模型形式，平差的误差方程为线性模型，因此在求解RFM参数过程中不需要初值。

步骤8.3：用求解的RFM参数来计算检查点对应的影像坐标，通过影像原有RFM及仿射变换补偿参数计算的检查点影像坐标的差值来评定求解RFM参数的精度。并输出新的RFM参数。

虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (5)

1.一种同轨激光测高点与三线阵立体影像联合测绘处理方法，其特征在于，所述方法包括以下步骤：

步骤1以卫星成像摄影轨道为单位，将同一轨多景三线阵立体影像构建一个平差区域网，并采用基于像方仿射变换补偿有理函数模型RFM的区域网平差模型；

步骤2在各三线阵立体影像内部布设连接点，在相邻的三线阵立体影像之间布设公共连接点；

步骤3根据激光测高点的经纬度坐标和高程值，通过正视影像的RFM，计算获取激光测高点在同轨获取的正视影像上的像点坐标；

步骤4分别将前视和后视影像与正视影像开展高精度的同名点匹配，获取正视影像上的激光测高点像点在前视和后视影像上的同名点，即获取了激光测高点在前视影像和后视影像上的像点坐标；

步骤5针对每个激光测高点像点，将影像RFM的像方仿射变换参数和激光测高点像点对应的地面经纬度坐标作为未知数，将激光测高点的高程值作为已知量，逐点构建误差方程；针对每个连接点，将影像RFM的像方仿射变换参数和连接点对应的地面三维坐标作为未知数，逐点构建误差方程；

步骤6对误差方程进行法化，得到法方程；采用最小二乘原理，进行整体平差解算，解算每张影像RFM的仿射变换参数、连接点对应的地面三维坐标改正数和激光测高像点对应的地面经纬度坐标改正数；

步骤7每完成一次平差解算，判断平差是否满足收敛条件，若满足，则平差结束；反之，利用平差解算获取的改正数更新影像RFM、连接点对应地面三维坐标和激光测高像点对应的地面经纬度坐标，重复步骤5-步骤6，再次进行平差迭代解算，不断更新影像RFM的仿射变换参数、连接点对应的地面三维坐标和激光测高像点对应的地面经纬度坐标，直到满足收敛条件为止；最后输出每张影像RFM的仿射变换参数；

步骤8利用最后一次平差迭代计算得到的每张影像RFM的仿射变换参数,对RFM模型进行像方仿射变换补偿，为每张影像生成得到新的RFM；

所述步骤4中，分别匹配正视影像上的激光测高点像点在前视和后视影像上的同名点具体包括以下步骤：

步骤4.1计算正视影像上激光测高点像点坐标位置的尺度不变特征变换SIFT特征，并保存特征描述符；

步骤4.2根据激光测高点的地面经纬度坐标和高程值，利用前视影像的RFM，计算激光测高点在前视影像上的大致像点位置，并以此为中心，向四周外扩一定像素范围作为匹配的搜索范围；计算前视影像搜索范围内每个像素的SIFT特征，并保存其特征描述符；

步骤4.3使用快速最近邻逼近搜索函数库搜索步骤4.2获取的多个特征描述符中与正视影像上激光测高点像点的SIFT特征描述符最近邻值，获取该特征描述符对应的前视影像上的像点坐标，亦即前视影像上激光测高点像点的整像素坐标值；

步骤4.4以前视影像上激光测高点像点的整像素坐标位置为中心，与正视影像上激光测高点像点开展最小二乘匹配，获取前视影像上激光测高点像点的精确坐标值；

步骤4.5重复执行步骤4.2-步骤4.4，同样获取后视影像上激光测高点像点的精确坐标值。

2.如权利要求1所述的同轨激光测高点与三线阵立体影像联合测绘处理方法，其特征在于，所述步骤1中，RFM指利用有理多项式建立影像的像素坐标与其对应的地面点的大地坐标之间的数学映射关系，其基本方程定义如下：

式中，(X_n,Y_n,Z_n)和(r_n,c_n)分别为地面点坐标(X,Y,Z)和影像坐标(r,c)经过平移和缩放后的归一化坐标，取值在[-1,1]之间；R₀,C₀为影像坐标正则化的平移参数、X₀,Y₀,Z₀为地面点坐标正则化的平移参数；R_s,C_s为影像坐标正则化的比例系数、X_s,Y_s,Z_s为地面点坐标正则化的比例系数；P_i(i＝1,2,3,4)分别表示一般多项式，式中各变量X_n,Y_n,Z_n的幂均不超过三次，所有变量的幂之和也不超过三次，其形式如下：

式中，a_ij(i＝1,2,3,4；j＝0,1,…,19)为有理多项式系数；

影像RFM的误差补偿采用像方仿射变换补充模型，则公式(1)中RFM描述的像点坐标(r,c)和地面点坐标(X,Y,Z)之间的关系被修正为：

式中，(Δr,Δc)为像点坐标(r,c)的系统误差像方补偿值，其值为：

式中，(a₀,a₁,a₂,b₀,b₁,b₂)表示补偿RFM系统误差的仿射变换参数。

3.如权利要求1所述的同轨激光测高点与三线阵立体影像联合测绘处理方法，其特征在于，所述步骤3是将激光测高点的地面经纬度坐标和高程值，代入同轨获取的正视影像RFM公式1，计算激光测高点在正视影像上的像点坐标。

4.如权利要求1所述的同轨激光测高点与三线阵立体影像联合测绘处理方法，其特征在于，所述步骤5中，连接点和激光测高点的误差方程如下：

式中，第一行误差方程为立体影像上连接点列出的误差方程，将影像RFM的仿射变换参数和连接点对应的地面三维坐标作为未知数；第二行误差方程为激光高点在立体影像上的像点坐标列出的误差方程，将影像RFM的仿射变换参数和激光点地面经纬度坐标作为未知数；V₁＝[v_R1 v_C1]^T和V₂＝[v_R2 v_C2]^T分别为连接点和激光点在影像上行、列坐标观测值残差向量；t＝[Δa₀ Δa₁ Δa₂ Δb₀ Δb₁ Δb₂]^T为影像RFM的像方坐标系统误差仿射变化补偿参数的改正数向量；x₁＝[ΔX_tie ΔY_tie ΔZ_tie]^T为连接点对应地面大地坐标的改正数向量；x₂＝[ΔX_laser ΔY_laser]^T为激光测高点地面经纬度的改正数向量；

是未知数t对应的偏导数系数矩阵；

和

分别是未知数x₁和x₂对应的偏导数系数矩阵；

为利用初值代入观测方程计算得到的常数项，(r⁰,c⁰)为利用未知数的近似值通过RFM计算出的连接点像点行、列坐标；P_i(i＝1,2)为各观测值的权矩阵，通常各类观测值的权值由他们的先验信息确定，一般可取观测值10倍先验标准差，并在每次平差迭代计算后重新计算各类观测值的权。

5.根据权利要求1所述的同轨激光测高点与三线阵立体影像联合测绘处理方法，其特征在于，所述步骤6具体包括以下步骤：

步骤6.1根据最小二乘平差原理，对误差方程进行法化，得到的法方程形式如下：

记为：

步骤6.2：由于影像上连接点数量较多，法方程中待求解的未知数个数过多，通过对误差方程进行变换消去X，只求解其中的仿射变换参数的改正数，然后通过再次前方交会的方式更新连接点对应的地面三维坐标，以提升解算的效率，将上式(8)变化为如下形式：

Nt＝G (9)

其中，

N中

为2阶对角矩阵，在解算时通过对每个小矩阵求逆的方式得到；这一步可以多线程并行处理，能够大大节省解算时间，提升计算效率；

步骤6.3：求解出

之后便能够分别得到N、G，针对公式(9)的方程，利用数学中的共轭梯度下降法进行迭代求解，在两次求解得到的t的差值小于设定的阈值，或者求解次数超过设定的次数之后结束迭代，输出得到最终的t，也就是仿射变换参数改正数。

Images