CN109166143A - 一种大区域网立体测绘卫星影像匹配方法 - Google Patents
一种大区域网立体测绘卫星影像匹配方法 Download PDFInfo
- Publication number
- CN109166143A CN109166143A CN201810735130.8A CN201810735130A CN109166143A CN 109166143 A CN109166143 A CN 109166143A CN 201810735130 A CN201810735130 A CN 201810735130A CN 109166143 A CN109166143 A CN 109166143A
- Authority
- CN
- China
- Prior art keywords
- image
- formula
- matched
- matching
- reference images
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 230000008569 process Effects 0.000 title claims abstract description 25
- 238000013507 mapping Methods 0.000 title claims abstract description 15
- 230000008030 elimination Effects 0.000 claims abstract description 4
- 238000003379 elimination reaction Methods 0.000 claims abstract description 4
- 238000004364 calculation method Methods 0.000 claims description 19
- 239000011159 matrix material Substances 0.000 claims description 12
- 230000009466 transformation Effects 0.000 claims description 12
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 9
- 241000668842 Lepidosaphes gloverii Species 0.000 claims description 6
- 238000003384 imaging method Methods 0.000 claims description 6
- 241000080795 Elaphoglossum affine Species 0.000 claims description 3
- 238000005314 correlation function Methods 0.000 claims description 3
- 235000013399 edible fruits Nutrition 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
- G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/02—Affine transformations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/80—Geometric correction
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种大区域网立体测绘卫星影像匹配方法，其整体步骤为：准备基准影像和待匹配影像；对基准影像进行规则格网化；利用影像自带的RPC参数和DEM数据，对待匹配影像进行遍历计算重叠区域；根据RPC参数和DEM数据将基准影像上获取的角点预测到待匹配的影像上；对基准影像与待匹配影像的分辨率进行一致性处理；对获取的基准影像和待匹配影像对应的角点集合进行相关系数匹配以及最小二乘匹配；记录最小二乘影像匹配的精确点位；对匹配点集进行RANSAC和匹配点前方交会粗差剔除；将基准影像标记为已匹配影像。本发明可以加快算法的运行速度和匹配的成功率和可靠性，对山地、高山等地形起伏大的地方也具有较高的可靠性。
Description
技术领域
本发明涉及一种方法，尤其涉及一种大区域网立体测绘卫星影像匹配方 法。
背景技术
长期以来，卫星影像的精确定位一直依赖于地面控制点，控制点数量与分 布直接影响影像对地目标定位的精度，然而，对地面特征不明显或者人员无法 到达的区域，获取足够数量的地面控制点通常非常困难，因此无控制点条件下 的精确定位已经成为一个发展趋势。大区域网立体测绘卫星匹配方法为空中三 角测量提供影像连接点，是实现高精度立体定位的前提条件。
在卫星遥感影像处理技术领域，影像匹配技术是数字摄影测量自动化的核 心技术，也是计算机视觉和图像分析等的关键技术。现有的影像匹配方法中应 用比较多的主要有核线匹配方法、金字塔匹配方法等。这些方法通常是对影像 进行全局的特征点提取，没有考虑到基准影像和待匹配影像的重叠类型，而且 在匹配的过程中采用全局遍历的方法，这样匹配不仅消耗了大量的时间，而且 获取的匹配点也不一定正确。此外，由于大区域中的多尺度遥感影像之间存在 巨大差异，因此利用现有的匹配方法进行影像匹配时，其匹配速度较慢，精度 较低，且对山地、高山等地形起伏较大的地方具有较低的可靠性。
发明内容
为了解决上述技术所存在的不足之处，本发明提供了一种大区域网立体测 绘卫星影像匹配方法。
为了解决以上技术问题，本发明采用的技术方案是：一种大区域网立体测 绘卫星影像匹配方法，其整体步骤为：
步骤一、获取参与匹配的n幅影像集合，将影像集合中的第一幅影像作为 基准影像，其它影像作为待匹配影像，影像准备就绪；
步骤二、对基准影像进行规则格网化，在格网化区域中进行Forstner角点 提取；
步骤三、利用影像自带的RPC参数和DEM数据，对待匹配影像进行遍历计算 重叠区域，如果基准影像与待匹配影像没有重叠则继续本步骤，如果有重叠则 进行下一步骤；
步骤四、根据RPC参数和DEM数据将基准影像上获取的角点预测到待匹配的 影像上；
步骤五、检查基准影像与待匹配影像的分辨率是否一致，如果不一致则进 行分辨率一致性处理；
步骤六、对获取的基准影像和待匹配影像对应的角点集合进行相关系数匹 配，如果相关系数大于阈值，则进行最小二乘匹配，否则对基准影像的下一个 角点重复该步骤；
步骤七、记录最小二乘影像匹配的精确点位，基准影像的下一个特征点重 复步骤六～七的过程，直至基准影像所有点集处理完毕；
步骤八、对匹配点集进行RANSAC和匹配点前方交会粗差剔除；
步骤九、将基准影像标记为已匹配影像，再从待匹配影像中选取基准影 像，进行步骤二～八的过程处理。
2、根据权利要求1所述的大区域网立体测绘卫星影像匹配方法，其特征在 于：所述步骤三中重叠区域的计算方法为：
根据两张影像附带的RPC参数以及辅助的全球DEM数据，参照公式1计算基 准影像和待匹配影像的四个角点(0,0),(width-1,0),(0,height-1),(width- 1,height-1)的经纬度坐标，进而求得投影到平均高程面h＝H上的影像四个角点 对应的经纬度坐标(P1i,L1j),(P2i,L2j)，将四个角点对应的经纬度坐标组成四边 形，则分别得到两个四边形S1和S2；
式中，F1、F2、F3、F4为一般多项式，计算方式如下：
式中，bijk(i,j,k＝0,1...20)为反解RPC参数；(P,L,H)为正则化的地面坐 标，(X,Y)为正则化的影像坐标，计算公式分别为：
式中，LAT_OFF、LAT_SCALE、LONE_OFF、LONG_SCALE、HEIGHT_OFF和 HEIGHT_SCALE为地面坐标的正则化参数；SMAP_OFF、SAMP_SCALE、LINE_OFF和 LINE_SCALE为影像坐标的正则化参数；
对两个四边形S1和S2进行求交运算，得到一个多边形，将多边形的外接矩 形确定为基准影像和待匹配影像的影像重叠区；再根据基准影像附带的RPC参 数以及影像重叠区的四个角点对应的经纬度坐标，按照公式4中所示的RPC反算 式从平均高程面上反算，从而获取基准影像和待匹配影像的同名区；
式中，F1、F2、F3、F4为一般多项式，计算方式如下：
其中，aijk(i,j,k＝0,1...20)为RPC系数；(P,L,H)为正则化的地面坐标， (X,Y)为正则化的影像坐标，计算公式分别为：
式中，LAT_OFF、LAT_SCALE、LONE_OFF、LONG_SCALE、HEIGHT_OFF和 HEIGHT_SCALE为地面坐标的正则化参数；SMAP_OFF、SAMP_SCALE、LINE_OFF和 LINE_SCALE为影像坐标的正则化参数；如果有交集则有重叠区域，否则没有。
3、根据权利要求2所述的大区域网立体测绘卫星影像匹配方法，其特征在 于：所述步骤六中相关系数的匹配原理为：
通过在目标窗口和搜索窗口之间统计相关系数ρ(c,r)，将相关系数最大的 位置作为同名点，计算公式如下所示：
其中，g和g′为目标窗口和搜索窗口中对应的像元灰度值，i和j为待匹 配窗口的行列坐标，r和c为对应搜索区域的行列坐标；和为目标窗口和 搜索窗口内的灰度的均值计算公式；
定义W1为基准影像中的目标窗口，大小为(2N+1)×(2N+1)(N＝1，2， 3...)，W2为待匹配影像的搜索窗口，大小为(2M+1)×(2M+1)(M＝1，2， 3...)，通常情况下M＞N，在W2窗口中按照行列顺序遍历目标窗口W1，计 算归一化互相关系数ρ(c,r)，ρ(c,r)最大时对应的待匹配影像中像素点为基准 影像中心点的同名点。
4、根据权利要求1所述的大区域网立体测绘卫星影像匹配方法，其特征在 于：所述步骤六中最小二乘的匹配原理为：
最小二乘影像匹配判断相似的标准是灰度差的平方和最小，将灰度差记为 v，则如公式10所示：
∑vv＝min 公式10
一般情况下两个二维影像之间存在几何变形，只有充分的考虑到影像的几 何变形，才能获得最佳的影像匹配结果，在实际计算中，由于影像匹配窗口所 选择的尺寸一般较小，所以用仿射变换关系来表示影像间的几何畸变，如公式11所示：
其中，a0、a1、a2和b0、b1、b2为左右影像间几何变形改正的参数，即 仿射变换参数，(x,y)为左影像窗口的像素坐标，(x2,y2)为经过仿射变换后的 右影像窗口的像素坐标；
同时考虑左右影像的线性灰度畸变，则得到公式12：
g(x,y)+n1(x,y)＝h0+h1g2(a0+a1x+a2y,b0+b1x+b2y)+n2(x,y) 公式12
其中，g(x,y)、g2(x2,y2)为灰度函数，h0、h1为线性畸变的参数，初始值 为h0＝0、h1＝1；n1(x,y)、n2(x,y)为随机噪声；
将公式12线性化后，可得最小二乘匹配误差方程式，如公式13所示：
v＝c1dh0+c2dh1+c3da0+c4da1+c5da2+c6db0+c7db1+c8db2-Δg公式13
其中，dh0、dh1、da0，...，db2为待定参数的改正值，Δg为相应像素的 灰度差，它们初值为h0＝0、h1＝1、a0＝0、a1＝1、a2＝0、b0＝0、b1＝0、b2＝1。 c1、c2、...、c8为误差方程式系数，如公式14所示，其中为影像灰 度的差分，g2为灰度值，x、y为像素坐标；
则误差方程式的矩阵形式如公式15所示：
V＝CX-L 公式15
其中，C表示误差方程式系数矩阵，X表示待定参数改正数矩阵，L表示 窗口内像素的灰度差，如公式16所示；T为矩阵的转置；
由于影像匹配中灰度为规则格网化的离散阵列，把采样间隔看作单位长 度，所以计算中的偏导数可用差分取代，如公式17所示：
最小二乘影像匹配是一个通过迭代进行匹配的算法，迭代开始时各个参数 的初始值为，h0＝a0＝a2＝b0＝b1＝0，h1＝a1＝b2＝1，设 为第i-1次变形时参数，为第i次 迭代时参数改正值，则第i次几何畸变改正参数值如公式18所示：
第i次辐射畸变改正参数值如公式19所示：
由于最小二乘匹配的精度与影像的灰度梯度有很大关系，所以左方影像可 通过梯度的平方权对坐标加权平均来精确点位，如公式20所示：
其中，和表示影像灰度的梯度；
根据最小二乘影像匹配仿射变换可得右影像同名点精确坐标，如公式21所 示。
本发明通过将RPC参数文件和数字高程模型相结合，获取基准影像和待匹 配影像的重叠区域，并在重叠区域内进行格网划分，在获取的重叠区域格网中 进行点集的提取和进行点集的匹配，这样不仅提高了匹配的速度，而且提高了 匹配的成功率和精确度。
附图说明
图1为本发明的整体流程示意图。
图2为基准影像和待匹配影像的相关统计过程展示图。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1所示的一种大区域网立体测绘卫星影像匹配方法，其整体流程为：
步骤一、获取参与匹配的n幅影像集合，将影像集合中的第一幅影像作为 基准影像，其它影像作为待匹配影像，影像准备就绪；
步骤二、对基准影像进行规则格网化，在格网化区域中进行Forstner角点 提取；
步骤三、利用影像自带的RPC参数和DEM数据，对待匹配影像进行遍历计算 重叠区域，如果基准影像与待匹配影像没有重叠则继续本步骤，如果有重叠则 进行下一步骤；
重叠区域判断方法如下：
由于卫星影像含有rpc参数文件，因此待匹配影像附带有理多项式系数(Rational Polynomial Coefficients，RPC)参数；
根据两张影像附带的RPC参数以及辅助的全球DEM(Digital Elevation Model，数字高程模型)数据，参照公式1计算基准影像和待匹配影像的四个角 点(0,0),(width-1,0),(0,height-1),(width-1,height-1)的经纬度坐标，进 而求得投影到平均高程面h＝H上的影像四个角点对应的经纬度坐标 (P1i,L1j),(P2i,L2j)，将四个角点对应的经纬度坐标组成四边形，即分别得到两 个四边形S1和S2。
式中，F1、F2、F3、F4为一般多项式，计算方式如下：
式中，bijk(i,j,k＝0,1...20)为反解RPC参数；(P,L,H)为正则化的地面坐 标，(X,Y)为正则化的影像坐标，计算公式分别为：
式中，LAT_OFF、LAT_SCALE、LONE_OFF、LONG_SCALE、HEIGHT_OFF和 HEIGHT_SCALE为地面坐标的正则化参数。SMAP_OFF、SAMP_SCALE、LINE_OFF和 LINE_SCALE为影像坐标的正则化参数。
对两个四边形S1和S2进行求交运算，得到一个多边形，将多边形的外接矩 形确定为基准影像和待匹配影像的影像重叠区。再根据基准影像附带的RPC参 数以及影像重叠区的四个角点对应的经纬度坐标，按照公式4中所示的RPC反算 式从平均高程面上反算，从而获取基准影像和待匹配影像的同名区。
式中，F1、F2、F3、F4为一般多项式，计算方式如下：
其中，aijk(i,j,k＝0,1...20)为RPC系数；(P,L,H)为正则化的地面坐标， (X,Y)为正则化的影像坐标，计算公式分别为：
式中，LAT_OFF、LAT_SCALE、LONE_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面坐标的正则化参数。SMAP_OFF、SAMP_SCALE、LINE_OFF和 LINE_SCALE为影像坐标的正则化参数。如果有交集则有重叠区域，否则没有。
步骤四、根据RPC参数和DEM数据将基准影像上获取的角点预测到待匹配的 影像上；
步骤五、检查基准影像与待匹配影像的分辨率是否一致，如果不一致则进 行分辨率一致性处理；
步骤六、对获取的基准影像和待匹配影像对应的角点集合进行相关系数匹 配，如果相关系数大于阈值，则进行最小二乘匹配，否则对基准影像的下一个 角点重复该步骤。
相关系数的匹配原理具体为：
通过在目标窗口和搜索窗口之间统计相关系数ρ(c,r)，将相关系数最大的 位置作为同名点，计算公式如下所示。
其中，g和g′为目标窗口和搜索窗口中对应的像元灰度值，i和j为待匹 配窗口的行列坐标，r和c为对应搜索区域的行列坐标。和为目标窗口和 搜索窗口内的灰度的均值计算公式。
具体实现过程如图2所示，左影像为基准影像，右影像为待匹配影像，W1 为基准影像中的目标窗口，大小为(2N+1)×(2N+1)(N＝1，2，3...)，W2为 待匹配影像的搜索窗口，大小为(2M+1)×(2M+1)(M＝1，2，3...)，通常情 况下M＞N，在W2窗口中按照行列顺序遍历目标窗口W1，计算归一化互相关 系数ρ(c,r)，ρ(c,r)最大时对应的待匹配影像中像素点为基准影像中心点的同 名点。
最小二乘的匹配原理具体为：
最小二乘影像匹配判断相似的标准是灰度差的平方和最小，将灰度差记为 v，如公式10所示。
∑vv＝min 公式10
一般情况下两个二维影像之间存在几何变形，只有充分的考虑到影像的几 何变形，才能获得最佳的影像匹配结果，在实际计算中，由于影像匹配窗口所 选择的尺寸一般较小，所以可以用仿射变换关系来表示影像间的几何畸变，如 公式11所示。
其中，a0、a1、a2和b0、b1、b2为左右影像间几何变形改正的参数，即 仿射变换参数，(x,y)为左影像窗口的像素坐标，(x2,y2)为经过仿射变换后的 右影像窗口的像素坐标。
如果同时考虑左右影像的线性灰度畸变则可得公式12。
g(x,y)+n1(x,y)＝h0+h1g2(a0+a1x+a2y,b0+b1x+b2y)+n2(x,y) 公式12
其中，g(x,y)、g2(x2,y2)为灰度函数，h0、h1为线性畸变的参数，初始值 为h0＝0、h1＝1；n1(x,y)、n2(x,y)为随机噪声；
将公式12线性化后，可得最小二乘匹配误差方程式，如公式13所示。
v＝c1dh0+c2dh1+c3da0+c4da1+c5da2+c6db0+c7db1+c8db2-Δg公式13
其中dh0、dh1、da0，...，db2为待定参数的改正值，Δg为相应像素的灰 度差，它们初值为h0＝0、h1＝1、a0＝0、a1＝1、a2＝0、b0＝0、b1＝0、b2＝1。c1、c2、...、c8为误差方程式系数，如公式14所示，其中为影像灰 度的差分，g2为灰度值，x、y为像素坐标。
则误差方程式的矩阵形式如公式15所示。
V＝CX-L 公式15
其中，C表示误差方程式系数矩阵，X表示待定参数改正数矩阵，L表示 窗口内像素的灰度差，如公式16所示。T为矩阵的转置。
由于影像匹配中灰度为规则格网化的离散阵列，可以把采样间隔看作单位 长度，所以计算中的偏导数可用差分取代，如公式17所示。
最小二乘影像匹配是一个通过迭代进行匹配的算法，迭代开始时各个参数 的初始值为，h0＝a0＝a2＝b0＝b1＝0，h1＝a1＝b2＝1，设 为第i-1次变形时参数，为第i次 迭代时参数改正值，则第i次几何畸变改正参数值如公式18所示。
第i次辐射畸变改正参数值如公式19。
由于最小二乘匹配的精度与影像的灰度梯度有很大关系，所以左方影像可 通过梯度的平方权对坐标加权平均来精确点位，如公式20所示。
其中，和表示影像灰度的梯度。
根据最小二乘影像匹配仿射变换可得右影像同名点精确坐标，如公式21所 示。
步骤七、记录最小二乘影像匹配的精确点位，基准影像的下一个特征点重 复步骤六、步骤七的过程，直至基准影像所有点集处理完毕；
步骤八、对匹配点集进行RANSAC和匹配点前方交会粗差剔除；
步骤九、将基准影像标记为已匹配影像，再从待匹配影像中选取基准影 像，进行步骤二至步骤八的过程处理。
本发明的关键创新点在于：
1.充分利用了影像自带的rpc参数文件，计算了影像之间的关系，判断影 像是否重叠，减小了冗余影像的处理过程；
2.根据rpc文件和DEM数据进行点位预测，缩小了匹配范围，加快了算法的 运行速度和匹配的成功率和可靠性。
本发明针对推扫式光学卫星遥感影像成像的特点，对输入的影像集，利用 遥感影像附带的RPC参数，并加入全球DEM数据作为辅助数据，实现了快速的同 名区位置预测，然后在同名区内提取角点进行匹配，有效减少了无效输入，加 快了匹配速度；另外，通过RPC参数和全球DEM数据进行对应角点预测，然后根 据相关系数匹配和最小二乘影像匹配，最终确定最佳的匹配点，可以解决传统 匹配方法在进行大区域影像匹配时存在的速度较慢、精度较低的问题，且对山 地、高山等地形起伏大的地方也具有较高的可靠性。本发明已经经过TH和ZY3 大区域网数据测试，匹配点数和匹配效率以及匹配精度都满足要求。
上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本 技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或 替换，也均属于本发明的保护范围。
Claims (4)
1.一种大区域网立体测绘卫星影像匹配方法，其特征在于：所述方法的整体步骤为：
步骤一、获取参与匹配的n幅影像集合，将影像集合中的第一幅影像作为基准影像，其它影像作为待匹配影像，影像准备就绪；
步骤二、对基准影像进行规则格网化，在格网化区域中进行Forstner角点提取；
步骤三、利用影像自带的RPC参数和DEM数据，对待匹配影像进行遍历计算重叠区域，如果基准影像与待匹配影像没有重叠则继续本步骤，如果有重叠则进行下一步骤；
步骤四、根据RPC参数和DEM数据将基准影像上获取的角点预测到待匹配的影像上；
步骤五、检查基准影像与待匹配影像的分辨率是否一致，如果不一致则进行分辨率一致性处理；
步骤六、对获取的基准影像和待匹配影像对应的角点集合进行相关系数匹配，如果相关系数大于阈值，则进行最小二乘匹配，否则对基准影像的下一个角点重复该步骤；
步骤七、记录最小二乘影像匹配的精确点位，基准影像的下一个特征点重复步骤六～七的过程，直至基准影像所有点集处理完毕；
步骤八、对匹配点集进行RANSAC和匹配点前方交会粗差剔除；
步骤九、将基准影像标记为已匹配影像，再从待匹配影像中选取基准影像，进行步骤二～八的过程处理。
2.根据权利要求1所述的大区域网立体测绘卫星影像匹配方法，其特征在于：所述步骤三中重叠区域的计算方法为：
根据两张影像附带的RPC参数以及辅助的全球DEM数据，参照公式1计算基准影像和待匹配影像的四个角点(0,0),(width-1,0),(0,height-1),(width- 1,height-1)的经纬度坐标，进而求得投影到平均高程面h＝H上的影像四个角点对应的经纬度坐标(P1i,L1j),(P2i,L2j)，将四个角点对应的经纬度坐标组成四边形，则分别得到两个四边形S1和S2；
式中，F1、F2、F3、F4为一般多项式，计算方式如下：
式中，bijk(i,j,k＝0,1...20)为反解RPC参数；(P,L,H)为正则化的地面坐标，(X,Y)为正则化的影像坐标，计算公式分别为：
式中，LAT_OFF、LAT_SCALE、LONE_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面坐标的正则化参数；SMAP_OFF、SAMP_SCALE、LINE_OFF和LINE_SCALE为影像坐标的正则化参数；
对两个四边形S1和S2进行求交运算，得到一个多边形，将多边形的外接矩形确定为基准影像和待匹配影像的影像重叠区；再根据基准影像附带的RPC参数以及影像重叠区的四个角点对应的经纬度坐标，按照公式4中所示的RPC反算式从平均高程面上反算，从而获取基准影像和待匹配影像的同名区；
式中，F1、F2、F3、F4为一般多项式，计算方式如下：
其中，aijk(i,j,k＝0,1...20)为RPC系数；(P,L,H)为正则化的地面坐标，(X,Y)为正则化的影像坐标，计算公式分别为：
式中，LAT_OFF、LAT_SCALE、LONE_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面坐标的正则化参数；SMAP_OFF、SAMP_SCALE、LINE_OFF和LINE_SCALE为影像坐标的正则化参数；如果有交集则有重叠区域，否则没有。
3.根据权利要求2所述的大区域网立体测绘卫星影像匹配方法，其特征在于：所述步骤六中相关系数的匹配原理为：
通过在目标窗口和搜索窗口之间统计相关系数ρ(c,r)，将相关系数最大的位置作为同名点，计算公式如下所示：
其中，g和g′为目标窗口和搜索窗口中对应的像元灰度值，i和j为待匹配窗口的行列坐标，r和c为对应搜索区域的行列坐标；和为目标窗口和搜索窗口内的灰度的均值计算公式；
定义W1为基准影像中的目标窗口，大小为(2N+1)×(2N+1)(N＝1，2，3...)，W2为待匹配影像的搜索窗口，大小为(2M+1)×(2M+1)(M＝1，2，3...)，通常情况下M＞N，在W2窗口中按照行列顺序遍历目标窗口W1，计算归一化互相关系数ρ(c,r)，ρ(c,r)最大时对应的待匹配影像中像素点为基准影像中心点的同名点。
4.根据权利要求1所述的大区域网立体测绘卫星影像匹配方法，其特征在于：所述步骤六中最小二乘的匹配原理为：
最小二乘影像匹配判断相似的标准是灰度差的平方和最小，将灰度差记为v，则如公式10所示：
∑vv＝min 公式10
一般情况下两个二维影像之间存在几何变形，只有充分的考虑到影像的几何变形，才能获得最佳的影像匹配结果，在实际计算中，由于影像匹配窗口所选择的尺寸一般较小，所以用仿射变换关系来表示影像间的几何畸变，如公式11所示：
其中，a0、a1、a2和b0、b1、b2为左右影像间几何变形改正的参数，即仿射变换参数，(x,y)为左影像窗口的像素坐标，(x2,y2)为经过仿射变换后的右影像窗口的像素坐标；
同时考虑左右影像的线性灰度畸变，则得到公式12：
g(x,y)+n1(x,y)＝h0+h1g2(a0+a1x+a2y,b0+b1x+b2y)+n2(x,y) 公式12
其中，g(x,y)、g2(x2,y2)为灰度函数，h0、h1为线性畸变的参数，初始值为h0＝0、h1＝1；n1(x,y)、n2(x,y)为随机噪声；
将公式12线性化后，可得最小二乘匹配误差方程式，如公式13所示：
v＝c1dh0+c2dh1+c3da0+c4da1+c5da2+c6db0+c7db1+c8db2-Δg公式13
其中，dh0、dh1、da0，...，db2为待定参数的改正值，Δg为相应像素的灰度差，它们初值为h0＝0、h1＝1、a0＝0、a1＝1、a2＝0、b0＝0、b1＝0、b2＝1。c1、c2、...、c8为误差方程式系数，如公式14所示，其中为影像灰度的差分，g2为灰度值，x、y为像素坐标；
则误差方程式的矩阵形式如公式15所示：
V＝CX-L 公式15
其中，C表示误差方程式系数矩阵，X表示待定参数改正数矩阵，L表示窗口内像素的灰度差，如公式16所示；T为矩阵的转置；
由于影像匹配中灰度为规则格网化的离散阵列，把采样间隔看作单位长度，所以计算中的偏导数可用差分取代，如公式17所示：
最小二乘影像匹配是一个通过迭代进行匹配的算法，迭代开始时各个参数的初始值为，h0＝a0＝a2＝b0＝b1＝0，h1＝a1＝b2＝1，设 为第i-1次变形时参数，为第i次迭代时参数改正值，则第i次几何畸变改正参数值如公式18所示：
第i次辐射畸变改正参数值如公式19所示：
由于最小二乘匹配的精度与影像的灰度梯度有很大关系，所以左方影像可通过梯度的平方权对坐标加权平均来精确点位，如公式20所示：
其中，和表示影像灰度的梯度；
根据最小二乘影像匹配仿射变换可得右影像同名点精确坐标，如公式21所示。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810735130.8A CN109166143A (zh) | 2018-07-06 | 2018-07-06 | 一种大区域网立体测绘卫星影像匹配方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810735130.8A CN109166143A (zh) | 2018-07-06 | 2018-07-06 | 一种大区域网立体测绘卫星影像匹配方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109166143A true CN109166143A (zh) | 2019-01-08 |
Family
ID=64897402
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810735130.8A Pending CN109166143A (zh) | 2018-07-06 | 2018-07-06 | 一种大区域网立体测绘卫星影像匹配方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109166143A (zh) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109711486A (zh) * | 2019-01-21 | 2019-05-03 | 湖北省国土资源研究院 | 基于相位相关的高重叠度遥感影像满度连接点匹配方法 |
| CN112560868A (zh) * | 2020-11-27 | 2021-03-26 | 西安中科星图空间数据技术有限公司 | 一种基于特征点库的多进程大区域网影像匹配方法及装置 |
| CN113313769A (zh) * | 2021-06-11 | 2021-08-27 | 湖北工业大学 | 一种光学卫星多面阵传感器片间无缝几何定标方法 |
| CN113884036A (zh) * | 2021-09-22 | 2022-01-04 | 浙江省水利河口研究院(浙江省海洋规划设计研究院) | 一种坐标转换方法及装置、电子设备、存储介质 |
| CN114596342A (zh) * | 2022-03-22 | 2022-06-07 | 感知天下(北京)信息科技有限公司 | 基于卫星影像的并行匹配方法、存储介质及计算机设备 |
| CN114998399A (zh) * | 2022-05-20 | 2022-09-02 | 中国人民解放军61540部队 | 一种异源光学遥感卫星影像立体像对预处理方法 |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101915913A (zh) * | 2010-07-30 | 2010-12-15 | 中交第二公路勘察设计研究院有限公司 | 一种稳健的高分辨率卫星影像连接点自动匹配方法 |
| CN102073874A (zh) * | 2010-12-29 | 2011-05-25 | 中国资源卫星应用中心 | 附加几何约束的航天三线阵ccd相机多影像立体匹配方法 |
| CN102693542A (zh) * | 2012-05-18 | 2012-09-26 | 中国人民解放军信息工程大学 | 一种影像特征匹配方法 |
| CN102968788A (zh) * | 2012-10-25 | 2013-03-13 | 武汉大学 | 一种基于规则格网面元的波段配准方法 |
| CN103822616A (zh) * | 2014-03-18 | 2014-05-28 | 武汉大学 | 一种图分割与地形起伏约束相结合的遥感影像匹配方法 |
| CN103886611A (zh) * | 2014-04-08 | 2014-06-25 | 西安煤航信息产业有限公司 | 一种适合于航空摄影飞行质量自动检查的影像匹配方法 |
| CN105913435A (zh) * | 2016-04-13 | 2016-08-31 | 西安航天天绘数据技术有限公司 | 一种适用于大区域的多尺度遥感影像匹配方法及系统 |
| CN106651927A (zh) * | 2016-12-30 | 2017-05-10 | 北京航天泰坦科技股份有限公司 | 一种正射影像镶嵌接边匹配的同名点粗差剔除方法 |
| US9740950B1 (en) * | 2013-03-15 | 2017-08-22 | Harris Corporation | Method and system for automatic registration of images |
| CN107705244A (zh) * | 2017-09-11 | 2018-02-16 | 中国国土资源航空物探遥感中心 | 一种适用于大区域多幅遥感影像的接边纠正方法 |
-
2018
- 2018-07-06 CN CN201810735130.8A patent/CN109166143A/zh active Pending
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101915913A (zh) * | 2010-07-30 | 2010-12-15 | 中交第二公路勘察设计研究院有限公司 | 一种稳健的高分辨率卫星影像连接点自动匹配方法 |
| CN102073874A (zh) * | 2010-12-29 | 2011-05-25 | 中国资源卫星应用中心 | 附加几何约束的航天三线阵ccd相机多影像立体匹配方法 |
| CN102693542A (zh) * | 2012-05-18 | 2012-09-26 | 中国人民解放军信息工程大学 | 一种影像特征匹配方法 |
| CN102968788A (zh) * | 2012-10-25 | 2013-03-13 | 武汉大学 | 一种基于规则格网面元的波段配准方法 |
| US9740950B1 (en) * | 2013-03-15 | 2017-08-22 | Harris Corporation | Method and system for automatic registration of images |
| CN103822616A (zh) * | 2014-03-18 | 2014-05-28 | 武汉大学 | 一种图分割与地形起伏约束相结合的遥感影像匹配方法 |
| CN103886611A (zh) * | 2014-04-08 | 2014-06-25 | 西安煤航信息产业有限公司 | 一种适合于航空摄影飞行质量自动检查的影像匹配方法 |
| CN105913435A (zh) * | 2016-04-13 | 2016-08-31 | 西安航天天绘数据技术有限公司 | 一种适用于大区域的多尺度遥感影像匹配方法及系统 |
| CN106651927A (zh) * | 2016-12-30 | 2017-05-10 | 北京航天泰坦科技股份有限公司 | 一种正射影像镶嵌接边匹配的同名点粗差剔除方法 |
| CN107705244A (zh) * | 2017-09-11 | 2018-02-16 | 中国国土资源航空物探遥感中心 | 一种适用于大区域多幅遥感影像的接边纠正方法 |
Non-Patent Citations (3)
| Title |
|---|
| 熊金鑫: "基于区域分割的多源、多时相卫星遥感影像联合匹配方法研究", 《中国博士学位论文全文数据库_基础科学辑》 * |
| 胡著智等: "《航天航空遥感技术与应用》", 31 March 2007, 南京大学出版社 * |
| 陈永明: "《航空摄影测量》", 31 July 2003, 《中国建筑工业出版社》 * |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109711486A (zh) * | 2019-01-21 | 2019-05-03 | 湖北省国土资源研究院 | 基于相位相关的高重叠度遥感影像满度连接点匹配方法 |
| CN109711486B (zh) * | 2019-01-21 | 2020-09-01 | 湖北省国土资源研究院 | 基于相位相关的高重叠度遥感影像满度连接点匹配方法 |
| CN112560868A (zh) * | 2020-11-27 | 2021-03-26 | 西安中科星图空间数据技术有限公司 | 一种基于特征点库的多进程大区域网影像匹配方法及装置 |
| CN113313769A (zh) * | 2021-06-11 | 2021-08-27 | 湖北工业大学 | 一种光学卫星多面阵传感器片间无缝几何定标方法 |
| CN113313769B (zh) * | 2021-06-11 | 2022-08-30 | 湖北工业大学 | 一种光学卫星多面阵传感器片间无缝几何定标方法 |
| CN113884036A (zh) * | 2021-09-22 | 2022-01-04 | 浙江省水利河口研究院(浙江省海洋规划设计研究院) | 一种坐标转换方法及装置、电子设备、存储介质 |
| CN113884036B (zh) * | 2021-09-22 | 2023-09-29 | 浙江省水利河口研究院(浙江省海洋规划设计研究院) | 一种坐标转换方法及装置、电子设备、存储介质 |
| CN114596342A (zh) * | 2022-03-22 | 2022-06-07 | 感知天下(北京)信息科技有限公司 | 基于卫星影像的并行匹配方法、存储介质及计算机设备 |
| CN114596342B (zh) * | 2022-03-22 | 2023-03-14 | 感知天下(北京)信息科技有限公司 | 基于卫星影像的并行匹配方法、存储介质及计算机设备 |
| CN114998399A (zh) * | 2022-05-20 | 2022-09-02 | 中国人民解放军61540部队 | 一种异源光学遥感卫星影像立体像对预处理方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109166143A (zh) | 一种大区域网立体测绘卫星影像匹配方法 | |
| Shean et al. | An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery | |
| Fahnestock et al. | Rapid large-area mapping of ice flow using Landsat 8 | |
| Bosch et al. | A multiple view stereo benchmark for satellite imagery | |
| Wolfe et al. | Achieving sub-pixel geolocation accuracy in support of MODIS land science | |
| CN105046251B (zh) | 一种基于环境一号卫星遥感影像的自动正射校正方法 | |
| CN107063296B (zh) | 一种卫星遥感传感器在轨辐射定标方法 | |
| CN108765488B (zh) | 一种基于阴影的高分辨率遥感影像建筑物高度估测方法 | |
| Gance et al. | Target Detection and Tracking of moving objects for characterizing landslide displacements from time-lapse terrestrial optical images | |
| CN106940887B (zh) | 一种gf-4卫星序列图像云与云下阴影检测方法 | |
| Ermolaev et al. | Automated construction of the boundaries of basin geosystems for the Volga Federal District | |
| CN111650579B (zh) | 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 | |
| CN111003214B (zh) | 基于云控制的国产陆地观测卫星姿轨精化方法 | |
| Crawford et al. | The 50-year Landsat collection 2 archive | |
| CN106887016A (zh) | 一种gf‑4卫星序列图像自动相对配准方法 | |
| CN116310883B (zh) | 基于遥感图像时空融合的农业灾害预测方法及相关设备 | |
| Xin et al. | High-precision co-registration of orbiter imagery and digital elevation model constrained by both geometric and photometric information | |
| Ren et al. | Automated SAR reference image preparation for navigation | |
| CN108876829B (zh) | 基于非线性尺度空间及径向基函数的sar高精度配准方法 | |
| Yan et al. | Topographic reconstruction of the “Tianwen-1” landing area on the Mars using high resolution imaging camera images | |
| CN111161186B (zh) | 一种推扫式遥感器通道配准方法及装置 | |
| CN106126839B (zh) | 一种三线阵立体测绘卫星成像仿真方法和系统 | |
| Restrepo et al. | Building Lunar Maps for Terrain Relative Navigation and Hazard Detection Applications | |
| Liu et al. | High-spatial-resolution nighttime light dataset acquisition based on volunteered passenger aircraft remote sensing | |
| Jovanovic et al. | Multi-angle geometric processing for globally geo-located and co-registered MISR image data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| CB02 | Change of applicant information | ||
| CB02 | Change of applicant information |
Address after: 101399 No. 2 East Airport Road, Shunyi Airport Economic Core Area, Beijing (1st, 5th and 7th floors of Industrial Park 1A-4) Applicant after: Zhongke Star Map Co.,Ltd. Address before: 101399 No. 2 East Airport Road, Shunyi Airport Economic Core Area, Beijing (1st, 5th and 7th floors of Industrial Park 1A-4) Applicant before: GEOVIS TECHNOLOGY (BEIJING) Co.,Ltd. |
|
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190108 |