激光雷达技术在应县木塔形变分析中的应用

郭 明，闫冰男，周腾飞，陈朝阳，张 晨，刘运明

(1.北京建筑大学 测绘与城市空间信息学院，北京 102616；2.北京建筑大学 代表性建筑与古建筑数据库教育部工程研究中心，北京 102616；3.北京建筑大学 现代城市测绘国家测绘地理信息局重点实验室，北京 102616；4.北京建筑大学 建筑遗产精细重构与健康监测北京市重点实验室，北京 102616；5.北京城建勘测设计研究院有限责任公司，北京 100101)

0 引 言

应县木塔又称佛宫寺释迦塔，位于山西省应县西北隅，建于辽代清宁二年(1056年)是中国现存最高最古老的木构塔式建筑，历史文化价值极高[1]。长期以来，应县木塔经历过多次自然和人为破坏，木材承载能力减弱，部分构件损坏，使木塔产生不同程度的倾斜、扭转等形变，并且随着时间的推移日趋严重[2]。应县木塔是中国木结构古建筑的典型代表，对应县木塔进行保护也是古建筑文化遗产保护领域的重要工作，通过测绘手段获取木塔现状数据，对其形变进行全面分析，以便制定合理的保护方案，科学保护木塔安全。

近年来大多学者从结构的角度采用有限元方法分析应县木塔的形变及结构受力情况[3]。陈特等[4]通过有限元模型评估应县木塔的稳定性，分析出木塔二层明层的破坏始于西南面外檐中柱的局部失稳，该柱因柱脚榫局部最大拉应力过大而发生损伤。陈平等[5]建立应县木塔二层明层柱网的有限元模型并进行纠偏分析，为木塔的进一步加固设计提供依据，并提出纠偏方案。杜雷鸣等[2，6]构建有限元分析模型，进行了重力和地震作用下的结构受力分析，对木塔结构抗震性能进行了安全评估。也有利用工程测量方法研究应县木塔，隋坤等[7]分析出木塔总的趋势是南高北低，各层间处于相对扭曲状态。由于古建筑木结构与现代混凝土建筑结构不同，单柱的局部失稳将引起木塔整体坍塌[8]，因此本文主要测量单柱的形变，从测绘的角度定量分析应县木塔的形变情况。

激光雷达扫描技术是一种全新的测绘手段，凭借其高效率、高精度、不接触古建筑本体的优势，在古建筑勘测中被越来越多地应用。通过激光雷达扫描技术可以获取古建筑精细的三维信息[9-10]，是古建筑文化遗产管理的重要手段和工具[11]；全站仪的测量精度可达到亚毫米级，在密集型遗产监测中常作为一种全局控制手段[12]。本文通过不同的地面激光雷达扫描仪获取木塔内、外的点云数据，同时在木塔周围布设标靶控制点，通过全站仪控制测量建立监测的绝对坐标系。采集到的点云数据首先要进行预处理，外部点云采用整体配准的算法，内部点云采用特征配准的粗配准和迭代最近点(ICP)算法精配准的两次配准，再通过控制点将内、外点云转换到绝对坐标系中，得到整体的点云模型。通过对点云模型进行多种剖切，测量出木塔各层柱子的偏移角度、距离以及整体的倾斜角度等，根据测量数据分析木塔单柱、单层到整体的偏移程度、倾斜方向以及扭转趋势等。

1 监测方案及数据获取

1.1 监测方案

根据应县木塔的地理位置、周围环境、建筑结构以及扫描精度要求，制定合理的监测方案，如图1所示。利用中远程、近程地面激光扫描仪以及全站仪等测量传感器，获得能反映木塔特征变化的高精度监测数据。通过中远程地面激光扫描仪采集木塔外部点云，但是木塔屋檐容易遮挡并且无法获取木塔内部点云，这时通过近程扫描仪采集内部点云；同时利用全站仪进行控制测量，根据已有监测控制点的分布在木塔周围布设标靶纸，每层至少布设3个标靶点。将塔内、塔外点云配准后，经过控制点进行坐标转换，得到整体的点云模型，实现传统测量与激光雷达扫描的坐标统一，便于进行形变分析。为了保证数据拼接质量，相邻站点之间数据重合度应达到30%以上，相邻站点之间至少有3个不共面的公共标靶。

图1 监测数据获取技术路线

1.2 塔外数据获取

中远程激光扫描仪的测程可以达到1 200 m，适用于远距离大范围的扫描。采集数据时在应县木塔外部设置10个测站(图2)，共布设25个标靶控制点。为了保证数据配准和坐标转换的精度，对25个标靶点控制进行高精度扫描，标靶纸及标靶点云如图3所示。

图2 外部扫描站点分布

图3 标靶纸及标靶点云

在单站扫描过程中，设置粗略扫描精度为0.04 cm·m-1，对木塔所在区域做高精度扫描，精度为0.01 cm·m-1，单站扫描数据如图4所示。

图4 塔外单站扫描数据

1.3 塔内数据获取

近程激光扫描仪的测程为0.5～130 m，扫描仪轻巧，便于携带，因此用于采集木塔内部数据。应县木塔每层均由8根内柱和24根外柱支撑，在扫描过程中对内中外三侧进行扫描(图5)，设置扫描分辨率为0.07 mm·m-1，保证每根柱子都能扫描清楚且没有遗漏部分。为方便数据配准，每两站之间通过标靶球连接(图6)，木塔上下2层通过楼梯上特征点、特征面连接。

图5 内部扫描站点分布

图6 标靶球

1.4 控制点数据获取

为了精确测定应县木塔的形变，需要建立较高精度的控制网[13]，通过全站仪进行控制测量，建立应县木塔的监测绝对坐标系。应县木塔院内已经埋设了固定的监测控制点，根据这些控制点的分布在木塔周围布设标靶纸，标靶纸易于粘贴和清理，不会对木塔造成伤害。标靶控制点布设应该精简，外侧主要分布在木塔东南面和西面，内侧布设在柱子上，每层至少布设3个标靶点。木塔东南面标靶纸分布如图7(a)所示，图7(b)，(c)分别为外部和内部布设的标靶纸。

图7 标靶纸控制点

2 监测数据预处理

监测数据预处理即为点云配准，点云配准是获得完整点云模型的必要步骤，配准质量也关系到后续成果的整体质量[14]。外部点云采用整体配准方法，内部点云首先通过特征进行粗配准再使用ICP算法精配准，最后由标靶控制点将内外点云均配准到绝对坐标系中。点云配准通过罗德里格矩阵求解空间转换参数，最终得到带有绝对坐标的整体点云模型，配准流程如图8所示。

图8 点云数据配准流程图

2.1 外部点云配准

整体配准将多站数据根据其中的约束关系一次性转换到统一的坐标系下，进行整体配准时必须以高精度的控制约束为基础[15]，同时点云中要具备足够的约束求解空间转换参数[16]。获取监测数据时，在应县木塔外部布设大量标靶，可以将标靶控制网作为基准，因此外部点云适用于整体配准。以标靶点构造点约束误差方程，则扫描点云中的标靶点X0(x0,y0,z0)与其观测值X(x,y,z)有如下关系[17]

X0-(-λRXt+ΔX)=0

(1)

式中：R为站点变换旋转矩阵；ΔX为平移参数，ΔX=[(Δx,Δy,Δz)]T；λ为尺度参数，一般在点云运算中取λ=1。

大旋角空间变换的角度参数求解通常用罗德里格矩阵，由反对称矩阵S[式(2)]构建，R=(I+S)·(I-S)-1，其中I为3阶单位阵。

(2)

式中：a,b,c为3个独立的参数。

结合罗德里格矩阵线性化展开得到点的误差方程为

V1=A1t+BX-L1

(3)

V1=[VxVyVz]

t=[dΔxdΔydΔzdadbdc]

B=R-1

X=[dx0dy0dz0]T

式中：V1为观测值改正数；A1为空间变换参数相关的系数矩阵；t为空间变换参数改正数；B为待定点系数矩阵；X为待定点改正值；L1为观测值残差。

配准后点云的平均误差在5 mm以下，如图9所示。

图9 应县木塔外部点云

2.2 内部点云数据配准

内部点云是形变分析的主要数据，对数据的精度要求高，因此进行两次配准。由于在采集木塔内部数据时布设了标靶球，其他的同名特征也比较明显，所以先利用特征进行粗配准，同一层相邻测站点云以标靶球为特征进行配准，上下层点云以楼梯上特征点、特征面进行配准。粗配准后点云大致重合但精度不高，此时利用ICP算法再次配准能提高配准精度。

(4)

重心坐标分别为Xm和Xn，重心化坐标为

旋转参数的误差方程可以表示为

V=At-L

(5)

求解罗德里格参数构造矩阵R后，平移参数的误差方程可表示为

(6)

Y=C(P,X)

(7)

式中：C为对应点的坐标变换，不停迭代直到满足精度为止。

粗配准后的点云误差在7 mm以下，经过ICP算法再次配准的误差在2 mm以下，二层至五层的内部点云模型如图10所示。

图10 应县木塔内部点云

2.3 点云数据坐标转换及内外部数据配准

由于内部点云和外部点云都完成了配准，因此内外部数据可看作两站点云，将内部点云和外部点云均配准到控制点坐标下，既实现了内外部点云的配准，也完成了点云坐标的转换。根据3个同名标靶点坐标可以求出旋转参数和平移参数，存在多余观测时,可根据最小二乘平差求出最优解，两站之间的空间转换模型为[21]

(9)

[XYZ]T和[xyz]T表示相邻两测站中的同名特征点坐标，结合罗德里格矩阵求出点误差方程，如式(5)所示，由于各控制点的权重相同，根据最小二乘准则得

X=(ATA)-1ATL

(9)

求出旋转参数的最优解，代入转换模型求出平移参数[ΔXΔYΔZ]T。配准后平均误差在5 mm以下，配准后的数据未出现分层情况，满足精度要求，整体点云如图11所示。

图11 应县木塔完整点云

3 形变分析

应县木塔形变主要是由柱子形变引起的，因此主要分析柱子形变。从不同视图的整体点云中截取柱子点云，测量出各柱子的倾斜角度、偏移距离；对模型整体同样进行剖分，测量整体的倾斜角度。根据测量数据分析单柱倾斜、单层位移、柱结构纵向倾斜、木塔整体倾斜和扭转等形变情况。

3.1 单柱倾斜分析

截取每根柱子的完整点云，在正视图和左视图中量测柱子的倾斜角度。将柱头、柱脚的点云数据拟合成圆，柱子边缘的点云数据拟合成直线，连接柱头、柱脚的圆心得到柱子的中心线，同时过柱脚圆心做出柱子垂线，垂线与中心线所形成的夹角视为柱子的倾斜角度。图12为应县木塔第2层部分柱子的倾斜角度。西南面柱子倾斜角度最大，向东最大倾斜3.43°，向北最大倾斜5.05°；东北面柱子倾斜角度最小，向东最大倾斜0.43°，向北最大倾斜0.26°；各面柱子向东北方向有不同程度的倾斜，西南面柱子倾斜情况最严重，东北面柱子倾斜角度小。

3.2 单层位移分析

截取每层柱头、柱脚的点云，根据点云拟合出圆，叠加2层数据量测出柱头圆心到柱脚圆心3个方向的位移，横向偏移表示沿东西方向的位移，纵向偏移表示沿南北方向的位移，综合偏移表示柱头、柱脚圆心的直线距离，图13为第2层柱子的偏移。

图13 应县木塔第2层位移

第2层西南面(M2W22～M2W24)外柱在3个方向的位移都很大，其中M2W23号外柱形变最大，向东偏移0.356 m，向北偏移0.326 m，综合偏移0.483 m，这面柱子严重向东北方向倾斜；南面、东南面、东面(M2W01～M2W10)外柱纵向偏移较大，柱子向北方向的倾斜程度比向东方向严重，并且倾斜程度逐渐减少；东北面(M2W10～M2W13)是平均倾斜程度最小的面，柱子变形程度最小，其中M2W12号外柱的南北位移是第2层中最小的，向北偏移0.002 m；北面(M2W13～M2W16)外柱横向偏移较大，柱子向东方向的倾斜程度较大；西北面和西面(M2W16～M2W22)外柱偏移量逐渐增大，柱子向东北方向的倾斜程度逐渐增大。内柱的位移量和外柱不同，但同方位内柱和外柱的偏移方向相同，东北角的M2N05号内柱变形最小，向东位移值最小为0.002 m，综合位移最小为0.005 m。其他层的偏移方向与第2层类似，但位移量偏小，分析得出应县木塔各层均向东北方向有不同程度的位移，每层西南面柱子位移最大，东北面柱子位移最小。

3.3 柱结构纵向倾斜分析

为了直观分析各层柱子的倾斜，从点云模型的俯视图中沿对角线方向截取整体的柱子点云，图14为北偏东67.5°方向上截取的点云。纵向来看，应县木塔各层的倾斜程度不同，二层和三层的倾斜程度远大于上部2层，第5层的倾斜程度相对最小，第2层的倾斜程度最大，第2层明层西南面柱子发生了严重的变形，现已对西南面内柱进行了加固。

图14 应县木塔柱结构纵向倾斜分析示意图

3.4 木塔整体倾斜分析

从整体模型中截出木塔的纵截面，在截面中将木塔第1层和第5层两侧的点云拟合出边缘线，连接两侧边缘线中点作为假定中心线，取中心线中点作为该层的中心点，连接2层中心点得到木塔的倾斜方向线，同时过底层中心做垂线，倾斜方向线和垂线的夹角即为倾斜角度。图15为东西方向和南北方向对角线的截面，可以看出木塔整体向东倾斜0.35°，向北倾斜0.54°。由此分析应县木塔整体的倾斜方向为东北方向，和单柱以及单层的倾斜方向一致。

图15 应县木塔整体倾斜示意图

3.5 木塔整体扭转分析

通过结构剖切，可以快速精确得出应县木塔的扭转姿态[22]。在木塔点云模型中沿着每一层最外沿进行剖切，拟合出每层的外侧边缘，连接对角点做出对角线，叠加各层数据，如图16所示。由于木塔每层相同位置的对角线不重合，因此木塔存在扭转变形。相对底层来说，各层对角线的位置和相对位移不同，表明木塔发生的扭转是复杂的，如图16箭头所指，木塔整体扭转趋势为西侧由南向北发生顺时针扭转，东侧由南向北发生逆时针扭转，并且上下2层的相对扭转程度各不相同。

图16 应县木塔扭转示意图

4 结语

(1)通过激光雷达扫描技术获取应县木塔内外部的点云数据，利用全站仪控制测量得到标靶控制点坐标，进行点云配准及坐标转换等处理后得到绝对坐标系下的整体点云模型。对应县木塔各层柱子的偏移量以及偏移角度进行了测量，由测量数据定量分析了目前的形变情况。分析得出应县木塔所有柱子都存在变形，第2层柱子的倾斜角度和位移都大于其他各层，其中M2W23号柱子形变最严重；单层来看西南面柱子偏移程度比其他面严重，各层向东北方向存在不同程度的倾斜；纵向看应县木塔二、三层的倾斜程度比四、五层要大；在纵截面中测量整体的倾斜角度，同样表明木塔的倾斜方向为东北方向；同时木塔存在复杂的扭转变形，扭转趋势为西侧由南向北顺时针扭转，东侧由南向北逆时针扭转。依据测量数据分析了应县木塔的形变，但没有分析各柱、各层的承载力大小和受力情况，还需进一步研究。

(2)通过三维激光雷达技术可以采集到相对完整的高精度点云数据，构建出精确的应县木塔点云模型，对点云模型测量可得出倾斜角度及位移量，方便对木塔的整体姿态、偏移、倾斜、扭转等情况进行全面分析，对研究应县木塔的形变有重要意义，本文的技术方案对研究中国其他物质文化遗产同样具有应用价值。