一种高精度小基高比立体测绘方法与流程

本发明属于航天光学遥感技术领域，涉及一种高精度小基高比立体测绘方法。

背景技术

小基高比立体测绘方法与传统大基高比立体测绘方法完全不同。大基高比是通过图像匹配在两幅图像中找到同名点，进而将同名控制点坐标代入共线方程求得共线方程参数，将待求点像方坐标代入已求解好的共线方程中，求得待求点的物方三维坐标。而小基高比立体测绘方法的核心计算公式为像方相对高程等于两幅图像的平面视差除以基高比数值。早在2002年，法国人就利用spot5卫星上的全色影像和多光谱影像组成立体像对，对小基高比情况下获取数字高程模型进行了试验，结果表明在立体交会角仅为0.02时，仍能获取一定精度的高程信息，从而验证了小基高比立体测绘技术的可行性。

传统的大基高比立体测绘方法是逐点代入共线方程，逐点获取相应点处的高程值；单台相机的大基高比立体测绘对平台的姿控提出了更高的要求，要卫星平台灵活、敏捷度高、稳定性强、姿态机动能力强，实现难度更大；大基高比条件下，拍摄城市区域(高大建筑物密集)时，丢失高程信息的遮挡区域将增大，不利于城市立体测绘。

2007年第28卷《journalofmathematicalimagingandvision》上，法国空间研究中心的juliedelon等人发表的《smallbaselinestereovision》论文通过理论推导和仿真分析，讨论了小基高比立体测绘方法的可行性及相关仿真实验验证，其不足之处在于：该方法仅通过两幅图像的平面视差除以基高比数值求得像方相对高程，没有对相机的内外方位元素进行校正，没有对原始采集到的图像进行合理的mtf补偿，不能计算出物方绝对高程。

技术实现要素：

本发明解决的技术问题是：克服现有技术的不足，提供了一种高精度小基高比立体测绘方法，该方法提出了必须对采集到的图像进行mtf补偿预处理，并校正相机的内外方位元素，最后通过控制点计算参考面的绝对高程来计算整个重叠区域的物方绝对高程值，实现最终的地物高精度立体测绘。

本发明的技术方案是：一种高精度小基高比立体测绘方法，步骤如下：

1)相机在轨道高度为h，基线长度为b的两个位置对地进行拍摄，采集获取两幅被噪声污染后的退化图像，并建立所采集图像的退化模型其中，g为采集到的被噪声污染后的退化图像，为傅里叶逆变换算子，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′为光学系统焦距，为物方地面采样距离；

7)通过约束k个控制点计算的绝对高程与真实高程的误差平方和最小，求解参考面高度h，即地面点c处的绝对高程为

进而可计算整幅图上所有点处的绝对高程值。

所述步骤2)进行基于频域晶胞的总变分正则化mtfc预处理的具体方法为：

21)采用总变分最小化模型对退化图像矩阵g进行去噪，具体公式为：

其中，α为正则化参数，β为可调参数，df是f的支持域，即成像系统的频域晶胞；表示散度，对于二维向量f＝[fx,fy]^t，

22)计算成像系统的频域晶胞，其模型为

其中，为点(ν,ω)处的系统调制传递函数数值，为系统调制传递函数矩阵，其大小为k×l，k、l均为正整数；comba(v,w)表示采样边带；θalias为阈值，其取值范围为[1,0)；

23)计算频域晶胞上的调制传递函数进行反卷积，得到复原后的图像f。

所述步骤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)计算获得校正外方位元素坐标值对应的原坐标值

其中，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′)：

其中，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)，即

其中，q为含有噪声的局部互相关功率谱，表示理论局部互相关功率谱，w是权重矩阵；wx、wy分别为x方向和y方向的权重系数，φ(δx,δy)局部互相关功率谱误差函数；

52)对相位相关矩阵进行奇异值分解，即

式中，∑1＝diag(σ1,σ2,…,σr)，r＝rank(q)；通过奇异值分解得到qx,qy，进而求解得到相应的相位px,py；由px＝kxx+bx,py＝kyy+by得平移量

53)通过初始q⁰(wx,wy)奇异值分解得到初始迭代平移量(δx⁰,δy⁰)，计算得到q¹(wx,wy)，即

54)迭代计算得到qⁱ⁺¹(wx,wy)和i＝1，2…………

判断是否小于条件阈值，如果小于某一阈值则迭代结束，得到不同迭代次数对应的平移量(δxⁱ,δyⁱ)，若不收敛，调整阈值再进行步骤54)迭代计算；

55)计算得到最终的平移量，即x、y方向的视差值为

所述步骤31)计算得到相机镜头x、y两个方向上的几何畸变量δx,δy的具体方法为：

其中，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)的像素值的具体方法为：

其中，δd为原始变形图像的采样间距，g′为原坐标(x,y)处待求像素灰度值，g1,g2,g3,g4分别为其周围4个像元的灰度值。

所述步骤42)经过双线性内插算法求得校正后坐标处(xk,yk)、(x′k,y′k)的像素值的具体方法为：

其中，δ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算法。首先，整个采样系统的模型可以表示为

其中，g为采集到的被噪声污染后的退化图像矩阵，为傅里叶逆变换算子，f为理想图像矩阵，n为噪声矩阵，h′为降晰函数傅里叶变换的模，*为卷积运算。基于频域晶胞模型的总变分正则化mtfc算法的技术途径如下：

a)采用总变分最小化模型进行去噪，即：

其中，β>0是可调参数。利用该总变分最小化模型进行去噪，该模型迭代计算公式为

b)计算成像系统的频域晶胞，其模型为

其中，为系统mtf，comba(v,w)表示采样边带，θalias为阈值，其取值范围为[1,0)。

c)小基高比立体测绘中需要利用频域晶胞进行h′的计算，由步骤b)中频域晶胞的计算公式可以得到与相应成像系统mtf最为匹配的频域晶胞，利用频域晶胞上的mtf进行图像复原，则有

基于频域晶胞的图像复原方法与传统方法的最大不同在于式(1)中的h′由式(3)计算求得。

由步骤a)的迭代计算公式可得到的值，将上式计算的h′代人，进行反卷积计算，求解复原图像f1、f1。

物理验证样机的试验数据在mtfc前最优高程精度为5.66个gsd，利用基于频域晶胞的总变分正则化mtfc方法的最优高程精度为4.76个gsd，基于频域晶胞的mtfc方法处理后的高程精度，相比于mtfc前提高了16.1％。

2)内方位元素校正

在小基高比立体测绘中，内方位元素校正主要是对由相机自身特性引起的几何变形进行校正。考虑改正项后的共线方程如式：

其中，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，经双线性内插求得：

3)核线重采样

核线在摄影测量中是一个重要概念。两个相机光心的连线称为基线，通过基线的平面(即核面)与像片的交线为核线，同一核面与左右像片的交线称为同名核线。由此可知任意物点一定位于通过该点的核面与像片的交线上，影像校正可以使匹配点搜索位于同名核线上，若如图所示基线平行于行方向，则平行核线重采样后，行就是核线，通过摄影基线oo'和物点a的平面为通过像点a的核平面，图4中i代表位于左方的航摄像片，k代表相应的平行于摄影基线的水平像片，ak为a点在左水平像片上的相应像点。设a、ak在各自的像片坐标系中的坐标分别为(x,y)和(xk,yk)，则有

其中，a1、a2、…、c3为第一幅图像的方向余弦，它们是这张像片相对于摄影基线的角方位元素的函数。

在相对水平像片上，当yk＝α(α为正整数)，则为核线。将yk＝α代入式(7)，且xk值取正整数，即求得一系列的像点坐标(x1,y1)、(x2,y2)、(x3,y3)、…，这些像点就位于倾斜像片i的核线上，将这些像点的灰度g′(x1y1)，g′(x2y2)，……直接赋给相对水平像片上相应的像点，就能获得相对水平像片上的核线影像。

由于在相对水平像片上，同名核线yt的坐标值相等，因此将同样的y′t＝α代入第二幅图像共线方程：

其中，a′1，b′1，…，c′3为第二幅图像的方向余弦，分别是第二幅图像相对于摄影基线的角方位元素的函数。由此即得第二幅图像上的同名核线。

图像是规则排列的灰度格网序列，为了获取核线的灰度序列，必须对原始图像灰度进行重采样。

4)亚像素级图像匹配算法

用小基高比(小于0.1)立体测绘可以克服大基高比遮挡和辐射差异大的缺点。在小基高比条件下，两幅图像的差异性小，图像匹配的相关性变大，可以显著提高同名点配准精度。但为了获得同样的高程精度，要求同名点配准精度有一定数量级的提高(一般要求优于1/20像元以上)。因此，研究小基高比成像条件下亚像素级的图像处理匹配方法是一项关键技术。本项目提出并采用了基于解最优化的相位相关法进行图像匹配，基本原理是构建优化目标函数

其中，q为含有噪声的局部互相关功率谱，表示理论局部互相关功率谱，w是权重矩阵。通过求解下面的最优化问题来得到目标区域中心点的相对平移量(δx,δy)，即

通过下面的步骤求解最优化问题：

a)有噪声干扰时，对相位相关矩阵奇异值分解，则

式中，∑1＝diag(σ1,σ2,…,σr)，r＝rank(q)。用最大奇异值和相应的奇异量来估计相位相关矩阵q，进而通过奇异值分解得到qx,qy，进而求解得到相应的相位px,py。用最小二乘法拟合px,py，得px＝kxx+bx,py＝kyy+by。则

b)由上面的方法通过初始q⁰(wx,wy)奇异值分解得到初始迭代平移量(δx⁰,δy⁰)，那么可以从下面的定义计算q¹(wx,wy)，即

c)由下面的迭代计算公式计算qⁱ⁺¹(wx,wy)和

判断是否小于某一阈值，如果小于某一阈值则迭代结束，否则返回步骤a)；

d)最终的平移量由下式计算

全局匹配的仿真试验结果表明，全局像差估计比较精确，可以达到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)

其中，为物方地面采样距离。

如图6所示为将实验所得视差图进行相对像方高程计算，得到的三维重建效果图。

6)基于控制点的绝对高程计算

本节主要阐述基于控制点的相对高程值计算绝对高程的方法。假设已知k个控制点的真实高程值相同参考点的高度为h，则计算绝对高程的模型为

z＝h+δz(12)

通过约束k个控制点计算的绝对高程与真实高程的误差平方和最小求解参考点高度，即

可得检测点的绝对高程计算公式

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。