土壤侵蚀形态演化数字摄影观测系统设计与实验

姜艳敏 郭明航 赵 军 温仲明 林 奇 史海静

(1.中国科学院水利部水土保持研究所， 陕西杨凌 712100； 2.中国科学院大学， 北京 100049；3.西北农林科技大学水土保持研究所， 陕西杨凌 712100； 4.西北农林科技大学草业与草原学院， 陕西杨凌 712100；5.西安敦瑞测量技术有限公司， 西安 710065)

0 引言

中国是世界上水土流失最严重的国家之一，全国土壤侵蚀面积为356万km2，约占国土面积的40%[1]，尤其是我国的黄土高原地区，水土流失情况十分严重，对生态环境及社会经济的可持续发展带来严重影响[2-3]。土壤侵蚀是侵蚀力与下垫面相互作用的过程，下垫面既是侵蚀作用的对象又是侵蚀作用的结果[4]。在侵蚀外力的作用下，土壤或其他地面组成物质被剥蚀、搬运以及沉积，从而产生侵蚀，而侵蚀作用的结果使得原来的下垫面形态及地表条件发生改变，从而引起新的侵蚀，如此反复，使得侵蚀作用不断发展。因此，土壤侵蚀过程是一个逐渐发展演化的过程。

长期以来，各国学者围绕土壤侵蚀量调查、侵蚀过程与机理，研发了很多的观测方法和技术手段[5-9]。从传统的人工手动测量法[10]、插钎法[11]、示踪法[12]、径流泥沙采样法[13]，到目前基于测绘、测量技术以及信息技术为基础的新的观测方法，如三维激光扫描法[14]、高精度全球定位系统(Global positioning system，GPS)法[15]、摄影测量法[16]等，为土壤侵蚀观测研究提供了多样化的选择。尽管在坡面侵蚀发生发展过程方法和手段的研究方面取得了一定进展，但无论哪种方法，尚不能解决连续降雨条件下土壤侵蚀形态的观测问题，目前大多数研究集中在一次性降雨条件下侵蚀形态发育过程的主观描述。如激光扫描观测法虽有较高的空间定位精度，但不能在降雨过程中观测，且在沟蚀阶段会在沟道出现漏测[17-18]。高精度GPS法虽然测量精度高，但在地形复杂情况下易受卫星信号的影响，稳定性差，且不能在连续降雨条件下观测。因此，目前土壤侵蚀坡面观测方法仍存在测量时空关系不一致、精度不高、实时性差等问题。研究更为先进的技术和手段对土壤坡面侵蚀演化过程进行观测，对于土壤侵蚀机理、土壤侵蚀演化过程等研究具有重要意义。

近年来，随着摄影观测技术的发展，数字化近景摄影观测技术逐渐应用于土壤侵蚀测量[19-21]。文献[22-26]利用数字摄影近景观测技术，在室内模拟降雨条件下对土壤侵蚀的演变过程进行观测，但无法实现在连续降雨过程中对土壤坡面的观测。GUO等[27]研发了一种手持式数字化近景摄影观测系统，实现了连续降雨条件下的观测。该系统采用高帧率的电荷耦合装置(Charge coupled device，CCD)工业相机，快速获取高度重叠的影像，在获取地面信息的整个过程中利用手持式的垂直扫描采集方式，巧妙避开了降雨过程中的大多数雨滴。但并未对影像中的雨滴噪声做直接处理，另外，手持式的扫描使得该系统影像获取的瞬时性低，限制了系统观测的时间分辨率，从而降低了系统的实用性。因此，在连续降雨条件下去除雨滴干扰，瞬时获取土壤侵蚀坡面形态变化的信息，是土壤侵蚀过程研究亟需解决的新问题。

基于此，本文耦合数字近景摄影观测技术和无线组网技术，设计一套能够在连续降雨条件下对土壤侵蚀下垫面形态演变过程进行观测的系统。通过组网技术，并行拍摄、解算降雨过程中拍摄的下垫面的数字影像，提取具有高时间分辨率和高空间分辨率的下垫面土壤侵蚀形态演化信息，从而为土壤侵蚀过程研究提供新的途径和技术手段。

1 工作原理与系统组成

1.1 工作原理

系统基于无线网络技术对若干相机进行组网。相机基于无线网络命令，并行采集数据，数据采集时将各组传感器单次采集的数字影像按时间排序，逐像素按其灰度值做二分类处理，进而实现雨滴噪声的剔除。系统基于摄影测量技术完成下垫面对象的高精度、高密度三维点云重建。一场降雨可以获得多个时间点的三维场景数据，以达到动态的观测效果。

1.2 系统组成

图1 数字近景摄影观测系统的逻辑结构设计Fig.1 Structural design of digital close-range photogrammetry system

数字近景摄影观测系统由影像采集、影像传输和影像解算3个功能子系统组成，且每个子系统由不同的软硬件单元组成(图1)。系统的各功能子系统都在一台运行环境为Windows 7的高容量PICO(Participant intervention comparison outcome)计算机控制下运行。 并针对各功能子系统开发了对应的软件系统，以z-map命名，其包括相机工作状况诊断、影像采集、影像解算及数字高程模型(Digital elevation model，DEM)生成等功能界面。

1.2.1影像采集系统

影像采集系统负责土壤侵蚀下垫面数字影像的采集、触发信号的接收、雨滴去除等工作。该系统的硬件部分主要由12台索尼CMOS(Complementary metal oxide semiconductor)相机和工控机组成的相机组、直流电源、标靶、防水转接件等部件构成。软件部分由总控制PICO计算机z-map软件的影像采集单元组成。影像采集信号触发后，12台相机组并行采集下垫面数字影像，通过工控机对各相机单次采集的多幅数字影像，按同一位置像素单元灰度大小排序，逐像素依据其灰度运用K-means算法进行聚类处理，去除雨滴在数字影像上所形成的噪声，获得去除雨滴后的下垫面的数字影像。

1.2.2影像解算系统

影像解算系统主要负责影像数字点云的匹配、三维重建、DEM生成以及土壤侵蚀量计算等。由超高容量的匹配机来实现数据的存储、匹配、三维重建等解算工作，与影像采集系统的软件部分一样，只需通过设置数据解算后存放的路径即可完成影像数据的解算。数字影像解算系统的软件有3个模块：并行计算管理模块、点云匹配和编辑模块、DEM生成和土壤侵蚀量计算模块。在各模块算法开发过程中采用python语言配合NumPy(Numerical python)计算库来做原型的研发，之后再采用C++语言重新实现。这样的流程减少了调试过程中的时间消耗，又能保证最终执行代码的效率。

1.2.3影像传输系统

影像传输系统在采集系统和解算系统之间起连接作用，主要负责控制命令的发出、信号接收、影像数据的传输。无线路由器、网络协议(Transmission control protocol/Internet protocol，TCP/IP)、千兆网硬件接口是影像传输系统的主要硬件单元。各子系统之间通过无线路由器组成一个局域网络，控制和计算单元通过无线网络发布并发采集命令，影像采集系统采集影像后并发作业，再通过影像传输系统把采集到的影像传输给影像解算系统。

2 关键技术与实现方法

2.1 多影像采集器的组网及空间坐标系的建立

2.1.1采集器的选取与组建

借助于无线路由器通过TCP/IP网络将若干组数字影像采集器进行组网，实现影像采集器的并发作业，获取同一时间节点下的土壤下垫面信息，每组数字影像采集器包括一个数码照相机和一个工业控制级别的计算机。本系统采用的相机是索尼CMOS相机，相机分辨率为3 264像素×2 448像素，配有12 mm镜头，实用光圈为F 1.2，相机帧率为15 f/s，为保证更大的拍摄视角，共选取了12台相机。12台相机共同组建在距离地面高18 m的钢筋板架上，且相机之间呈均匀排列，与地面土槽呈垂直方向布设。与每个相机相匹配工作的硬件单元是电源和工控机，电源负责给相机和工控机供电，工控机控制相机的影像采集、雨滴去除等工作。

2.1.2空间坐标系的建立

空间坐标系的建立是将所有影像采集器的坐标进行统一。影像采集器所记录的是空间物体信息的二维图像，为了获取实际空间物体表面某点的三维几何位置，必须建立物体的三维空间坐标和对应的二维图像坐标之间的对应关系，影像采集器的几何成像模型决定了目标物体表面点的坐标与其在二维图像中的像素坐标的对应关系，而解算相机参数是建立几何成像模型的前提。参数的求解通过相机标定来完成[28]，相机的参数包括内部参数和外部参数。内部参数包括相机的焦距f，图像主点的x、y坐标(cx,cy)，畸变参数K1、K2和K3；外部参数是拍摄图像时相机的投影中心点坐标(x0，y0，z0)和3个旋转角(ψ，ω，κ)[29]。本文以针孔模型为相机标定的理论基础，借助棋盘格和控制点作为相机标定参照物，以Microsoft Visual Studio为开发平台，采用开源计算机视觉库(Open source computer vision library，OpenCV)编译相机标定程序，求解相机的参数。

图2 相机标定Fig.2 Calibration of camera

标定实验在黄土高原土壤侵蚀与旱地农业国家重点实验室人工模拟降雨实验大厅进行，实验场地为10 m(长)×4 m(宽)的液压式可调坡度钢制土槽。相机内部参数依据张正友棋盘格标定算法获取[30]。首先，将小网格长和宽均为0.05 m的棋盘格面板平整地放置于标定实验场土槽平面(图2a)。采用所有相机拍摄多幅图像，改变棋盘格面板的方位及倾斜度，再次拍摄多幅图像。将图像导入算法中求解内部参数。相机的外姿态通过在标定实验场内布设控制点来获取。选取102个小方块标志作为标定控制点，将所有控制点均匀地布设于土槽表面(图2b)。采用所有相机拍摄多幅图像，调整控制点的位置和距离，再次拍摄多幅图像。在相机内部参数已知的基础上，通过编制解算程序代码求解相机外部参数，从而得到每个控制点的x、y、z坐标。

经过标定算法及编译程序的迭代运算，得到各参数的近似值，若近似值在一定的容许范围内收敛，则结束迭代运算，得到最终的参数值。实验中各相机内部参数的标定结果如表1所示。

表1 各相机的标定参数Tab.1 Calibration parameters of each camera

2.2 数字影像中雨滴的去除

坡面数字影像的采集是在降雨条件下进行的，降雨过程中，雨滴在空间的场分布近似于随机场，相机拍摄得到的影像混合了雨滴和下垫面对象的两类信息。雨滴的去除是获取坡面物点精确信息的前提和必要工作。在短暂的时间段内，比如几秒的时间段内，下垫面对象可以认为是一个稳定的空间对象，主要变化的是随机性很高的雨场数据。根据这一思路，对各组传感器单次采集的数字影像按时间排序，逐像素按其灰度值做二分类处理，并通过K-means[31-32]算法去除雨滴在数字影像上所形成的噪声。K-means算法是一种基于形心划分的聚类算法[33]，它以数据到形心的距离作为目标函数，并以误差平方和准则函数作为聚类质量的度量函数，不断进行迭代计算求极值优化聚类结果[34]。具体算法过程如下：

(1)

(2)

min(Pix(i,j)1,Pix(i,j)2,…,Pix(i,j)n))/2

(3)

式中i、j——像素行、列数

n——迭代次数

(2)以像素的灰度距离聚类，并且构造选择集

(4)

(5)

(6)

Pix|Pix=min(dis0,dis1,dis2)

(7)

(3)以三分类选择集内元素聚合，平均得到新的类中心

(8)

(9)

(10)

(11)

对每个工控机采集的60幅原始影像像素灰度进行方差分析。结果显示，60幅原始影像的方差在28.96～29.95之间，而去除雨滴后的影像方差为24.14，由此可知K-means算法能够较好地去除雨滴噪声。

2.3 点云匹配

数字点云匹配是将多幅影像进行匹配得到同名点的过程。在匹配过程中首先提取影像上的SIFT(Scale invariant feature transform)特征[35]，SIFT特征的提取主要是依据多幅影像上具有明显特征的点，包括边缘点，具有明显特征的点、孤立的点[36]，以SIFT算子作为特征提取的工具，获取多幅影像的SIFT特征。通过SIFT特征匹配影像之间的同名点，同名点匹配是基于特征点及其描述子的相似性来进行的。获得同名点后，以光束法区域网平差原理[37]进行平差处理，以单幅影像的光线束为平差单元，以中心投影的共线方程作为平差的基础方程，将相机摄影点、相机成像点及其相应地面物点坐标作为一个整体，组建共线误差方程，通过迭代计算求出该匹配点的坐标。再通过空间前方交绘计算同名点的三维坐标，所有匹配成功的同名点都计算完毕后，可获得土壤侵蚀坡面的三维点云。

2.4 点云的修补和DEM生成

当土壤侵蚀出现沟道后，径流就会沿沟道汇集，数字影像难以拍摄到水流下方物点的影像，即难以得到水流下方物点对应的像点信息，这给水流下方数字点云的匹配带来了困难。通过对沟道存在水流情形下数字影像的分析以及匹配解算发现，由于沟底凸凹不平，总能匹配出一些稀疏的点云。基于这种客观存在，根据地形变化的连续性，利用地学普遍采用的反距离权重法[38](Inverse distance weighted，IDW)，依据水流周边和水流区域稀疏的点云坐标，便可内插得到水淹区域空缺物点的三维坐标，从而拟合出其数字点云(图3)。修补后的点云，利用同样的距离平方反比法，在系统软件模块下插值便可生成下垫面DEM。反距离权重法计算式为

(12)

式中a——参与计算的像素点数

m——像素点总数

Va——控制点的属性值

da——控制点与当前计算点间的距离

V——计算所得当前点的属性值

图3 沟道底部数字点云修补Fig.3 Digital point cloud reparation at bottom of flow channel

3 实例分析

3.1 实验布设

应用实验在黄土高原土壤侵蚀与旱地农业国家重点实验室人工模拟降雨实验大厅进行。实验小区为可调坡度钢制土槽，小区规格：长×宽×深为10 m×1.0 m×0.5 m，土槽坡度为20°(图4a)，土槽下端设置集流装置，用于收集径流泥沙。供试土壤为黄绵土，装填土容重为1.3 g/cm3。小区布设后，正式降雨前在实验土槽表面均匀撒水，再静置24 h，使土槽内部土壤水分的再分配达到应力均匀、土壤结构稳定。降雨强度设置为120 mm/h，降雨历时为150 min。至坡面产流后开始收集径流含沙量全样，每隔5 min采样1次(图4b)。降雨结束后，采用便携式径流泥沙仪测量泥沙含量。在收集径流泥沙的同时，对土壤侵蚀坡面进行全覆盖的数字影像采集，降雨开始前采集第1次坡面的数字影像，降雨开始后每隔5 min采集1次，直至降雨结束。所采集图像的重叠度至少为4°，采集帧率不少于15 f/s。

图4 实验布设Fig.4 Parallel experiments and layout

3.2 精度检测

精度是指对同一对象多次测量值的稳定程度，选用标准差来衡量。为检测数字近景摄影观测系统的测量精度，将具有标准尺寸的标尺均匀地布设于土槽表面，并在土槽表面任意位置布设2个标靶以进行长度约束。采用数字近景摄影观测法，在相同的光照和纹理条件下重复拍照60次，单独对每次的图像集合进行匹配计算，测量每把标尺的尺寸，对60次测量的结果进行统计分析并选取中误差作为衡量测量精度的指标。

通过SPSS 18软件对标尺的60次测量结果统计分析可知，标尺测量的平均长度为309.270 3 mm，测量的标准差为1.711 3 mm，单次最小测量误差为0.006 2 mm，说明该观测的精度达到毫米级。经过K-S(Kolmogorov-Smirnov)检验，得出该标尺尺寸测量结果的Z值为0.392，P值为0.999，大于0.05，由此可知该标尺尺寸的测量结果均服从正态分布，测量数据的分布如图5所示。

图5 标尺测量数据分布图Fig.5 Histogram of measured results of ruler

3.3 准确度检测

3.3.1凹槽尺寸观测法

准确度是测量值与实际值之间的偏差，以相对误差来衡量。在人工模拟降雨条件下，通过不同雨强(30、60、90、120 mm/h)、土槽坡度(0°、5°、10°、15°)共20种组合条件下，采用数字近景摄影观测系统获取土槽坡面3个已知凹槽长、宽、深的尺寸，并计算其与实际值之间的相对误差，从而评估数字近景摄影观测系统对土壤侵蚀坡面几何尺寸的观测准确度。

(1)同雨强不同坡度观测

对数字近景摄影观测系统在雨强60 mm/h，0°、5°、10°、15°共4个不同坡度条件下的观测准确度进行检测(表2)。测量值与实际值两者之间的最大相对误差为-2.556 2%，最高精度可达到99.996 8%；4个不同坡度下的平均相对误差分别为0.005 0%、-0.251 3%、-0.353 9%、-0.396 5%。对相对误差进行频率分布分析(图6a)，发现大部分观测的相对误差都较小且在0附近分布，相对误差为-0.5%～1%约占85%。以上结果表明数字化摄影观测系统对土壤侵蚀坡面的几何尺寸的观测是准确的，且坡度对该系统的观测准确度无显著影响。

表2 同雨强不同坡度观测结果Tab.2 Observation results of same rain intensity and different slopes

图6 测量值与实际值相对误差分布Fig.6 Distributions of relative error between measured and actual values

(2)同坡度不同雨强观测

对数字近景摄影观测系统在坡度为10°，雨强为30、60、90、120 mm/h条件下的观测准确度进行检测(表3)。数字近景摄影观测系统得测量值与实际值两者之间的最大相对误差为-2.968 3%，最高精度可达到99.990 1%；不同雨强条件下观测的平均相对误差分别为-0.495 8%、-0.353 9%、-0.475 1%、-0.637 6%。对相对误差进行频率分布分析(图6b)，发现大部分观测的相对误差在0附近分布，相对误差为-0.5%～1%约占75%。以上结果表明数字化摄影观测系统对土壤侵蚀坡面的几何尺寸的观测是准确的，且雨强对该系统的观测准确度无显著影响。

表3 同坡度不同雨强观测结果Tab.3 Observation results of same slope and different rain intensities

3.3.2三维激光扫描法

采用激光扫描仪观测降雨前后侵蚀坡面的三维数字地形，解算观测所得到的土壤侵蚀量(表4)。结果表明，数字近景摄影测量的土壤侵蚀总量为452 180 cm3，激光扫描观测的土壤侵蚀总量为407 971.36 cm3，由此可得数字近景摄影测量和激光扫描相对于径流泥沙法的土壤侵蚀总量观测误差分别为3.87%和6.28%，数字近景摄影测量系统可以更加精确地量化侵蚀坡面的土壤侵蚀量。激光扫描仪由于受扫描视角的限制，测量时存在扫描盲区、漏测的问题；而数字摄影观测系统由于采用多影像采集器组网技术，增大了影像采集的视野范围，可采集到足够数量沟道底部、沟壁的数字影像，从而弥补了激光扫描法在沟道观测时存在的数据缺失等缺陷。

表4 数字摄影测量与激光扫描观测对比Tab.4 Comparison of digital photogrammetry and laser scanning observation

3.3.3径流泥沙观测法

采集径流泥沙全样，是观测土壤流失量最为可靠的方法，在本次检测中以径流泥沙观测法得到的结果作为实际值。采用数字摄影观测系统将降雨过程中不同时间点的数字影像进行体积解算得到土壤流失量结果，并与相同时间段收集的径流泥沙含量进行对比(表5)。结果表明，两种观测方法的平均相对误差为-1.73%；从不同侵蚀阶段两种观测方法的相对误差来看，在降雨初期，两种观测方法相对误差比较大，在降雨历时达到50 min时，观测精度开始变高，说明此时正是坡面沟道快速发育的明显分界点，沟道快速发育前后两种观测方法观测的相对误差分别为20.85%和3.47%；随着降雨历时的延长，土壤坡面侵蚀沟发育形态的变化越来越明显，数字近景摄影观测系统的观测精度逐渐提高，观测精度最高可达99.26%。

数字点云密度代表着数字点云对坡面地表形态表达的准确程度，高密度的数字点云可以更加精确地表达坡面侵蚀沟的形态信息。通过不同时间节点侵蚀坡面数字点云的数量和坡面观测面积计算数字点云密度。结果表明，数字点云的平均数量为1.335 1×106个，平均点云密度为0.134个/mm2。表5中列出了由不同时间节点侵蚀坡面的高密度数字点云转换生成的DEM，其空间分辨率可达到2 mm，可准确表达侵蚀形态的空间变化，实现了对土壤侵蚀坡面形态变化过程的动态监测。

表5 不同时间点土壤流失量的观测结果Tab.5 Observation results of soil loss at different times

4 结论

(1)提出了一种在连续降雨过程中对土壤侵蚀坡面动态监测的数字近景摄影观测方法，通过无线组网技术，并行拍摄、解算降雨过程中下垫面的数字影像，提取坡面精细地貌动态变化的信息。该系统时间观测分辨率可达到分钟级别，空间分辨率达到2 mm。可从时间和空间尺度上更加准确地描述土壤侵蚀过程，解决了侵蚀观测中时空不一致的问题。

(2)与传统径流泥沙法的平行观测结果比较表明，数字化摄影观测法在坡面土壤侵蚀过程的不同阶段其准确度不同，随降雨历时的延长，数字近景摄影观测系统的观测精度逐渐提高，土壤流失量估算平均相对误差为-1.73%，单次观测精度最高可达99.26%。

(3)与激光扫描法平行观测结果比较表明，数字近景摄影测量系统对坡面土壤侵蚀量的观测精度高于激光扫描仪。数字近景摄影观测法采用多影像采集器组网技术，增大了影像采集的视野范围，克服了激光扫描仪观测时沟道底部激光线不能投射到位而造成的漏测现象，实现了全覆盖的数字影像采集。