CN109029379B - 一种高精度小基高比立体测绘方法 - Google Patents

一种高精度小基高比立体测绘方法 Download PDF

Info

Publication number
CN109029379B
CN109029379B CN201810589576.4A CN201810589576A CN109029379B CN 109029379 B CN109029379 B CN 109029379B CN 201810589576 A CN201810589576 A CN 201810589576A CN 109029379 B CN109029379 B CN 109029379B
Authority
CN
China
Prior art keywords
image
images
calculating
height
corrected
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
CN201810589576.4A
Other languages
English (en)
Other versions
CN109029379A (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.)
Beijing Institute of Space Research Mechanical and Electricity
Original Assignee
Beijing Institute of Space Research Mechanical and Electricity
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 Beijing Institute of Space Research Mechanical and Electricity filed Critical Beijing Institute of Space Research Mechanical and Electricity
Priority to CN201810589576.4A priority Critical patent/CN109029379B/zh
Publication of CN109029379A publication Critical patent/CN109029379A/zh
Application granted granted Critical
Publication of CN109029379B publication Critical patent/CN109029379B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C11/00Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
    • G01C11/04Interpretation of pictures
    • G01C11/06Interpretation of pictures by comparison of two or more pictures of the same area
    • G01C11/08Interpretation of pictures by comparison of two or more pictures of the same area the pictures not being supported in the same relative position as when they were taken

Landscapes

  • Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Image Processing (AREA)

Abstract

本发明一种高精度小基高比立体测绘方法,步骤为:(1)一台相机在轨道高度为H,基线长度为B的两个位置对地进行拍摄,采集获取两幅图像;(2)对两幅图像分别进行基于频域晶胞的总变分正则化MTFC预处理;(3)计算得到校正畸变后的两幅图像;(4)计算得到两幅核线重采样后的图像;(5)计算得到一幅匹配后的视差图;(6)求得整幅图各个像素点处的相对高程值;(7)通过约束K个控制点计算的绝对高程与真实高程的误差平方和最小,来求解参考面高度,通过参考面高程加上各个点处的相对高程值,求得整幅图像各个点处的绝对高程。

Description

一种高精度小基高比立体测绘方法
技术领域
本发明属于航天光学遥感技术领域,涉及一种高精度小基高比立体测绘方法。
背景技术
小基高比立体测绘方法与传统大基高比立体测绘方法完全不同。大基高比是通过图像匹配在两幅图像中找到同名点,进而将同名控制点坐标代入共线方程求得共线方程参数,将待求点像方坐标代入已求解好的共线方程中,求得待求点的物方三维坐标。而小基高比立体测绘方法的核心计算公式为像方相对高程等于两幅图像的平面视差除以基高比数值。早在2002年,法国人就利用SPOT5卫星上的全色影像和多光谱影像组成立体像对,对小基高比情况下获取数字高程模型进行了试验,结果表明在立体交会角仅为0.02时,仍能获取一定精度的高程信息,从而验证了小基高比立体测绘技术的可行性。
传统的大基高比立体测绘方法是逐点代入共线方程,逐点获取相应点处的高程值;单台相机的大基高比立体测绘对平台的姿控提出了更高的要求,要卫星平台灵活、敏捷度高、稳定性强、姿态机动能力强,实现难度更大;大基高比条件下,拍摄城市区域(高大建筑物密集)时,丢失高程信息的遮挡区域将增大,不利于城市立体测绘。
2007年第28卷《Journal of Mathematical Imaging and Vision》上,法国空间研究中心的Julie Delon等人发表的《Small Baseline Stereovision》论文通过理论推导和仿真分析,讨论了小基高比立体测绘方法的可行性及相关仿真实验验证,其不足之处在于:该方法仅通过两幅图像的平面视差除以基高比数值求得像方相对高程,没有对相机的内外方位元素进行校正,没有对原始采集到的图像进行合理的MTF补偿,不能计算出物方绝对高程。
发明内容
本发明解决的技术问题是:克服现有技术的不足,提供了一种高精度小基高比立体测绘方法,该方法提出了必须对采集到的图像进行MTF补偿预处理,并校正相机的内外方位元素,最后通过控制点计算参考面的绝对高程来计算整个重叠区域的物方绝对高程值,实现最终的地物高精度立体测绘。
本发明的技术方案是:一种高精度小基高比立体测绘方法,步骤如下:
1)相机在轨道高度为H,基线长度为B的两个位置对地进行拍摄,采集获取两幅被噪声污染后的退化图像,并建立所采集图像的退化模型
Figure GDA0002610363280000021
其中,g为采集到的被噪声污染后的退化图像,
Figure GDA0002610363280000022
为傅里叶逆变换算子,f为理想图像,n为噪声,H′为频域晶胞上的调制传递函数,*为卷积运算;
2)对两幅被噪声污染后的退化图像分别进行基于频域晶胞的总变分正则化MTFC预处理,得到放大后的两幅复原图像f1、f2
3)计算得到校正内方位元素畸变后的两幅图像fcal1、fcal2
4)计算得到校正外方位元素畸变后的两幅图像fext1、fext2
5)采用基于解最优化的相位相关法对校正了外方位元素畸变的两幅图像fext1、fext2进行亚像素级匹配,得到匹配后的视差图V;
6)利用视差图V计算得到整幅图各个像素点处的相对高程值ΔZ;
ΔZ=(s(a1-s(c1)))·GSD/(B/H)
其中,s(a1)、s(c1)分别为a1、c1点视差值,a1、c1为的地面任意两点A和C在第一幅图像上的成像点,P为探测器像元大小,f′为光学系统焦距,
Figure GDA0002610363280000031
为物方地面采样距离;
7)通过约束K个控制点计算的绝对高程与真实高程的误差平方和最小,求解参考面高度h,即地面点C处的绝对高程为
Figure GDA0002610363280000032
进而可计算整幅图上所有点处的绝对高程值。
所述2)步骤进行基于频域晶胞的总变分正则化MTFC预处理的具体方法为:
21)采用总变分最小化模型对退化图像矩阵g进行去噪,具体公式为:
Figure GDA0002610363280000033
Figure GDA0002610363280000034
Figure GDA0002610363280000035
其中,α为正则化参数,β为可调参数,
Figure GDA0002610363280000036
Df是f的支持域,即成像系统的频域晶胞;
Figure GDA0002610363280000037
表示散度,对于理想图像f=[fx,fy]T
Figure GDA0002610363280000038
22)计算成像系统的频域晶胞,其模型为
Figure GDA0002610363280000039
其中,
Figure GDA0002610363280000041
为点(ν,ω)处的系统调制传递函数数值,
Figure GDA0002610363280000042
为系统调制传递函数矩阵,其大小为K×L,K、L均为正整数;comba(v,w)表示采样边带;θalias为阈值,其取值范围为[1,0);
23)计算频域晶胞上的调制传递函数
Figure GDA0002610363280000043
进行反卷积,得到复原后的图像。
所述3)步骤计算得到校正内方位元素畸变后的两幅图像fcal1、fcal2的具体过程为:
31)计算得到相机镜头x、y两个方向上的几何畸变量Δx,Δy;
32)对两幅复原后的图像f1、f2进行校正,计算每幅图像校正后的图像坐标(xcal,ycal)对应的原坐标(x,y),即
x=xcal-abs(Δx)
y=ycal-abs(Δy)
其中,abs(·)为取绝对值操作;
33)经过双线性内插算法求得校正后坐标处(xcal,ycal)的像素值;
34)由上面的过程最终实现对复原后图像f1、f2的内方位元素校正,得到校正畸变后的两幅图像fcal1、fcal2
所述4)步骤计算得到校正外方位元素畸变后的两幅图像fext1、fext2的具体方法为:
41)计算获得校正外方位元素坐标值对应的原坐标值
Figure GDA0002610363280000051
其中,x、y分别为内方位元素校正后图像fcal1上像点的横、纵坐标;a1,a2,a3,b1,b2,b3,c1,c2,c3为第一幅图像的方向余弦;取xk、yk为整数,求得一系列的像点坐标(x,y);
对于图像fcal2,将y′k=yk代入下面的第一幅图像共线方程,求得对应外方位元素校正后坐标(x′k,y′k)的原坐标值(x′,y′):
Figure GDA0002610363280000052
其中,x′k为整数,x′、y′为内方位元素校正后图像fcal2上像点的横、纵坐标;a′1,b′1,…,c′3为第一幅图像的方向余弦,分别是第二幅图像相对于摄影基线的角方位元素的函数;
42)经过双线性内插算法求得校正后坐标处(xk,yk)、(x′k,y′k)的像素值;
43)将校正过内方位元素的图像fcal1、fcal2上的所有点进行41)步骤、42)步骤的核线重采样操作,得到校正了外方位元素畸变的两幅图像fext1、fext2
所述5)步骤中基于解最优化的相位相关法的具体过程为:
51)通过求解下面的最优化问题来得到目标区域中心点的相对平移量(Δx,Δy),即
Figure GDA0002610363280000061
其中,Q为含有噪声的局部互相关功率谱,
Figure GDA0002610363280000062
表示理论局部互相关功率谱,W是权重矩阵;wx、wy分别为x方向和y方向的权重系数,φ(Δx,Δy)为局部互相关功率谱误差函数;
52)对相位相关矩阵进行奇异值分解,即
Figure GDA0002610363280000063
式中,
Figure GDA0002610363280000064
1=diag(σ12,…,σr),r=rank(Q);通过奇异值分解得到qx,qy,进而求解得到相应的相位px,py;由px=kxx+bx,py=kyy+by得平移量
Figure GDA0002610363280000065
53)通过初始Q0(wx,wy)奇异值分解得到初始迭代平移量(Δx0,Δy0),计算得到Q1(wx,wy),即
Figure GDA0002610363280000066
54)迭代计算得到Qi+1(wx,wy)和
Figure GDA0002610363280000067
i=1,2…………
Figure GDA0002610363280000068
判断
Figure GDA0002610363280000069
是否小于条件阈值,如果
Figure GDA00026103632800000610
小于某一阈值则迭代结束,得到不同迭代次数对应的平移量(Δxi,Δyi),若不收敛,调整阈值再进行54)步骤迭代计算;
55)计算得到最终的平移量,即x、y方向的视差值为
Figure GDA0002610363280000071
所述31)步骤计算得到相机镜头x、y两个方向上的几何畸变量Δx,Δy的具体方法为:
Figure GDA0002610363280000072
Figure GDA0002610363280000073
其中,x、y为退化图像上像点的横、纵坐标;X、Y、Z为地面坐标系,a1,a2,a3,b1,b2,b3,c1,c2,c3为通过相机的姿态构建的旋转矩阵,x0、y0为x、y两个方向上的主点偏移,XS,YS,ZS为相机的轨道参数;Δx、Δy为x、y两个方向上的几何畸变量。
所述33)步骤经过双线性内插算法求得校正后坐标处(xcal,ycal)的像素值的具体方法为:
Figure GDA0002610363280000074
其中,Δd为原始变形图像的采样间距,g′为原坐标(x,y)处待求像素灰度值,g1,g2,g3,g4分别为其周围4个像元的灰度值。
所述42)步骤经过双线性内插算法求得校正后坐标处(xk,yk)、(x′k,y′k)的像素值的具体方法为:
Figure GDA0002610363280000081
其中,Δd′为原始变形图像的采样间距,g″为原坐标(x′,y′)处待求像素灰度值,g′1,g′2,g′3,g′4分别为原坐标(x′,y′)周围4个像元的灰度值;同理求出第二幅图像坐标处(x′k,y′k)的像素值。
所述1)步骤中相机在轨道高度为H,基线长度为B的两个位置对地进行拍摄,两个位置关于光轴方向对称。
B/H≤0.12。
本发明与现有技术相比的优点在于:
现有方法只是核线重采样后,通过两幅图像的平面视差除以基高比数值求得像方相对高程。而本发明通过限制复原后的混叠效应,构建了频域晶胞模型,利用频域晶胞上降晰函数进行反卷积,实现了采集图像无虚假添加的MTF补偿,弥补相机自身的不足;通过内外方位元素重采样进行系统几何误差校正,提高小基高比立体测绘精度;通过基于解最优化的相位相关匹配法,实现亚像素级的平面视差计算,提高最终的测绘精度;通过计算控制点计算参考面的绝对高程来计算整个重叠区域的物方绝对高程值,实现最终的地物高精度立体测绘。
附图说明
图1为本发明方法的流程框图;
图2为相机几何畸变图(50×50网格);
图3为双线性内插重采样;
图4为过A点的核平面;
图5(a)为原参考图;图5(b)为匹配后的视差图;
图6为三维重建效果图。
具体实施方式
如图1所示,本发明的实现过程为:
1)小基高比图像预处理方法
小基高比立体测绘图像预处理主要进行了MTF补偿,本发明提出了基于频域晶胞模型的总变分正则化MTFC算法。首先,整个采样系统的模型可以表示为
Figure GDA0002610363280000091
其中,g为采集到的被噪声污染后的退化图像矩阵,
Figure GDA0002610363280000092
为傅里叶逆变换算子,f为理想图像,n为噪声矩阵,H′为降晰函数傅里叶变换的模,*为卷积运算。基于频域晶胞模型的总变分正则化MTFC算法的技术途径如下:
a)采用总变分最小化模型进行去噪,即:
Figure GDA0002610363280000093
Figure GDA0002610363280000094
其中,
Figure GDA0002610363280000095
β>0是可调参数。利用该总变分最小化模型进行去噪,该模型迭代计算公式为
Figure GDA0002610363280000096
Figure GDA0002610363280000097
b)计算成像系统的频域晶胞,其模型为
Figure GDA0002610363280000098
其中,
Figure GDA0002610363280000099
为系统MTF,comba(v,w)表示采样边带,θalias为阈值,其取值范围为[1,0)。
c)小基高比立体测绘中需要利用频域晶胞进行H′的计算,由步骤b)中频域晶胞的计算公式可以得到与相应成像系统MTF最为匹配的频域晶胞,利用频域晶胞上的MTF进行图像复原,则有
Figure GDA0002610363280000101
基于频域晶胞的图像复原方法与传统方法的最大不同在于式(1)中的H′由式(3)计算求得。
由步骤a)的迭代计算公式可得到的
Figure GDA0002610363280000102
值,将上式计算的H′代人,进行反卷积计算,求解复原图像f1、f2
物理验证样机的试验数据在MTFC前最优高程精度为5.66个GSD,利用基于频域晶胞的总变分正则化MTFC方法的最优高程精度为4.76个GSD,基于频域晶胞的MTFC方法处理后的高程精度,相比于MTFC前提高了16.1%。
2)内方位元素校正
在小基高比立体测绘中,内方位元素校正主要是对由相机自身特性引起的几何变形进行校正。考虑改正项后的共线方程如式:
Figure GDA0002610363280000103
Figure GDA0002610363280000104
其中,x,y为像点坐标,X,Y,Z为地面点坐标,a1,a2,a3,b1,b2,b3,c1,c2,c3为通过相机的姿态构建的旋转矩阵,x0,y0为主点偏移,f′为光学系统焦距,XS,YS,ZS为相机的轨道参数。Δx,Δy为附加参数,主要由相机镜径向畸变系数和切向畸变系数、CCD安装误差系数等组成,主要依据实际定标过程中考虑哪些误差项。
通过精确测量墙面214个控制点的物方精确坐标,与图像上像方坐标匹配成控制点对,M个控制点对就有关于共线方程系数的M对方程组,通过对方程组求解得到修正后的共线方程参数,进而得到内方位元素。如图所示,为通过计算墙面靶标得到的相机几何畸变图。
由通过共线方程计算出的几何畸变量Δx,Δy,对两幅复原后的图像f1、f2进行校正,计算每幅图像校正后的图像坐标(xcal,ycal)对应的原坐标(x,y),即
x=xcal-abs(Δx)
y=ycal-abs(Δy) (5)
其中,abs(·)为取绝对值操作。
变形图形坐标不一定为整数,必须进行灰度重采样得到变形坐标点(x,y)的灰度值作为对应校正后坐标的灰度值。如图3所示,Δd为原始变形图像的采样间距,待求像素灰度值为g′,可由坐标点(x,y)周围4个像元的灰度值g1,g2,g3,g4,经双线性内插求得:
Figure GDA0002610363280000111
3)核线重采样
核线在摄影测量中是一个重要概念。两个相机光心的连线称为基线,通过基线的平面(即核面)与像片的交线为核线,同一核面与左右像片的交线称为同名核线。由此可知任意物点一定位于通过该点的核面与像片的交线上,影像校正可以使匹配点搜索位于同名核线上,若如图所示基线平行于行方向,则平行核线重采样后,行就是核线,通过摄影基线OO'和物点A的平面为通过像点a的核平面,图4中I代表位于左方的航摄像片,k代表相应的平行于摄影基线的水平像片,ak为A点在左水平像片上的相应像点。设a、ak在各自的像片坐标系中的坐标分别为(x,y)和(xk,yk),则有
Figure GDA0002610363280000112
其中,a1、a2、…、c3为第一幅图像的方向余弦,它们是这张像片相对于摄影基线的角方位元素的函数。
在相对水平像片上,当yk=α(α为正整数),则为核线。将yk=α代入式(7),且xk值取正整数,即求得一系列的像点坐标(x1,y1)、(x2,y2)、(x3,y3)、…,这些像点就位于倾斜像片I的核线上,将这些像点的灰度g′(x1y1),g′(x2y2),……直接赋给相对水平像片上相应的像点,就能获得相对水平像片上的核线影像。
由于在相对水平像片上,同名核线yt的坐标值相等,因此将同样的y′t=α代入第二幅图像共线方程:
Figure GDA0002610363280000121
其中,a′1,b′1,…,c′3为第二幅图像的方向余弦,分别是第二幅图像相对于摄影基线的角方位元素的函数。由此即得第二幅图像上的同名核线。
图像是规则排列的灰度格网序列,为了获取核线的灰度序列,必须对原始图像灰度进行重采样。
4)亚像素级图像匹配算法
用小基高比(小于0.1)立体测绘可以克服大基高比遮挡和辐射差异大的缺点。在小基高比条件下,两幅图像的差异性小,图像匹配的相关性变大,可以显著提高同名点配准精度。但为了获得同样的高程精度,要求同名点配准精度有一定数量级的提高(一般要求优于1/20像元以上)。因此,研究小基高比成像条件下亚像素级的图像处理匹配方法是一项关键技术。本项目提出并采用了基于解最优化的相位相关法进行图像匹配,基本原理是构建优化目标函数
Figure GDA0002610363280000122
其中,Q为含有噪声的局部互相关功率谱,
Figure GDA0002610363280000123
表示理论局部互相关功率谱,W是权重矩阵。通过求解下面的最优化问题来得到目标区域中心点的相对平移量(Δx,Δy),即
Figure GDA0002610363280000124
通过下面的步骤求解最优化问题:
a)有噪声干扰时,对相位相关矩阵奇异值分解,则
Figure GDA0002610363280000125
式中,
Figure GDA0002610363280000126
1=diag(σ12,…,σr),r=rank(Q)。用最大奇异值和相应的奇异量来估计相位相关矩阵Q,进而通过奇异值分解得到qx,qy,进而求解得到相应的相位px,py。用最小二乘法拟合px,py,得px=kxx+bx,py=kyy+by。则
Figure GDA0002610363280000131
Figure GDA0002610363280000132
b)由上面的方法通过初始Q0(wx,wy)奇异值分解得到初始迭代平移量(Δx0,Δy0),那么可以从下面的定义计算Q1(wx,wy),即
Figure GDA0002610363280000133
c)由下面的迭代计算公式计算Qi+1(wx,wy)和
Figure GDA0002610363280000134
Figure GDA0002610363280000135
判断
Figure GDA0002610363280000136
是否小于某一阈值,如果
Figure GDA0002610363280000137
小于某一阈值则迭代结束,否则返回步骤a);
d)最终的平移量由下式计算
Figure GDA0002610363280000138
全局匹配的仿真试验结果表明,全局像差估计比较精确,可以达到1/100个像素精度。但用于局部像素点平移量预估时会受到局部噪声的干扰比较大,因此在应用该方法进行亚像素级匹配时,一方面要对遥感器的噪声进行控制,另一方面还需对匹配数据进行抑噪预处理。
由上述算法进行图像匹配的结果如图5(b)所示。视差图中颜色越浅代表高度越高。图5(a)为原参考图。
上述基于解最优化的相位相关匹配算法,可以实现小基高比成像条件下的亚像素级图像匹配。然而还需要针对遥感图像进行亚像素级匹配验证,加快算法的处理速度,使其适用于航天立体测绘。
5)相对高程的计算
本节主要工作是基于上节匹配的视差图求取相对高程值。视差图上的像素点a1表征参考图上点a1与待匹配图像上对应匹配点a2的像方像素差。假设a1点视差值为s(a1),探测器像元大小为P,A点像方坐标差为
pa=a1-a2=s(a1)·P (9)
A点和C点的像方视差为
Δp=(a1-a2)-(c1-c2)=(s(a1)-s(c1))·P (10)
可计算物方相对高程
ΔZ=(s(a1-s(c1)))·GSD/(B/H) (11)
其中,
Figure GDA0002610363280000141
为物方地面采样距离。
如图6所示为将实验所得视差图进行相对像方高程计算,得到的三维重建效果图。
6)基于控制点的绝对高程计算
本节主要阐述基于控制点的相对高程值计算绝对高程的方法。假设已知K个控制点的真实高程值
Figure GDA0002610363280000142
相同参考点的高度为h,则计算绝对高程的模型为
Z=h+ΔZ (12)
通过约束K个控制点计算的绝对高程与真实高程的误差平方和最小求解参考点高度,即
Figure GDA0002610363280000143
可得检测点的绝对高程计算公式
Figure GDA0002610363280000144
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (10)

1.一种高精度小基高比立体测绘方法,其特征在于步骤如下:
1)相机在轨道高度为H,基线长度为B的两个位置对地进行拍摄,采集获取两幅被噪声污染后的退化图像,并建立所采集图像的退化模型
Figure FDA0002610363270000011
其中,g为采集到的被噪声污染后的退化图像,
Figure FDA0002610363270000012
为傅里叶逆变换算子,f为理想图像,n为噪声,H′为频域晶胞上的调制传递函数,*为卷积运算;
2)对两幅被噪声污染后的退化图像分别进行基于频域晶胞的总变分正则化MTFC预处理,得到放大后的两幅复原图像f1、f2
3)计算得到校正内方位元素畸变后的两幅图像fcal1、fcal2
4)计算得到校正外方位元素畸变后的两幅图像fext1、fext2
5)采用基于解最优化的相位相关法对校正了外方位元素畸变的两幅图像fext1、fext2进行亚像素级匹配,得到匹配后的视差图V;
6)利用视差图V计算得到整幅图各个像素点处的相对高程值ΔZ;
ΔZ=(s(a1-s(c1)))·GSD/(B/H)
其中,s(a1)、s(c1)分别为a1、c1点视差值,a1、c1为的地面任意两点A和C在第一幅图像上的成像点,P为探测器像元大小,f′为光学系统焦距,
Figure FDA0002610363270000013
为物方地面采样距离;
7)通过约束K个控制点计算的绝对高程与真实高程的误差平方和最小,求解参考面高度h,即地面点C处的绝对高程为
Figure FDA0002610363270000021
进而可计算整幅图上所有点处的绝对高程值。
2.根据权利要求1所述的一种高精度小基高比立体测绘方法,其特征在于:所述2)步骤进行基于频域晶胞的总变分正则化MTFC预处理的具体方法为:
21)采用总变分最小化模型对退化图像矩阵g进行去噪,具体公式为:
Figure FDA0002610363270000022
Figure FDA0002610363270000023
Figure FDA0002610363270000024
其中,α为正则化参数,β为可调参数,
Figure FDA0002610363270000025
Df是f的支持域,即成像系统的频域晶胞;
Figure FDA0002610363270000026
表示散度,对于理想图像f=[fx,fy]T
Figure FDA0002610363270000027
22)计算成像系统的频域晶胞,其模型为
Figure FDA0002610363270000028
其中,
Figure FDA0002610363270000029
为点(ν,ω)处的系统调制传递函数数值,
Figure FDA00026103632700000210
为系统调制传递函数矩阵,其大小为K×L,K、L均为正整数;comba(v,w)表示采样边带;θalias为阈值,其取值范围为[1,0);
23)计算频域晶胞上的调制传递函数
Figure FDA00026103632700000211
进行反卷积,得到复原后的图像。
3.根据权利要求1所述的一种高精度小基高比立体测绘方法,其特征在于:所述3)步骤计算得到校正内方位元素畸变后的两幅图像fcal1、fcal2的具体过程为:
31)计算得到相机镜头x、y两个方向上的几何畸变量Δx,Δy;
32)对两幅复原后的图像f1、f2进行校正,计算每幅图像校正后的图像坐标(xcal,ycal)对应的原坐标(x,y),即
x=xcal-abs(Δx)
y=ycal-abs(Δy)
其中,abs(·)为取绝对值操作;
33)经过双线性内插算法求得校正后坐标处(xcal,ycal)的像素值;
34)由上面的过程最终实现对复原后图像f1、f2的内方位元素校正,得到校正畸变后的两幅图像fcal1、fcal2
4.根据权利要求1所述一种高精度小基高比立体测绘方法,其特征在于:所述4)步骤计算得到校正外方位元素畸变后的两幅图像fext1、fext2的具体方法为:
41)计算获得校正外方位元素坐标值对应的原坐标值
Figure FDA0002610363270000031
其中,x、y分别为内方位元素校正后图像fcal1上像点的横、纵坐标;a1,a2,a3,b1,b2,b3,c1,c2,c3为第一幅图像的方向余弦;取xk、yk为整数,求得一系列的像点坐标(x,y);
对于图像fcal2,将y′k=yk代入下面的第一幅图像共线方程,求得对应外方位元素校正后坐标(x′k,y′k)的原坐标值(x′,y′):
Figure FDA0002610363270000041
其中,x′k为整数,x′、y′为内方位元素校正后图像fcal2上像点的横、纵坐标;a′1,b′1,…,c′3为第一幅图像的方向余弦,分别是第二幅图像相对于摄影基线的角方位元素的函数;
42)经过双线性内插算法求得校正后坐标处(xk,yk)、(x′k,y′k)的像素值;
43)将校正过内方位元素的图像fcal1、fcal2上的所有点进行41)步骤、42)步骤的核线重采样操作,得到校正了外方位元素畸变的两幅图像fext1、fext2
5.根据权利要求1所述的一种高精度小基高比立体测绘方法,其特征在于:所述5)步骤中基于解最优化的相位相关法的具体过程为:
51)通过求解下面的最优化问题来得到目标区域中心点的相对平移量(Δx,Δy),即
Figure FDA0002610363270000042
其中,Q为含有噪声的局部互相关功率谱,
Figure FDA0002610363270000043
表示理论局部互相关功率谱,W是权重矩阵;wx、wy分别为x方向和y方向的权重系数,φ(Δx,Δy)为局部互相关功率谱误差函数;
52)对相位相关矩阵进行奇异值分解,即
Figure FDA0002610363270000044
式中,
Figure FDA0002610363270000045
1=diag(σ12,…,σr),r=rank(Q);通过奇异值分解得到qx,qy,进而求解得到相应的相位px,py;由px=kxx+bx,py=kyy+by得平移量
Figure FDA0002610363270000051
53)通过初始Q0(wx,wy)奇异值分解得到初始迭代平移量(Δx0,Δy0),计算得到Q1(wx,wy),即
Figure FDA0002610363270000052
54)迭代计算得到Qi+1(wx,wy)和
Figure FDA0002610363270000053
Figure FDA0002610363270000054
判断
Figure FDA0002610363270000055
是否小于条件阈值,如果
Figure FDA0002610363270000056
小于某一阈值则迭代结束,得到不同迭代次数对应的平移量(Δxi,Δyi),若不收敛,调整阈值再进行54)步骤迭代计算;
55)计算得到最终的平移量,即x、y方向的视差值为
Figure FDA0002610363270000057
6.根据权利要求3所述的一种高精度小基高比立体测绘方法,其特征在于:所述31)步骤计算得到相机镜头x、y两个方向上的几何畸变量Δx,Δy的具体方法为:
Figure FDA0002610363270000058
Figure FDA0002610363270000059
其中,x、y为退化图像上像点的横、纵坐标;X、Y、Z为地面坐标系,a1,a2,a3,b1,b2,b3,c1,c2,c3为通过相机的姿态构建的旋转矩阵,x0、y0为x、y两个方向上的主点偏移,XS,YS,ZS为相机的轨道参数;Δx、Δy为x、y两个方向上的几何畸变量。
7.根据权利要求3所述的一种高精度小基高比立体测绘方法,其特征在于:所述33)步骤经过双线性内插算法求得校正后坐标处(xcal,ycal)的像素值的具体方法为:
Figure FDA0002610363270000061
其中,Δd为原始变形图像的采样间距,g′为原坐标(x,y)处待求像素灰度值,g1,g2,g3,g4分别为其周围4个像元的灰度值。
8.根据权利要求4所述的一种高精度小基高比立体测绘方法,其特征在于:所述42)步骤经过双线性内插算法求得校正后坐标处(xk,yk)、(x′k,y′k)的像素值的具体方法为:
Figure FDA0002610363270000062
其中,Δd′为原始变形图像的采样间距,g″为原坐标(x′,y′)处待求像素灰度值,g′1,g′2,g′3,g′4分别为原坐标(x′,y′)周围4个像元的灰度值;同理求出第二幅图像坐标处(x′k,y′k)的像素值。
9.根据权利要求1-8任一所述的一种高精度小基高比立体测绘方法,其特征在于:所述1)步骤中相机在轨道高度为H,基线长度为B的两个位置对地进行拍摄,两个位置关于光轴方向对称。
10.根据权利要求1-8任一所述的一种高精度小基高比立体测绘方法,其特征在于:B/H≤0.12。
CN201810589576.4A 2018-06-08 2018-06-08 一种高精度小基高比立体测绘方法 Active CN109029379B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810589576.4A CN109029379B (zh) 2018-06-08 2018-06-08 一种高精度小基高比立体测绘方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810589576.4A CN109029379B (zh) 2018-06-08 2018-06-08 一种高精度小基高比立体测绘方法

Publications (2)

Publication Number Publication Date
CN109029379A CN109029379A (zh) 2018-12-18
CN109029379B true CN109029379B (zh) 2020-10-20

Family

ID=64612560

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810589576.4A Active CN109029379B (zh) 2018-06-08 2018-06-08 一种高精度小基高比立体测绘方法

Country Status (1)

Country Link
CN (1) CN109029379B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109711486B (zh) * 2019-01-21 2020-09-01 湖北省国土资源研究院 基于相位相关的高重叠度遥感影像满度连接点匹配方法
CN110307858B (zh) * 2019-06-26 2021-09-03 中国空间技术研究院 一种测绘卫星交汇误差的自适应校正方法
CN113469896B (zh) * 2021-05-21 2022-04-15 贵州师范学院 一种提高地球同步轨道卫星对地观测影像几何校正精度的方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5463397A (en) * 1993-10-25 1995-10-31 Hughes Aircraft Company Hyper-precision SAR interferometry using a dual-antenna multi-pass SAR system
KR100912715B1 (ko) * 2007-12-17 2009-08-19 한국전자통신연구원 이종 센서 통합 모델링에 의한 수치 사진 측량 방법 및장치
CN101718550B (zh) * 2009-12-18 2011-04-27 北京空间机电研究所 小基高比立体测绘光学系统
CN102184540B (zh) * 2011-05-03 2013-03-20 哈尔滨工程大学 基于尺度空间的亚像素级立体匹配方法
CN103630120B (zh) * 2013-07-16 2016-09-14 中国人民解放军信息工程大学 基于严密几何模型的火星表面线阵影像核线重采样方法
CN107144293A (zh) * 2017-04-07 2017-09-08 武汉大学 一种视频卫星面阵相机的几何定标方法

Also Published As

Publication number Publication date
CN109029379A (zh) 2018-12-18

Similar Documents

Publication Publication Date Title
CN109903352B (zh) 一种卫星遥感影像大区域无缝正射影像制作方法
De Franchis et al. An automatic and modular stereo pipeline for pushbroom images
Li et al. Rigorous photogrammetric processing of HiRISE stereo imagery for Mars topographic mapping
CN104537707B (zh) 像方型立体视觉在线移动实时测量系统
CN110500995A (zh) 利用rpc参数建立高分辨率卫星影像等效几何成像模型的方法
Li et al. A new analytical method for estimating Antarctic ice flow in the 1960s from historical optical satellite imagery
CN109696182A (zh) 一种星载推扫式光学传感器内方位元素定标方法
CN109029379B (zh) 一种高精度小基高比立体测绘方法
CN113514829B (zh) 面向InSAR的初始DSM的区域网平差方法
CN111044037B (zh) 一种光学卫星影像的几何定位方法及装置
CN103697864B (zh) 一种基于大虚拟相机的窄视场双相机影像拼接方法
Wu et al. Co-registration of lunar topographic models derived from Chang’E-1, SELENE, and LRO laser altimeter data based on a novel surface matchingmethod
CN106885585A (zh) 一种基于光束法平差的星载摄影测量系统一体化检校方法
CN110853140A (zh) 一种dem辅助的光学视频卫星影像稳像方法
CN111156969A (zh) 一种宽幅遥感影像立体测绘方法及系统
Pi et al. On-orbit geometric calibration using a cross-image pair for the linear sensor aboard the agile optical satellite
Chen et al. Large-scale block bundle adjustment of LROC NAC images for Lunar South Pole mapping based on topographic constraint
Li et al. Photogrammetric processing of Tianwen-1 HiRIC imagery for precision topographic mapping on Mars
Li et al. The generalized-line-based iterative transformation model for imagery registration and rectification
CN117830543A (zh) 基于星载双站InSAR和激光雷达数据估计DEM的方法、装置、设备及介质
CN110503604B (zh) 一种基于高精度pos的航空面阵影像实时正射拼接方法
Hu et al. Planetary3D: A photogrammetric tool for 3D topographic mapping of planetary bodies
CN113920046B (zh) 一种多分片卫星影像拼接和几何模型构建方法
Ma Rational function model in processing historical aerial photographs
Awange et al. Fundamentals of photogrammetry

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
GR01 Patent grant
GR01 Patent grant