CN113674415B - 利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法 - Google Patents

利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法 Download PDF

Info

Publication number
CN113674415B
CN113674415B CN202110986336.XA CN202110986336A CN113674415B CN 113674415 B CN113674415 B CN 113674415B CN 202110986336 A CN202110986336 A CN 202110986336A CN 113674415 B CN113674415 B CN 113674415B
Authority
CN
China
Prior art keywords
image
parallax
point
images
quasi
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110986336.XA
Other languages
English (en)
Other versions
CN113674415A (zh
Inventor
岳庆兴
唐新明
王霞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Original Assignee
Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ministry Of Natural Resources Land Satellite Remote Sensing Application Center filed Critical Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Priority to CN202110986336.XA priority Critical patent/CN113674415B/zh
Publication of CN113674415A publication Critical patent/CN113674415A/zh
Application granted granted Critical
Publication of CN113674415B publication Critical patent/CN113674415B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/97—Determining parameters from multiple pictures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Remote Sensing (AREA)
  • Computer Graphics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法，该方法包括：基于高分七号和资源三号立体影像及其RPC参数生成多视准核线影像；生成高分七号立体影像的多视准核线影像的金字塔影像，对金字塔影像进行匹配，以得到高分七号立体影像的视差图；生成资源三号立体影像的多视准核线影像的金字塔影像，对金字塔影像进行匹配，以得到资源三号立体影像的视差图；利用资源三号立体影像的视差图对高分七号立体影像的视差图的漏洞进行修补，并利用修补后的高分七号立体影像的视差图生成连续、无空洞DSM。本发明实现了不同准核线像对视差图的快速变换，保障了视差图漏洞修补的精度和速度，实现了高质量的连续、无空洞DSM制作。

Description

利用高分七号和资源三号影像联合制作连续、无空洞DSM的 方法
技术领域
本发明涉及图像处理技术领域，特别涉及一种利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法。
背景技术
摄影测量和计算机视觉领域对多视立体影像匹配方面的研究主要针对三视立体影像，重点往往集中在如何同时利用三个视角的匹配代价以及如何综合利用多种匹配代价方面。物方几何约束法是连接三视立体影像的核心方法，这种方法难以兼顾多视几何约束和准核线影像的“同名像点位于同一影像行”的几何约束,卫星影像多于三视条件下的高效几何约束模型构建问题还未解决。
另外，目前常见的多视角匹配都是基于分辨率接近、波段范围基本一致的“同源”立体影像。通过独立的两两构建立体像对的方法匹配各自的点云或DSM，然后把点云或DSM融合到一起。匹配核心算法方面，全局/半全局优化匹配是应用最多的匹配方法，衍生了很多改进的匹配技术。但这些算法往往难以兼顾高精度和高效率，在卫星影像匹配方面需要进一步的自适应改进。
高分七号和资源三号立体影像都属于沿轨方向的立体影像，相对几何变形主要在卫星飞行方向，目前，尚没有利用高分七号和资源三号立体影像联合制作连续、无空洞DSM的方法。
发明内容
针对上述问题，本发明提供了一种利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法，可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。
本发明提供了一种利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法，其特征在于，所述方法具体包括以下步骤：
步骤1，基于高分七号两线阵立体影像及其RPC参数和资源三号三线阵立体影像及其RPC参数生成多视准核线影像；
步骤2，生成所述高分七号两线阵立体影像的多视准核线影像的金字塔影像，对所述金字塔影像进行匹配，以得到所述高分七号两线阵立体影像的视差图；
步骤3，生成所述资源三号三线阵立体影像的多视准核线影像的金字塔影像，对所述金字塔影像进行匹配，以得到所述资源三号三线阵立体影像的视差图；
步骤4，利用所述资源三号三线阵立体影像的视差图对所述高分七号两线阵立体影像的视差图的漏洞进行修补，并利用修补后的所述高分七号两线阵立体影像的视差图生成连续、无空洞DSM。
优选的，步骤1具体包括以下子步骤：
步骤1.1，将高分七号后视影像作为主影像，高分七号前视影像作为第一辅影像，将资源三号正视影像、后视影像、前视影像分别作为第二辅影像、第三辅影像、第四辅影像，分别生成所述主影像与所述第一辅影像、第二辅影像、第三辅影像、第四辅影像间的投影连接点；
步骤1.2，建立准核线影像与对应的原始影像间的坐标转换模型，其中，所述原始影像指所述主影像、所述第一辅影像、所述第二辅影像、所述第三辅影像、所述第四辅影像，每幅原始影像与所述多视准核线影像中的一幅准核线影像相对应。
步骤1.3，基于所述多视准核线影像与对应的原始影像间的坐标转换模型，将所述准核线影像的像面坐标转换到所述原始影像的像面坐标，并基于所述原始影像的像面坐标、原始影像的RPC参数和物方高程，计算所述多视准核线影像的RPC参数；
步骤1.4，根据所述多视准核线影像与对应的原始影像间的坐标转换模型，逐点计算多视准核线影像像点对应的原始影像像点坐标，通过双线性插值获取原始影像像点的灰度作为多视准核线影像的灰度，以生成所述多视准核线影像。
优选的，步骤1.1具体包括以下子步骤：
步骤1.1.1，将所述第一辅影像作为目标辅影像，确定所述主影像与所述第一辅影像的导入顺序，并计算平均核线倾角；
步骤1.1.2，根据所述平均核线倾角确定多视准核线影像的范围；
步骤1.1.3，生成覆盖所述主影像准核线影像范围的所述主影像和所述第一辅影像间的投影连接点；
步骤1.1.4，依次将所述第二辅影像、所述第三辅影像和所述第四辅影像作为目标辅影像，计算得到所述主影像与所述第二辅影像、所述第三辅影像和所述第四辅影像的投影连接点。
优选的，步骤1.2具体包括以下子步骤：
步骤1.2.1，分别计算所述第一辅影像、第二辅影像、第三辅影像、第四辅影像沿所述准核线影像方向和垂直于所述准核线影像方向的偏移比例系数。
步骤1.2.2，计算投影轨迹倾斜参数，其中，所述投影轨迹倾斜参数表示投影轨迹线与原始影像行方向的夹角。
步骤1.2.3，根据所述偏移比例系数和所述投影轨迹倾斜参数，建立准核线影像与对应的原始影像间的坐标转换模型，具体包括：
首先，计算投影连接点行坐标与每条投影轨迹起算点坐标、影像行坐标和影像列坐标为变量的函数：
y(i,j)＝y(i,j0)+a0+a1·x(i,j)+a2·x(i,j)·x(i,j)+a3·y(i,j0)+a4·y(i,j0)·y(i,j0)+a5·x(i,j)·y(i,j0)
式中，i＝1，2，3，4，5为主影像和第一至第四辅影像的序号，y(i,j)为影像i的第j个投影连接点行坐标，x(i,j)为影像i的第j个投影连接点列坐标，j0表示第j个投影连接点所在投影轨迹的起算点点号，y(i,j0)为影像i的第j个投影连接点所在投影轨迹的起始点行坐标，a0,a1…a5为每幅原始影像待求的变换系数，由每幅影像上的投影连接点通过最小二乘法求得。
将每幅原始影像第一个投影连接点列、行坐标用c0i、r0i表示，i＝1，2，3，4，5为原始影像序号。设准核线影像任一点的行、列坐标为r、c，则对应原始影像的行、列坐标ro、co计算方法为：
ro＝r·sclci·slp+c0i-c·sclri
co＝oc+a0+a1·ro+a2·ro·ro+a3·oc+a4·oc·oc+a5·ro·oc
其中，oc＝c0i+r·sclci，sclci和sclri(i＝1,2,3,4,5)分别表示第i个原始影像沿投影轨迹方向和垂直投影轨迹方向的平均比例系数，slp为投影轨迹倾斜参数，r、c分别为准核线影像的任一点的行、列坐标，ro、co分别为与所述准核线影像的任一点在对应的原始影像上的行、列坐标，c0i为第i个原始影像的第1个投影连接点列坐标，a0,a1…a5为待求的变换系数，a0,a1…a5为待求的变换系数。
优选的，在步骤1.1.3中，根据下述公式计算所述主影像和所述第一辅影像间的投影连接点：
P0＝n·epioff+ΔW
P1＝H-n·epioff·tan(θ)+ΔH
其中，ΔH为主影像准核线影像相对主影像在行方向的外扩范围，ΔW为主影像准核线影像相对主影像在列方向的外扩范围，θ为平均核线倾角，P0为第n条投影轨迹在主影像上的起算点列坐标，P1为第n条投影轨迹在主影像上的起算点行坐标,H为主影像的长度，n为投影轨迹序号。
优选的，步骤2具体包括以下子步骤：
步骤2.1，对高分七号两线阵立体影像的准核线影像的顶层金字塔影像进行匹配，以得到顶层金字塔影像的视差图；
步骤2.2，对高分七号两线阵立体影像的多视准核线影像的顶层和底层之外的其他层的金字塔影像进行匹配，以得到所述其他层的金字塔影像的视差图；
步骤2.3，对高分七号两线阵立体影像的准核线影像的底层金字塔影像进行匹配，以得到顶层金字塔影像的视差图。
优选的，金字塔的层数优选的为5层，其中，在金字塔第3层，执行大面积可疑区域检测步骤，具体包括：
以每个无效视差为中心在一个搜索范围ssd内向8个邻域方向搜索有效视差，如果某个方向没有搜到有效视差，则将该方向的搜索距离设为最大搜索范围，8个邻域方向依次为上、下、左、右、左上、右下、左下、右上，在所述8个领域方向的搜素距离依次为d1、d2、d3、d4、d5、d6、d7、d8，根据以下公式分别求取d1和d2，d3和d4，d5和d6，d7和d8共四个对角距离：
d12＝d1+d2
d34＝d3+d4
d56＝d5+d6
d78＝d7+d8
将求取的四个对角距离d12、d34、d56、d78从小到大进行排序，如果最小的对角距离大于搜索范围ssd，则认为该无效视差属于大面积可疑区域。
优选的，在金字塔第4层，执行遮挡检测步骤，具体包括：
如果高分七号准核线影像像点的视差为无效视差时，则向该像点左右两侧搜索有效视差，设右侧有效视差为dR，右侧搜索距离为ssdR，左侧有效视差为dL，左侧搜索距离为ssdL。若同时满足：
dR-dL>3
1.2>(ssdL+ssdR)/(dR-dL)>0.8
则该无效视差对应的像点为位于遮挡区域的点。
优选的，步骤4具体包括以下子步骤：
步骤4.1，对高分七号主影像的视差图和高分七号辅影像的视差图进行修补，其中，高分七号两线阵立体影像的视差图包括高分七号主影像的视差图和高分七号辅影像的视差图；
步骤4.2，利用修补后的高分七号主影像的视差图和高分七号辅影像的视差图计算点云；
步骤4.3，通过点云网格化生成所述连续、无空洞DSM。
优选的，步骤4.3具体包括以下子步骤：
步骤4.3.1，计算所有点云的平面坐标范围，所述平面坐标范围包括最小经度Lmin、最大经度Lmax、最小纬度Bmin、最大纬度Bmax，并基于所述平面坐标范围，计算要生成的连续、无空洞DSM的长和宽；
步骤4.3.2，基于所述要生成的连续、无空洞DSM的长和宽和每个点云的三维坐标，计算所述每个点云对应的DSM网格点的行列坐标，根据每个DSM网格点的隶属点云计算DSM网格点的经纬度坐标和高程，由此生成所述连续、无空洞DSM。
本发明由于采用了上述的技术方案，使之与现有技术相比，具有以下的优点和积极效果:通过生成多视准核线影像，使得高分七号和资源三号的同名像点位于同一行影像上，并且通过多视准核线影像的坐标快速映射方法，实现不同准核线像对视差图的快速变换，保障了视差图漏洞修补的精度和速度，实现了高质量的连续、无空洞DSM制作。
附图说明
图1为根据本发明实施例的利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法的流程图；
图2为根据本发明实施例的根据双投影面法获取投影连接点的过程示意图；
图3为根据本发明实施例的用于逐像素CENSUS变换的像点；
图4为根据本发明实施例的用于间隔1像素CENSUS变换的像点。
具体实施方式
下面对本发明进行更全面的描述，其中说明本发明的示例性实施例。
如图1所示，本发明实施例提供的利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法，具体包括以下步骤：
步骤1，基于高分七号两线阵立体影像及其RPC参数和资源三号三线阵立体影像及其RPC参数生成多视准核线影像。
核线影像是将框幅式立体像对基于核线几何关系重新排列后的影像对。核线影像一般是针对两幅影像而言的，而同一轨道获取的二视以上卫星立体影像相对几何变形主要位于影像行方向，列方向具备天然的近似平行特点，可以通过多项式拟合方法将影像重新排列，保证同名像点位于同一影像行上，其作用同框幅式核线影像一样，称为准核线影像。基于本发明的方法，在多于二视的条件下，可以保证同名像点位于多幅准核线影像的同一行上，称多于二视的准核线影像为多视准核线影像。
准核线影像可以保证同名像点位于同一影像行上，因此，对于主影像的一个像点，只需在辅影像的同一行进行一维搜索即可找到同名像点，因而可以提高影像匹配的效率、精度和成功率。
本发明的所述高分七号两线阵立体影像包括后视影像和前视影像，所述资源三号三线阵立体影像包括正视影像、后视影像和前视影像，即，原始影像共有5幅影像，每幅影像相应的生成1幅准核线影像，所以共生成5幅准核线影像。一幅原始影像对应一幅准核线影像，总的准核线影像称为多视准核线影像。本发明中，有了两幅原始影像及其RPC参数，才能确定准核线影像的参数，即，本发明将2视扩展到多视，这是本发明的主要改进点之一。
步骤1具体包括以下子步骤：
步骤1.1，将高分七号后视影像作为主影像，高分七号前视影像作为第一辅影像，将资源三号正视影像、后视影像、前视影像分别作为第二辅影像、第三辅影像、第四辅影像，分别生成所述主影像与所述第一辅影像、第二辅影像、第三辅影像、第四辅影像间的投影连接点。
根据本发明的优选实施例，通过两个投影基准面获取主影像与第一辅影像的投影连接点和一个高程基准面上对应的地面点。与现有算法相比，本发明将算法扩展到多视模式。如图2所示，主影像像点P1通过主影像RPC(RPCN)和高程面H1计算地面点Pg1，将Pg1通过第一辅影像RPC(RPCB)投射到第一辅影像上得第一辅影像像点，将该第一辅影像像点通过RPCB和高程面H0计算地面点Pg2，将Pg2通过RPCN计算主影像像点P2，通过P2和RPCN及高程面H1计算地面点Pg3，将Pg3通过RPCB计算第一辅影像像点，通过这种“追赶”法得到主影像与第一辅影像间的投影连接点，即高程面H1上的对应地面点Pgn(n＝1,3,5,……)。将获取的地面点Pgn(n＝1,3,5,……)通过第二辅影像RPC(RPCF)投射到第二辅影像上，得到第二辅影像的投影连接点。主影像与第三、第四辅影像间的投影连接点的生成方法与主影像和第一辅影像间的投影连接点生成方法是相同的，在此不再赘述。
步骤1.1具体包括以下子步骤：
步骤1.1.1，将所述第一辅影像作为目标辅影像，确定所述主影像与所述第一辅影像的导入顺序，并计算平均核线倾角。
由于实际导入高分七号和资源三号原始影像的顺序可能并不确定，因此需要确定“左右”片(左右片分别指主影像和第一辅影像)的顺序。目的在于确定高程基准面的高程变化与投影坐标变化的正相关或负相关关系，用一个参数S为1或-1表示。方法为从影像的一端(本发明为下端)开始获取一个投影连接点，若该点行坐标小于初始坐标则令S＝1，否则令S＝-1。根据参数S重新设置起算高程面获取一条投影轨迹线上的投影连接点。用多组相邻连接点计算多个倾角取平均值得到平均核线倾角θ。
步骤1.1.2，根据所述平均核线倾角确定多视准核线影像的范围。
在计算投影连接点前需要首先计算几个参数，包括高程分辨率参数dh，两个投影面间距dH，主影像的准核线影像相对主影像在行列两个方向的外扩范围ΔH、ΔW。
将主影像中心坐标根据主影像RPC的HIEGHT_OFF参数投射到第一辅影像得投影坐标p0，将主影像中心坐标根据高程HIEGHT_OFF+1投射到第一辅影像得投影坐标p1，p0、p1间的像面距离为dp，则高程分辨率为dh＝1.0/dp(m)。两个投影基准面H0与H1的高程间隔设为hoff＝H/M*dh(M为每条投影轨迹线上的点数)。
主影像准核线影像的范围为主影像准核线影像相对主影像在行列两个方向的外扩范围，其中，根据以下公式计算主影像准核线影像相对主影像在列、行方向的外扩范围：
ΔW＝H·cos(θ)·sin(θ)
ΔH＝-H·sin(θ)·sin(θ)
其中，ΔW为主影像准核线影像相对主影像在列方向的外扩范围，ΔH为主影像准核线影像相对主影像在行方向的外扩范围，θ为平均核线倾角，H为主影像的长度，其中，所述主影像指高分七号后视影像。
投影起算点在主影像准核线影像行方向的基准点(起算点)数设为N，因此相邻投影起算点在主影像准核线影像行方向的距离为2·ΔW/N。
步骤1.1.3，生成覆盖所述主影像准核线影像范围的所述主影像和所述第一辅影像间的投影连接点。
如上所述，主影像与第一辅影像间需要生成N条投影轨迹线。其中，根据下述公式计算所述主影像和所述第一辅影像间的投影连接点(也称第n条投影轨迹在主影像上的起算点)：
P0＝n·epioff+ΔW
P1＝H-n·epioff·tan(θ)+ΔH
其中，epioff为垂直投影轨迹方向的相邻两条投影轨迹的平均间隔(取决于投影轨迹线的数目，优选的，本发明将总轨迹线数目N设为10，因此epioff＝W/10)，ΔH为主影像准核线影像相对主影像在行方向的外扩范围，ΔW为主影像准核线影像相对主影像在列方向的外扩范围，θ为平均核线倾角，P0为第n条投影轨迹在主影像上的起算点列坐标，P1为第n条投影轨迹在主影像上的起算点行坐标。H为主影像的长度，n为投影轨迹序号。
每条投影轨迹上的下一个点只需调整H0和H1的使用顺序即可依次获得。主影像和第一辅影像间的投影连接点是通过上述方法逐点按顺序计算的。计算的结果是均匀布满主、辅影像的投影连接点，而第二、三、四辅影像的投影连接点只需要将主影像上所有投影连接点经一个高程面全部投到第二、三、四辅影像得到。
步骤1.1.4，依次将所述第二辅影像、所述第三辅影像和所述第四辅影像作为目标辅影像，计算得到所述主影像与所述第二辅影像、所述第三辅影像和所述第四辅影像的投影连接点。
通过步骤1.1.3得到均匀布满主、辅影像的投影连接点。将主影像上所有投影连接点经一个高程面H0全部投到第二辅影像上得到第二辅影像的投影连接点；将主影像上所有投影连接点经一个高程面H0全部投到第三辅影像上得到第三辅影像的投影连接点；将主影像上所有投影连接点经一个高程面H0全部投到第四辅影像上得到第四辅影像的投影连接点。
步骤1.2，建立准核线影像与对应的原始影像间的坐标转换模型，其中，所述原始影像指所述主影像、所述第一辅影像、所述第二辅影像、所述第三辅影像、所述第四辅影像，每幅原始影像与所述多视准核线影像中的一幅准核线影像相对应。
下面对该步骤进行详细说明。获取主影像与第一辅影像、第二辅影像、第三辅影像、第四辅影像各自的投影连接点后，如何建立准核线影像与对应原始影像间的坐标变换模型仍然是一个待解决的关键问题。本发明根据每幅影像(指主影像和第一辅影像、第二辅影像、第三辅影像、第四辅影像)各自的投影连接点计算各自准核线影像的平均倾角、范围、起算点及准核线影像像点到原始影像像点的坐标变换系数。
步骤1.2具体包括以下步骤：
步骤1.2.1，分别计算所述第一辅影像、第二辅影像、第三辅影像、第四辅影像沿所述准核线影像方向和垂直于所述准核线影像方向的偏移比例系数。
本发明用辅影像第一条投影轨迹上首末两投影点的距离dci与主影像对应投影点的距离dcm的比值表示该列比例系数，用辅影像第一条投影轨迹上首投影点和最后一条投影轨迹上首投影点的距离dri与主影像第一条投影轨迹上首投影点和最后一条投影轨迹上首投影点的距离drm的比值表示该行比例系数。也可分段计算偏移比例系数并将偏移比例系数表示为影像平面坐标的函数。用符号sclci和sclri(i＝1,2,3,4,5表示主影像和所有辅影像序号)表示沿投影轨迹方向和垂直投影轨迹方向的平均比例，其中sclc0＝1，sclr0＝1。计算公式如下：
sclci＝dci/dcm，sclri＝dri/drm。
步骤1.2.2，计算投影轨迹倾斜参数，其中，所述投影轨迹倾斜参数表示投影轨迹线与原始影像行方向的夹角。
用首、末两条投影轨迹起始点坐标计算，假设首条投影轨迹起始点坐标为(rs,cs)，最后一条投影轨迹起始点坐标为(re,ce)。计算方法为两点的行坐标差除以列坐标差，用符号slp表示：
slp＝(re-rs)/(ce-cs)。
步骤1.2.3，根据所述偏移比例系数和所述投影轨迹倾斜参数，建立准核线影像与对应的原始影像间的坐标转换模型。
首先，计算投影连接点行坐标与每条投影轨迹起算点坐标、影像行坐标和影像列坐标为变量的函数：
y(i,j)＝y(i,j0)+a0+a1·x(i,j)+a2·x(i,j)·x(i,j)+a3·y(i,j0)+a4·y(i,j0)·y(i,j0)+a5·x(i,j)·y(i,j0)
式中，i＝1，2，3，4，5为主影像和第一至第四辅影像的序号，y(i,j)为影像i的第j个投影连接点行坐标，x(i,j)为影像i的第j个投影连接点列坐标，j0表示第j个投影连接点所在投影轨迹的起算点点号，y(i,j0)为影像i的第j个投影连接点所在投影轨迹的起始点行坐标，a0,a1…a5为每幅原始影像待求的变换系数，由每幅影像上的投影连接点通过最小二乘法求得。
将每幅原始影像第一个投影连接点列、行坐标用c0i、r0i表示，i＝1，2，3，4，5为原始影像序号。设准核线影像任一点的行、列坐标为r、c，则对应原始影像的行、列坐标ro、co计算方法为：
ro＝r·sclci·slp+c0i-c·sclri
co＝oc+a0+a1·ro+a2·ro·ro+a3·oc+a4·oc·oc+a5·ro·oc
其中，oc＝c0i+r·sclci，sclci和sclri(i＝1,2,3,4,5)分别表示第i个原始影像沿投影轨迹方向和垂直投影轨迹方向的平均比例系数，slp为投影轨迹倾斜参数，r、c分别为准核线影像的任一点的行、列坐标，ro、co分别为与所述准核线影像的任一点在对应的原始影像上的行、列坐标，c0i为第i个原始影像的第1个投影连接点列坐标，a0,a1…a5为待求的变换系数。
假设有N条投影轨迹，每条投影轨迹上有M个投影连接点，则每幅影像有N*M个投影连接点，可以建立N*M个ro、co的关系式，通过最小二乘法计算a0-a5。
步骤1.3，基于所述多视准核线影像与对应的原始影像间的坐标转换模型，将所述准核线影像的像面坐标转换到所述原始影像的像面坐标，并基于所述原始影像的像面坐标、原始影像的RPC参数和物方高程，计算所述多视准核线影像的RPC参数。
计算准核线影像的长Hepi、宽Wepi
Hepi＝W·cos(θ)+H·sin(θ)
Wepi＝W·sin(θ)+H·cos(θ)
其中，θ为平均核线倾角，W为主影像的宽度，H为主影像的长度。
需要说明的是，所有准核线影像的长、宽均相同。
要生成准核线影像RPC需要首先计算准核线影像像点坐标-地面坐标映射格网。具体的，本发明通过从准核线影像像面坐标(c,r)(c和r的范围分别是0-Wepi-1、0-Hepi-1)转换到原始影像像面坐标(co,ro)，再通过原始影像RPC坐标和物方高程H计算物方大地坐标经纬度(L,B)。(c,r,L,B,H)构成计算核线影像RPC的一个控制点。像面坐标密度为21*21，物方高程面分5层，中间层高程设为原始影像RPC的HEIGHT_OFF参数。相邻两层间隔为HEIGHT_SCALE参数的二分之一。在计算控制点时并不是直接采用规则像方网格坐标和物方规则高程面，而是在规则位置的基础上加入一个随机数。像方坐标加入随机数的范围为0到像方格网间距的二分之一，物方高程加入随机数的范围则为0到两个规则高程面间距的二分之一。除了控制格网外，同时生成检查点格网用与检查RPC精度。检查格网像方坐标计算方法与控制格网一致，物方高程在中间高程面的基础上做-HEIGHT_SCALE到HEIGHT_SCALE的随机浮动。
步骤1.4，根据所述多视准核线影像与对应的原始影像间的坐标转换模型，逐点计算多视准核线影像像点对应的原始影像像点坐标，通过双线性插值获取原始影像像点的灰度作为多视准核线影像的灰度，以生成所述多视准核线影像。
所述多视准核线影像的长度、宽度分别为Hepi、Wepi，根据以下公式确定所述多视准核线影像的长Hepi、宽Wepi
Hepi＝W·cos(θ)+H·sin(θ)
Wepi＝W·sin(θ)+H·cos(θ)
其中，θ为平均核线倾角，W为主影像的宽度，H为主影像的长度。
步骤2，生成所述高分七号两线阵立体影像的多视准核线影像的金字塔影像，对所述金字塔影像进行匹配，以得到所述高分七号两线阵立体影像的视差图。
本发明通过金字塔逐层精化的方式获取高分七号视差图，根据本发明的优选实施例，金字塔的层数为5层，其中，分辨率最低的一层称为顶层，原始分辨率层称为底层。假设高分七号两线阵立体影像的多视准核线影像的底层金字塔影像的宽度、高度分别为w0_g7、h0_g7，则第i层(i＝1,…,5)金字塔影像的宽、高分别为：
wi_g7＝w0_g7/2i,hi_g7＝h0_g7/2i
也就是金字塔影像底层宽、高与第i层金字塔影像的宽、高是比例关系，比例系数为2i。
第i层金字塔影像的RPC与底层金字塔影像的RPC间的关系是：
LINE_OFF_i＝(LINE_OFF_0-2i/2.0+0.5)/2i
SAMP_OFF_i＝(SAMP_OFF_0-2i/2.0+0.5)/2i
LINE_SCALE_i＝LINE_SCALE_0/2i
SAMP_SCALE_i＝SAMP_SCALE_0/2i
其中，LINE_OFF_i、LINE_OFF_0分别为第i层、底层金字塔影像的行偏移量，SAMP_OFF_i、SAMP_OFF_0分别为第i层、底层金字塔影像的列偏移量，LINE_SCALE_i、LINE_SCALE_0分别为第i层、底层金字塔影像的行缩放量，SAMP_SCALE_i、SAMP_SCALE_0分别为第i层、底层金字塔影像的列缩放量。第i层金字塔影像的每个像素灰度由原始分辨率准核线影像每2i*2i个像素取平均值得到。
步骤2具体包括以下子步骤：
步骤2.1，对高分七号两线阵立体影像的准核线影像的顶层金字塔影像进行匹配，以得到顶层金字塔影像的视差图。
顶层匹配采用类似经典SGM匹配的方法，但在匹配代价聚合步骤采用一种SGM的改进算法。
高分七号两线阵立体影像包括后视影像和前视影像，高分七号的每层金字塔影像都将高分七号后视影像作为主影像，前视影像作为辅影像，主影像的平均视差根据RPC的HEIGHT_OFF(高程偏移量)参数确定：
利用顶层主影像RPC计算(0,0)像点和高程为HEIGHT_OFF的椭球面的交点坐标，将该坐标用顶层辅影像RPC计算辅影像坐标(r,c)。则平均视差为c。
假设顶层视差搜索范围为ssh，首先计算每个像点在-ssh/2到ssh/2范围内的匹配代价立方体，每个匹配代价通过9×7窗口的CENSUS变换及汉明距离计算。
匹配代价聚合采用一种改进的消息传递方法实现。具体来说，顶层主影像的每个像素在辅影像相对主影像同名像点的列坐标差处，同时独立获取该像素的8邻域像素中的两个相邻像素的匹配代价，并取所述两个相邻像素的匹配代价的平均值作为最终的匹配代价，两个顶层主影像邻接像素的传递策略采用SGM方法。这两个相邻像素根据传递方向决定：
以左上角像素为总起始点，从左往右偶数列从上往下传播时，接收左、上邻接像素传递的匹配代价；
以左上角像素为总起始点，从左往右奇数列从下往上传播时，接收左、下邻接像素传递的匹配代价；
以右下角像素为总起始点，从右往左偶数列从下往上传播时，接收右、下邻接像素传递的匹配代价；
以左上角像素为总起始点，从右往左奇数列从上往下传播时，接收右、上邻接像素传递的匹配代价；
以左上角像素为总起始点，从上往下偶数行从左往右传播时，接收左、上邻接像素传递的匹配代价；
以左上角像素为总起始点，从上往下奇数行从右往左传播时，接收右、上邻接像素传递的匹配代价；
以右下角像素为总起始点，从下往上偶数行从右往左传播时，接收右、下邻接像素传递的匹配代价；
以左上角像素为总起始点，从下往上奇数行从左往右传播时，接收左、下邻接像素传递的匹配代价；
起始点在左侧从上往下偶数行，向右上传播时，接收左下、上邻接像素传递的匹配代价；
起始点在上侧从左往右奇数列，向左下传播时，接收右上、左邻接像素传递的匹配代价；
起始点在右侧从下往上偶数行，向左下传播时，接收右上、下邻接像素传递的匹配代价；
起始点在下侧从右往左奇数列，向右上传播时，接收左下、右邻接像素传递的匹配代价；
起始点在左侧从下往上偶数行，向左上传播时，接收右下、左邻接像素传递的匹配代价；
起始点在下侧从左往右奇数列，向右下传播时，接收左上、下邻接像素传递的匹配代价；
起始点在右侧从上往下偶数行，向右下传播时，接收左上、右邻接像素传递的匹配代价；
起始点在上侧从左往右奇数列，向左上传播时，接收右下、上邻接像素传递的匹配代价；
完成8个方向的消息传播后，将每个像素传递增强后的八个方向的匹配代价相加，得到每个像素最终的聚合匹配代价向量。
最优匹配代价的获取包括两个步骤：首先选取每个像点ssh个聚合匹配代价中最小的一个，记录其序号d。如果该序号为0或ssh-1，则视差取该序号整数值d，否则利用序号为d-1、d、d+1的匹配代价通过抛物线拟合方法计算亚像素级视差df。完成所有像素视差获取，得到主影像视差图。
交换主、辅影像顺序，同样的方法计算辅影像视差图，通过左右视差一致性检测剔除误差点。假设匹配获取的主影像第r行、c列像元的视差为dr,c，辅影像第r行、INT(c+dr,c+0.5)列(INT表示取整)像元的视差为d1，则INT(c+dr,c+0.5)+d1-c即为左右视差一致性度量，如果该值超出一个门限，即为无效值。本发明将门限设为0.8，即左右视差一致性度量超过0.8的视差为无效视差。
完成左右视差一致性检测后，还要将“被无效视差包围的面积较小的视差”设为无效值。以每个有效视差为中心在一个搜索范围ssd内(如15个像素)向8个邻域方向搜索无效视差，如果某个方向没有搜到无效视差，就该方向的搜索距离设为最大搜索范围。假设8个邻域方向依次为上、下、左、右、左上、右下、左下、右上，搜素距离依次为d1、d2、d3、d4、d5、d6、d7、d8。四个对角距离为：
d12＝d1+d2
d34＝d3+d4
d56＝d5+d6
d78＝d7+d8
将d12、d34、d56、d78从小到大排序，假设依次为d11、d21、d31、d41，如果d41大于ssd-2，则该视差设为无效值，以下称这个过程为孤立碎片过滤。
利用视差有效的影像点计算互信息参数，采用CENSUS变换后的汉明距离和互信息加权的方式重新计算匹配代价立方体并获取主辅影像的视差图，经左右一致性检测后获得有效视差图。
步骤2.2，对高分七号两线阵立体影像的多视准核线影像的顶层和底层之外的其他层的金字塔影像进行匹配，以得到所述其他层的金字塔影像的视差图。
根据本发明的优选实施例，金字塔影像共有5层，在步骤2.2，对第2、3、4层金字塔影像进行匹配。
获取上一层金字塔影像的主辅影像视差图以后，需要将视差图转移到当前层。每层视差图的长、宽与该层影像的长、宽一致，而当前层的影像和视差图的长、宽是上一层影像、视差图长、宽的两倍。用disp_crt表示当前层视差，用disp_lst表示上一层视差。当前层视差图每个偶数行、列的点与其右邻接点、下邻接点、右下邻接点构成一个“计算单元”。也就是说，“计算单元”的行、列数与上一层视差图的长、宽相等。第i行j列“计算单元”对应当前层的四个视差图像素，行列分别为：(i*2,j*2)、(i*2+1,j*2)、(i*2,j*2+1)、(i*2+1,j*2+1)。
首先将上一层视差数值乘以2：disp_lsti,j＝disp_lsti,j*2
第i行j列“计算单元”四个像素视差的计算方法为：
disp_crti*2,j*2＝disp_lsti,j
disp_crti*2+1,j*2＝(disp_lsti,j+disp_lsti+1,j)/2
disp_crti*2,j*2+1＝(disp_lsti,j+disp_lsti,j+1)/2
disp_crti*2+1,j*2＝(disp_crti*2+1,j*2+disp_crti*2,j*2+1)/2
如果等号右侧的disp_lsti,j、disp_lsti+1,j、disp_lsti,j+1中任何一项为无效值，则对应的等号左侧也为无效值。这样求得的当前层视差称为初始视差。
假设当前层主影像(0，0)像素投射到一个高程面H上得到地面坐标(L,B,H)，将(L,B,H)投射到辅影像上得到像点坐标(r,c)。主影像(0，0)像素投射到一个高程面H+dH上得到地面坐标(L1,B1,H+dH)，将(L1,B1,H+dH)投射到辅影像上得到像点坐标(r1,c1)。如果c1>c，表示“视差越大高程越大”，反之表示“视差越大高程越小”。
当前层每个像素的视差范围根据该像素初始视差(也就是由上一层视差计算的结果)及周围像素视差决定，具体包括三种情况：
一是该像素视差为有效值disp，但在某个方向或某几个方向距离无效值较近，判断依据是在8邻域方向一定搜索距离(如5个像素)之内有没有无效值存在。找出这几个方向无效值邻接的有效视差值以及沿着这几个方向跨过无效区域找到的有效视差值，统计其中的最大最小视差值(dmax和dmin)。统计以该像素为中心的5×5像素范围内的有效视差最大最小值(dmax1和dmin1)，最后统计dmax、dmin、dmax1和dmin1中的最大最小值(dmax2和dmin2)，则该像素的视差范围为(dmin2-dadd1，dmax2+dadd2)。
二是该像素视差为有效值，且不存在距离无效值较近的情况。统计以该像素为中心的5×5像素范围内的有效视差最大最小值(dmax3和dmin3)，该像素的视差范围为(dmin3-dadd3，dmax3+dadd4)。
三是该像素视差为无效值。在8邻域方向搜索有效视差，统计其中的最大最小视差值(dmax4和dmin4)，该像素的视差范围为(dmin2-dadd5，dmax2+dadd6)。以上三种条件下计算视差范围时用到的dadd1、dadd2、dadd3、dadd4、dadd5、dadd6称为视差外扩范围。在“视差越大高程越大”条件下，增大dadd2、dadd4或dadd6表示向高处搜索的范围更大，而在“视差越大高程越小”条件下，增大dadd1、dadd3或dadd5表示向高处搜索的范围更大。本发明在金字塔影像的第1、2、3层匹配时设置的所有外扩范围均为2，在金字塔影像的第4层匹配时设置的所有外扩范围均为3。
将搜索范围的最小值直接取整数，假设为dminI，最大值加0.5后取整数，假设为dmaxI，则主影像行、列坐标为(r,c)的像元在辅影像的搜索起始坐标为(r,c+dminIr,c),搜索距离为dmaxIr,c-dminIr,c+1。
根据当前层初始视差计算互信息参数，采用CENSUS变换后的汉明距离和互信息加权的方式计算每个像素的匹配代价向量。
匹配代价聚合与步骤2.1的聚合方法类似，也是采用类似的消息传递方法实现，不同的是因为相邻两个像素的搜索范围最小值可能不同，在进行消息传播时首先要确定匹配代价向量与视差的对应关系。假设当前像元(其行列为(r,c))的视差最小值为dminIr,c，邻接像元(其行列为(r1,c1))视差最小值为为dminIr1,c1，则当前像元匹配代价向量的第i个值与邻接像元匹配代价向量(用L数组表示，用整数下标表示向量中的元素，如L0，表示最小视差对应的匹配代价)的第i-dminIr1,c1+dminIr,c个值对应视差相同。视差惩罚参数为P1、P2(P1<P2)。假设邻接像元匹配代价最小值为Lmin，i1＝i-dminIr1,c1+dminIr,c，则当前像元匹配代价向量的第i个值接收的匹配代价是Lmin+P2、Li1-1+P1、Li1、Li1+1+P1中的最小值。
完成8个方向的消息传播后，将每个像素传递增强后的八个方向的匹配代价相加，得到每个像素最终的聚合匹配代价向量。与2.1的方法一样，通过抛物线拟合法得到每个像素的视差。
交换主、辅影像顺序得到辅影像视差图，主、辅影像视差图通过左右视差一致性检测和孤立碎片过滤得到剔除粗差的主、辅影像视差图。
在金字塔第3层，执行大面积可疑区域检测步骤，这是本发明增加的步骤，是本发明的主要改进点之一。方法与步骤2.1类似，具体的，以每个无效视差为中心在一个搜索范围ssd内(如95个像素)向8个邻域方向搜索有效视差，如果某个方向没有搜到有效视差，则将该方向的搜索距离设为最大搜索范围(表示搜索到了视差图边界，不知道具体数值，故设为一个较大的系数)。根据本发明的优选实施例，8个邻域方向依次为上、下、左、右、左上、右下、左下、右上，在所述8个领域方向的搜索距离依次为d1、d2、d3、d4、d5、d6、d7、d8，根据以下公式分别求取d1和d2，d3和d4，d5和d6，d7和d8共四个对角距离：
d12＝d1+d2
d34＝d3+d4
d56＝d5+d6
d78＝d7+d8
将求取的四个对角距离d12、d34、d56、d78从小到大进行排序，如果最小的对角距离大于搜索范围ssd，则认为该无效视差属于大面积可疑区域。
具体的，将d12、d34、d56、d78从小到大排序，假设依次为d11、d21、d31、d41，如果d11大于搜素范围ssd，则该无效视差属于大面积可疑区域。
经这一步确定为大面积可疑区域的点在第4层和底层不再进行处理。
在金字塔第4层，执行遮挡检测步骤，这是本发明增加的步骤，是本发明的主要改进点之一。具体的，遮挡检测步骤如下：
若当前主、辅影像为“视差越大高程越大”情况，如果高分七号准核线影像像点的视差为无效视差时，则向该像点左右两侧搜索有效视差，设右侧有效视差为dR，右侧搜索距离为ssdR，左侧有效视差为dL，左侧搜索距离为ssdL。若同时满足：
dR-dL>3
1.2>(ssdL+ssdR)/(dR-dL)>0.8
则该无效视差对应的像点为位于遮挡区域的点，也就是说该像点两侧的视差或高程差异较大，且遮挡的宽度与视差的差异接近，则有效视差中间的无效视差点为遮挡造成的。在完成金字塔影像整幅视差图的遮挡检测的基础上进行第二次遮挡检测，若某个无效视差往该无效视差点上下方向在一定搜索范围(如5个像素)内可以搜索到遮挡点，则该无效视差点也位于遮挡区域内。
确定为遮挡区域的点(以下简称遮挡点)，在确定下一层的搜索范围时，与其它类型的点计算方法不同，详见步骤2.3。
步骤2.3，对高分七号两线阵立体影像的准核线影像的底层金字塔影像进行匹配，以得到顶层金字塔影像的视差图。
底层匹配的主要步骤与顶层之外的其他层类似，不同之处在于无效初始视差及距离无效初始视差较近的像元搜索范围计算方法不同，假设当前主、辅影像为步骤2.2所述的“视差越大高程越大”情况，首先排除大面积可疑区域，若不是大面积可疑区域，则视差分以下几种类型：
类型1，与无效区域距离较近的有效点，“距离较近”的判断方法与步骤2.2的方法一样，无效区域的点包含遮挡点，下限外扩3，上限外扩5；
类型2，与无效区域距离较近的有效点，“距离较近”的判断方法与步骤2.2的方法一样，无效区域的点不包含遮挡点，下限外扩3，上限外扩25；
类型3，有效点，且不属于“与无效区域距离较近”的情况，下限外扩3，上限外扩3；
类型4，无效点且非遮挡点，下限外扩3，上限外扩25；
类型5，无效点且为遮挡点，下限外扩3，上限外扩5；
若主、辅影像为“视差越大高程越小”的情况，外扩范围设置为：
类型1，与无效区域距离较近的有效点，“距离较近”的判断方法与步骤2.2的方法一样，无效区域的点包含遮挡点，下限外扩5，上限外扩3；
类型2，与无效区域距离较近的有效点，“距离较近”的判断方法与步骤2.2的方法一样，无效区域的点不包含遮挡点，下限外扩25，上限外扩3；
类型3，有效点，且不属于“与无效区域距离较近”的情况，下限外扩3，上限外扩3；
类型4，无效点且非遮挡点，下限外扩25，上限外扩3；
类型5，无效点且为遮挡点，下限外扩5，上限外扩3；
步骤3，生成所述资源三号三线阵立体影像的多视准核线影像的金字塔影像，对所述金字塔影像进行匹配，以得到所述资源三号三线阵立体影像的视差图。
以资源三号的正视准核线影像和前视准核线影像为主、辅影像，采用步骤2的方法获取视差图，称为正-前视差图。
与步骤2唯一不同之处在于底层CENSUS变换的计算方法，步骤2采用的是图3所示的9×7像素窗口的CENSUS变换，白色像素为中心点，像素之间没有间隔，称为逐像素CENSUS变换。步骤3采用的则是图4所示的间隔1个像素的9x7个像点进行CENSUS变换，白色像素为中心点，称为间隔1像素CENSUS变换，这样做的原因是资源三号底层分辨率被采样成与高分七号一致，窗口内像素过少则特征不显著，影响匹配质量。
获取正-前视差图和正-后视差图后，需要将两个视差图进行融合。方法为：首先计算正-后视差dn-b到正-前视差dn-f的变换系数，在小范围内，两者之间呈比例关系，即
dn-f＝a×dn-b
比例系数a主要由正视相机-前视相机夹角、正视相机-后视相机夹角的相对大小决定，在小范围内接近常数。严密的计算方法为：
根据正视准核线影像像点坐标Pn和正-后视差dn-b，计算后视准核线影像像点Pb，Pn和Pb通过前方交会得到交会地面点坐标G(L,B,H)；
将交会地面点G通过前视影像RPC计算前视准核线影像像点Pf。通过Pn和Pf的列坐标计算正-前视差dn-f；
计算比例系数a＝dn-f/dn-b。
本发明每隔一定的像素间隔(如t行、t列像素，称为“节点像素”，本发明将t设为100，称为节点间隔)计算一个比例系数，以该“节点像素”为中心，长宽为节点间隔t的矩形范围内均采用该“节点像素”的比例系数。
假设正-后视差图的任一像素(r,c)(r和c分别代表行、列坐标)，如果其视差dn-b为有效值，通过“节点像素”的比例系数计算正-前视差dn-f。相同的方法将正-后视差图每个像点的视差变换成一幅新的正-前视差图(以下称原本匹配获取的正-前视差图为“原正-前视差图”，通过比例系数计算的正-前视差图为“新正-前视差图”，以便区分)。
对比“原正-前视差图”和“新正-前视差图”的每个像素的视差(用dorg和dnew区分)，有以下几种类型：
类型1，dorg有效，dnew无效；
类型2，dorg有效，dnew有效，dorg与dnew差异较小；
类型3，dorg有效，dnew有效，dorg与dnew差异较大；
类型4，dorg无效，dnew有效。
如果为类型3，则将dorg设为无效值，如果为类型4，将dorg替换为dnew。完成这一步骤，“原正-前视差图”的遮挡区域得到了修补，但也剔除了一些原本判断为有效的视差，因此会形成一些新的孤立碎片。通过步骤2.1的孤立碎片过滤方法剔除这些孤立碎片，得到最终的“原正-前视差图”。
步骤4，利用所述资源三号三线阵立体影像的视差图对所述高分七号两线阵立体影像的视差图的漏洞进行修补，并利用修补后的所述高分七号两线阵立体影像的视差图生成连续、无空洞DSM。
步骤4具体包括以下子步骤：
步骤4.1，对高分七号主影像的视差图和高分七号辅影像的视差图进行修补，其中，高分七号两线阵立体影像的视差图包括高分七号主影像的视差图和高分七号辅影像的视差图。
经过步骤2获取了高分七号两线阵立体影像的视差图，经过步骤3获取了资源三号三线阵立体影像的“原正-前视差图”(以下简称资源三号三线阵立体影像的视差图)。高分七号两线阵立体影像的视差图经过遮挡检测等粗差剔除步骤剔除了地物遮挡等因素引起的错误匹配。资源三号三线阵立体影像的视差图不仅剔除了错误匹配，而且实现了视差图的修补。接下来，需要通过分辨率较低但连续性更好的资源三号三线阵立体影像的视差图来修复分辨率高但遮挡漏洞较多的高分七号两线阵立体影像的视差图。
采用与步骤3类似的方法分两次将资源三号视差图分别变换为高分七号的主影像视差图和辅影像视差图：
第一次是以高分七号主影像视差图对应的主影像-辅影像代替步骤3中资源三号正视影像-后视影像的角色。得到修补后新的高分七号主影像视差图，具体方法为：
将资源三号正视准核线影像每个像点对应地面点的坐标(经度L、纬度B、高程H)均表示成以资源三号正视准核线影像平面坐标(行坐标r、列坐标c)和资源三号正视-前视视差d为变量的多项式形式：
L＝a0+a1·d+a2·r+a3·c+a4·r·c+a5·r·r+a6·c·c
B＝b0+…+b6·c·c
H＝c0+…+c6·c·c
其中，ai、bi、ci(i＝0～6)为前方交会拟合系数，L为地面点经度，B为地面点纬度，H为地面点高程，r为资源三号正视准核线影像行坐标，c为正视准核线影像列坐标，d为正视准核线影像像点(r,c)的视差。
具体的，在资源三号正视准核线影像范围内均匀选择n×n个(本发明将n设为11)像点坐标和通过严密前方交会计算的地面坐标，然后通过最小二乘法计算拟合系数。
计算获取资源三号正视-前视准核线影像前方交会拟合系数后，逐个像点点计算资源三号正视-前视准核线影像每个有效视差对应的地面点，假设资源三号正视准核线影像像点行、列坐标为(r,c)，则对应的地面坐标经度、纬度、高程(L,B,H)计算方法为：
L＝a0+a1·d+a2·r+a3·c+a4·r·c+a5·r·r+a6·c·c
B＝b0+…+b6·c·c
H＝c0+…+c6·c·c
将每个地面坐标经度、纬度、高程(L,B,H)由高分七号主影像准核线影像RPC计算高分七号主影像准核线影像像点行、列坐标(rb，cb)。
同时，将每个地面坐标经度、纬度、高程(L,B,H)由高分七号辅影像准核线影像RPC计算高分七号辅影像准核线影像像点行、列坐标(rf，cf)。
计算cb的整数部分cbi和小数部分cbf，cf的整数部分cfi和cff：
cbi＝INT(cb)，cbf＝cb-cbi；
cfi＝INT(cf)，cff＝cf-cfi；
其中，INT表示取整。
将rb和rf取整得到rbi和rfi：
rbi＝INT(rb+0.5),rfi＝INT(rf+0.5)。
则高分七号主影像像点(rbi,cbi)的新视差为cff-cbf，高分七号主影像像点(rfi,cfi)的新视差为cbf-cff。
按上述方法逐点计算资源三号正视准核线影像像点对应的高分七号主影像和辅影像对应整数点视差，生成新的高分七号主影像视差图和辅影像视差图(以下称原本匹配获取的高分七号主影像视差图为“原主影像视差图”，通过资源三号主影像视差图计算的高分七号主影像视差图为“新主影像视差图”；称原本匹配获取的高分七号辅影像视差图为“原辅影像视差图”，通过资源三号主影像视差图计算的高分七号辅影像视差图为“新辅影像视差图”，以便区分)。
对比“原主影像视差图”和“新主影像视差图”的每个像素的视差(用dorg和dnew区分)，有以下几种类型：
类型1，dorg有效，dnew无效；
类型2，dorg有效，dnew有效，dorg与dnew差异较小；
类型3，dorg有效，dnew有效，dorg与dnew差异较大；
类型4，dorg无效，dnew有效。
如果为类型3，则将dorg设为无效值，如果为类型4，将dorg替换为dnew。
完成这一步骤，“原主影像视差图”的遮挡区域得到了修补，但也剔除了一些原本判断为有效的视差，因此会形成一些新的孤立碎片。通过步骤2.1的孤立碎片过滤方法剔除这些孤立碎片，得到最终的“原主影像视差图”。
对比“原辅影像视差图”和“新辅影像视差图”的每个像素的视差(用dorg和dnew区分)，同样有上述4种类型。
如果为类型3，则将dorg设为无效值，如果为类型4，将dorg替换为dnew。完成这一步骤，“原辅影像视差图”的遮挡区域得到了修补，但也剔除了一些原本判断为有效的视差，因此会形成一些新的孤立碎片。通过步骤2.1的孤立碎片过滤方法剔除这些孤立碎片，得到最终的“原辅影像视差图”。
步骤4.2，利用修补后的高分七号主影像的视差图和高分七号辅影像的视差图计算点云。
获取完高分七号主、辅影像“视差图”后，要获得影像像点对应的地面坐标仍然需要通过同名像点的前方交会获取地面坐标。为了提高计算效率，本发明采用一种局部拟合算法，将每个地面点的坐标(经度L、纬度B、高程H)均表示成以高分七号准核线影像平面坐标(行坐标r、列坐标c)和视差d为变量的多项式形式：
L＝a0+a1·d+a2·r+a3·c+a4·r·c+a5·r·r+a6·c·c
B＝b0+…+b6·c·c
H＝c0+…+c6·c·c
其中，ai、bi、ci(i＝0～6)为前方交会拟合系数，L为地面点经度，B为地面点纬度，H为地面点高程，r为主影像准核线影像行坐标，c为主影像准核线影像列坐标，d为主影像准核线影像像点(r,c)的视差。
具体的，在高分七号主影像范围内均匀选择n×n个(本发明将n设为11)像点坐标和通过严密前方交会计算的地面坐标，然后通过最小二乘法计算拟合系数。
计算获取高分七号主、辅准核线影像前方交会拟合系数后，逐个像点点计算主、辅准核线影像每个有效视差对应的地面点，假设像点行、列坐标为(r,c)，则对应的地面坐标经度、纬度、高程(L,B,H)计算方法为：
L＝a0+a1·d+a2·r+a3·c+a4·r·c+a5·r·r+a6·c·c
B＝b0+…+b6·c·c
H＝c0+…+c6·c·c
最后，所有地面点一起形成用于计算最终DSM的点云。
步骤4.3，通过点云网格化生成所述连续、无空洞DSM
步骤4.3具体包括以下子步骤：
步骤4.3.1，计算所有点云的平面坐标范围，所述平面坐标范围包括最小经度Lmin、最大经度Lmax、最小纬度Bmin、最大纬度Bmax，并基于所述平面坐标范围，计算要生成的连续、无空洞DSM的长和宽。
具体的，根据以下公式计算所述连续、无空洞DSM的长HDSM、宽WDSM：
HDSM＝(Bmax-Bmin)/res
WDSM＝(Lmax-Lmin)/res
其中，res为所述连续、无空洞DSM分辨率，单位为度。
步骤4.3.2，基于所述要生成的连续、无空洞DSM的长和宽和每个点云的三维坐标，计算所述每个点云对应的DSM网格点的行列坐标，根据每个DSM网格点的隶属点云计算DSM网格点的经纬度坐标和高程，由此生成所述连续、无空洞DSM。
在该步骤，将每个点云规划到对应的规则网格位置，具体的，假设某个点云的三维坐标为(L,B,H)，它隶属于DSM的第r行、第c列(r、c的范围分别为0-HDSM-1、0-WDSM-1)(也称(L,B,H)为(r,c)网格点的“附属点”)，即，根据以下公式计算三维坐标为(L,B,H)的点云隶属的DSM的行列坐标：
r＝INT((L-Lmin)/res+0.5)
c＝INT((B-Bmin)/res+0.5)
其中，r为该点云隶属的DSM的行坐标，c为该点云隶属的DSM的列坐标，L为该点云的经度坐标，B为该点云的纬度坐标，Lmin为所有点云经度坐标的最小值，Bmin为所有点云纬度坐标的最小值，res为所述连续、无空洞DSM分辨率。
根据以下公式计算DSM的(r,c)网格点对应的经纬度(Lr,c，Br,c)：
Lr,c＝Lmin+c×res
Br,c＝Bmin+r×res
其中，Lr,c为该点云隶属的DSM的(r,c)网格点经度，Br,c为该点云隶属的DSM的(r,c)网格点经度，Lmin为所有点云经度坐标的最小值，Bmin为所有点云纬度坐标的最小值，res为所述连续、无空洞DSM分辨率，c为该点云隶属的DSM的列坐标。
DSM的(r,c)网格点的高程通过所有点云中隶属于以DSM像面坐标(r,c)为中心，s×s个DSM网格点的所有地面坐标通过距离倒数加权法计算得到。s取大于等于1的奇数，如s等于3表示取包括(r,c)及其8邻域网格点在内的9个网格点的所有“附属点”参与(r,c)网格点的高程计算。
假设附属点个数为q，附属点地面坐标表示为(Li,Bi,Hi)(i＝0,…,q-1)，Li为附属点经度坐标，Bi为附属点纬度坐标，Hi为附属点高程。根据以下公式计算DSM的(r,c)网格点的高程Hr,c：
其中，Qi为(Li,Bi,Hi)的权，计算方法为：
a称为缓冲系数，a越大，“附属点”与网格点的距离对权值Qi的影响越不明显，最终的DSM会更加平滑。反之，距离网格点较近的附属点对网格点高程的“贡献”越大，越容易呈现细节。本发明取a等于res的二分之一。
以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (3)

1.一种利用高分七号和资源三号影像联合制作连续、无空洞DSM的方法，其特征在于，所述方法具体包括以下步骤：
步骤1，基于高分七号两线阵立体影像及其RPC参数和资源三号三线阵立体影像及其RPC参数生成多视准核线影像；
步骤2，生成所述高分七号两线阵立体影像的多视准核线影像的金字塔影像，对所述金字塔影像进行匹配，以得到所述高分七号两线阵立体影像的视差图；
步骤3，生成所述资源三号三线阵立体影像的多视准核线影像的金字塔影像，对所述金字塔影像进行匹配，以得到所述资源三号三线阵立体影像的视差图；
步骤4，利用所述资源三号三线阵立体影像的视差图对所述高分七号两线阵立体影像的视差图的漏洞进行修补，并利用修补后的所述高分七号两线阵立体影像的视差图生成连续、无空洞DSM；其中，步骤1具体包括以下子步骤：
步骤1.1，将高分七号后视影像作为主影像，高分七号前视影像作为第一辅影像，将资源三号正视影像、后视影像、前视影像分别作为第二辅影像、第三辅影像、第四辅影像，分别生成所述主影像与所述第一辅影像、第二辅影像、第三辅影像、第四辅影像间的投影连接点；
步骤1.2，建立准核线影像与对应的原始影像间的坐标转换模型，其中，所述原始影像指所述主影像、所述第一辅影像、所述第二辅影像、所述第三辅影像、所述第四辅影像，每幅原始影像与所述多视准核线影像中的一幅准核线影像相对应；
步骤1.3，基于所述多视准核线影像与对应的原始影像间的坐标转换模型，将所述准核线影像的像面坐标转换到所述原始影像的像面坐标，并基于所述原始影像的像面坐标、原始影像的RPC参数和物方高程，计算所述多视准核线影像的RPC参数；
步骤1.4，根据所述多视准核线影像与对应的原始影像间的坐标转换模型，逐点计算多视准核线影像像点对应的原始影像像点坐标，通过双线性插值获取原始影像像点的灰度作为多视准核线影像的灰度，以生成所述多视准核线影像；其中，
步骤1.1具体包括以下子步骤：
步骤1.1.1，将所述第一辅影像作为目标辅影像，确定所述主影像与所述第一辅影像的导入顺序，并计算平均核线倾角；
步骤1.1.2，根据所述平均核线倾角确定多视准核线影像的范围；
步骤1.1.3，生成覆盖所述主影像准核线影像范围的所述主影像和所述第一辅影像间的投影连接点；
步骤1.1.4，依次将所述第二辅影像、所述第三辅影像和所述第四辅影像作为目标辅影像，计算得到所述主影像与所述第二辅影像、所述第三辅影像和所述第四辅影像的投影连接点；
步骤1.2具体包括以下子步骤：
步骤1.2.1，分别计算所述第一辅影像、第二辅影像、第三辅影像、第四辅影像沿所述准核线影像方向和垂直于所述准核线影像方向的偏移比例系数；
步骤1.2.2，计算投影轨迹倾斜参数，其中，所述投影轨迹倾斜参数表示投影轨迹线与原始影像行方向的夹角；
步骤1.2.3，根据所述偏移比例系数和所述投影轨迹倾斜参数，建立准核线影像与对应的原始影像间的坐标转换模型，具体包括：
首先，计算投影连接点行坐标与每条投影轨迹起算点坐标、影像行坐标和影像列坐标为变量的函数：
y(i,j)＝y(i,j0)+a0+a1·x(i,j)+a2·x(i,j)·x(i,j)+a3·y(i,j0)+a4·y(i,j0)·y(i,j0)+a5·x(i,j)·y(i,j0)
式中，i＝1，2，3，4，5为主影像和第一至第四辅影像的序号，y(i,j)为影像i的第j个投影连接点行坐标，x(i,j)为影像i的第j个投影连接点列坐标，j0表示第j个投影连接点所在投影轨迹的起算点点号，y(i,j0)为影像i的第j个投影连接点所在投影轨迹的起始点行坐标，a0,a1…a5为每幅原始影像待求的变换系数，由每幅影像上的投影连接点通过最小二乘法求得；
将每幅原始影像第一个投影连接点列、行坐标用c0i、r0i表示，i＝1，2，3，4，5为原始影像序号；设准核线影像任一点的行、列坐标为r、c，则对应原始影像的行、列坐标ro、co计算方法为：
ro＝r·sclci·slp+c0i-c·sclri
co＝oc+a0+a1·ro+a2·ro·ro+a3·oc+a4·oc·oc+a5·ro·oc
其中，oc＝c0i+r·sclci，sclci和sclri(i＝1,2,3,4,5)分别表示第i个原始影像沿投影轨迹方向和垂直投影轨迹方向的平均比例系数，slp为投影轨迹倾斜参数，r、c分别为准核线影像的任一点的行、列坐标，ro、co分别为与所述准核线影像的任一点在对应的原始影像上的行、列坐标，c0i为第i个原始影像的第1个投影连接点列坐标，a0,a1…a5为待求的变换系数；
在步骤1.1.3中，根据下述公式计算所述主影像和所述第一辅影像间的投影连接点：
P0＝n·epioff+△W
P1＝H-n·epioff·tan(θ)+△H
其中，epioff为垂直投影轨迹方向的相邻两条投影轨迹的平均间隔，△H为主影像准核线影像相对主影像在行方向的外扩范围，△W为主影像准核线影像相对主影像在列方向的外扩范围，θ为平均核线倾角，P0为第n条投影轨迹在主影像上的起算点列坐标，P1为第n条投影轨迹在主影像上的起算点行坐标,H为主影像的长度，n为投影轨迹序号；
步骤2具体包括以下子步骤：
步骤2.1，对高分七号两线阵立体影像的准核线影像的顶层金字塔影像进行匹配，以得到顶层金字塔影像的视差图；
步骤2.2，对高分七号两线阵立体影像的多视准核线影像的顶层和底层之外的其他层的金字塔影像进行匹配，以得到所述其他层的金字塔影像的视差图；
步骤2.3，对高分七号两线阵立体影像的准核线影像的底层金字塔影像进行匹配，以得到顶层金字塔影像的视差图；
步骤3具体包括：
以资源三号的正视准核线影像和前视准核线影像为主、辅影像，采用步骤2的方法获取视差图，称为正-前视差图；获取正-前视差图和正-后视差图后，再将两个视差图进行融合；方法为：
首先计算正-后视差dn-b到正-前视差dn-f的变换系数，严密的计算方法为：
根据正视准核线影像像点坐标Pn和正-后视差dn-b，计算后视准核线影像像点Pb，Pn和Pb通过前方交会得到交会地面点坐标G(L,B,H)；
将交会地面点G通过前视影像RPC计算前视准核线影像像点Pf；通过Pn和Pf的列坐标计算正-前视差dn-f；
计算比例系数a＝dn-f/dn-b；
假设正-后视差图的任一像素(r,c)，其中，r和c分别代表行、列坐标，如果其视差dn-b为有效值，通过“节点像素”的比例系数计算正-前视差dn-f；相同的方法将正-后视差图每个像点的视差变换成一幅新的正-前视差图，其中，将原本匹配获取的正-前视差图称为“原正-前视差图”，通过比例系数计算的正-前视差图称为“新正-前视差图”；
对比“原正-前视差图”和“新正-前视差图”的每个像素的视差，其中用dorg和dnew区分，有以下几种类型：
类型1，dorg有效，dnew无效；
类型2，dorg有效，dnew有效，dorg与dnew差异较小；
类型3，dorg有效，dnew有效，dorg与dnew差异较大；
类型4，dorg无效，dnew有效；
如果为类型3，则将dorg设为无效值，如果为类型4，将dorg替换为dnew；
步骤4具体包括以下子步骤：
步骤4.1，对高分七号主影像的视差图和高分七号辅影像的视差图进行修补，其中，高分七号两线阵立体影像的视差图包括高分七号主影像的视差图和高分七号辅影像的视差图；
步骤4.2，利用修补后的高分七号主影像的视差图和高分七号辅影像的视差图计算点云；
步骤4.3，通过点云网格化生成所述连续、无空洞DSM；
步骤4.3具体包括以下子步骤：
步骤4.3.1，计算所有点云的平面坐标范围，所述平面坐标范围包括最小经度Lmin、最大经度Lmax、最小纬度Bmin、最大纬度Bmax，并基于所述平面坐标范围，计算要生成的连续、无空洞DSM的长和宽；
步骤4.3.2，基于所述要生成的连续、无空洞DSM的长和宽和每个点云的三维坐标，计算所述每个点云对应的DSM网格点的行列坐标，根据每个DSM网格点的隶属点云计算DSM网格点的经纬度坐标和高程，由此生成所述连续、无空洞DSM。
2.根据权利要求1所述的方法，其特征在于，金字塔的层数为5层，其中，在金字塔第3层，执行大面积可疑区域检测步骤，具体包括：
以每个无效视差为中心在一个搜索范围ssd内向8个邻域方向搜索有效视差，如果某个方向没有搜到有效视差，则将该方向的搜索距离设为最大搜索范围，8个邻域方向依次为上、下、左、右、左上、右下、左下、右上，在所述8个领域方向的搜索距离依次为d1、d2、d3、d4、d5、d6、d7、d8，根据以下公式分别求取d1和d2，d3和d4，d5和d6，d7和d8共四个对角距离：
d12＝d1+d2
d34＝d3+d4
d56＝d5+d6
d78＝d7+d8
将求取的四个对角距离d12、d34、d56、d78从小到大进行排序，如果最小的对角距离大于搜索范围ssd，则认为该无效视差属于大面积可疑区域。
3.根据权利要求1所述的方法，其特征在于，在金字塔第4层，执行遮挡检测步骤，具体包括：
如果高分七号准核线影像像点的视差为无效视差时，则向该像点左右两侧搜索有效视差，设右侧有效视差为dR，右侧搜索距离为ssdR，左侧有效视差为dL，左侧搜索距离为ssdL；若同时满足：
dR-dL>3
1.2>(ssdL+ssdR)/(dR-dL)>0.8
则该无效视差对应的像点为位于遮挡区域的点。
CN202110986336.XA 2021-08-26 2021-08-26 利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法 Active CN113674415B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110986336.XA CN113674415B (zh) 2021-08-26 2021-08-26 利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110986336.XA CN113674415B (zh) 2021-08-26 2021-08-26 利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法

Publications (2)

Publication Number Publication Date
CN113674415A CN113674415A (zh) 2021-11-19
CN113674415B true CN113674415B (zh) 2024-01-09

Family

ID=78546492

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110986336.XA Active CN113674415B (zh) 2021-08-26 2021-08-26 利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法

Country Status (1)

Country Link
CN (1) CN113674415B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117036627A (zh) * 2023-05-25 2023-11-10 北京道达天际科技股份有限公司 一种异源立体像对影像生成数据表面模型的方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105160702A (zh) * 2015-08-20 2015-12-16 武汉大学 基于LiDAR点云辅助的立体影像密集匹配方法及系统
CN108415871A (zh) * 2017-02-10 2018-08-17 北京吉威时代软件股份有限公司 基于物方半全局多视影像匹配的密集dsm生成方法
CN112434709A (zh) * 2020-11-20 2021-03-02 西安视野慧图智能科技有限公司 基于无人机实时稠密三维点云和dsm的航测方法及系统
CN112991525A (zh) * 2021-05-07 2021-06-18 北京道达天际科技有限公司 像方和物方混合匹配基元的数字表面模型生成方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105160702A (zh) * 2015-08-20 2015-12-16 武汉大学 基于LiDAR点云辅助的立体影像密集匹配方法及系统
CN108415871A (zh) * 2017-02-10 2018-08-17 北京吉威时代软件股份有限公司 基于物方半全局多视影像匹配的密集dsm生成方法
CN112434709A (zh) * 2020-11-20 2021-03-02 西安视野慧图智能科技有限公司 基于无人机实时稠密三维点云和dsm的航测方法及系统
CN112991525A (zh) * 2021-05-07 2021-06-18 北京道达天际科技有限公司 像方和物方混合匹配基元的数字表面模型生成方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
一种适用于低空影像的数字表面模型自动获取方法;金姣;王树根;罗超;陈奇;孙明伟;;测绘与空间地理信息(05);全文 *
基于半全局优化的资源三号卫星影像DSM提取方法;岳庆兴;高小明;唐新明;;武汉大学学报(信息科学版)(10);全文 *
高分辨率遥感影像DSM的改进半全局匹配生成方法;杨幸彬;吕京国;江珊;张丹璐;;测绘学报(10);全文 *

Also Published As

Publication number Publication date
CN113674415A (zh) 2021-11-19

Similar Documents

Publication Publication Date Title
CN107977997B (zh) 一种结合激光雷达三维点云数据的相机自标定方法
CN106780590B (zh) 一种深度图的获取方法及系统
CN112434709B (zh) 基于无人机实时稠密三维点云和dsm的航测方法及系统
CN113034568B (zh) 一种机器视觉深度估计方法、装置、系统
WO2021138989A1 (zh) 多波段立体相机的深度估计加速方法
CN104820991B (zh) 一种基于代价矩阵的多重软约束立体匹配方法
CN108074218A (zh) 基于光场采集装置的图像超分辨率方法及装置
CN106023230B (zh) 一种适合变形图像的稠密匹配方法
CN110060283B (zh) 一种多测度半全局密集匹配方法
KR101681095B1 (ko) 컬러 영상과 시점 및 해상도가 동일한 깊이 영상 생성 방법 및 장치
CN108399631B (zh) 一种尺度不变性的倾斜影像多视密集匹配方法
CN113358091B (zh) 一种利用三线阵立体卫星影像生产数字高程模型dem的方法
CN108876861B (zh) 一种地外天体巡视器的立体匹配方法
US8867826B2 (en) Disparity estimation for misaligned stereo image pairs
CN103996202A (zh) 一种基于混合匹配代价和自适应窗口的立体匹配方法
CN103996201A (zh) 一种基于改进梯度和自适应窗口的立体匹配方法
CN112991420A (zh) 一种视差图的立体匹配特征提取及后处理方法
CN105005964A (zh) 基于视频序列影像的地理场景全景图快速生成方法
CN107103610B (zh) 立体测绘卫星影像匹配可疑区域自动检测方法
CN104331890B (zh) 一种全局视差估计方法和系统
CN113674415B (zh) 利用高分七号和资源三号影像联合制作连续、无空洞dsm的方法
CN111126418A (zh) 一种基于平面透视投影的倾斜影像匹配方法
CN114663298A (zh) 基于半监督深度学习的视差图修补方法及系统
CN107578429B (zh) 基于动态规划和全局代价累积路径的立体影像密集匹配方法
CN112991421A (zh) 一种机器人视觉立体匹配方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information
CB02 Change of applicant information

Address after: 1 Baisheng Village, Zizhuyuan, Haidian District, Beijing, 100089

Applicant after: Ministry of Natural Resources Land Satellite Remote Sensing Application Center

Address before: 100048 No.1 Baisheng village, Zizhuyuan, Haidian District, Beijing

Applicant before: Ministry of Natural Resources Land Satellite Remote Sensing Application Center

GR01 Patent grant
GR01 Patent grant