CN102073874A - 附加几何约束的航天三线阵ccd相机多影像立体匹配方法 - Google Patents

附加几何约束的航天三线阵ccd相机多影像立体匹配方法 Download PDF

Info

Publication number
CN102073874A
CN102073874A CN 201010622773 CN201010622773A CN102073874A CN 102073874 A CN102073874 A CN 102073874A CN 201010622773 CN201010622773 CN 201010622773 CN 201010622773 A CN201010622773 A CN 201010622773A CN 102073874 A CN102073874 A CN 102073874A
Authority
CN
China
Prior art keywords
image
space
matching
matched
linear array
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.)
Granted
Application number
CN 201010622773
Other languages
English (en)
Other versions
CN102073874B (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.)
China Center for Resource Satellite Data and Applications CRESDA
Original Assignee
China Center for Resource Satellite Data and Applications CRESDA
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 China Center for Resource Satellite Data and Applications CRESDA filed Critical China Center for Resource Satellite Data and Applications CRESDA
Priority to CN 201010622773 priority Critical patent/CN102073874B/zh
Publication of CN102073874A publication Critical patent/CN102073874A/zh
Application granted granted Critical
Publication of CN102073874B publication Critical patent/CN102073874B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Processing (AREA)

Abstract

附加几何约束的航天三线阵CCD相机多影像立体匹配方法,(1)将原始三视影像的灰度均值和方差映射至给定的灰度均值和方差值;(2)将经过步骤(1)处理后的三视影像分别生成金字塔影像;(3)在生成的金字塔影像中的正视影像上提取特征点和边缘点作为待匹配点;(4)对生成的金字塔影像进行几何畸变预改正;(5)利用附加几何约束的多影像匹配方法匹配出待匹配点在前视、后视影像上的同名点;将上述匹配出的同名点作为观测值,利用最小二乘匹配方法提高同名点的匹配精度;(6)得到数字表面模型DSM;(7)对数字表面模型DSM中的点云进行滤波处理,滤除高频率信号或者噪声信号,生成数字高程模型DEM,完成航天三线阵CCD相机多影像立体匹配。

Description

附加几何约束的航天三线阵CCD相机多影像立体匹配方法
技术领域
本发明涉及一种航天三线阵CCD相机多影像立体匹配方法,属于对地观测与遥感技术。
背景技术
三线阵CCD(Charge Coupled Device)相机是实现航天遥感立体成像的重要方式之一,其稳定灵活的立体构成(同轨立体、三度重叠)、理想的基高比、立体影像的同时相等优点,为高精度目标点定位、高质量DEM(Digital ElevationModel)生成等提供了有效技术途径,显示出其独特的优势。
数字摄影测量与遥感、计算机视觉的核心问题是在不同的影像中高精度、自动化地提取相应同名点即影像匹配。因此,影像匹配技术已发展成为当今数字摄影测量与遥感、数字图像处理、计算机视觉、机器视觉的热门研究内容之一。在DEM(数字高程模型)自动生成、DOM(数字正射影像)自动生成、城市建筑物变形观测、城市三维重建、工业生产自动化等方面具有巨大的应用潜力。
立体影像匹配是高精度、自动化提取DEM高程数据的必经步骤,是影响DEM高程几何精度的重要环节,也是自动生成DEM的基础。在提高DEM高程几何精度、提高产品线作业效率等方面有重要现实意义。
在高可靠性的立体匹配的基础上自动生成数字高程模型数据,还需要解决匹配结果中房屋、植被等的滤除问题,即从DSM到DEM的自动化处理。如何在确保地形特征细节完整保留的基础上,完全滤除房屋、植被等对象,具有一定的技术难度。自动影像匹配效果一般由匹配精度、匹配成功率、误匹配率和微细地形地貌的表现能力等几个指标进行评价,使匹配结果可以更好的反映地形地貌特征,减少所需的人工编辑工作量。
影像匹配作为遥感影像数据处理的核心,近年来仍处在不断发展与完善之中。从当前的发展来看,其发展趋势可归纳如下:
(1)向多影像处理方向发展
近年来,大面阵航空数码相机、ADS-40三线阵航空数码相机、三线阵卫星摄影机、SPOT卫星传感器等大大地增加了数据量,同时也为影像自动匹配提供了多余的观测数据。因此,多影像匹配成为研究的重点内容之一。欧空局与德国Inpho公司开发了的处理大重叠度影像的高性能航空数字摄影测量处理软件,都是使用多像匹配算法。另外,大重叠度影像在近景摄影测量、计算机视觉领域也大量出现,得到了应用广泛。
(2)强调场景与影像成像几何约束
针对不同的应用目的,不同的传感器特点,不同的物方景观,探索其内在的几何规律和约束关系,以增加匹配的可靠性,在当前的匹配算法中已逐步开始应用。
发明内容
本发明的技术解决问题是:克服现有技术的不足,提供一种附加几何约束的航天三线阵CCD相机多影像立体匹配方法,该方法充分利用三线阵影像的三度多视角重复观测,能够大幅提高特征点的匹配数量、匹配精度和可靠性,达到改善航天三线阵前/后视影像由于视角差异而导致影像变形过大、实现影像自动匹配的目的。
本发明的技术解决方案是:附加几何约束的航天三线阵CCD相机多影像立体匹配方法,其特征在于步骤如下:
(1)对航天三线阵CCD相机原始三视影像进行预处理,将原始三视影像的灰度均值和方差映射至给定的灰度均值和方差值;
(2)将经过步骤(1)处理后的三视影像分别生成金字塔影像;
(3)在生成的金字塔影像中的正视影像上提取特征点和边缘点作为待匹配点;
(4)对生成的金字塔影像进行几何畸变预改正,消除金字塔影像中的几何畸变;
(5)利用附加几何约束的多影像匹配方法匹配出待匹配点在前视、后视影像上的同名点;将上述匹配出的同名点作为观测值,利用最小二乘匹配方法提高同名点的匹配精度;
(6)利用摄影测量中的前方交会关系及步骤(5)得到的三视影像上的同名点,求出地面点的坐标,得到数字表面模型DSM;
(7)对数字表面模型DSM中的点云进行滤波处理,滤除高频率信号或者噪声信号,生成数字高程模型DEM,完成航天三线阵CCD相机多影像立体匹配。
所述步骤(1)采用Wallis滤波器对航天三线阵CCD相机原始三视影像进行预处理。
所述步骤(1)中的给定的灰度均值和方差分别为127和40~70,其中方差的选取随Wallis滤波中矩形区域的减小而减小。
所述步骤(4)中的几何畸变预改正步骤如下:
(4.1)对于步骤(3)中提取的每个待匹配点p0,p0在空间中的对应点为P(X,Y,Z),分别在空间中取一个高程为Z的面元Γw,利用影像的定向模型,将面元Γw投影到正视影像、前视或后视影像上的两个四边形Γ和Γ′;
(4.2)利用四边形Γ和Γ′在正视影像、前视或后视影像之间定义一个局域仿射变换;
(4.3)利用局域仿射变换将前视或后视影像中的四边形Γ′重采样为新的影像窗口Γ″;
(4.4)使用基准影像窗口Γ和重采样的影像窗口Γ″进行相关匹配。
所述步骤(5)中附加核线约束的多影像匹配方法步骤为:
(5.1)以正视影像为参考影像,前视、后视影像为搜索影像,对正视影像上的待匹配点及其近似高度值,由已知的影像外方位元素并通过直线拟合的方法确定待匹配点在搜索影像上的近似核线;
(5.2)计算参考影像与搜索影像的相应的相关窗口之间的相关系数值ρ(p0,Z),并求出参考影像上待匹配点p0处的相关系数平均值Sρ
(5.3)寻找Z值,使Sρ最大的Z值为p0对应的正确高度值,进而得到参考影像p0处的唯一正确匹配。
本发明与现有技术相比有益效果为:
(1)本发明针对现有的航天遥感立体影像匹配方法成功率不高、匹配条件单一的情况,充分利用三线阵影像的三度多视角重复观测,大幅提高特征点的匹配数量、匹配精度和可靠性,改善航天三线阵前/后视影像由于视角差异而导致的影像变形过大、难以实现自动匹配的技术难题。
(2)本发明采用采用Wallis滤波器对航天三线阵CCD相机原始三视影像进行预处理,在提取影像中的特征点时可提高特征点的数量和精度,而在影像匹配中则提高了匹配结果的可靠性和稳定性。
(3)本发明在互相关匹配之前根据影像的定向参数和必要的地形信息消除影像中由于传感器方位元素不同和地形坡度所产生的影像几何畸变的干扰。
(4)本发明采用附加核线(几何)约束的多影像匹配方法,可以同时匹配三视影像,在匹配过程中间接地运用了核线(几何)条件,匹配成功率较传统单独像对的方法有很大的提高,可有效地解决错误匹配的问题,提高匹配精度,而且能有效解决影像遮挡问题。
附图说明
图1为本发明流程图;
图2为本发明附加核线约束的多影像匹配示意图;
图3为本发明点线混合松弛算法示意图;
图4本发明为金字塔匹配示意图;
图5为本发明金字塔分层示意图,其中a为四像元平均法示意图,b为九像元平均法示意图;
图6为本发明每层影像中的匹配示意图;
图7为本发明几何畸变预改正方法示意图。
具体实施方式
下面结合附图详细介绍本发明的实现过程。如图1所示,本发明具体步骤如下:
(1)对航天三线阵CCD相机原始三视影像进行预处理(即影像增强),将原始三视影像的灰度均值和方差映射至给定的灰度均值和方差值;
影像增强:通常情况下获取的原始影像在灰度和纹理等方面存在一定程度的差别,例如由成像时间、成像角度、以及光照条件等方面的不同引起影像的辐射畸变。另外原始卫星影像还普遍存在色调偏暗,对比度偏低等问题。这些因素对于特征提取和影像匹配都会造成严重影响,因此需要在影像匹配前对影像进行增强,提高匹配的可靠性和成功率。
Wallis滤波作为一种比较特殊的滤波器,可增强原始影像的反差并同时抑制噪声,特别是可以大大增强影像中不同尺度的影像纹理模式,所以在提取影像中的特征点时可提高特征点的数量和精度,而在影像匹配中则提高了匹配结果的可靠性和精度。该滤波器的目的是将影像的灰度均值和方差(即影像的动态范围)映射至给定的灰度均值和方差。这是一种局部影像变换,它使在影像不同位置处的灰度方差和灰度均值都具有近似相等的数值,即影像反差小的区域反差增强,影像反差大的区域反差减小,使得影像中灰度的微小变化信息得到增强。上述特性使Wallis滤波器对低反差影像和反差不均匀的影像有特殊作用。由于Wallis滤波器在计算影像的局部灰度方差和均值时使用了平滑算子,所以在增强影像有用信息的同时抑制了噪声,提高了影像的信噪比,使影像中存在的极为模糊的纹理模式得到增强。因此,处理后的影像虽然在视觉效果上有些像一幅噪声影像,但在进行特征点提取或立体匹配时,其效果要理想得多。
Wallis滤波作为一种局部的影像变换,使影像反差小的区域反差增大,影像反差大的区域反差减小,使得影像中灰度的微小变换信息得到增强。本发明方法在计算影像的局部灰度均值和方差时引入了平滑算子,能抑制噪声,改善影像质量,提高影像匹配成功率。
Wallis滤波器表现形式为:
f(x,y)=[g(x,y)-mg]·(csf)/[csg+(1-c)sf]+bmf+(1-b)mg      (1)
或f(x,y)=g(x,y)r1+r0
r1=(csf)/(csg+sf/c),r0=bmf+(1-b-r1)mg
其中,g(x,y)为原始影像的灰度值,f(x,y)为Wallis变换后结果影像的灰度值;mg为原始影像的局部灰度均值,sg为原始影像的局部灰度方差值;mf为结果影像局部灰度均值的目标值,Sf为结果影像的局部灰度方差值的目标值;c∈[0,1]为影像方差的扩展常数;b∈[0,1]为影像的亮度系数。
Wallis滤波实现过程如下:
a.将原始影像分为互不重叠的矩形区域,矩形区域的尺度对应于要增强的纹理模式的尺度;
b.计算各矩形区域的灰度均值与方差;
c.给定的灰度均值和方差分别设定为127和40~70之间的数值,其中后者应随矩形区域尺度的减小而减小,以防止大量像素的灰度值被饱和(即落于[0,255]之外),然后计算出各区域的Wallis滤波器乘性系数r1和加性系数r0
d.由于各矩形区域互不重叠,所以数字影像的任一像素的系数r1,r0采用双线性内插得到,并计算出所有像素新的灰度值。
r1=(csf)/(csg+sf/c),r0=bmf+(1-b-r1)mg
(2)将经过步骤(1)处理后的三视影像分别生成金字塔影像;
根据影像灰度相关的谱分析结果可知,当信号中的高频成分较少时,相关函数曲线较平滑,但相关的拉入范围较大;反之,当高频成分较多时,相关函数曲线较陡,相关精度较高,但相关拉入范围较小。此外,当信号中存在高频窄带随机噪声或信号中存在较强的高频信号时,相关函数出现多峰值,因此容易导致误匹配。综合考虑相关结果的正确性(或称为可靠性)与精度(准确性),可以得出从粗到精的相关策略:即先通过低通滤波进行粗相关,找到同名点的粗略位置,然后利用高频信息进行精确相关。
对原始影像进行低通滤波进行粗相关,将其结果作为预测值,逐步加入较高的频率成分,在逐渐变小的搜索区中进行相关,最后用原始信号,以得到最好的精度。对于二维影像逐次进行低通滤波,并增大采样间隔,得到一个像元总数逐渐变小的影像序列,依次在这些影像对中进行相关计算,并将前一级相关的结果作为下一级相关的初值,直至原始影像,实现对影像的分频道相关匹配。由于将这些影像叠置起来很像一座金字塔,因而通常称之为金字塔影像或分层结构影像,相应的匹配策略称之为金字塔影像匹配策略。图4为金字塔匹配示意图。
金字塔影像的建立可按l×l像元变换成一个像元逐层形成。一般取l=2(如图5a)的较多,但取l=3(如图5b)是计算量最小的方法,匹配结果从上一层传递到下一层时正好与3×3个像元的中心像元相对应,而l=2时上一层的结果与下一层2×2个像元的公共角点相对应。将原始影像称为第零层,则第一层影像的每一个像素相当于零层的l×l个像素,第k层影像的每一个像素相当于零层的(l×l)k个像素。
在一层影像中的匹配情况可用图6表示,模板f遍历被搜索的g影像,利用相关性比较、寻找最大相关系数的像元,可以看到被搜索图越大,匹配速度越慢;模板越小,匹配速度越快,匹配速度影响取决于模板的尺寸和搜索窗口的大小。另外,匹配的可靠性是随着搜索窗口的增大而增强,但是,如果两个匹配影像存在明显差异,大的搜索窗口反而会使匹配精度降低,甚至匹配错误。太小的搜索窗口情况下,如果待匹配区域没有足够的特征信息也会出现误匹配。所以窗口大小对匹配成功率有很大的影响。窗口大小一般根据经验值可基本确定,一般在正视影像上以特征点为中心,大小为7×7或者9×9等奇数个像素,在前视或者后视影像上以特征点对应的同名点为中心,大小与正视影像上窗口大小相同或者略大于。
(3)在生成的金字塔影像中的正视影像上提取特征点、格网点以及边缘点作为待匹配点;
特征点提取:匹配窗口内的灰度影像信噪比越大,匹配的精度越高。为了提高影像匹配的精度和鲁棒性,影像匹配时应基于点特征进行匹配。点特征主要指影像上的明显点,如角点、圆点等。提取点特征的算子称为兴趣算子或有利算子,即运用某种算法从影像中提取出我们感性的,即有利于某种目的的点。
本发明利用Harris算子在影像上提取特征点。Harris算子可给出与自相关函数相联系的矩阵M。M阵的特征值是自相关函数的一阶曲率,如果两个曲率值都高,则认为该点是特征点。
Harris算子只涉及图像的一阶导数:
M = G ( s ~ ) ⊗ g x 2 g x g y g x g y g y 2 - - - ( 2 )
但是解特征向量需要比较多的计算量,且两个特征值的和等于矩阵M的迹,两个特征值的积等于矩阵M的行列式。所以用下式来判定角点质量:(k常取0.04-0.06)
I=det(M)-k·tr2(M),k=0.04                     (3)
式中,gx是x方向的梯度,gy是y方向的梯度;
Figure BSA00000411784200082
为高斯模板,
Figure BSA00000411784200083
为卷积操作;I为每点的兴趣值,det是矩阵的行列式,tr是矩阵的迹,k为默认常数。
(4)对生成的金字塔影像进行几何畸变预改正,消除金字塔影像中的几何畸变;
基于灰度互相关的影像匹配算法经常受到由于摄影机方位不同和由于地形坡度所产生的影像几何畸变等的严重干扰。因此为了解决影像中几何畸变对灰度互相关匹配的影响,在互相关匹配之前根据影像的定向参数和必要的地形信息采用如图7所示的几何算法在影像匹配前消除影像中的几何畸变。具体为:
(4.1)对于步骤(3)中提取的每个待匹配点p0,p0在空间中的对应点为P(X,Y,Z),分别在空间中取一个高程为Z的面元Γw,利用影像的定向模型,将面元Γw投影到正视影像、前视或后视影像上的两个四边形Γ和Γ′;
(4.2)利用四边形Γ和Γ′在正视影像、前视或后视影像之间定义一个局域仿射变换;
x′=a0+a1x+a2y
y′=b0+b1x+b2y
其中,x′和y′为地面控制点投影到影像面的坐标,x和y为控制点在影像上的量测坐标,ai和bi(i=0,1,2)为仿射变换系数。
(4.3)利用局域仿射变换将前视或后视影像中的四边形Γ′重采样为新的影像窗口Γ″;
(4.4)使用基准影像窗口Γ和重采样的影像窗口Γ″进行相关匹配,即完成对生成的金字塔影像的几何畸变预改正。
(5)利用附加核线约束的多影像匹配方法匹配出待匹配点在前视、后视影像上的同名点;将上述匹配出的同名点作为观测值,利用最小二乘匹配方法提高同名点的匹配精度;
多影像匹配方法如图2所示。图中匹配方法的搜索空间定义为一定高程范围内基准点所对应的成像光线,高程的搜索步长设定为使得待匹配影像上沿同名核线的最小步长为1个像元。与此同时,将传统的特征点整体松弛匹配也扩展为点、线特征的混合松弛匹配如图3所示。最终的匹配结果由归一化相关系数(是高程的函数)峰值所对应的高程确定(确定峰值时要进行二次曲线拟合,从而达到子像元级的匹配精度)。因此算法实际上是在核线约束条件下同时匹配所有影像,而不是先进行各个立体像对的匹配再把匹配结果综合起来考虑。
这种方法主要有三个优点:一是匹配点的确定是同时综合考虑了所有影像间的互相关结果,因此可以有效的减少误匹配;二是匹配的最终结果为基准点对应成像光线上的空间点;三是由于空间点坐标的确定实际上等价于多影像的空间前方交会,因此相对于单个立体模型(两张影像)而言具有更高的高程精度。
附加核线约束的多影像匹配方法步骤为:
(5.1)以正视影像为参考影像,前视、后视影像为搜索影像,对正视影像上的待匹配点及其近似高度值,由已知的影像外方位元素(空中三角测量提供)并通过直线拟合的方法确定待匹配点在搜索影像上的近似核线;
(5.2)计算参考影像与搜索影像的相应的相关窗口之间的相关系数值ρ(p0,Z),并求出参考影像上待匹配点p0处的相关系数平均值Sρ
ρ ( p 0 , Z ) = Σ s ∈ W [ I 0 ( s ) - I 0 ‾ ] × { I i [ s i ( Z ) ] - I i ‾ } Σ s ∈ W [ I 0 ( s ) - I 0 ‾ ] 2 · Σ s ∈ W { I i [ s i ( Z ) ] - I i ‾ } 2 - - - ( 4 )
式中
I 0 ‾ = 1 M × N Σ s ∈ W I 0 ( s ) - - - ( 5 )
I ‾ i = 1 M × N Σ s ∈ W I 0 [ s i ( Z ) ] - - - ( 6 )
其中,W和s分别代表参考影像上点p0处的相关窗口和窗口内的一个像素;M与N代表了相关窗口W的大小;si(Z)代表第i张搜索影像上对应点s的像点。s0(·)和si(·)分别是参考影像I0与搜索影像I1的传感器模型。运用近似DSM以及影像的定向参数,通过窗口变换过程可以计算出si(Z)。通过(4)式可以看到,ρ是高度值Z的函数,Z的取值区间为[Z0-ΔZ,Z0+ΔZ]。因此,如果给定参考影像上一个像点以及其近似高度值Z0,在一个可能高程误差ΔZ的情况下,所有单独立体像对的ρ函数可以定义到一个唯一的框架中。在此,定义点p0处关于Z值的ρ值的总和的平均值为Sρ,即:
S ρ = 1 2 Σ i = 1 2 ρ i ( p 0 , Z ) - - - ( 7 )
(5.3)于是通过寻找Z值Z∈[Z0-ΔZ,Z0+ΔZ],使Sρ最大的Z值为p0对应的正确高度值,进而得到参考影像p0处的唯一正确匹配。这里高度误差ΔZ决定了在搜索影像上沿着相应核线的搜索距离。
最小二乘匹配方法
由前所述,附加核线约束条件的多影像匹配可以获得子像元精度的匹配结果。为了进一步提高匹配精度,可以在原始影像匹配后,增加一级最小二乘匹配。在本算法中,最小二乘匹配解求的未知数只有匹配点的空间坐标Xk,Yk,Zk(对于特征点而言k=1,对于特征线而言k=1…N,N为特征线中的节点数)和待匹配影像(前视或后视影像)与基准影像(正视影像)间的灰度辐射畸变参数(gi′=ai+bi*g0)。另外除了根据影像窗口中的灰度差列出误差方程外,还需要根据共线约束条件列出相应的误差观测方程,如果是直线则需要添加直线约束误差方程,曲线则添加三维Cardinal样条的误差观测方程。由于该方法充分利用了影像窗口内的信息以及各种空间几何约束条件进行平差计算,因此可以使影像匹配可以达到1/3甚至1/10像元的精度,确保影像匹配精度可达到子像元级。
现假定两幅影像只存在辐射变形,即右片灰度分布g2相对于左片灰度分布g1存在线形变形,可用下式表示:
g1(x1,y1)+n1=h0+h1g2(x2,y2)+n2                   (11)
其中h0、h1为发生畸变改正参数;
n1、n2分别为g1,g2中存在的随机噪声。
此时仅考虑辐射线性畸变的最小二乘匹配的数学模型为:
v=h0+h1g2-(g1-g2)                                  (12)
但在实际应用中,两幅影像不仅存在辐射畸变,还存在几何变形,只有充分考虑影像的几何变形和辐射畸变,才能获得最佳的匹配结果。由于影像匹配的窗口很小,所以几何变形只考虑一次畸变:
x 2 = a 0 + a 1 x 1 + a 2 y 1 y 2 = b 0 + b 1 x 1 + b 2 y 1 - - - ( 13 )
式中,x1,y1为左片像点坐标,x2,y2为右片像点坐标。
ai,bi(i=0,1,2)为仿射变形系数,将式(13)代入(11)中,得:
g1(x1,y1)+n1=h0+h1g2(a0+a1x1+a2y1,b0+b1x1+b2y1)+n2          (14)
线性化后,即得最小二乘影像匹配的误差方程式:
v=dh0+c2dh1+c3da0+c4da1+c5da2+c6db0+c7db1+c8db2-Δg          (15)
式中dh0,dh1,da0,da1,da2,db0,db1,db2为待定参数的改正数,Δg是相应像素的灰度差,ci(i=1…8)为方程系数。
c 1 = 1 c 2 = g 2 c 3 = ∂ g 2 ∂ x 2 ∂ x 2 ∂ a 0 = g x ‾ c 4 = ∂ g 2 ∂ x 2 ∂ x 2 ∂ a 1 = x g x ‾ c 5 = ∂ g 2 ∂ x 2 ∂ x 2 ∂ a 2 = y g x ‾ c 6 = ∂ g 2 ∂ y 2 ∂ y 2 ∂ b 0 = g y ‾ c 7 = ∂ g 2 ∂ y 2 ∂ y 2 ∂ b 1 = x g y ‾ c 8 = ∂ g 2 ∂ y 2 ∂ y 2 ∂ b 2 = y g y ‾ - - - ( 16 )
在数字影像中,由于灰度是按规则格网排列的离散阵列,所以偏导数可用差分代替,即:
g x ‾ = 1 2 [ g 2 ( I + 1 , J ) - g 2 ( I - 1 , J ) ] g y ‾ = 1 2 [ g 2 ( 1 , J + I ) - g 2 ( 1 , J - I ) ] - - - ( 17 )
按式(5)建立误差方程式,其矩阵形式为:
V=CX-L                                 (18)
然后对式(18)进行最小二乘计算。
最小二乘匹配方法的具体步骤为步:
1)设(x1,y1),(x2,y2)分别是参考影像、搜索影像上对应点的坐标,根据几何变形参数ai,bi(i=0,1,2),对搜索影像窗口内坐标进行变换;一般取初值:h0=0,h1=1,a0=0,a1=1,b0=0,b1=0,b2=1;
2)重采样。由于计算所得的坐标x2,y2,一般不可能刚好是改为搜索影像阵列中的整数行列号,因而利用双线性内插进行重采样获取g2(x2,y2);
3)辐射畸变改正。h0,h1,是辐射畸变改正参数,g1,g2分别是参考影像、搜索影像的灰度,并有g1(x1,y1)=h0+h1g2(x2,y2);
4)根据参考影像、搜索影像灰度阵列g1(x1,y1)和h0+h1g2(x2,y2),计算参考影像、搜索影像匹配窗口中的相关系数P,如果P>1,转到步骤7);
5)采用最小二乘影像匹配(最小均方误差原则),求解变形参数改正值:dh0,dh1,da0,…db2
6)计算变形参数h0,h1,a0,a1,a2,b0,b1,b2,转到步骤1);
7)计算最佳匹配点位置:以参考影像窗口中的灰度重心为目标点,利用步骤1)、2)中的h0,…,b2对其坐标作校正后得右影像的重采样点作为匹配点。
一般来说,实现高精度最小二乘影像匹配需要比较准确的初值才能获得较好的结果。
(6)利用摄影测量中的空间前方交会关系及步骤(5)得到的同名点,求出地面点的坐标,得到数字表面模型DSM;
摄影测量中的空间前方交会关系具体可以参照《摄影测量学》教材,此处不进行详细的描述。
(7)对数字表面模型DSM中的点云进行滤波处理,滤除高频率信号或者噪声信号,生成数字高程模型DEM,完成航天三线阵CCD相机多影像立体匹配。
(8)采用基于Wallis滤波器的方法进行多景影像的色彩一致性等处理,使处理后的正射影像具有全局一致的视觉效果,同时最大程度地保持原始影像的色彩和亮度等信息。
(9)正射影像纠正主要利用三线阵影像及其辅助参数,根据生成的DEM数据进行纠正,获得大范围正射影像图。根据纠正后的多景正射影像,确定其重叠区域和镶嵌线的具体位置,并将区域正射影像拼接为一个特殊格式的巨型文件,便于存储管理及交互式编辑和浏览。
上述步骤(8)、(9)为采用本发明立体匹配后对影像的后续处理步骤,该步骤为本领域的惯常做法,不进行详细的描述。
本发明未详细说明的部分属于本领域技术人员公知常识。

Claims (5)

1.附加几何约束的航天三线阵CCD相机多影像立体匹配方法,其特征在于步骤如下:
(1)对航天三线阵CCD相机原始三视影像进行预处理,将原始三视影像的灰度均值和方差映射至给定的灰度均值和方差值;
(2)将经过步骤(1)处理后的三视影像分别生成金字塔影像;
(3)在生成的金字塔影像中的正视影像上提取特征点和边缘点作为待匹配点;
(4)对生成的金字塔影像进行几何畸变预改正,消除金字塔影像中的几何畸变;
(5)利用附加几何约束的多影像匹配方法匹配出待匹配点在前视、后视影像上的同名点;将上述匹配出的同名点作为观测值,利用最小二乘匹配方法提高同名点的匹配精度;
(6)利用摄影测量中的前方交会关系及步骤(5)得到的三视影像上的同名点,求出地面点的坐标,得到数字表面模型DSM;
(7)对数字表面模型DSM中的点云进行滤波处理,滤除高频率信号或者噪声信号,生成数字高程模型DEM,完成航天三线阵CCD相机多影像立体匹配。
2.根据权利要求1所述的附加几何约束的航天三线阵CCD相机多影像立体匹配方法,其特征在于:所述步骤(1)采用Wallis滤波器对航天三线阵CCD相机原始三视影像进行预处理。
3.根据权利要求1所述的附加几何约束的航天三线阵CCD相机多影像立体匹配方法,其特征在于:所述步骤(1)中的给定的灰度均值和方差分别为127和40~70,其中方差的选取随Wallis滤波中矩形区域的减小而减小。
4.根据权利要求1所述的附加几何约束的航天三线阵CCD相机多影像立体匹配方法,其特征在于:所述步骤(4)中的几何畸变预改正步骤如下:
(4.1)对于步骤(3)中提取的每个待匹配点p0,p0在空间中的对应点为P(X,Y,Z),分别在空间中取一个高程为Z的面元Γw,利用影像的定向模型,将面元Γw投影到正视影像、前视或后视影像上的两个四边形Γ和Γ′;
(4.2)利用四边形Γ和Γ′在正视影像、前视或后视影像之间定义一个局域仿射变换;
(4.3)利用局域仿射变换将前视或后视影像中的四边形Γ′重采样为新的影像窗口Γ″;
(4.4)使用基准影像窗口Γ和重采样的影像窗口Γ″进行相关匹配。
5.根据权利要求1所述的附加几何约束的航天三线阵CCD相机多影像立体匹配方法,其特征在于:所述步骤(5)中附加核线约束的多影像匹配方法步骤为:
(5.1)以正视影像为参考影像,前视、后视影像为搜索影像,对正视影像上的待匹配点及其近似高度值,由已知的影像外方位元素并通过直线拟合的方法确定待匹配点在搜索影像上的近似核线;
(5.2)计算参考影像与搜索影像的相应的相关窗口之间的相关系数值ρ(p0,Z),并求出参考影像上待匹配点p0处的相关系数平均值Sρ
(5.3)寻找Z值,使Sρ最大的Z值为p0对应的正确高度值,进而得到参考影像p0处的唯一正确匹配。
CN 201010622773 2010-12-29 2010-12-29 附加几何约束的航天三线阵ccd相机多影像立体匹配方法 Active CN102073874B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN 201010622773 CN102073874B (zh) 2010-12-29 2010-12-29 附加几何约束的航天三线阵ccd相机多影像立体匹配方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN 201010622773 CN102073874B (zh) 2010-12-29 2010-12-29 附加几何约束的航天三线阵ccd相机多影像立体匹配方法

Publications (2)

Publication Number Publication Date
CN102073874A true CN102073874A (zh) 2011-05-25
CN102073874B CN102073874B (zh) 2013-04-24

Family

ID=44032408

Family Applications (1)

Application Number Title Priority Date Filing Date
CN 201010622773 Active CN102073874B (zh) 2010-12-29 2010-12-29 附加几何约束的航天三线阵ccd相机多影像立体匹配方法

Country Status (1)

Country Link
CN (1) CN102073874B (zh)

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102506824A (zh) * 2011-10-14 2012-06-20 航天恒星科技有限公司 一种城市低空无人机系统生成数字正射影像图的方法
CN103076004A (zh) * 2011-12-27 2013-05-01 范大昭 多基线立体匹配系统的断裂特征自适应窗口变化方法
CN103218787A (zh) * 2013-04-23 2013-07-24 国家测绘地理信息局卫星测绘应用中心 多源异构遥感影像控制点自动采集方法
CN103646389A (zh) * 2013-03-26 2014-03-19 中国科学院电子学研究所 一种基于几何模型的sar斜距图像同名点自动提取方法
CN105069843A (zh) * 2015-08-22 2015-11-18 浙江中测新图地理信息技术有限公司 一种面向城市三维建模的密集点云的快速提取方法
CN106846384A (zh) * 2016-12-30 2017-06-13 中国人民解放军61540部队 一种多视角大倾斜线阵影像匹配方法及装置
CN106886794A (zh) * 2017-02-14 2017-06-23 湖北工业大学 顾及高阶结构特征的异源遥感影像同名点匹配方法
CN108171731A (zh) * 2017-09-28 2018-06-15 中国矿业大学(北京) 一种顾及拓扑几何多要素约束的最小影像集自动优选方法
CN108846436A (zh) * 2018-06-13 2018-11-20 武汉朗视软件有限公司 一种多影像立体匹配方法
CN109166143A (zh) * 2018-07-06 2019-01-08 航天星图科技(北京)有限公司 一种大区域网立体测绘卫星影像匹配方法
CN111046906A (zh) * 2019-10-31 2020-04-21 中国资源卫星应用中心 一种面状特征点可靠加密匹配方法和系统
CN111583119A (zh) * 2020-05-19 2020-08-25 北京数字绿土科技有限公司 一种正射影像拼接方法、设备以及计算机可读介质
CN111693028A (zh) * 2020-06-23 2020-09-22 上海海洋大学 一种基于投影影像获取数字水深模型的方法
CN112729130A (zh) * 2020-12-29 2021-04-30 四川天奥空天信息技术有限公司 卫星遥感测量树木冠层高度的方法
CN112765095A (zh) * 2020-12-24 2021-05-07 山东省国土测绘院 一种立体测绘卫星影像数据归档方法和系统
CN113159103A (zh) * 2021-02-24 2021-07-23 广东拓斯达科技股份有限公司 图像匹配方法、装置、电子设备以及存储介质
CN114612631A (zh) * 2022-03-02 2022-06-10 自然资源部重庆测绘院 一种基于InSAR技术的高精度无漏洞DSM提取方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004233138A (ja) * 2003-01-29 2004-08-19 Kurabo Ind Ltd 写真測量における計測点の対応付け方法及び装置
JP2007033157A (ja) * 2005-07-25 2007-02-08 Ntt Data Corp 画像解析装置、画像解析方法及びプログラム
CN101464149A (zh) * 2008-12-31 2009-06-24 武汉大学 Pos辅助航空影像匹配方法
CN101604018A (zh) * 2009-07-24 2009-12-16 中国测绘科学研究院 高分辨率遥感影像数据处理方法及其系统
CN101915913A (zh) * 2010-07-30 2010-12-15 中交第二公路勘察设计研究院有限公司 一种稳健的高分辨率卫星影像连接点自动匹配方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004233138A (ja) * 2003-01-29 2004-08-19 Kurabo Ind Ltd 写真測量における計測点の対応付け方法及び装置
JP2007033157A (ja) * 2005-07-25 2007-02-08 Ntt Data Corp 画像解析装置、画像解析方法及びプログラム
CN101464149A (zh) * 2008-12-31 2009-06-24 武汉大学 Pos辅助航空影像匹配方法
CN101604018A (zh) * 2009-07-24 2009-12-16 中国测绘科学研究院 高分辨率遥感影像数据处理方法及其系统
CN101915913A (zh) * 2010-07-30 2010-12-15 中交第二公路勘察设计研究院有限公司 一种稳健的高分辨率卫星影像连接点自动匹配方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
《ISPRS Journal of Photogrammetry & Remote Sensing》 20060228 Li Zhang,Armin Gruen Multi-image matching for DSM generation from IKONOS imagery 195-211 1-5 , *
《武汉大学学报·信息科学版》 20080930 张力等 基于多基线影像匹配的高分辨率遥感影像DEM的自动生成 943-946 1-5 第33卷, 第9期 *
《测绘学报》 20090630 袁修孝等 一种综合利用像方和物方信息的多影像匹配方法 220-222 1-5 第38卷, 第3期 *

Cited By (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102506824B (zh) * 2011-10-14 2014-08-27 航天恒星科技有限公司 一种城市低空无人机系统生成数字正射影像图的方法
CN102506824A (zh) * 2011-10-14 2012-06-20 航天恒星科技有限公司 一种城市低空无人机系统生成数字正射影像图的方法
CN103076004A (zh) * 2011-12-27 2013-05-01 范大昭 多基线立体匹配系统的断裂特征自适应窗口变化方法
CN103076004B (zh) * 2011-12-27 2015-09-09 范大昭 多基线立体匹配系统的断裂特征自适应窗口变化方法
CN103646389A (zh) * 2013-03-26 2014-03-19 中国科学院电子学研究所 一种基于几何模型的sar斜距图像同名点自动提取方法
CN103218787B (zh) * 2013-04-23 2015-08-19 国家测绘地理信息局卫星测绘应用中心 多源异构遥感影像控制点自动采集方法
CN103218787A (zh) * 2013-04-23 2013-07-24 国家测绘地理信息局卫星测绘应用中心 多源异构遥感影像控制点自动采集方法
CN105069843A (zh) * 2015-08-22 2015-11-18 浙江中测新图地理信息技术有限公司 一种面向城市三维建模的密集点云的快速提取方法
CN106846384A (zh) * 2016-12-30 2017-06-13 中国人民解放军61540部队 一种多视角大倾斜线阵影像匹配方法及装置
CN106886794A (zh) * 2017-02-14 2017-06-23 湖北工业大学 顾及高阶结构特征的异源遥感影像同名点匹配方法
CN108171731B (zh) * 2017-09-28 2021-12-14 中国矿业大学(北京) 一种顾及拓扑几何多要素约束的最小影像集自动优选方法
CN108171731A (zh) * 2017-09-28 2018-06-15 中国矿业大学(北京) 一种顾及拓扑几何多要素约束的最小影像集自动优选方法
CN108846436A (zh) * 2018-06-13 2018-11-20 武汉朗视软件有限公司 一种多影像立体匹配方法
CN109166143A (zh) * 2018-07-06 2019-01-08 航天星图科技(北京)有限公司 一种大区域网立体测绘卫星影像匹配方法
CN111046906A (zh) * 2019-10-31 2020-04-21 中国资源卫星应用中心 一种面状特征点可靠加密匹配方法和系统
CN111046906B (zh) * 2019-10-31 2023-10-31 中国资源卫星应用中心 一种面状特征点可靠加密匹配方法和系统
CN111583119A (zh) * 2020-05-19 2020-08-25 北京数字绿土科技有限公司 一种正射影像拼接方法、设备以及计算机可读介质
CN111693028A (zh) * 2020-06-23 2020-09-22 上海海洋大学 一种基于投影影像获取数字水深模型的方法
CN112765095A (zh) * 2020-12-24 2021-05-07 山东省国土测绘院 一种立体测绘卫星影像数据归档方法和系统
CN112729130A (zh) * 2020-12-29 2021-04-30 四川天奥空天信息技术有限公司 卫星遥感测量树木冠层高度的方法
CN113159103A (zh) * 2021-02-24 2021-07-23 广东拓斯达科技股份有限公司 图像匹配方法、装置、电子设备以及存储介质
CN113159103B (zh) * 2021-02-24 2023-12-05 广东拓斯达科技股份有限公司 图像匹配方法、装置、电子设备以及存储介质
CN114612631A (zh) * 2022-03-02 2022-06-10 自然资源部重庆测绘院 一种基于InSAR技术的高精度无漏洞DSM提取方法
CN114612631B (zh) * 2022-03-02 2023-06-09 自然资源部重庆测绘院 一种基于InSAR技术的高精度无漏洞DSM提取方法

Also Published As

Publication number Publication date
CN102073874B (zh) 2013-04-24

Similar Documents

Publication Publication Date Title
CN102073874B (zh) 附加几何约束的航天三线阵ccd相机多影像立体匹配方法
CN106780590B (zh) 一种深度图的获取方法及系统
CN104156968B (zh) 大面积地形复杂区域无人机序列影像快速无缝拼接方法
US5309522A (en) Stereoscopic determination of terrain elevation
CN103106688B (zh) 基于双层配准方法的室内三维场景重建方法
CN104156957B (zh) 一种稳定高效的高分辨率立体匹配方法
CN104820991B (zh) 一种基于代价矩阵的多重软约束立体匹配方法
CN104299228B (zh) 一种基于精确点位预测模型的遥感影像密集匹配方法
CN102184540B (zh) 基于尺度空间的亚像素级立体匹配方法
CN104237887B (zh) 一种sar遥感影像匹配方法
CN103559737A (zh) 一种对象全景建模方法
CN106023230B (zh) 一种适合变形图像的稠密匹配方法
CN104318583B (zh) 一种可见光宽带光谱图像配准方法
CN103093459A (zh) 利用机载LiDAR点云数据辅助影像匹配的方法
CN105139355A (zh) 一种深度图像的增强方法
CN102855628A (zh) 多源多时相高分辨率卫星遥感影像自动匹配方法
CN104166992A (zh) 基于网格变形的内容感知双目图像缩放方法
Li et al. Research on multiview stereo mapping based on satellite video images
Liao et al. A linear pushbroom satellite image epipolar resampling method for digital surface model generation
CN102944308B (zh) 一种时空联合调制干涉成像光谱仪姿态误差校正方法
CN105137431B (zh) 一种sar立体模型构建与量测方法
Zhao et al. Alignment of continuous video onto 3D point clouds
CN115719320B (zh) 基于遥感影像的倾斜校正稠密匹配方法
Kim et al. Automatic registration of LiDAR and optical imagery using depth map stereo
Fayard et al. Matching stereoscopic SAR images for radargrammetric applications

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant