一种水土流失情况测定的方法

摘要

本发明涉及一种水土流失情况测定的方法。方法包括：图像采集、对采集的图像进行分块处理、对各区块分别进行三维重建、拼接三维重建后的各区块得到该时刻和根据不同时刻的三维点云计算水土流失量。本发明的方法能够在降雨过程中观测。通过拍摄、解算降雨过程中下垫面的数字影像，提取土壤侵蚀过程中下垫面变化的时空信息，更加准确地描述土壤侵蚀过程。同时本发明也能够应用于类似的体积变化的量测。

权利要求

1.一种基于立体视觉的图像三维重建方法，其特征在于，方法包括：步骤一，采集图像

采集目标区域某一时刻时间窗口的图像；

步骤二，区域分块

按采集时间先后对所有图像进行序列标号；每连续多帧图像为一任务区块；相邻区块与区块之间有重叠图像，重叠部分的图像有相同序列标号；

步骤三，三维重建

对步骤二得到的各区块分别采用立体视觉技术进行三维重建，得到各区块的三维点云；

对各区块三维点云的各三维点进行编号，得到各三维点的编号特征；任意三维点A的编号特征为[(A1，a1),(A2，a2),…,(An，an),…,(AN-1，aN-1),(AN，aN)],其中：n＝1,2,3,…,N；N为三维点A对应的原始图像的张数，N取正整数；An为三维点A对应的所有原始图像中某一图像n的ID号；an为三维点A在ID号为An的原始图像上的SIFT特征点的ID号,(An，an)为编号特征的编号元组；

步骤四，拼接

(1)确定邻接区块

①寻找有重叠图像的相邻区块与区块之间的同名点:相邻区块与区块之间的某两个三维点的编号特征交集包括至少一个编号元组时，两个三维点为同名点；

②两个区块的同名点的个数占其中一个区块中总三维点数的10％以上时，这两个区块为邻接区块；

③同名点个数占区块中三维点个数比例最高的区块为母站；

(2)邻接区块拼接

以母站为基准对其他区块进行拼接：

(2.1)确定稳定的同名点

①构建同名边

对于两个邻接区块，分别在各区块中利用区块中的同名点组合构建边网，并计算每条边的长度；其中一个区块中有n个同名点，则该区块中组合构建的边网有 条边，两个邻接区块之间具有 组同名边，同名边对应同名点；

②确定稳定的同名点

邻接区块之间的各组同名边中两条边长度的比值构成比例数组，比例数组包括 个数值；采用基于拉依达准则(3σ准则)的奇异数据滤波法对 个数值进行检验，得到可靠数值点，可靠数值点对应的同名点作为稳定的同名点；

③基于稳定的同名点做两个邻接区块之间的旋转、平移、缩放参数的计算，从而完成两个区块之间的拼接；

完成所有邻接区块之间的拼接后得到该时刻下测量区域的点云集合，去除点云集合中观测区域之外的点云，得到该时刻下测量区域的三维点云。

2.一种水土流失情况测定方法，其特征在于，方法包括：步骤A，构建参考面：采集目标区域降雨前的图像，执行权利要求1所述步骤二至步骤四，利用所得三维点云中的多个三维点拟合得到的平面为参考面；

步骤B，获取不同时刻目标区域下雨时的三维点云：每一时刻的图像采集和处理执行权利要求1所述步骤一至步骤四；

步骤C，统一坐标系：将降雨时各时刻三维点云的坐标统一到同一坐标系 下；

步骤D，计算各时刻参考面以下的DEM体积量:将某时刻的三维点云内插得到对应的DEM数据，计算该时刻参考面以下的DEM体积量；

步骤E，相邻时刻参考面以下的DEM体积量差值为相邻时刻时间段内的土壤侵蚀量。

3.权利要求2所述的水土流失情况测定方法，其特征在于，步骤A和步骤B的图像采集包括：采集目标区域某一时间窗口的图像，所采集图像的重叠度至少为5度；单位平方米内，在

10-12秒内采集图像的数量为120-240帧，且单位平方米至少采集120帧。

4.权利要求2所述的水土流失情况测定方法，其特征在于，所述目标区域设有多个标靶。

5.如权利要求3所述的水土流失情况测定方法，其特征在于，所述时间窗口时长：15-

180s。

6.如权利要求3所述的水土流失情况测定方法，其特征在于，步骤二按采集时间先后对所有图像进行序列标号，每连续50-230帧图像为一任务区块；相邻区块与区块之间有图像重叠，重叠图像的帧数占两个相邻区块中每个区块内图像帧数的30％-45％。

说明书

一种水土流失情况测定的方法

技术领域

[0001] 本发明涉及图像三维重建技术和水土流失情况测定方法。

背景技术

[0002] 土壤侵蚀过程是侵蚀力与下垫面相互作用的过程，下垫面既是土壤侵蚀作用的对象，也是土壤侵蚀作用的结果，在土壤侵蚀过程中的每一时刻,都使得下垫面的土体发生位移,下垫面的形态发生变化。围绕观测下垫面土体的位移量和下垫面形态的变化，发展了多种多样的观测技术和方法：

[0003] ⑴接触式观测法

[0004] 通过皮尺、卷尺、游标卡尺等常规的测量工具，对下垫面的形态变化进行简单测量。优点：简单易行；缺点：精度偏低，且只能观测次降雨条件下土壤侵蚀的总量，而不能观测土壤侵蚀过程中在不同时段或不同空间的变化量；

[0005] ⑵示踪观测法

[0006] 利用137Cs、7Be等天然放射性核素、稳定性稀土元素等按照指定要求添加到土壤里面，对下垫面的形态进行监测。优点：可观测土壤侵蚀速率、空间分布特征；缺点：观测的时间间隔太长、空间定位精度一般在米级以上；

[0007] ⑵维激光扫描观测法

[0008] 通过三维激光扫描仪对下垫面形态进行测量；原理：利用激光测距技术，通过记录被测物体表面大量的密集的点的三维坐标、反射率和纹理等信息，可快速测量被测物表面的三维点数据。优点：空间定位精度高；缺点：仪器不防水，不能在降雨过程中观测，且在沟蚀阶段会在沟道出现漏测，最大缺陷是测量效率太慢，不能监测下垫面形态动态演变过程；

[0009] ⑷传统摄影测量观测法

[0010] 用于监测河道侵蚀和重力侵蚀，其观测时间间隔分别为次降雨、天、年甚或更长，适合长时段侵蚀总量的观测，而不能观测短时段内相应分量的变化；

[0011] ⑸径流泥沙采样法

[0012] 在降雨过程当中，利用容器对下垫面侵蚀产生的泥沙按照不同时间段分别采样，对下垫面的土壤侵蚀量监测。优点：简单易行；缺点：泥沙的空间来源无法区分；

[0013] 那么，如何在降雨条件下,从时空两个维度上观测下垫面形态演化以及表征土壤侵蚀各物理量。亦即，如何随着土壤侵蚀的发生发展，能够在设定的时间结点上，观测并量化表达下垫面的空间形态，进而根据下垫面空间形态计算土壤侵蚀量、输沙率、侵蚀速率等，以上的观测技术和方法都未能很好地解决这个问题。

发明内容

[0014] 针对现有技术的缺陷或不足，本发明提供了一种基于立体视觉的图像三维重建方法。

[0015] 本发明提供的基于立体视觉的图像三维重建方法包括：

[0016] 步骤一，采集图像

[0017] 采集目标区域某一时间窗口的图像；

[0018] 步骤二，区域分块

[0019] 按采集时间先后对所有图像进行序列标号；每连续多帧图像为一任务区块；相邻区块与区块之间有重叠图像，重叠部分的图像有相同序列标号；

[0020] 步骤三，三维重建

[0021] 对步骤二得到的各区块分别采用立体视觉技术进行三维重建，得到各区块的三维点云；

[0022] 对各区块三维点云的各三维点进行编号，得到各三维点的编号特征；任意三维点A的编号特征为[(A1，a1),(A2，a2),…,(An，an),…,(AN-1，aN-1),(AN，aN)],其中：n＝1,2,3,…,N；N为三维点A对应的原始图像的张数，N取正整数；An为三维点A对应的所有原始图像中某一图像n的ID号；an为三维点A在ID号为An的原始图像上的SIFT特征点的ID号,(An，an)为编号特征的编号元组；

[0023] 步骤四，拼接

[0024] (1)确定邻接区块

[0025] ①寻找有重叠图像的相邻块与块之间的同名点，相邻块与块之间的两个三维点的编号特征交集包括至少一个编号元组时，两个三维点为同名点；

[0026] ②两个块的同名点的个数占其中一个块中总三维点数的10％以上时，这两个区块为邻接区块；

[0027] ③同名点个数占区块中三维点个数比例最高的区块为母站；

[0028] (2)邻接区块拼接

[0029] 以母站为基准对其他区块进行拼接：

[0030] (2.1)确定稳定的同名点

[0031] ①构建同名边

[0032] 对于两个邻接区块，分别在各区块中利用该区块中的同名点组合构建边网，并计算每条边的长度；其中一个区块中有n个同名点，则该区块中组合构建的边网有 条边；

[0033] ②确定稳定的同名点

[0034] 邻接区块之间的各组同名边中两条边长度的比值构成比例数组，比例数组包括个数值；采用基于拉依达准则(3σ准则)的奇异数据滤波法对 个数值进行检验，得到可靠数值点，，可靠数值点对应的同名点作为稳定的同名点；

[0035] ③基于稳定的同名点做两个三维区块之间的旋转、平移、缩放参数的计算，从而完成两个区块之间的拼接；

[0036] 完成所有邻接区块之间的拼接后得到该时刻下测量区域的点云集合，去除点云集合中观测区域之外的点云，得到该时刻下测量区域的三维点云。

[0037] 本发明还提供了一种水土流失情况测定方法，方法包括：

[0038] 步骤A，构建参考面：采集目标区域下雨前的图像，执行上述步骤二至步骤四，利用所得三维点云中的多个三维点拟合得到的平面为参考面；

[0039] 步骤B，获取不同时刻目标区域下雨时的三维点云：每一时刻的图像采集和处理执行上述步骤一至步骤四；

[0040] 步骤C，统一坐标系：将下雨时各时刻三维点云的坐标统一到同一坐标系下；

[0041] 步骤D，计算参考面各时刻参考面以下的DEM体积量:将某时刻的三维点云内插得到对应的DEM数据，计算该时刻参考面以下的DEM体积量；

[0042] 步骤E，相邻时刻参考面以下的DEM体积量差值为相邻时刻时间段内的土壤侵蚀量。

[0043] 步骤A和步骤B的图像采集包括：采集目标区域某一时间窗口的土壤图像，所采集图像的重叠度至少为5度；单位平方米内，在10-12秒内采集图像的数量为120-240帧，且单位平方米至少采集120帧。

[0044] 所述目标区域设有多个标靶。

[0045] 所述时间窗口时长：15-180s。

[0046] 步骤二按采集时间先后对所有图像进行序列标号，每连续50-230帧图像为一任务区块；相邻区块与区块之间有图像重叠，重叠图像的帧数占两个相邻区块中每个区块内图像帧数的10％-20％。

[0047] 本发明的技术效果为：

[0048] 1、本发明的方法能够在降雨过程中观测。通过拍摄、解算降雨过程中下垫面的数字影像，提取土壤侵蚀过程中下垫面变化的时空信息，更加准确地描述土壤侵蚀过程。同时本发明也能够应用于类似的体积变化的量测。

[0049] 2、本发明的方法观测精度高，方法可观测一定时间段内自然界土壤侵蚀过程，对土壤侵蚀观测的空间和时间分布率更高。

[0050] 3、本发明的方法对计算机硬件单元适用范围宽，可以根据实际的硬件环境，计算过程中的多个区块可通过联网进行并行计算或通过单机分批次计算。

[0051] 4、本发明在确定邻接区块时采用SIFT特征算子的ID号构造编号组集合，通过区块之间的集合交集计算同名点，较之传统特征匹配算法计算速度至少提高了4倍。

附图说明

[0052] 图1A为实施例中降雨前下垫面的示意图；图1B为实施例中降雨后下垫面的示意图；

[0053] 图2为实施例中步骤一所采集图像的示意图；

[0054] 图3为实施例中步骤三各区块的点云示意图；

[0055] 图4为步骤四中三维重建后的示意图；

[0056] 图5为步骤五中的参考面及计算示意图；

[0057] 图6A为激光扫描仪观测得到的数字点云图；图6B为本发明的方法得到的数字点云图。

具体实施方式

[0058] 以下是提供的具体实施例，以对本发明的技术效果作进一步解释说明。

[0059] 实施例：

[0060] 该实施例于2014年7月9日在黄土高原土壤侵蚀与旱地农业国家重点实验室人工模拟降雨厅进行。

[0061] 实验中人工模拟降雨的强度为90毫米每小时，所构造的下垫面如图1A、图1B所示，其中图1A为未发生侵蚀的(降雨前拍摄)下垫面图片，图1B为发生侵蚀的(降雨后拍摄)下垫面图片，钢槽底部为斜面，钢槽内部放置有待测定土壤，土壤类型为黄绵土，土壤装填容重为1.3g/cm3，下垫面装置的土槽尺寸为1米宽，5米长的立方体容器，下垫面底面与水平面构成15度倾斜角度。

[0062] 分时段进行摄影并进行数字影像的解算，计算机视觉的下垫面土壤侵蚀观测法：从降雨开始至降雨结束对实验下垫面进行全覆盖的数字影像采集，采集时间分别在降雨开始后50min、70min、90min、100min，第100min后每隔10min采集一次，采集器距离下垫面的平均高度为80±5cm，数字影像的采集密度为150-170帧/m2；

[0063] 收集水沙全样并换算流失量：在降雨开始后按照计算机视觉的下垫面土壤侵蚀观测的相同时间间隔收集实验过程中的水沙全样,计算样品中所含的土壤重量，并按土槽装填时1.3g/cm3的容重换算相应的体积量。

[0064] 图像采集和图像处理过程如下：

[0065] 一、获取不同时刻下测量区域的三维点云

[0066] 步骤一，采集图像

[0067] 区域采集面积：5平方米，本次采集得到1013张照片，采集过程历时62秒，如图2所示；

[0068] 步骤二，对步骤一所采集的所有图像进行区域任务分块

[0069] 按采集时间先后对步骤一所有图像进行序列标号，将步骤一中得到的1013张数字影像，以每块180张、相邻区块重叠65张的任务分割模式分割为9个任务子区，其中第9子区具有93张照片；

[0070] 步骤三，三维重建

[0071] ①对步骤二得到的各区块分别采用立体视觉技术进行三维重建，得到各区块的数字点云。

[0072] ②特征点编号

[0073] 对步骤三所得到的各个区块中三维点进行特征编号：任意三维点A的编号特征为[(A1，a1),(A2，a2),…,(An，an),…,(AN-1，aN-1),(AN，aN)],其中：n＝1,2,3,…,N；，N为三维点A对应的原始图像的张数，N取正整数；An为三维点A对应的所有原始图像中某一图像n的ID号；an为三维点A在ID号为An的原始图像上的SIFT特征点的ID号,(An，an)为编号特征的编号元组；比如点云中的P点，该点在4张原始图像成像，其四张图像的ID号码是100,101,104,106，分别在其上的SIFT特征点的ID号码是：52,78,99,48.则编号特征表.[(100，52),(101，
78),(104，99),(106，48)]；如图3所示。步骤四，拼接

[0074] (1)确定邻接区块

[0075] ①有同名图像的区块为邻接区块；寻找有同名图像块与块之间的同名点，三维点与三维点的编号组集合的交集包括至少两个元素时为同名点；

[0076] ②采用sift描述子匹配技术确定没有同名序列区块之间的同名特征点，两个区块的同名点的个数占其中一个区块中总三维点数的10％以上时，这两个区块为邻接区块；

[0077] 定义同名点个数占区块中点云个数比例高的区块为母站；

[0078] (2)邻接区块拼接

[0079] 以母站为基准对其他区块进行拼接：

[0080] (2.1)确定稳定的同名点

[0081] ①构建同名边

[0082] 对于两个邻接区块，分别在各区块中利用该区块中的同名点组合构建边网，并计算每条边的长度；例如：设其中一个区块中有n个同名点，则该区块中组合构建的边网有条边；

[0083] ②确定稳定的同名点

[0084] 邻接区块之间的各组同名边长度的比值构成比例数组，包括 个数值；

[0085] 采用基于拉依达准则(3σ准则)的奇异数据滤波法对 个数值进行检验，得到可靠数值点，可靠数值点对应的同名点作为稳定的同名点：

[0086] 对比例数组中的数值按大小排序，接着去掉排序后数组中的前、后各占总数值个数3％-5％的数值，得到中间数组，然后求取中间数组的平均值和方差，以3倍方差作为置信区间，中间数组中的数值落在置信区间以内的为可靠数值点，可靠数值点对应的同名点作为稳定的同名点；

[0087] 基于稳定的同名点做两个三维区块之间的旋转、平移、缩放参数的计算，从而完成两个区块之间的拼接；

[0088] 完成所有邻接区块之间的拼接后得到该时刻下测量区域的点云集合，去除点云集合中观测区域之外的点云，得到该时刻下测量区域的三维点云，如图4所示；

[0089] 以两个标靶之间的实际尺寸为参考，确定三维点云的实际尺寸；

[0090] 重复步骤一至步骤四，得到其他时刻下测量区域的三维点云；

[0091] 二，计算土壤侵蚀量

[0092] (1)定义参考面：以未发生侵蚀的(降雨前)三维点云为基础拟合参考平面，具体为：选取150到200个三维点，这些三维点经过拟合得到一张平面，该平面即定义为整个测量过程中的参考面。

[0093] (2)统一降雨时各时刻三维点云的坐标系：某一时刻的三维点云作为一个基准，统一按照均匀分布在下垫面的周边的4-6个标靶的空间位置进行空间相似变换，由此，三维点云的坐标系统纳入到标靶坐标系统中，标靶在整个降雨过程中空间位置是稳定不变的，因此所有的点云数据都具有了统一的坐标系。

[0094] (3)参考面以下体积的计算：将某一时刻三维点云内插得到对应的DEM数据。人机交互在DEM数据上指定方向参考点，如图5中十字交叉点，计算参考面以下的DEM体积量。

[0095] (4)相邻时刻体积量的较差即为该时段的侵蚀量。

[0096] 该实施例参考面相关数据为：图5中格网为参考平面；十字为方向参考点位置；参考面左顶点:(1.319944,160.862244)；参考面右底点:(59.419609,-15.33840)；DEM步距：0.1毫米；方向参考点坐标：x＝20.219984,y＝79.760979,z＝-21.923674。

[0097] 表1本发明的方法与径流泥沙采集法观测结果对比表

[0098]

[0099]

[0100] 与传统的径流泥沙采集法平行观测结果表明全程平均准确度为1.05％。

[0101] 将本发明的检测方法与三维激光扫描仪检测观测方法进行了对比：

[0102] 采用实施例最后一次的观测结果与三维激光扫描仪观测结果进行对比，在实验观测降雨结束后，用三维激光扫描仪和本方法观测并计算下垫面的三维数字点云，结果如表2和图6所示。

[0103] 表2本发明的方法与激光扫描观测结果对比

[0104]

[0105] 与激光扫描仪平行观测表明：本发明的方法克服了激光扫描仪不能够在降雨过程中观测下垫面演化的缺陷；二是下垫面形态的数字点云的密度提高了27.8％，即本发明的方法测量每个点云所代表的物点尺寸要更小，对下垫面三维形态的表达要更准确；三是本发明的方法克服了激光扫描仪观测时沟道底部激光线不能投射到位而造成的漏测现象(图6A中黑色部分区域是漏测区域)，实现了全覆盖的数字影像采集。