CN113358091B - 一种利用三线阵立体卫星影像生产数字高程模型dem的方法 - Google Patents

Abstract

本发明公开了一种利用三线阵立体卫星影像生产数字高程模型方法，包括：将生产区域内多度重叠覆盖的遥感卫星三线阵立体影像构建平差区域网并完成区域网平差，实现所有影像的高精度定向；随后根据任意两张影像之间重叠率和同名光线交会角等指标值，将多度重叠覆盖的三线阵立体影像构建多对双像立体影像，并分别采用影像密集匹配算法提取数字高程模型DEM；最后根据双像立体影像的交会角度、相对精度，以及构成影像的时相差异等指标，计算各个双像立体影像的权重值，并基于该权重值，将区域内所有的双像立体影像提取的DEM进行加权融合，生成最终DEM产品。本发明充分利用三线阵立体像的信息优势，能够有效提升DEM的匹配成功率、分辨率和产品质量。

Description

一种利用三线阵立体卫星影像生产数字高程模型DEM的方法

技术领域

本发明涉及测绘地理信息产品中数字高程模型产品生产技术领域，尤其涉及一种基于多度重叠覆盖的遥感卫星三线阵立体卫星影像的数字高程模型生产方法。

背景技术

测绘地理信息资源是一个国家基础性、战略性信息资源，在政府管理决策、产业发展、人民生活等方面发挥着越来越重要的作用。数字高程模型(Digital Elevation Model，DEM)是通过有限的地形高程数据实现对地面地形的数字化模拟，是用一组有序数值阵列形式表示地面高程的一种实体地面模型。DEM做为一种重要的测绘地理信息成果形式，广泛应用于测绘、水文、气象、地貌、地质、土壤、工程建设、通讯、军事等国民经济和国防建设以及人文和自然科学领域。

DEM产品的生产方式多样。近年来，采用遥感卫星光学立体影像，通过航天摄影测量手段获取DEM，已经成为大范围高质量DEM产品的重要生产方法。该方法的主要步骤按流程依次为：立体影像区域网平差、立体影像立体匹配提取DEM、DEM编辑等。所谓卫星立体影像是由两张或多张从不同成像视角拍摄获取的同一地区的遥感卫星平面影像构建而成。目前利用由两张影像构建的双像立体影像开展DEM生产的理论、技术和方法已经逐渐趋于成熟。然而近年来，能够获取三视立体影像(一般指由垂直指向天底方向相机获取的“正视影像”，以及由指向前进方向的前方和后方相机获取的“前视影像”和“后视影像”构建的立体影像)的测绘遥感卫星相继发射(如日本的ALOS卫星、中国的资源三号系列卫星、天绘一号系列卫星等)，三视立体影像已成为重要的立体测绘应用影像源，且拥有比双像立体更多的信息优势。

如何充分利用三视立体影像相比双像立体影像多出一张影像的信息优势，并且当三视立体影像多度重叠覆盖时，如何充分利用影像信息冗余，实现更高质量的DEM产品生产，已经成研究热点。相比于利用双像立体影像生产DEM，利用多度重叠覆盖的三视立体影像生产DEM的主要差异体现在“立体影像立体匹配提取DEM”环节，目前国内外在这方面的研究和应用中，一种方法是直接挑选基高比最大的两张影像构建一个双像立体影像，开展DEM产品生产，放弃并浪费了多度重叠覆盖的三视立体影像的信息优势。另一种方法是直接针对三视(或多视)立体影像采用多基线、多光线、多视点影像匹配方法来实现DEM提取，该方法的优点是能将双像立体(两度重叠)影像匹配方法中的“病态解”转化为“确定解”，较好解决了影像中相似性特征、地表断裂特征的误匹配等问题。但是该方法的算法复杂性较高、成熟度较低，不容易检测遮挡等可疑区域，当构成立体影像的不同影像之间时相、辐射和几何差异较大时，对匹配效果影响极大，匹配生成的DEM质量不高。总之，目前利用多度重叠覆盖的三视(或多视)立体影像生成DEM的方法均存在较大局限性，急需寻找其他新的技术方法，以便能充分利用三视(或多视)立体像的信息优势，更加有效实现高质量DEM产品生产。

发明内容

为解决上述技术问题，本发明的目的是提供一种利用三线阵立体卫星影像生产数字高程模型DEM的方法，该方法利用多度重叠覆盖的三线阵立体影像的冗余信息优势，提出了一种利用遥感卫星三线阵立体影像加权融合生产数字高程模型产品方法，能够有效减少因地形遮挡或阴影等导致局部影像上DEM无法成功匹配情况，并大幅提升DEM分辨率以及地形地貌细节复杂等区域的DEM匹配质量。

本发明的目的通过以下的技术方案来实现：

一种利用三线阵立体卫星影像生产数字高程模型DEM的方法，包括：

步骤A将生产任务区域内多度重叠覆盖的遥感卫星三线阵立体影像构建平差区域网并完成区域网平差，实现区域内所有卫星影像的高精度相对和绝对定向；

步骤B计算所有拥有重叠关系的任意两张影像之间的重叠率和同名光线交会角，当影像重叠率和同名光线交会角同时符合设定的门槛阈值时，将其构建为一个双像立体影像，则多度重叠覆盖的三线阵立体影像可拆分构建为多个多度重叠覆盖的双像立体影像；

步骤C针对各个双像立体影像，首先采用投影轨迹法生成近似核线立体影像，然后采用半全局密集匹配算法生成视差图，最后利用近似核线立体影像成像几何模型进行交会获取每个像素的物方坐标，生成任务区域内的物方三维点云，并通过点云栅格化生成规则格网DEM产品；

步骤D综合考虑双像立体影像的同名光线交会角度、定向精度、时相差异等因素，设置各个双像立体影像提取的DEM在最终加权融合时的权重值w；

步骤E利用获取的权重值，对所有双像立体影像提取的DEM进行加权融合，生成整个任务区域的融合后DEM产品；

步骤F对融合后DEM产品，进行编辑，修正影像自动匹配的错误区域、剔除人工建筑和人工地物高程，形成最终的DEM产品。

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

三视立体影像可以获取区域内更多视角的地表信息，本发明采用三线阵立体影像开展DEM生成，可以有效降低了因地形起伏遮挡等导致的局部区域无法成功匹配DEM的风险，提升了整个区域的DEM匹配成功率。

充分利用多度重叠覆盖立体影像拥有的信息冗余优势，通过将任务区域内的同轨和异轨三线阵立体影像拆分为多个双像立体影像(如同轨影像构建的前-后视立体影像，前-正视立体影像，正-后视立体影像，以及异轨影像构建的前-后视立体影像，前-正视立体影像，正-后视立体影像等)后，分别匹配DEM，再进行加权融合生成一个总的DEM，其效果类似于“图像超分辨率技术”，因此可以有效提升DEM产品的分辨率，也可以提高地形地貌细节复杂区域的DEM匹配质量。

相对于直接针对多视立体影像采用多基线、多光线、多视点匹配方法提取DEM，本发明的方法更加简单可靠，匹配效果更好，可以有效避免云斑、阴影等对匹配的影响，并且方法的适用性和扩展性更强，如可以轻易扩展到多度重叠的二、三视立体影像混合的DEM生产，影像重叠的度数不会影响算法的复杂性。

附图说明

图1是利用三线阵立体卫星影像生产数字高程模型DEM的方法流程图。

具体实施方式

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

如图1所示，为利用三线阵立体卫星影像生产数字高程模型DEM的方法，包括以下步骤：

步骤1)将生产任务区域内多度重叠的遥感卫星三线阵立体影像构建平差区域网，量测密集连接点并根据需要量测适量像控点，选取合适的区域网平差模型，在全区域内统一进行平差迭代计算，直至满足平差精度要求，获取所有影像成像几何模型误差补偿参数，输出新的影像成像几何模型，完成区域内卫星影像的高精度相对和绝对定向。

上述步骤1)进一步包括：

步骤1.1)根据DEM生产任务区域范围，收集获取区域内遥感卫星三线阵立体影像，为保证区域网的连接性，区域内不应该有影像覆盖漏洞区域，各个三线阵立体影像之间的重叠区域不应小于20％影像幅宽。为了更好地展现本发明方法的效果，可以收集区域内不同时间获取的三线阵立体影像，使区域内的三线阵立体影像形成多度重叠覆盖。最后将这些影像构成平差区域网。

步骤1.2)在区域网内的立体影像上量测连接点(注：连接点指能够构建立体模型或建立相邻立体模型之间连接关系的同名像点)，连接点可采用影像自动匹配和人工判读/量测两种方法获取。具体为：

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

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

步骤1.3)根据区域网平差精度需要，在平差区域网内的立体影像上布设适量的高精度像控点(注：像控点是指位于影像特定位置和特定目标上，具有成图坐标系中坐标信息的控制点，其中用于解算影像定向参数而布设在特定位置的必要数量的像控点称为基本定向点；用于检查成果正确性的像控点称为检查点)。像控点布设时以一个三线阵立体影像为基本布点单元，布设数量和布设位置取决于最终的平差精度需求、控制点收集情况等，通常情况下，在每个基本布点单元的四个角点分别布设1个控制点、中央位置布设1个检查点。控制点应尽量布设在相邻立体影像的重叠范围内同一目标点上，控制点尽量共用，利于增加重复观测条件，提高平差精度。具体为：

(1)当采用已有高精度正射纠正影像(DOM)、数字表面模型(DSM)或数字高程模型(DEM)数据作为控制资料时，采用影像匹配提取像控点，即首先将三线阵立体影像与高精度DOM进行影像同名点匹配，经过粗差剔除后得到的同名目标作为像控点，并通过高精度DOM、DSM或DEM获取该像控点的物方平面位置和高程值。

(2)当采用外业GNSS点作为控制资料，或者前述采用影像匹配提取像控点分布不均匀、数量不足时，采取人工目视判读法补充获取像控点，即在三线阵立体影像需要布设像控点位置，通过人工目视识别的方法，选取清晰、明显的像点目标作为像控点，并从已有DOM、DSM、DEM或外业GNSS点获取像控点的物方平面坐标和高程值。

步骤1.4)选取像方仿射变换模型补偿的基于RFM的区域网平差模型，针对所有连接点和像控点分别构建误差方程，在全区域内统一进行平差迭代计算，直至满足平差精度要求，解算出每张影像RFM的像方补偿参数，并基于此更新并输出所有影像的RFM。具体为：

步骤1.4.1)选择RFM作为立体影像区域网平差时的成像几何模型，选取像方仿射变换模型作为立体影像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次。

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

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

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

步骤1.4.2)给定补偿每张影像RFM系统误差的仿射变换参数初值和所有连接点地面坐标的初值。

本实施例中，所有连接点地面坐标初值通过RFM对连接点进行直接前方交会得到地面点坐标(lat,lon,h)作为初值。

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

步骤1.4.3)将补偿影像RFM的仿射变换模型参数和连接点对应的地面点坐标一并求解，并对其求偏导，即可逐点构建连接点和控制点的误差方程式。

V＝Bt+AX-l,P (5)

其中，V＝[v_x v_y]^T表示像点行和列坐标观测值的残差向量；t＝[Δa₀ Δa₁ Δa₂Δb₀ Δb₁ Δb₂]^T为六个仿射变换参数的改正数向量，X＝[Δlat Δlon Δh]^T为连接点地面坐标的改正数向量；

为未知数t的系数矩阵；

为未知数X的系数矩阵；

为初值计算得到的常数项；P为权矩阵。

控制点可认为其地面坐标是准确的，不需要对其地面坐标进行改正，因此其误差方程中X为零向量。

步骤1.4.4)根据最小二平差原理，对误差方程进行法化，得到法方程；采用边法化边消元的循环分块法解求改化法方程，确定出补偿每张影像RFM的仿射变换模型参数的改正数和每个连接点的地面坐标的改正数。进一步描述如下：

步骤1.4.4.1)列法方程形式如下：

记为：

由于影像上连接点数量较多，如果直接对上式进行求解的话，求解的未知数个数过多，因此通过对误差方程进行变换消去X，只求解其中的仿射变换参数的改正数，然后通过再次前方交会的方式更新连接点地面点坐标，提升解算的效率。

步骤1.4.4.2)将上式(7)变化为如下形式：

Nt＝G (8)

其中，

N中

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

步骤1.4.4.3)求解出

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

步骤1.4.5)用这些改正数修正每张影像RFM的仿射变换参数和连接点的地面坐标，并重复步骤1.4.3—1.4.4，通过迭代过程不断更新连接点地面坐标和影像RFM的仿射变换参数，当仿射变换参数中的平移参数a₀,b₀小于阈值时(本实施例为0.1个像元pixel)时，平差迭代结束；当不满足预设条件时，继续迭代至满足条件为止。

步骤1.4.6)平差解算完成后，利用最终的影像RFM仿射变换参数改正数，对每张影像的RFM进行像方仿射变换补偿，生成得到新的RFM参数文件。此时立体影像的几何精度已经大幅度提升，可以进行后续的DEM提取步骤。

步骤2)计算所有拥有重叠关系的任意两张影像之间的重叠率和同名光线交会角度，当影像重叠率和同名光线交会角度同时符合设定门槛阈值时，将其构建为一个双像立体影像，则多度重叠覆盖的三线阵立体影像可以构建多个双像立体影像。

上述步骤2)进一步包括：

步骤2.1，利用卫星影像RFM的反解公式(见公式9)，分别计算任务区域内每张影像四个角点的经纬度坐标(lon_n,lat_n)，进而获取每张影像的地理范围多边形

利用多边形面积计算公式(见公式10)求解其面积S_Pi；再利用多边形求交算法依次计算区域内所有存在相交关系的任意两张影像之间的相交多边形

同样利用多边形面积计算公式(见公式10)求解其面积S_U；最后利用公式11计算区域内所有存在相交关系的任意两张影像之间的重叠率M_n，(n＝1,2,...,k)。

其中，各符号含义参见公式1和公式2，Z_n取值为区域内的平均高，可由外部获取。

其中lon_i+1＝lon₁,lat_i+1＝lat₁。

其中S_P1，S_P2分别为两个多边形的面积，S_U为这前述两个多边形相交多边形的面积。

步骤2.2)依次针对任务区域内所有存在重叠关系的任意两张影像，在其上任意选取1个同名像点(也可以直接从步骤1.2中获取的连接点中选取)，计算该同名像点对应地面目标的成像交会角β(亦即这两张影像之间的同名光线空间交会角度)，具体如下：

步骤2.2.1)假定Z_n＝0和Z_n＝100，采用公式(9)，分别计算影像1上同名像点对应的地面坐标(lon₀,lat₀,0)和(lon₂,lat₂,100)，并分别将他们转换到地心直角坐标系中，记为P₀(X₀,Y₀,Z₀)和P₂(X₂,Y₂,Z₂)，可得影像1上该同名像点对应地面目标的入射向量为

同理可得影像2上该同名像点对应地面目标的入射向量为α₂。

步骤2.2.2)采用下式计算这两张影像之间的同名光线空间交会角度。

步骤2.2.3)根据影像源情况，设置一个合理的影像重叠率指标值M⁰(本实施例为40％，但是不限于此)和影像同名光线空间交会角度指标值β⁰(本实施例为10°，但是不限于此)，依次对区域内所有存在重叠关系的任意两张影像进行资格判断，当满足下式时，即将其构建为一个双像立体影像。

M≥M⁰&β≥β⁰ (13)

其中M为任意两张影像的重叠率，β为这两张影像的交会角。

此时，同轨获取的前、后视影像，前、正视影像，正、后视影像，以及不同轨获取的前、后视影像，前、正视影像，正、后视影像都将可能构建双像立体。

步骤3)分别针对各个双像立体影像，首先采用投影轨迹法生成近似核线立体影像，然后采用半全局密集匹配算法生成视差图，最后利用近似核线立体影像成像几何模型进行交会获取每个像素的物方坐标，生成任务区域内的物方三维点云，并通过点云栅格化生成规则格网DEM产品。

上述步骤3)进一步包括：

步骤3.1)将构成双像立体的两张影像分别采用投影轨迹法生成近似核线影像,消除了立体影像上下视差，将二维匹配搜索问题转化为一维搜索问题，从而降低影像匹配算法复杂性和计算量并提升匹配的效率和正确率。具体包括：

步骤3.1.1)根据双像立体影像的相互位置关系，选取一张影像作为左影像，另一张作为右影像，假设影像行数为L,列数为S。

步骤3.1.2)左影像第i行(i＝1,2,3,···,S)上选取一点m(S_i,L_i)，m和左摄影中心S(X_S,Y_S,Z_S)形成一条成像光线，在该光线上任意选取十个点a_j(j＝1,2,3,···,10)，并根据区域平均高程值赋予其高程值h_j，由RFM的正反算公式计算a_j在右影像上的坐标(S_aj,L_aj)。

步骤3.1.3)根据a_j在右影像上的坐标(S_aj,L_aj)采用最小二乘线性拟合法拟合出右核线。

步骤3.1.4)在右核线上选取若干点，按照步骤3.1.2)和步骤3.1.3)可以拟合出左核线，则m(S_i,L_i)必然位于左核线上。

步骤3.1.5)通过上述方法得到左右影像近似同名核线，并利用一维灰度线性内插进行重采样，得到左右近似核线影像，最后还需分别计算生成左右近似核线影像的RFM。

步骤3.2)针对立体核线影像，采用半全局密集匹配算法，开展密集匹配，获取视差图。具体包括：

步骤3.2.1)基于Census变换法计算匹配代价影像各像素的匹配代价。

采用Census变换将图像像素的灰度值编码成二进制码流，以此来获取邻域像素灰度值相对于中心像素灰度值的大小关系，变换公式如下：

式中，p为某一像点位置；q为以像点p为中心的邻域窗口中的其他像素，N_p表示中心像素p的邻域窗口；

为比特位的逐位连接运算；I(*)表示像素点*处的灰度值；ξ运算定义如下：

通过计算左右影像对应的两个像素的Census变换值的汉明(Hamming)距离，获取该像素的匹配代价：

C(p,d):＝Hamming(C_sl(p),C_sr(d)) (16)

式中，d为其他影像上与像素p对应的像素。

步骤3.2.2)采用动态规划的方法，合并至少8个方向的一维动态规划结果实现二维的代价积聚。

计算像素p沿某一路径r的匹配代价聚合值：

式中，p为像素位置；d为像素p在立体影像上的左右视差值；L_r(p,d)为像素p路径r上在视差d处的代价积聚值；C(p,d)为像点p视差为d时的匹配代价；右边等式第2项表示像素p在r方向上的前一个像素代价累计的最小值；最后一项是为了防止数值过大，对最优路径的选择没有影响；P₁和P₂分别是当前像素p在路径r方向上前一个像素的视差差异为1的惩罚系数和视差差异大于1的惩罚系数，且P₂>P₁。

将像素p所有方向的匹配代价聚合值相加得到像素的最终的匹配代价，亦即聚合代价立方体：

步骤3.2.3)对于每一个像素，从上式(18)选择聚合代价最小的视差作为最优视差，从而生成最终的视差图D：

D＝min_dS(p,d) (19)

步骤3.2.4)为了获得更加精细平滑的视差图，对计算得到的视差图进行进一步的优化，包括利用左右一致性检查剔除粗差，使用二次曲线拟合相邻匹配代价进行亚像素插值，基于快速双边算子进行平滑处理等。

步骤3.3)根据视差图，获取核线影像同名点，通过立体核线影像的空间前方交会得到密集的地面三维点云数据，剔除粗差后，根据DEM要求的分辨率规格，将高程点云数据进行规则内插，生成该双像立体影像提取的DEM数据。

步骤4)综合考虑双像立体影像的交会角度、定向精度、时相差异等因素，通过公式20，获取各个双像立体影像提取的DEM在构建最终融合DEM时的权重值w。

其中β为双像立体影像同名光线交会角度，β₀为交会角度基准(一般设置为90°)；λ_time为双像立体影像是否为同轨成像的权重分，t_left为第一张影像成像时间，t_right为第二张影像成像时间，均精确到小时；σ为双像立体影像的绝对高程中误差(单位为米)，σ₀为立体影像绝对高程中误差基准，可根据实际情况设置，本实施例中，由于采用的是区域网平差，可以认为所有的双像立体影像的绝对高程中误差一致，则上式可以简化为：

步骤5)利用上一步骤获取的权重值，针对所有双像立体影像提取的DEM，进行加权融合，生成任务区域的融合后DEM产品。

上述步骤5)进一步包括：

步骤5.1，计算所有双像立体影像提取DEM的地理范围，进而获取融合后DEM的地理范围。根据融合后DEM的分辨率需求，计算出融合后DEM像元阵列的行列数量，以及各个像元的地理坐标。具体为：

步骤5.1.1，根据公式(22)，分别计算各个双像立体影像提取DEM的四个角点的地理坐标，求取所有双像立体影像提取DEM地理覆盖范围求并集，亦即获得了融合后DEM的地理范围。

其中(D_C,D_R)为地面坐标(可以是经纬度，也可以是投影坐标，取决于DEM采用的平面坐标系统)，(c,r)为DEM上像元的坐标，a，b为偏移参数，d_c，d_r为DEM上行和列方向上的分辨率。

步骤5.1.2)当给定融合后DEM的行和列分辨率为d_r和d_c(一般d_r＝d_c)时，通过公式(23)可计算融合后DEM像元阵列的总行数R和总列数C：

其中D_{max_c}和D_{min_c}是融合后DEM地理范围在东西方向最大和最小值，D_{max_r}和D_{min_r}是融合后DEM地理范围在南北方向最大和最小值。

此外，可知融合后DEM的偏移参数a＝D_{min_c}，b＝D_{min_r}。

步骤5.2)针对融合后DEM上的每一个像元(r′,c′)，利用公式(22)计算其对应的地面坐标(D_C,D_R)，搜寻覆盖地面坐标(D_C,D_R)的所有双像立体影像提取的DEM，获取(D_C,D_R)在各个双像立体提取DEM上的像素位置，通过内插获取该像素位置的高程值，最后将这些高程值进行加权相加，作为融合后DEM上像元(r′,c′)的高程值。具体为：

步骤5.2.1)采用公式(24)获取地面坐标(D_C,D_R)在各个双像立体提取DEM上的像元位置(r_i,c_i)：

步骤5.2.2)获取双像立体提取DEM上像元位置(r_i,c_i)处对应高程值h_i，若(r_i,c_i)不是整像素值，可采用双线性内插获取高程值。

步骤5.2.3)针对各个双像立体提取DEM上获取的高程值h_i，利用公式(21)获取其对应的权重值w_i，采用下式进行加权求和，获得最终高程值H：

其中n为覆盖地面坐标(D_C,D_R)的双像立体提取DEM的数量。

注：当h_i是空值或无效值时，不参与上式的加权融合计算。

步骤5.2.4)将H作为融合后DEM上像元(r′,c′)的高程值。

步骤5.3)按照从上到下、从左到右的顺序，按照步骤5.2方法逐像元获取融合后DEM各个像元的高程值，即可生成了最终的融合后DEM。

步骤6，针对融合后DEM产品，进行必要的编辑，形成最终的DEM产品。

DEM编辑内容主要包括修正影像自动匹配的错误区域(主要是水体区域、影像辐射较差区域)、优化匹配效果不理想区域，同时还需剔除地表人工建筑和人工地物高度。编辑细则取决于最终的DEM产品质量和细节要求。

DEM编辑方式主要包括自动方式编辑和人工交互编辑，自动编辑包括DEM自动内插、DEM滤波处理、DEM自动平滑等；人工交互编辑一般需要借助匹配生产DEM时的立体影像，恢复构建立体观测及量测环境并叠加DEM数据，开展人工交互编辑。

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

Claims (6)

1.一种利用三线阵立体卫星影像生产数字高程模型DEM的方法，其特征在于，所述方法包括以下步骤：

步骤A将生产任务区域内多度重叠覆盖的遥感卫星三线阵立体影像构建平差区域网并完成区域网平差，实现区域内所有卫星影像的高精度相对和绝对定向；

步骤B计算所有拥有重叠关系的任意两张影像之间的重叠率和同名光线交会角，当影像重叠率和同名光线交会角同时符合设定的门槛阈值时，将其构建为一个双像立体影像，则多度重叠覆盖的三线阵立体影像可拆分构建为多个多度重叠覆盖的双像立体影像；

步骤C针对各个双像立体影像，首先采用投影轨迹法生成近似核线立体影像，然后采用半全局密集匹配算法生成视差图，最后利用近似核线立体影像成像几何模型进行交会获取每个像素的物方坐标，生成任务区域内的物方三维点云，并通过点云栅格化生成规则格网DEM产品；

步骤D设置各个双像立体影像提取的DEM在最终加权融合时的权重值；

步骤E利用获取的权重值，对所有双像立体影像提取的DEM进行加权融合，生成整个任务区域的融合后DEM产品；

步骤F对融合后DEM产品，进行编辑，修正影像自动匹配的错误区域、剔除人工建筑和人工地物高程，形成最终的DEM产品；

所述步骤B中：将多度重叠覆盖的三线阵立体影像拆分构建为多个多度重叠覆盖的双像立体影像，具体包括以下步骤：

步骤B1利用卫星影像有理函数模型RFM，分别计算任务区域内每张影像四个角点的经纬度坐标，进而获取每张影像的地理范围多边形，利用多边形面积计算公式，求解其面积S_Pi；再利用多边形求交算法依次计算区域内所有存在相交关系的任意两张影像之间的相交多边形，同样利用多边形面积计算公式求解其面积S_U；最后利用下述公式，计算区域内所有存在相交关系的任意两张影像之间的重叠率；

其中S_P1，S_P2分别为两个多边形的面积，S_U为这前述两个多边形相交多边形的面积；

步骤B2对任务区域内所有存在重叠关系的任意两张影像，在其上任意选取1个同名像点，假定高程值分别为0和100，利用卫星影像RFM公式，分别计算影像1上该同名像点对应的地面坐标，转换到地心直角坐标系后可得影像1上该同名像点对应地面目标的入射向量为α₁，同理可得影像2上该同名像点对应地面目标的入射向量为α₂；采用下述公式计算该同名像点对应地面目标的成像交会角β，即两张影像之间的同名光线空间交会角度；

步骤B3设置一个合理的影像重叠率指标值M⁰和影像同名光线空间交会角度指标β⁰作为判断条件，依次对区域内所有存在重叠关系的任意两张影像进行资格判断，当满足M≥M⁰&β≥β⁰时，则将该两张影像构建为一个双像立体影像。

2.如权利要求1所述的利用三线阵立体卫星影像生产数字高程模型DEM的方法，其特征在于，所述步骤A中：首先量测平差区域网内密集连接点并根据需要量测适量像控点；其次选取合适的区域网平差模型，在全区域内统一进行平差迭代计算，直至满足平差精度要求，获取所有影像成像几何模型误差补偿参数；最后输出新的影像成像几何模型，完成区域内卫星影像的高精度相对和绝对定向。

3.如权利要求1所述的利用三线阵立体卫星影像生产数字高程模型DEM的方法，其特征在于，所述步骤B3中：同轨获取的前、后视影像，前、正视影像，正、后视影像，以及不同轨获取的前、后视影像，前、正视影像，正、后视影像都将能构建双像立体影像。

4.如权利要求1所述的利用三线阵立体卫星影像生产数字高程模型DEM的方法，其特征在于，所述步骤C中：规则格网DEM产品的生成具体包括如下步骤：

步骤C1将构成双像立体的两张影像分别采用投影轨迹法生成近似核线影像；

步骤C2针对立体核线影像，采用半全局密集匹配算法，开展密集匹配，获取视差图；

步骤C3根据视差图，获取核线影像同名点，通过立体核线影像的空间前方交会得到密集的地面三维点云数据，剔除粗差后，根据DEM要求的分辨率规格，将高程点云数据进行规则内插，生成DEM。

5.如权利要求1所述的利用三线阵立体卫星影像生产数字高程模型DEM的方法，其特征在于，所述步骤D中：获取各个双像立体影像提取的DEM在最终加权融合时的权重值w的计算公式为：

其中，β为双像立体影像同名光线交会角度，β₀为交会角度基准，σ为双像立体影像的绝对高程中误差，σ₀为立体影像绝对高程中误差基准，可根据实际情况设置；λ_time为双像立体影像是否为同轨成像的权重分，t_left为第一张影像成像时间，t_right为第二张影像成像时间，均精确到小时。

6.如权利要求1所述的利用三线阵立体卫星影像生产数字高程模型DEM的方法，其特征在于，所述步骤E具体包括以下步骤：

步骤E1计算所有双像立体影像提取DEM的地理范围，进而获取融合后DEM的地理范围；根据融合后DEM的分辨率需求，计算出融合后DEM像元阵列的行列数量，以及各个像元的地理坐标；

步骤E2对融合后DEM上的每一个像元(r′,c′)，计算对应地面坐标(D_C,D_R)，搜寻覆盖(D_C,D_R)的所有双像立体影像提取的DEM，然后计算(D_C,D_R)在各个双像立体影像提取的DEM上的像素位置，通过双线性内插获取该像素位置的高程值h_i，最后通过下式将这些高程值进行加权相加，作为像元(r′,c′)的高程值H：

其中w_i为h_i在加权融合时的权重值，n为覆盖地面坐标(D_C,D_R)位置的双像立体提取DEM的数量；

步骤E3按照从上到下、从左到右的顺序，按照步骤E2方法逐像元获取融合后DEM各个像元的高程值，即生成了最终的融合后DEM。

Images