一种三峡库区消落带岸坡变化检测方法

王勋，蒋钦伟，姜正容，黄波林，顾东明

(1.重庆市地质矿产勘查开发局107地质队，重庆 401120；2.重庆一零七市政建设工程有限公司，重庆 401120；3.三峡大学 防灾减灾湖北省重点实验室，湖北 宜昌 443002；4.中国地质大学(武汉)工程学院，武汉 430074)

0 引言

三峡库区消落带不稳定岸坡的监测一直是库区地质灾害防治工作的重点和难点。受库水位周期性变化的影响和当前监测技术的限制，监测数据采集往往只针对175 m水位线以上岸坡区域，对于普遍存在土体掏蚀、局部崩塌、岩体劣化掉块等现象的消落带[1-3]，无法周期性地进行有效的岸坡形变数据采集，继而无法实现对消落带岸坡形变的长期定量对比分析，三峡库区不稳定岸坡监测存在监测盲区。

由于消落带地貌的变化机理复杂多样、变化区域狭长以及作业环境艰苦，传统的RTK、全站仪、经纬仪等仪器的单点数据采集模式无法满足地表形变监测的数据需求且工作量巨大。另外，由于岸坡坡度近乎垂直，一般的垂直摄影测量获取的三维模型无法准确细致地反映出岸坡的真实地貌[4-6]。而机载激光雷达测量系统虽然能够用于地貌扫描，但在顾及运行安全的前提下，获取的点云密度较小，对于地貌特征的描述不够精细，且三维激光扫描仪成本高昂，很难在库区大面积使用[7-8]。

贴近摄影测量技术[9]的应用优势在于能够对岸坡近距离成像，从而获取岸坡超高分辨率的影像，利用贴近物体表面(摄影距离为5～50 m)摄影的高清影像获取被摄物体精确坐标、精细形状，生产出精细的地理信息成果[10-12]。基于此，本文将“贴近摄影测量”技术引入到消落带岸坡变形调查工作中，以获取长江沿岸消落带的精细三维模型，并在精细三维模型的基础上采样得到密集三维点云数据，通过两期点云的对比，定量分析消落带岸坡的变化情况。

针对两期点云的对比，国内外学者进行了大量的研究，点云的对比方法主要分为点到点、点到模型两大类。Girardeau-Montaut 等[13]基于八叉树搜索提出了两期点云点到点(cloud-to-cloud，C2C)的形变量计算方法，直接搜索两期点云之间的最邻近点对，计算最邻近点对之间的距离作为形变量；后续其他学者提出了建立局部高程模型或最小二乘曲面拟合等方法减少点密度、噪声点以及表面粗糙度对计算结果的影响[14]。Huang等[15]提出了一种适用于复杂地形形变计算的自适应局部点云到网格的方法(adaptive local cloud-to-mesh，ALC2M)，在进行全局配准过后，利用法向量与Z轴的夹角筛选出陡坡点，建立局部参考坐标系与表面模型，计算点到模型的距离作为变形量，经过与传统点到模型的计算方法相比，能够消除差匹配对计算结果的影响。Lague等[16]提出了一种多尺度模型到模型的点云比较方法(multi-scale model to model cloud comparison，M3C2)，通过在三维点云中使用与局部表面粗糙度一致的尺度对核心点(core point)进行曲面法线估计和定向，以局部置信区间和法线方向构建空间圆柱体，计算两期点云与空间圆柱体的局部相交区域，通过局部相交区域计算两期点云模型的平均变化量，此方法能够减小噪声点以及点云模型表面粗糙度对变形量计算的影响，但是法向量的计算对结果的影响较大。

已有的算法虽然都经过严密的论证与工程应用检验，但是消落带岸坡的变化主要是掏蚀、冲蚀、崩塌以及剥落等造成，只利用最近邻点或者法向量等作为特征依据并不能较好地检测出变化结果。本文针对消落带岸坡变化的一般类型，提出了一种结合最近邻点搜索与法向圆柱搜索的检测算法。在已经完成配准的基础上，该算法首先计算出源点云与对比点云的法向量，通过KD树(K dimensional tree)搜寻两期点云之间的最近邻点，通过判断近邻点对之间法向量夹角的关系计算变形量，若近邻点对的法向量相互垂直，则变形量即为近邻点对的距离；否则以当前检测点的法向量方向为轴线构建空间圆柱体，使圆柱体与原始点云相交，计算相交处点云的重心点，则变化量即为检测点到重心点与近邻点连线线段的距离。

1 贴近摄影测量技术

1.1 岸坡消落带变化模式

自然岸坡按照物质组成包括土质岸坡、岩质岸坡、岩-土质岸坡，根据前期调查资料及现场核实，长江干流库岸再造模式主要以冲蚀、剥蚀、掏蚀为主，沿线有涉水滑坡的地段主要以滑移为主，少量段存在崩塌型库岸再造。土质岸坡受水位变动影响，整体变化表现为水对土体的掏蚀和冲蚀作用，形成大量的冲蚀沟与掏蚀沟，如图1(a)、图1(b)所示。岩质岸坡以及岩-土质岸坡整体变化表现为水对土体的掏蚀、浪蚀，以及对岩石的溶蚀作用，在周期性水位变动条件下，岩体呈现局部的强烈劣化崩解、溶蚀、滑移等现象，如图1(c)、图1(d)所示。

图1 消落带岸坡变化类型

1.2 贴近摄影测量技术

贴近摄影测量是由武汉大学张祖勋院士团队针对精细化测量需求所提出的一种全新的摄影测量技术[17]，通过贴近飞行能够达到亚厘米甚至毫米级分辨率的影像数据采集。贴近摄影测量主要有三个特点：近距离摄影、相机正对物体表面以及已知测量目标初始形态。

贴近摄影测量总体技术流程主要遵循“从无到有”“由粗到精”的基本思想。首先，采用无人机按照常规摄影测量方式获取测量目标的低分辨率影像，解算测量目标的初始地形信息；然后，根据初始地形信息、贴近飞行距离、高精度影像分辨率、影像重叠度进行三维航线规划，使无人机沿航线实现自动贴近飞行，通过航摄三维建模软件得到测量目标的高精度三维实景模型。

图2 贴近摄影测量三维航线设计示意图

贴近摄影测量技术流程最为关键的步骤为三维航线规划[18]。如图2所示，常规摄影测量飞行平面垂直于竖直方向，垂直航高为h，在获取测量目标初始地形信息后，根据真实地貌表面拟合一个空间平面即仿地平面，根据贴近摄影测量影像分辨率计算相应的倾斜航高H，从而确定贴近摄影测量飞行平面，最后根据影像重叠度规划三维航线，通过设置相机的倾斜角度参数保证相机始终垂直仿地平面进行拍摄。

2 消落带岸坡变化识别方法

通过贴近摄影测量获取到测量目标的高精度三维模型后，对三维模型进行精细化采样获得高密度的三维点云模型数据，对不同时期的三维点云模型进行对比分析，从而实现形变区域的变化识别。本文提出的变化检测算法具体步骤包括法向量计算、网格化采样、变化检测、检测结果可视化。

2.1 法向量计算

点云模型由空间中离散的点组成，为使岸坡几何特征更具真实性，需通过法向量构建可视化的三维结构特征。霍夫变换[19]可以快速检测出复杂环境中的平面，从而计算出模型的法向量，并且对具有噪声和异常值的点云模型鲁棒性较强。

2.2 网格化采样

由贴近摄影测量技术生成的三维实景模型可以高精度地还原岸坡真实地貌，但是点云数据在空间上的分布并不连续，点与点之间存在一定的间距，变化检测的结果受点间距的影响，并不能保证点与点的正确对应。此外，点云数据过于密集，会严重降低数据处理效率。原始点云越密集，采样网格尺寸越小，点与点之间的对应关系越准确，同时数据量会大幅增加，数据处理的效率会降低。所以，如要在保证检测结果准确性的同时提高效率，需要对点云数据进行网格化采样，具体算法如下。

1)根据给定的检测范围三维坐标，对点云数据进行裁剪，根据后期数据处理的要求，确定格网尺寸l，计算各坐标轴方向的格网数量。

(1)

式中：ceil(·)为向上取整标识符，返回不小于给定数字的最小整数值；Nx为X轴上网格数量，Ny为Y轴上网格数量，Nz为Z轴上网格数量；xmax为X方向上的最大值，xmin为X方向上的最小值，ymax为Y方向上的最大值，ymin为Y方向上的最小值，zmax为Z方向上的最大值，zmin为Z方向上的最小值；l为格网尺寸。

2)根据各坐标轴方向的格网数量和格网尺寸，按照式(2)遍历计算出每一个格网中心点的坐标。

(2)

式中：(centerx，centery，centerz)为每一个格网的中心点坐标，i∈(0，Nx)，j∈(0，Ny)，k∈(0，Nz)。

3)以格网中心点为原点，采用最小近邻搜索的方式查找最邻近点Pi，计算两点之间的距离d，若距离d小于格网尺寸l，则Pi即为当前格网中的唯一有效点得以保留。

4)遍历完成所有的格网，使用格网内离点云中心点距离最近的点作为格网内唯一的点云数据，以便建立点云数据的对应关系与精简数据。

2.3 变化检测

对于目标点云cloudt和源点云clouds，在进行法向量计算与网格化采样后，两期点云之间的变化量计算示意图如图3所示，具体算法如下。

1)构建clouds的KD树检索模型，对于cloudt中的每一个点P，查找其在clouds中的最邻近对应点Pne；通过三维空间中两向量之间的夹角公式计算两点法向量Dp、Dne之间的夹角α的余弦值，即有式(3)。

(3)

式中：‖‖表示模长运算。对于给定的阈值m0，若|cosα|≤m0，则认为法向量相互垂直，则P点形变量即为P点到Pne的距离。

3)按照上述算法步骤遍历计算完目标点云，即可获取整个区域内的变形情况。

图3 变化量计算示意图

2.4 检测结果可视化

当计算完区域内全部变形量后，将包含变化量的点云加载到原始三维实景模型中，根据变化量的大小对映色值建立形变变化图，以不同的颜色对应不同的变化量数值大小，故发生不同程度的变化时，其颜色分布各有不同。

3 实验与分析

3.1 实验设计

为证明本文提出的消落带岸坡变化检测方法的有效性，选取了一段消落带实地场景进行实验，设置了不同尺寸不同形态的标志物，如图4所示。同时为了消除其他因素如天气等对建模精度的影响，在同一地点前一天完成常规摄影测量影像采集并建模规划三维航线，在第二天1 h内采用已规划好的贴近摄影测量航线，完成两期影像数据的采集。另外，为验证两期影像的建模稳定性，设置了4个数据检查点，通过两期建模采样的点云坐标进行稳定性评估。

图4 变形标志物示意图

3.2 影像数据采集与建模

实验数据采集于长江沿岸巫山段消落带某一典型地貌，如图5所示。所使用的采集仪器是大疆精灵Phantom 4 RTK。Phantom 4 RTK采用双备份GNSS系统，高精度的GNSS系统采用实时差分定位技术。GPS/北斗/GLONASS三系统6频点RTK为飞行器提供厘米级定位，同时能够通过云台控制载荷相机的俯仰角和旋转角，使之能够对测量目标多角度拍摄。数据采集时，首先采用常规摄影测量方式，通过云台控制相机垂直向下拍摄，在距离崖壁最高点70 m的水平面作业，获得该区域的低分辨率影像并进行初始三维建模，然后根据初始模型规划贴近摄影测量航线，如图6所示。为保证毫米级精度的地貌获取，实验中倾斜航高设置为20 m，相片航向重叠度为80%、旁向重叠度为70%，单次采集影像41张。两期影像采集完成，进行岸坡形态三维重建，建模结果如图7、图8所示。不同尺寸变形标志物重建效果如图9所示。

图5 实验区段消落带地貌

图6 贴近摄影测量航线规划示意图

图7 第一期岸坡三维重建效果图

图8 第二期岸坡三维重建效果图

图9 标志物三维重建效果图

注：图中Point#0、Point#1、Point#2、Point#3为检查点。图10 高密度采样点云效果图

完成三维重建之后，精细化网格采样得到三维点云数据，网格大小为1 cm，对点云按照同一边界裁剪之后，单期点云数据包含约4 600万个点，采样后点云模型如图10所示。

虽然两期数据采集方式、采集环境等条件都相同，但是整体建模结果依旧会有较小的误差。通过计算两期检查点的坐标差进行建模稳定性评价，实验中设置了4个检查点，分布位置如图10所示。检查点两期点云坐标对比结果如表1所示。两期检查点x坐标差值均值为2.15 mm，y坐标差值均值为20.05 mm，z坐标差值均值为26.57 mm，可以认为两期数据的稳定性差异在30 mm以内。

表1 检查点两期点云坐标对比结果统计表 mm

3.3 岸坡变化检测

将裁剪好的两期数据配准后按照所提出的算法进行变化检测，同时为了对实验结果进行评估，利用传统最近邻点到面点云对比方法和M3C2点云对比算法进行计算，并进行对比。由检测结果可以看出，使用本文的计算方法，可以检测出放置的各个尺寸的变形标志物，且变形量符合变形标志物的实际尺寸，如图11(a)所示，变化区域明显，图11(b)中，固定变化量分布区间较小。采用最邻近点到面的计算方法，虽然可以检测出变形标志物，但是变形量与标志物实际尺寸相差较大，如图12所示。采用M3C2点云对比算法的计算结果十分杂乱，如图13所示，其中图13(b)变化量直方统计图中变现为固定变化量区间大，可视为无固定变化量，无法检测出标志物的实际形状，更无法量测尺寸。对比之下，本文方法能够较为准确地计算出变化区域的细小变形量，并对点云法向量估计的结果具有一定的鲁棒性。对检测出的变形标志物进行量测，与实际尺寸比较，结果对照表如表2所示，检测结果的比例误差均在10%以内。

表2 标志物检测结果尺寸对照表

图11 本文方法

图12 最近邻点到面计算方法

图13 M3C2点云对比算法

3.4 变化分类

由于消落带岸坡形态极其复杂，通过贴近摄影测量方式获取测量数据，依旧存在观测死角，无法对全场景内每一个目标进行准确建模，同时建模过程中，测量场景内的移动目标(如被风吹动的植被)也会导致建模结果的变化，需要人为的去判别出这类变化并予以剔除。另外，导致变化的因素有多种多样，如岩质岸坡崩解坍塌、土质岸坡冲蚀堆积等，需要人为进行分类、标注，以便后期进行有针对性的治理。将检测经过可视化渲染后结果叠加到模型上，进行人为判别并分类、标注。

4 结束语

本文基于贴近摄影测量技术，提出了一种三峡库区消落带岸坡变化检测算法，通过两期点云的对比，检测变化区域并计算出形变量。实验表明，该方法能够检测出较为细小的变化，其点云法向量的计算结果具有一定的鲁棒性，具有良好的消落带岸坡变化检测能力。