一种基于二维图像信息的三维地形测量
翻译:杜雷
班级：测绘一班
学号：[1**********]0
【摘要】 研究目的:利用数字图像测量技术对河流模型实验中的河床地形测量研究。创新要点:以高质量的图像径向畸变校正为基础,依据多幅图像间映射换算以及单幅图像灰度信息变换两方面的研究,来提取三维地形信息。研究方法:使用考虑残差的一阶误差传播的最大似然估计方法求解图像与地形坐标的映射变换参数;借助漫反射光照模型,通过适当假设,推导出地形高程与图像灰度、灰度梯度及太阳方向的关系式。重要结论:1.以Devernay基于直线的图像畸变校正方法为基础,通过使用更高精度的Canny边缘检测算子提取图像边缘,黄金分割法和二次插值法相结合的最优化方法求取畸变系数等,可以得到较好的校正效果,求算的畸变系数精度更高,求算过程更快;2.利用图像与地面间的透视投影关系方法求解地形高程坐标,均匀并尽可能多地将控制点布置在控制区域内,使计算值更接近实测值,并满足河流演变实验的精度要求;3.通过推导关系式,理论上可以根据图像灰度、灰度梯度、太阳方向向量算得地形高程。
【关键词】 河床地形测量；数字图像测量技术；径向畸变矫正；投影变换；灰度信息变换； 1 介绍
河床地形测量是一个在河工模型试验中最复杂的任务。传统的地形测量工具，如电子水平仪，电子中转仪，以及全站仪，很难在河床地形测量之中应用。目前，电阻式测距仪是通常使用的装置。它具有良好的适应性，可以测量多种地形。它适用于不同的含沙量和沙质条件。但是，它的低测量的效率和河床干扰问题是它的主要缺点。
为了满足动态测量的需求，各种各样的地形测量方法已经进化。地面激光扫描，空中光探测和测距（LIDAR），多波束声纳，实时动态（RTK），全球定位系统，和全站仪测量不断被开发以满足增加地貌科学研究的兴趣和工程应用。然而，这样的方法更适合于大规模河床模型，甚至是实际的河流。此外，还有很多软件及硬件资源。在最近几年，随着允许的精确校准方法的发展，非量测相机变得越来
越可靠，摄影的自动化变得可访问到广泛的用户群。使用在线结构从运动（SFM）计划等。2013年Bouratsis等。利用自由和开源软件的优势从照片制作中创建的河流环境模型的高分辨率数字高程；2013有人使用的一对商业摄像机来记录河床的演化，并且由一组计算机视觉的一个计算方法和图像处理算法被用来分析视频和重建的瞬时3D床的表面上。詹姆斯和罗布森（2012）集成SFM和多视点立体（MVS）算法研究在方法测量地形需要很少的专业知识和实现自动化处理。阿斯特吕克等在2012年研究了立体技术测量沙床抬高斜在区波时间尺度。此方法是非侵入性，所以有很高的精度，时间分辨率也很容易让人抓获。此外，鲁等人在2008年提出了一些有关平原和高原地形的专题制图相关的研究
图象数据是简单地2D数据。然而，许多线索可能在单个图像或多个图像中找到，其中可以通过图像场景进行三维重建。最典型的重建方法包括从校正方法的形状（Zhang等，1999年），从质地法的形状（福赛思和庞塞，2002），手动交互方法（Shashua，1997）。也有根据一些方法在多个图像进行重建，例如立体视觉的方法，运动图像序列的方法，和光度立体法（李，1991;张Y.J.，2000;Pollefeys和GOOL，2002）。还值得一提的是，三维重建应用广泛的还有地形测量。还有在其他领域的一些显著的应用。 Gómez等人，（2013年）进行了新的方法用于重建，从多个3D + T彩色多普勒3D + T速度场图像。 Mikes等，（2013年）集中在扫描光束偏折表面重建，并提出和分析的新的数学方法从而解决了三维表面形状重建问题。柳（2013）提出了一种新的深度提取方法技术与计算积分成像小透镜阵列。此外，图像分析技术也已用于从CCD空间和时间上检索水面仰角（Benetazzo，2006; Gallego等，2011）。在一定程度上，这些技术是类似于地形测量。他们用重叠的视频图像恢复地形信息。
2 摄像机失真校正
图像的非线性失真校正是图像测量及三维重建必不可少的步骤。它是CCD照相机的几何模型用于确定三维几何的关系对象表面的信息和相应的图像匹配点。由于光学透镜的结构和理想的线性与小孔成像模型之间的不同，从CCD拍摄的图像会产生非线性几何失真。这严重影响了3D信息的准确性通过2D图像获得
的。
2.1数学模型
在设计和制造透镜的过程中不可避免地产生在图像畸变。径向变形，偏心变形和薄：一般效果可以通过数学模型来表示棱镜失真。镜头畸变通常可表示为：
其中，xd和yd是原始失真图像坐标，xu和yu是理想的失真图像坐标。xd和yd是在x和y方向上的失真，分别描述为：
其中第一项是径向失真，所述第二术语是偏心失真，而第三项是薄棱镜变形; K1，K2，P1，P2，S1和S2为非线性失真参数（翁等人，1992;Ahmed和法拉格，2005; Wang等人，2008）。以前的研究（张Z.Y.，2000;江等人，2001;个人和Kovacic，2002; Miks和诺瓦克，2012）指出，透镜径向失真是多比其他像差更加严重。切向失真和畸变棱镜不必考虑。太多的非线性参数不会提高解的精度，而是将导致其不稳定。
即最常用的描述径向失真多项式模型（PM）可以被写入为：
其中，rd和ru是从各距离扭曲的点（xd，yd）和失真点（xu，yu），以扭曲的中心P;ki（ⅰ= 1，2，...）是径向失真参数。 Devernay和Faugeras（1995）表明所述第一阶径向对称畸变参数，K1，可以实现一定的准确性。 为了实现更高的精度，我们使用了一阶和二阶畸变参数，k1和k2，测量失真的图像失真程度。由于图像中心是失真（Tordoff和Murray，2000）中心的一个很好的近似，
我们随便拿图像作为失真中心P的中心（xc，yc）。图像校正公式坐标表示为：
然后，改写后的未失真坐标和扭曲的坐标由下式给出：
2.2校准过程
为了找到畸变参数，我们使用Devernay和Faugeras（1995年）的直线法和一些新的改进。这种方法使用基本属性（Devernay和Faugeras，2001年）是在空间中的每个直线上的投影是一条直线。如果我们可以发现径向失真图像上的变换，这样在空间中的每个直线被认为是
转换后的图像中的直线，那么我们就可以估计的图像的畸变的参数。 通过使用此属性，采用一个反复的过程来估计失真参数k1和k2。
2.2.1边缘检测
图像边缘具有两个属性：方向和大小。边缘像素沿着边缘方向平缓变化;然而在垂直于边缘的方向存在一个巨大的变化。首先，用一个新的算术Canny边缘检测器（穆罕默德等人，2013）被用来获得的幅度和边缘的方向，然后将数据由高斯低通滤波器进行处理。其次，非极大值抑制和滞后阈值被用于的精确位置具有子像素精度的优势。这种方法可以提高检测的边缘的分辨率，保护低强度边缘，增强抗干扰能力。
2.2.2提取物扭曲线段
边缘检测之后，我们需要提取扭曲线段，这是在三维空间中最有可能的直线。因为某些部分可以被打破边缘检测，当两端之间的距离小于一个0.5像素时，我们加入碎段。我们也可以把最低门槛上段长度定为2个像素，理由是较短的片段可以包含比关于失真的有用信息更多的噪音。
2.2.3测量失真错误
为了得到失真参数，计算测量一个扭曲线段的曲率，通过它多少扭曲来分割我们使用的扭曲线段，通过最小二乘逼近以形成一条直线。在一般情况下，直线的方程式接合端点（x1，y2）和（x2，y2），表示为：
为一个线段的任何点（x，y）的，我们使用作为距离，从点到直线，然后畸变误差是所有的距离的平方的总和。其结果是，该畸变误差为零，当失真段是一个直线。扭曲段的曲率越大，畸变误差将是较大的。
2.2.4优化方法
我们结合黄金分割搜索方法GSSM）和二次插值法（QIM）（龚和王，2009年）来计算最佳失真参数。 QIM是准确的，但不非常有效的。所以，我们首先得到一个小的GSSM值范围，然后使用QIM获取，满足我们的精度要求。
整个优化过程是：首先，我们使用该范围的金色部分作为初始点，然后执行由GSSM迭代，直到有一点，从获得的函数f（x）值的差相邻的两步骤迭代。然后一个由迭代QIM完成。如果没有太大的区别，那么从两种方法求出的f（x）的值，另一迭代由QIM完成。再次，如果没有求出的f（x）的值相差无几从由QIM相邻两步骤迭代，并且自变量之间的差别非常小，迭代继续进行，直至精度要求得到满足。否则，我们改为GSSM。我们设定精度为11010。
由于参数数据减少，并且包含许多异常值，有一个从所述最优参数进一步移
动的危险。因为这个原因，我们在K1第一只优化，直到它给出了一个稳定的解决方案。然后，K2加入，并最终全面优化失真参数下进行。
2.3标定实验结果
要验证校准效果，我们进行两个图像校准实验。该原始图像坐标是不可或缺的。 然而，校准后，图像的坐标非整数，且需要估计的整数坐标的新的灰度值。考虑这个问题，我们使用的三次样条插值法灰度级内插。
图像边缘轮廓如图1所示，该两个图像都产生桶形失真。通过Devernay和Faugeras实验结果（1995）的方法示于图2，计算时间是65.6秒和32.9秒，失
777.9415108.691410真参数是K1=和K1=，残余误差分别为24和1.8。该通
过该方法实验结果如图3，计算时间是55.3秒和30.2秒，失真参数K1=6.9464×10-7，K2=4.6261×10-12和K1=7.9267×10-7，K2=6.2357×10-13，残余误差分别为1和0.8。结果看起来很合理，表明Devernay ，Faugeras（1995年）提出的方法是有效的。我们也提高了精度和失真参数的计算效率，因此，我们可以使用校准图像信息做一个更精确的图像调查。
图1 两个用于校正失真图像图
（a）从我们的CCD照相机;（二）从http://www.21tx.com/
DC/ 2007/11/05 /11645.html（五月，2013访问）
图2通过Devernay和Faugeras（1995年）的原始方法使用参数校准来计算图像
（a）：经校准图1（a）中的结果;
（b）：校准图1（b）中的结果。
图3 通过该方法使用参数校准来计算图像
（a）:经校准图1(a)中的结果;
（b）：校准图（b）中结果。
3 2D图像间的映射转换和三维空间坐标
3.1透视变换
考虑到校准的图像可以看作地面中央突起，则存在图像与地面之间的透视投影变换的关系。在摄影测量中，常用的数学模型包括直接线性变换（DLT）（Chen等，1994）和空间后方交会的方法（王，2007）。空间后方交会的方法需要较少的控制点，但它具有较高的要求,在实际迭代过程的初始值要有较高的准确度。所以我们使用DLT来计算转换参数。
DLT是一个直接描述关系角度图像和地面之间的解析表达式。它包含11个未知参数（L1-L11），并能提高图像的线性错误。让（U，V）是图像坐标，（X，Y，Z）是相应的三维空间坐标，然后直接线性关系可以表示为：
当给定的n（n≥6）控制点，我们通常通过最小二乘法来计算这些11个参数。传统的方法是：
当
残留方格的最小总和是在传统方法的目标函数。然而，坐标图像之间的统计属性和相应的三维空间中的坐标被致命地忽略。该溶液的精度不足够大，而且它只能用于低精度摄影。在这项研究中，我们使用一个近似最大似然估计法（AMLE），并考虑残差的一阶误差传播（周和邓，2011年）来计算转换参数。另外，更简单的迭代算法被应用。
3.2变换参数解决方案基于角法
由于每个点带来了两个方程，用i来代表点数，让：
那么我们的控制方程：
是观测值。假定观测是彼此独立的
，
和
据误差传播法，A1的方差，i可以表示为：
是观测协方差矩阵。根
其中，带S是由相机拍摄规模确定的等效误差的因素。因为在摄影的空间坐标误差以毫米和图像测量坐标通常表示是像素，所需要的像素错误转化为毫米。
这是一个半线性错误 - 在变量（EIV）模型参数估计问题（Muhlich 和 Meste，2001年）
。在传统的问题可以表示为
即，忽略剩余的权重，是单位矩阵，该AMLE问题是一个普遍的最小二乘问题，并要确定的。考虑各种方法相对复杂，估计的参数是对应于矩阵的最小特征值的特征向量
我们使用一种迭代加权总最小二乘法（WTLS）。其精度可以比传统的方法大10％，而且，将有一个显著影响时的控制点的观测精度相差很大。
改进后的表达式如下：
其中：
问题的参数的解决方案是对应于矩阵S的下值分解的最小特征值，而S被参数确定。具体迭代过程描述如下：
1.让解决方案通过传统的方法来计算初始值;
2.使用初始值根据公式（15）计算S；
3.执行在S奇异值分解;该解决方案是对应于最小特征值的特征向量;
4.与欧几里得范数进行比较的初始值的的解决方案中，如果差小于1105，则计算停止阈值;如果不继续做迭代。
用这种方法我们可以得到所有公式（10）的11个参数。我们希望有的点的三维坐标，可以从不同的角度同一个地方的照片与两个不同的CCD摄像机，然后使用该方法来计算两个图像的11个参数，从而获得通过耦合方程的三维坐标。如显示在下面的原始和派生方程：
其中，（ui，vi）是从一个CCD摄像头，Lj目标点的图像坐标（J =1，2，...，
11）是第i个图象的解决透视变换参数。
从式（16）明显的看出，根据图像坐标的对应点，可以容易地获得的三维坐标。此外，我们可以从三个或更多个摄像机获得较高的测量精度。如果我们从两个或多个位置捕捉显著图像，并至少有六个非共面的控制点，一个单一的相机也可以具有类似的效果。
3.3实验结果
为了验证结果，两个模型实验如图4。他们在不同的位置，由Microvision公司（MV-1394）提供的一台相机捕获。该图像是640×400像素。河流模型中的河流是由河床地形演变冲积蜿蜒的河流实验研究而建造的。该实验是在小规模的水盆中进行的尺寸为4米×1.5米×0.3米（长×宽×深）。盆充满沙子上至厚度为10厘米。这些图像显示了稳定的弧形范围由流形。这个保护范围面积约2平方米，图像已与失真校正和歪斜校正完成。控制点是图像的热点。
图4 在不同位置拍摄的两幅图像
（a）从角度1，（b）从角度2
图5示出所测量的地形轮廓（小时）。一种光纤液位测量仪（IH型）在测量中被使用。通过该方法计算出的地形轮廓示于图6.双线性插值算法所示。 两幅图像在不同位置从1角捕获;从2角以分散的数据插值并应用。从两个图像，我们可以看到，这两个图像的地形变化趋势是基本相同。它们几乎完全匹配于所述控制区域，只有该区域以外的一些错误。
我们从控制区域中选择了10个控制点，即（xi，yi）（i = 1,2，...，10）。让hi第i个点插值测量高程值，h1i是用10个控制点的第i个点计算的高程值，h2i可以使用八个控制点的第i个点计算的高程值，并h3i是使用六个控制点的第i个点计算的高程值。表1表示的是对于仰角从不同的方向的坐标混淆矩阵。
图5 测量地形等高线
图6 计算地形等高线
表1 海拔混淆矩阵从不同角度的数据处理方法的坐标
从表1，我们可以看到，随着更多的控制点，计算值和实验值更接近。测量和计算的数据之间，即偏置结果，（J = 1，2，3），分别是3.1101，43.8393和47.8101。考虑到hi来自于测量数据插值，并且差别不大，从精确的实际值，
然而，对于三次的最大差异是1.049，4.422，h1i能满足河工模型试验的精度要求。
和3.316，h2i和h3i明显是不实际的。在实际的实验中，许多因素可能影响计算，精度的透视变换参数，数目和控制点的位置，目标点中心坐标的精度，并依此类推。周和邓（2011）比较了经典最小二乘法和通过实验的AMLE方法叠加在图像坐标高斯噪声，并表明，如果在计算的点具有不相等的精度，由AMLE驱动的结果会更好。在这种情况下，我们可以减少在图6更所示的噪声。 此外，它也有必要进行进一步的优化和改进，并建立改进的错误校正和补偿的机制。
4 2D图像灰度间的转换值和三维空间坐标
4.1 理论假设和数学模型
灰度特性是2D图像最重要特征之一。从CCD图像传感器获得的图像灰度值受到各种条件的影响。主要影响因素包括照明，大气环境（取决于大气特性：大气衰减，悬浮散射大气中的颗粒等），物体表面（取决于物体温度：对象表面表征，观察方向，发生率等的角度），光学成像系统。
假设理想的光学透镜，一个理想的光学成像系统，无辐射的场景区域源，以及该区域的照明与从大气中的相同，然后该照明的图像的灰度值仅与每个点的位置和气候条件有关。
灰度值与照明相应对象的表面是成比例的。假设没有外界的影响，在给定时刻的理想灰值和所述地形高程值具有正相关，描述如下：
其中，IP0是点P的理想的灰度值，HP是对应海拔坐标，f是正相关函数。 通常利用照明模型研究来对象的表面反射率的效果。现有的光照模型包括整体和局部光照模型。整体光照模型使用光线追踪技术，并有强烈的真实感，但是需要大量的计算。当地的光照模型主要包括漫反射朗伯模型和进一步加镜面反射的Phong模型。在这项研究中，我们使用的漫反射光照模型。设Sp是扩散反射灰色系数
其中Cp在P点的一个特定的反射系数光波长，a的入射角，Ke是漫反射衰减系数。
然后实际的图像灰度值可以表示成如下的公式：
其中，k是在考虑气候因素和室内因素情况下的衰减系数，将式（17）和式
（18）代入式（19），我们得到：
其中是可以根据归一化的衍生物体表面法线的内积[-p，-q，1]和光源（太阳）方向矢量[-ps，-qs，1]：
其中，分别代表在梯度的x和y方向。在这项研究中，我们使用h = h（x，y）来表示物体表面方程。
4.2太阳方向矢量计算模型
太阳方向矢量由太阳位置确定：
其中是太阳方位角，是太阳仰角。
太阳赤纬太阳高度角的变化与本地时间有关。我们使用代表太阳赤纬，并使用代表观察位置的纬度，我们规定北纬为正，且南纬为负。太阳时角由Ω表示，然后太阳仰角的计算公式为：
太阳方位角可以有以下计算：
对于上述等式中，我们可以参考“的地面气象观测规范“为太阳赤纬和太阳小时，但角度以避免耗时的操作，我们用于实时溶液方法计算它们。
太阳时角：
其中TT是实际太阳时，
其中S和M分别是各自小时值和观察时间最小值，LC是经度修改后的值，而且EQ时间等式。
其中LC是观察东经度值，
其中：
其中t由两部分组成：
其中，N是与世界时间同步的，
其中，Y是年份，floor是INT类型。
太阳赤纬是：
因此，如果我们了解经度，纬度和图像采集时间，我们可以计算出太阳 偏角和太阳时角根据公式（25）和（31），然后根据方程（22）-（24）推导出太阳的方向向量。
4.3公式改造和推导
在具体的实验中，太阳方向随着时间的推移矢量的变化。如果我们让图像采集时间是T，那么太阳的方向向量会[-ps（T），-qs（T），1]。对于任何图像点（x，y）时，根据方程（20）和（21），我们有：
我们做如下的改造：
然后我们将上边的代入公式（32），然后我们会有：
其中
可以表示为：
是I到h的偏导，和Gx和Gy是x和y方向的灰度梯度，分别
假设f是线性方程，表示为：，根据公式（33），我们可以得到：
使，并且，然后通过化简我们可以得到：
解方程（36），我们可以得到四阶方程的解：
其中，
根据式（33）和式（37），我们可以得到：
其中，r是指示特征量灰度梯度和太阳之间的关系的方向向量。式（38）解释了海拔地形坐标，灰度值，灰梯度，和太阳方向矢量之间的关系，这具有十分
重要的意义，如果我们可以得到对式（38）的精确分析和表达。
4.4 初步拟合实验
我们对初步实验进行了上面的推导，样品图像如图7（相同的弯曲触及如图
4）所示：
图7 由CCD相机拍摄的实验样品图像
根据谷歌地球，我们得到了实验地点坐标的经度和纬度：39°06′37″N,117°10′03″ E。根据图像采集的时间在2013年6月1号下午4:00
，我们得到：
，。在大量初步简化（合理或不合理的可能）后，我们得到如下相符合的表述通过使用部分测量数据：
通过拟合公式（39）计算出的地形等高线如图8所示，与图5和图6相比，我们只能观察地形变化的基本趋势，并有相当大的数值误差。很显然，初步拟合方程（39）是远远不够的。这种试探性的拟合实验不能证明非常成功。在后来的研究中，我们作出广泛的尝试以获得式（38）更精确的表达。不幸的是，我们还没有得到一个非常满意的结果。主要的原因是：
（1），我们使用其中仅考虑漫反射的局部光照模型，但没有考虑隐面消除;
（2）我们假定在推导过程的未知的正相关性函数f是线性的，可能不是这种情况;
（3）对应的三维空间中的点的通过投影变换获得的图像的灰度值和灰度梯度是错误的;
（4）解算的太阳方位角和太阳高度角是错误的;
（5）有实际图像，其影响是难以消除一些高反射离散噪声点;
（6）由于该物体的材料，环境和气候的复杂性，各种反射系数和衰减系数难以确定。此外，一个方程如方程（39）只能通过使用特定图象数据来解算，即，图7。
在本次研究中。我们未来的工作重点将放在如何获得更精确和更广义的公式。
5 总结
在这项研究中，我们提出了两个3D断层扫描重建方法，并尝试将其应用到河床模型试验。得出以下结论：
1.由于光学透镜的结构和理想的线性小孔成像模式的区别，在从CCD摄像机捕获的图像中产生的非线性几何失真。为了从2D图像信息获得更精确的三维地形信息，我们使用了基于（Devernay和Faugeras，1995年）线性化方法的一些显著的改善，如用Matlab支持的新算术Canny边缘检测线性方法的失真校正方法（2013a），获得边缘大小和方向，兼有黄金分割搜索方法和二次插值法计算出最优畸变参数等。实验结果表明，我们的方法像Devernay和Faugeras（1995）的方法一样有效，并且提高了精度和增加了失真参数的计算效率。
2.图像和地面之间的透视投影关系是2D图像和3D空间之间的最重要的关系之一。我们观察到的校准图像作为接地中央凸起，并且使用直接线性变换来计算变换参数，以提高参数的精度。根据传统的方法，我们使用了AMLE考虑的残余误差的第一阶的错误传播，来计算变换参数，并用于迭代加权总最小二乘法做一个迭代计算。我们得出一些结论，在控制区域的控制点应该分布尽可能均匀，最好能有更多的控制点。通过这些手段，计算出的值更接近实验值，能满足一个河床模型实验的精度要求。实验结果表明，该方法可以具有一定的实用价值。
3.灰度特征是二维图像的最重要的特征之一。本研究假定图像的灰度值仅与
每个点的实际空间位置和气候条件（时间）有关;在该条件下，无任何外部影响，在给定时刻存在有理想的灰度值和所述地形高程值时，与实际的灰度值之间的有正相关性，是理想的灰度值的各个灰度因素的影响下的变形的值。我们得出有关地形标高实际图像灰度值，灰度梯度，和太阳方向矢量一个封闭的公式，具有一定的中间变量。但是，由于种种复杂的原因，有限的条件下所获得的公式不是很满意，虽然我们已经做了广泛的尝试。今后的工作重点是在该领域中一定要了解如何获得更精确和更广义的公式。
参考文献
21