三维激光扫描是一种重要的隧道测量手段，具有效率高、数据量大、精度高等优点。对于扫描所得点云的计算与处理，如何通过点云数据来体现隧道的变化，已有许多学者进行相关的研究。刘燕萍等^[1]使用了站扫方式提取出断面的圆心和半径，计算出隧道的收敛变形情况和预估变化趋势。刘绍堂等^[2]提出了用基于三维激光扫描技术的多点整体监测方法，引入了椭圆拟合法描述隧道的整体变形。陈勇^[3]通过分析野外轨道点云数据，重构了既有线曲线参数，对现有轨道进行了整正计算。许磊等^[4]提出了一种基于射线相交的隧道断面自动测量方法，满足多种结构类型的隧道断面测量要求。王令文等^[5]通过采用全站仪对点云数据进行绝对定位来提高拼接精度，以设计轴线为基准提取隧道断面，分析了隧道的断面收敛情况。

隧道三维激光扫描变形测量大多用于计算隧道直径收敛，研究隧道内壁断面形状是否发生了变形^[6]，未见用此方法测量分析轨道的变形。地铁点云数据密度大、精度高，可提取任意里程处轨道面的几何信息，本文利用隧道扫描点云提取出的轨道面高程数据，给出了局部相对变形波动值的概念，提出了一种基于小波变换和维格纳-威尔分布的变形提取方法，结合理论分析、模拟数据和实例计算验证了该方法的有效性与可行性。

1 线状结构局部相对变形波动指标

线状结构的地铁变形通常采用周期性重复观测的方法，获取特征处的多期坐标，求取周期间的变形量，以期掌握地铁结构的整体变形情况。由于地铁隧道内基准的建立和维护困难，作业效率低下，掌握地铁全线隧道变形的代价很高，受地铁隧道线状分布的影响，其横向和垂向传递误差较大，绝对精度较低。从线状横向和垂向变形对隧道结构的危害程度来看，分布在长区段相对变化平缓的变形，即使累积变形量较大，结构损伤的可能性也较小；而局部较为剧烈的相对变形对线状地铁结构的影响更大，因此，对多期测量成果在局部的相对变形进行分析和研究更有实际意义。

为提取两期测量成果中隐含的线状结构相对波动信息，本文参考了高铁轨道不平顺性检测原理与方法^[7]，定义了线状结构局部相对变形波动性指标，如图 1所示，实线为首次测量成果，虚线为第$k$次测量成果，$\mathrm{2}s$为检测弦长。

图  1  线状结构局部相对变形波动性指标示意图

Figure  1.  Schematic of Linear Structure Partial Relative Deformation Fluctuation Index

下载: 全尺寸图片 幻灯片

定义线状结构局部相对变形波动值为：

式中，$d_i^{\mathrm{0}}$为首次测量成果中里程$i$测点到里程i－s与里程$i+s$测点连线的垂直距离；$d_i^k$为第$k$次测量的相应距离；${\delta }_i$表示首期与第$k$期之间的相对变形波动值。首期与第$k$期之间的垂直距离作差表示累积变形，相邻两期之间作差表示本期变形。

线状结构局部相对变形波动值的实质为两期检测弦长中点偏距之差，形式上近似于两期中点弦偏距变化值对里程的导数，与两期中点弦偏距变化曲线的导数密切相关。

2 相对变形波动提取与定位

2.1   小波变换与小波分解

小波变换利用系数和规范小波基的积作为基础进行信号分析，从而可以在时、频两域表征信号的局部特性^[8-9]。将满足容许条件^[10]的基本小波函数$\psi \left(t\right)$进行伸缩和平移得到连续小波函数族$\left\{{\psi }_{a,b}(t)\right\}$^[11]。将连续小波与信号$f(t)$进行内积，得到信号的连续小波变换：

式中，${\psi }^{*}$表示$\psi$的共轭函数；$W_f(a,b)$为小波系数^[10]，其中，$a$是信号频率相关的伸缩（尺度）因子，$b$是与信号时间相关的平移因子。

用Mallet算法^[9]将原始信号$f(t)$分解为小波系数低频部分$A_{\mathrm{1}}$和小波系数高频部分$D_{\mathrm{1}}$，再对$A_{\mathrm{1}}$进行分解得到$A_{\mathrm{2}}$和$D_{\mathrm{2}}$，逐层分解最终得到$A_n$和$D_n$。对$A_i$和$D_i$进行重构，从而得到原始信号在第$i$层的近似信号$a_i$和细节信号$d_i$^[12]。

2.2   维格纳-威尔分布

信号$f(t)$的能量在时域和小波域是相同的^[13]，$|W_f\left(a,b\right){|}^{\mathrm{2}}$可以看作关于尺度$a$和平移$b$的二元能量分布。维格纳-威尔分布（Wigner-Ville distribution，WVD）用于描述信号能量在时频两域的分布，并且具有平移不变性、不损失信号幅值和相位信息等优点^[14]，比$|W_f\left(a,b\right){|}^{\mathrm{2}}$更适合能量型时-频联合小波分析。对信号$f(t)$的瞬时自相关函数进行傅里叶变换，得到了WVD的一般表达式为：

式中，$f^{*}$表示$f$的共轭函数；t表示时间；ω表示频率；τ表示时间延迟。本文计算方法采用加窗后的平滑伪维格纳-威尔分布（smoothed pseudo Wigner-Ville distribution，SPWVD）。

2.3   相对变形波形提取

如果原始信号中存在异常相对变形波形，则在该波形的特征频率下，位于异常波形处的小波能量将产生突变。由于轨道不平顺幅值沿里程基本满足正态分布^[14]，根据拉依达准则（$\mathrm{3}\sigma$准则）设定能量阈值用于确定相对变形波形发生的位置。具体提取算法分为确定变形波长和确定变形位置两大步骤。

2.3.1   确定相对变形波长

使用小波分解对原始变形波动值进行处理，采用的小波基函数由原始数据的特性决定，由数据量和小波基函数共同确定分解的层数。（1）对原始信号进行小波分解，得到不同层数下的近似信号$a_i$和细节信号$d_i$；（2）对所有$a_i$和$d_i$进行快速傅里叶变换（fast Fourier transform，FFT），分析后确定特征波长。在$a$中确定特征波长后，选择对该波长分辨率最大的$d$作为下一步分析的原始数据。

2.3.2   确定相对变形位置

和原始信号相比，经过小波分解再重构后的近似信号和细节信号成分更加单一，能量分布更为集中，可以更好地在SPWVD中提取相对变形波形的位置。提取步骤为：（1）对确定的细节信号$d$进行SPWVD计算，得到频率-里程-能量三维分布图；（2）在频率轴中找到特征波长所对应的频率，提取出该频率下不同里程处的SPWVD能量值$P\left(i\right)$，即特征频率对应的里程-能量分布图；（3）从里程的第$L$米开始，根据$\mathrm{3}\sigma$准则设定检测阈值$U(i)$，$U(i)$是由从$P(\mathrm{0})$到P（i-1）的均值$E_{P(\mathrm{0},i-\mathrm{1})}$与标准差${\sigma }_{P\left(\mathrm{0},i-\mathrm{1}\right)}$决定的，计算公式为：

将检测阈值$U(i)$和SPWVD能量值$P(i)$绘制到同一张图中，能量值满足$P\left(i\right)>U(i)$的波形最高点对应的里程数即为相对变形波形的近似位置。

3 算法验证与实例分析

3.1   模拟数据计算与分析

模拟高程值与计算相对变形波动值沿里程的分布图如图 2所示，其中，$h_{\mathrm{1}}$、$h_{\mathrm{2}}$为模拟的两期轨道中心点高程值。高程值采样间隔为0.5 m，中点弦的弦长设为10 m，分别加入了中误差为0.33 mm、幅值为1 mm的正态随机噪声。两期模拟高程值并不重合，在最高点处相差5 cm，但相对关系保持相近。在$h_{\mathrm{2}}$中加入3个沉降幅值为1 mm、波长为10 m的标准下凹全波余弦型沉降曲线模型，沉降最大处的位置分别为140 m、205 m和340 m。

图  2  模拟高程值与计算相对变形波动值

Figure  2.  Simulation Elevation Data and Calculated Deformation Fluctuations

下载: 全尺寸图片 幻灯片

采用db6（Daubechies6）作为小波基函数，对原始变形波动值进行小波分解。经过查阅文献和实际试验，分解的层数定为6层。各层近似信号的波长-FFT频谱图如图 3所示，细节信号频谱图如图 4所示，标出每张图幅值最大点的横坐标，即特征波长，由于$a_{\mathrm{5}}$与$a_{\mathrm{6}}$幅值最大点对应的特征频率为0，波长为无穷大，所以不标出峰值。

图  3  近似信号$a$的FFT频谱

Figure  3.  Amplitude Spectrum of Approximate Signal $a$

下载: 全尺寸图片 幻灯片

图  4  细节信号$d$的FFT频谱

Figure  4.  Amplitude Spectrum of Detail Signal $d$

下载: 全尺寸图片 幻灯片

分析FFT频谱图可得：（1）特征波长10.427 m在近似信号$a_{\mathrm{1}}$、$a_{\mathrm{2}}$、$a_{\mathrm{3}}$中都有最大幅值，并且与中点弦的弦长接近，因此，将其定为待提取相对变形波形的特征波长；（2）特征波长10.427 m在各层细节信号$d$中都有对应的波峰，$d_{\mathrm{4}}$中最为显著，可知$d_{\mathrm{4}}$对于待提取特征波长10.427 m，具有最大的时空分辨率，换算成特征频率为0.096 /m，将$d_{\mathrm{4}}$定为下一步进行波形位置确定的数据。

计算$d_{\mathrm{4}}$的SPWVD得到能量三维分布图见图 5，反映了能量在里程和频率两域的分布。

图  5  细节信号$d_{\mathrm{4}}$的SPWVD能量分布

Figure  5.  SPWVD Energy Distribution of the Reconstructed Detail Signal $d_{\mathrm{4}}$

下载: 全尺寸图片 幻灯片

提取出特征波长10.427 m对应的不同里程处的SPWVD能量值$P(i)$，如图 6中实线所示。$L$取50 m，设定检测阈值$U(i)$，将$U(i)$和$P(i)$绘制到图 6中，满足$P\left(i\right)>U(i)$的波形最高点对应的里程数即为相对变形波形的近似位置，可知存在3处波形，位置分别位于里程136.25 m、200.25 m、336.75 m处附近。

图  6  特征波长对应的里程-SPWVD分布曲线

Figure  6.  Mileage-SPWVD Energy Distribution of Characteristic Wavelengths

下载: 全尺寸图片 幻灯片

从原始变形波动值中的相应位置提取出的波形形状非常不规律，这是由于原始变形波动值含有许多高频噪声，不利于幅值的确定。多层近似信号可以看作是原始数据逐步去噪、简化后得到的信号，在近似信号中，提取相对变形波形可以得到更好的结果。

图 7是从第2层近似信号$a_{\mathrm{2}}$的相应位置提取出的局部相对变形波形，标注出沉降最大处的里程和幅值。这些波形和标准沉降曲线波形相似，幅值分别为0.91 mm、0.94 mm、0.98 mm，与预设值1 mm相近。由图 7可知，提取的相对变形波形的里程范围为136.25~146.68 m、200.25~210.43 m、336.75~347.18 m，与预设相对变形波形位置接近。将提取局部相对变形波形信息与预设值汇总到表 1中，可以发现，通过本文算法检测得到的结果与预设值接近。

图  7  从$a_{\mathrm{2}}$中提取的相对变形波形

Figure  7.  Extraction Results from the 2 Reconstructed Approximate Signal $a_{\mathrm{2}}$

下载: 全尺寸图片 幻灯片

表  1  提取波形与预设波形数据对比

Table  1.  Comparison of Extraction Results and Default Values

数据来源	波形数量	波长/m	里程范围/m			幅值/mm
算法提取	3	10.427	136.25~146.68	200.25~210.43	336.75~347.18	0.91	0.94	0.98
预设值	3	10	135~145	200~210	335~345	1	1	1

下载: 导出CSV 

| 显示表格

3.2   实测数据计算与分析

采用以下两个实测数据进行计算与分析：（1）将轨道设计值加上随机误差作为首期观测值，实际扫描的数据作为当前期值；（2）同一轨道短时间内扫描两次所得点云。

3.2.1   实测数据1

数据来源于武汉地铁某区间三维激光扫描点云，使用的仪器为徕卡ScanStation P30，精度指标为测距误差1.2 mm+10 mm/km，角度误差‍8″，每隔0.5 m截取一个断面^[15]，从断面点云中提取出左右轨面处的点云，取平均值得轨道中心点高程，截取里程为42 613.656~43 113.156 m，总长为499.5 m。在高程设计值中加入中误差为0.85 mm的随机噪声。对高程的测量值与设计值进行上述相对变形波形提取算法的处理，将中点弦的弦长分别设为5 m、10 m和15 m，分别计算不同中点弦的弦长对应的相对变形波形。

对10 m弦长变形波动值进行小波分解、重构和FFT分析，近似信号$a$的FFT结果如图 8所示，波长9.615 m表现最显著，在$a_{\mathrm{1}}$、$a_{\mathrm{2}}$、$a_{\mathrm{3}}$中均表现为最大幅值，因此，定特征波长为9.615 m。

图  8  近似信号$a$的FFT频谱

Figure  8.  Amplitude Spectrum of Approximate Signal $a$

下载: 全尺寸图片 幻灯片

特征波长在第4层细节信号$d_{\mathrm{4}}$中幅值最大，具有最高的分辨率，故选择$d_{\mathrm{4}}$计算其SPWVD。提取特征波长9.615 m所对应的里程-能量分布曲线，设置阈值，提取能量分布曲线超过阈值的波峰顶点对应里程，即为相对变形位置，如图 9所示。由图 9可知，存在4个相对变形波形，位置分别为42 700.656 m、42 940.656 m、43 024.906 m、43 054.906 m。在近似信号$a_{\mathrm{2}}$的相应位置提取4个相对变形波形，将波形的位置、幅值、范围等信息汇总如表 2所示。

图  9  相对变形波形位置的提取

Figure  9.  Extraction Locations of Abnormal Deformation

下载: 全尺寸图片 幻灯片

表  2  10 m弦长下的相对变形波形

Table  2.  Information of Abnormal Deformation Waves in 10 m Midpoint Chord

近似位置/m	幅值/mm	里程范围/m	形变最大位置/m
42 700.656	2.99	42 696.156~42 705.156	42 700.156
42 940.656	3.18	42 936.656~42 946.156	42 941.656
43 024.906	2.57	43 024.650~43 034.156	43 030.156
43 054.906	5.55	43 051.656~43 062.156	43 057.156

下载: 导出CSV 

| 显示表格

使用相同方法对弦长为5 m和15 m的情况提取相对变形，结果如表 3所示。该算法可以提取出波长与事前给定的中点弦的弦长近似的相对变形波形，不同弦长所提取出的波形互不相同。

表  3  5 m和15 m弦长下的相对变形波形

Table  3.  Information of Abnormal Deformation Waves in 5 m and 15 m Midpoint Chord

弦长/m	近似位置/m	幅值/mm	里程范围/m	形变最大位置/m
5	43 008.656	2.79	43 006.656~43 011.156	43 008.656
5	43 046.906	1.92	43 042.656~43 047.156	43 044.656
15	43 011.406	3.30	42 997.156~43 011.656	43 004.156

下载: 导出CSV 

| 显示表格

3.2.2   实测数据2

数据2为武汉地铁某区段同一晚上推行两趟获取的数据，小车上使用的扫描仪为徕卡ScanStation P30，截取里程范围为13 633.36~14 134.36 m，总长为501.00 m。将中点弦的弦长分别设为8 m和20 m，分别提取出变形波形。对于8 m弦长，对变形波动值进行小波分解、重构和FFT分析，结果显示波长8.358 m表现最为显著，将此波长定为特征波长。特征波长在$d_{\mathrm{4}}$中幅值最大，具有最高的分辨率，提取$d_{\mathrm{4}}$的SPWVD中特征波长8.358 m所对应的能量值，组成里程-能量分布曲线，设置阈值，结果如图 10所示，由图 10可知，原始相对变形波动波形中不存在波长在8 m附近的变形波形。

图  10  8 m中点弦下的相对变形波形的位置提取

Figure  10.  Extraction Locations of Abnormal Deformation Under 8 m Midpoint Chord

下载: 全尺寸图片 幻灯片

图 11为中点弦的弦长为20 m（特征波长21.804 m）的SPWVD-阈值图，从图 11中不能提取出相应波长的相对变形波动波形。

图  11  20 m中点弦下的相对变形波形的位置提取

Figure  11.  Extraction Locations of Abnormal Deformation Under 20 m Midpoint Chord

下载: 全尺寸图片 幻灯片

4 结语

本文在定义线状结构局部相对变形波动值的基础上，对地铁三维激光扫描点云中提取的轨道数据进行小波分析与SPWVD分析，可以识别相应波长的相对变形波形，并确定其影响范围和最大变形幅值。对原始的变形波动值进行多层小波分解与重构，使信号中的待提取成分表现更为显著，分辨率更高，从而确定相对变形波形的波长。SPWVD具有明确的物理意义，可以用于能量检测，从而确定相对变形波形的位置。实例分析的研究结果表明，本文提出的相对变形提取算法具有一定的有效性和准确性，可以提取波长与中点弦的弦长接近的相对变形波形，能够达到较好的相对变形波形信息检测的效果。

后续研究将采用此方法分析地铁隧道结构侧壁的横向水平和顶壁的垂向线状相对变形，期望从多期地铁扫描数据中检测出更多有效变形信息。