基于LiDAR点云的多分裂导线自动提取与精细建模方法与流程

文档序号:12064791阅读:722来源:国知局
基于LiDAR点云的多分裂导线自动提取与精细建模方法与流程

本发明属于激光雷达点云数据信息提取技术领域,涉及一种基于LiDAR点云的多分裂导线自动提取与精细建模方法。



背景技术:

电力工业是国民经济基础工业之一,是国家的重要支柱产业。随着我国经济高速发展,超高压大容量输电线路越建越多,线路走廊穿越地理环境更复杂,给线路维护带来很多困难。

机载激光扫描(Light Detection and Ranging,LiDAR)技术作为近年来快速发展的一项新技术,可以快速获取高精度的三维信息。该技术具有全天候工作等优点,可以弥补传统航空摄影测量难以快速进行电力线测量的缺憾,提高电力巡线的效率,减少输电事故的发生。因此,研究机载激光扫描数据中电力线的自动提取和精细重建技术对电力巡线工作有着非常重要的现实意义。

在高压、超高压、特高压输电线路架设中,为抑制电晕放电和线路电抗,一般采用多分裂导线进行架设,且每隔一段距离采用间隔棒对分裂导线进行固定。但现有文献中的LiDAR电力线建模研究绝大部分局限于非分裂导线,或者把多分裂导线当成单导线来进行处理,且现有方法只适用于数据质量较好的情况,对于存在较多噪声的点云数据,提取与重建结果并不理想。例如,澳大利亚的McLaughlin根据LiDAR点云的维数特征自动提取电力线,通过点的聚集方法,得到同一条导线上的点数据,再通过悬链线方程拟合,但该方法在检测电力线点的过程中易受噪声的影响,提取的各条电力线有间断,并不连续;余洁采用滤波的方法滤除地面和植被点,采用二维Hough变换分离各条电力线,根据双曲余弦函数拟合单条电力线,但树木和电力线混合区域的提取结果会相对变差,当电力线出现垂直排列时,用二维Hough变换很难检测出多条电力线;陈驰采用同一条电力线的局部相邻点的高程比较接近这一原则,统计每条电力线在XOY投影面上的邻近点,判断其高程差是否超过一定的阈值的方法解决重投影的问题,但对于多分裂导线仅利用高程值特征很难区分每一根分裂导线。

在多分裂导线提取方面,王成、夏少波首次提出一种多分裂导线全自动三维重建方法,采用空间聚类识别单股电力线点云,然后分段进行主成分分析并投影到与主方向垂直的平面上,对第一段数据采用全排列的方法进行k-means聚类确定分裂数,然后对每段进行聚类,最后基于悬链线进行最小二乘拟合。该方法处理过程复杂且只适用于数据质量较好的情况,对于存在较多噪声的点云数据,提取与重建结果并不理想。主要原因是该算法对噪声敏感,抗噪性不强,提取的电力线不连续、不完整,拟合精度的不够高,且当分裂导线间存在间隔棒时,会出现聚类误差过大甚至错误。

因此,随着机载LiDAR电力巡线的推广应用,电力部门急需一种能够克服连续噪声影响的多分裂电力线提取与精细重建的稳健方法。



技术实现要素:

本发明目的在于提供一种更稳健的基于LiDAR点云的多分裂导线自动提取与精细建模方法。

本发明的技术方案如下:

一、基于LiDAR点云的单相电力线自动提取方法,包括:

步骤1,从LiDAR点云中提取电力线和电力塔点云;

步骤2,根据电力线点在局部区域内高程基本相同、而电力塔点在局部区域内高程变化大的特点,从电力线和电力塔点云中进一步提取电力塔点云,剔除电力塔点云得电力线点云;

步骤3,对电力线点云进行空间分割,得到各相电力线点云;采用随机一致性检测法检测各相电力线点云中的噪声点,并剔除。

步骤1进一步包括:

1.1基于点云密度和地形坡度从LiDAR点云中滤除地面点,获得非地面点云;

1.2基于点云密度和高差从非地面点云中滤除植被点;

1.3基于三维空间k-d树在滤除了植被点的非地面点云中搜索剩余植被点,并滤除;

1.4基于种子填充法从子步骤1.3所得的剩余点云中提取电力线和电力塔点云。

子步骤1.4具体为:

(1)将剩余点云记为数据集B,对数据集B建立k-d树,初始化计数label为0;

(2)扫描数据集B,若当前点B(i)未被处理,把B(i)加入栈Q,然后,执行如下步骤:

(2a)栈Q中弹出一未处理点作为种子点q,标记为已处理并存入点集C;在数据集B中搜索种子点q的n’邻近点,得n’邻近点集N,将n’邻近点集N中未被处理、且与种子点q的距离和高差均小于对应的第二距离阈值Td2和第二高差阈值H2的点加入栈Q;第二距离阈值Td2和第二高差阈值H2为与点云质量和点云密度有关的经验值;

(2b)重复子步骤(2a)直至栈Q为空,此时,点集C中点即构成一个连通区域,该连通区域内点标记为label;

(2c)获取该连通区域的矩形范围,计算对角线长度,若对角线长度满足长度阈值,则该连通区域为电力线和电力塔点云,长度阈值根据所处理电力线的长度进行取值;

(3)令label+1,重复步骤(2),直至数据集B中所有点均被处理。

步骤2进一步包括:

2.1根据电力线点在局部区域内高程基本相同、而电力塔点在局部区域内高程变化大的特点,通过对电力线和电力塔点云进行张量分析,提取潜在电力塔点云;

2.2基于种子填充法对潜在电力塔点云进行空间分割,获取电力塔点云,剔除电力塔点云得电力线点云。

子步骤2.1具体为:

(1)对电力线和电力塔点云建立k-d树;

(2)对k-d树中各数据点p分别获取以p为中心的m邻近点,对p进行张量分析,获得mxx、myy、mzz分别用来描述m邻近点在X方向、Y方向、Z方向的拟合误差;其中,(xj,yj,zj)为m邻近点中第j个邻近点的三维坐标,(xp,yp,zp)为数据点p的三维坐标;

(3)根据①电力线在Z方向的拟合误差小、在X和Y方向的拟合误差大,以及②电力塔在Z方向的拟合误差大、在X和Y方向的拟合误差小的特点,根据各数据点p的mxx、myy、mzz值提取潜在电力塔点。

子步骤2.2中,基于种子填充法对潜在电力塔点云进行空间分割,获取电力塔点云,具体为:

采用种子填充法对潜在电力塔点云进行空间分割,获取多个连通区域,高差大于第三高差阈值H3的连通区域则认为是电力塔连通区域;

对各电力塔连通区域分别执行:获取电力塔连通区域的中心坐标,在步骤1所提取的电力线和电力塔点云中,以中心坐标为圆形,预设半径的圆形范围内点云即电力塔点云;

上述,第三高差阈值H3根据电力塔长度进行取值,预设半径设为电力塔长度的一半。

步骤3中,采用随机一致性检测法检测各相电力线点云中的噪声点,具体为:

对各相电力线点云分别进行如下步骤:

(1)从该相电力线点云中随机产生两个点;

(2)计算两点在XOY水平面上的直线方程,统计到该直线距离小于第四距离阈值Td4的点数;

(3)重复步骤(1)~(2),直至达到预设的迭代次数;

(4)取统计点数最多的直线方程为该相电力线的平面直线方程,该相电力线点云中到该平面直线方程距离小于第四距离阈值Td4的所有点均为该相电力线点,其他点为噪声点;

上述,第四距离阈值Td4为经验值,其取值要确保区分出电力线点云中任意两条电力线。

二、基于LiDAR点云的多分裂导线自动提取与精细建模方法,包括:

对前一技术方案所得各相电力线点云分别进行如下:

(1)对该相电力线点云进行误差分析,判断电力线类型,若该相电力线点云为单相地线点云,执行步骤(2);若该相电力线点云为单相导线点云,执行步骤(3);

(2)对单相地线点云进行悬链线拟合并输出;

(3)基于二分法从单相导线点云中提取分裂子导线点云,并对分裂子导线点云进行悬链线拟合,本步骤进一步包括:

3.1将单相导线点云投影到XOZ竖直面,将投影等分为若干导线段,对各导线段分别进行线性拟合,得到各导线段的直线方程Ax+Bz+C=0;

3.2对各导线段上所有点(xk,zk)分别计算所有导线段上d>DIS的点构成上分裂子导线点云S(0),所有导线段上d<-DIS的点构成下分裂子导线点云S(2),所有导线段上-DIS≤d≤DIS的点构成中分裂子导线点S(1);DIS为与多分裂导线规格相关的距离阈值,其值不小于0;

3.3分别判断各分裂子导线点云S(s)中点数是否大于点数阈值NUM,s=0、1、2,点数阈值NUM用来判断粗差点,根据经验设置;剔除点数不大于NUM的分裂子导线点云S(s),对点数小于NUM的分裂子导线点云S(s)执行步骤3.4;

3.4对分裂子导线点云进行悬链线拟合,计算拟合残差,输出拟合残差小于残差阈值的分裂子导线点云拟合结果;对拟合残差不小于残差阈值的分裂子导线点云分别执行步骤3.5;残差阈值为与点云质量有关的经验值;

3.5将分裂子导线点云投影到XOY水平面,对投影进行线性拟合,得到直线方程ax+by+c=0,对分裂子导线点云上所有点(x'k,y'k)分别计算分裂子导线点云上t>DIS的点构成左分裂子导线点云T(0),分裂子导线点云上t<-DIS的点构成右分裂子导线点云T(2),分裂子导线点云上-DIS≤t≤DIS的点构成中分裂子导线点云T(1);

3.6分别判断各分裂子导线点云T(u)中点数是否大于点数阈值NUM,u=0、1、2;剔除点数不大于NUM的分裂子导线点云T(u);对点数小于NUM的分裂子导线点云T(u)进行悬链线拟合,并输出拟合结果。

步骤(1)具体为:

对该相电力线点云进行拟合,计算拟合残差,若拟合残差小于残差阈值,则该相电力线点云为地线点云;否则,为多分裂的导线点云;同时,计算拟合残差过程中,将误差大于三倍残差阈值的点作为粗差点剔除;残差阈值为与点云质量有关的经验值。

和现有技术相比,本发明具有如下优点和有益效果:

本发明具有抗噪性强、拟合精度高、人工干预少等特点,提高了三维巡线系统的效率与精度,可获得更精准的三维坐标,降低电力线走廊巡线成本,具有一定的科研意义和重要的使用价值。

附图说明

图1为本发明的整体流程图;

图2为自动提取电力线点云和电力塔点云的流程图;

图3为电力塔点云的提取与定位的流程图;

图4为单相导线或地线的提取与拟合的流程图;

图5为多分裂导线的自动I去与精细建模的流程图。

具体实施方式

下面结合附图对本发明做进一步详细描述,本发明整体流程见图1,概括起来,本发明包括四个步骤:

一、从LiDAR点云中自动提取电力线和电力塔点云。

本步骤流程见图2。根据点云的分布特性,可基于点云密度和地形坡度从LiDAR点云中滤除地面点,基于点云密度和高差从LiDAR点云中初步分离出植被点和非植被点,然后在三维空间基于k-d树搜索植被同类点。其中,密度和高差阈值根据直方图统计结果进行设置。基于种子填充法进行点云分割得到多个连通区域,根据连通区域长度提取电力线和电力塔点云。

1.1基于点云密度和地形坡度从LiDAR点云中滤除地面点,获得非地面点云。

对原始LiDAR点云在XOY水平面建立规则矩形格网,规则矩形格网大小记为gridsize,找到各格网的最低点,并计算各格网的点云密度。若格网的点云密度大于预设的密度阈值S1,则对格网内所有点,计算其与最低点的坡度,坡度小于预设的坡度阈值T1的点则为地面点。密度阈值S1的设值与所采用扫描仪的频率有关,一般设为扫描仪的平均点云密度。坡度阈值T1的设置与当地地形有关,当该区域地势变化较大时,可减小规则矩形格网大小gridsize或增大坡度阈值T1;当该区域地势较平坦时,可增大规则矩形格网大小gridsize或减小坡度阈值T1。本实施例中,坡度阈值T1取0.3。

1.2基于点云密度和高差从非地面点云中滤除植被点。

对非地面点云,找到各格网的最低点,计算各格网点云密度,若点云密度大于预设的密度阈值S1,则该格网中与最低点高差小于预设的第一高差阈值H1的所有点为植被点。第一高差阈值H1的设值与植被平均高度有关,当该区域植被较高大时,可通过增大第一高差阈值H1来尽可能多的提取植被点;当该区域大多为低矮灌木时,可通过减小第一高差阈值H1以避免过多的剔除电力线或电力塔点。本实施例中,第一高差阈值H1取15。

1.3在三维空间基于k-d树搜索滤除了植被点的非地面点云中的剩余植被点,并滤除剩余植被点。

对已经分类的植被点按三维坐标建立k-d树,在k-d树中搜索未分类的非地面点云中各点的n邻近点,若该未分类点到n邻近点的平均距离小于预设的第一距离阈值Td1,则认为该未分类点是剩余植被点。这里,未分类的非地面点云即子步骤1.2所得的滤除了植被点的非地面点云。n和第一距离阈值Td为经验值,本实施例中,n取3,第一距离阈值Td也取3。n取值越大或Td1取值越小,则说明判断为剩余植被点的条件越苛刻,提取的剩余植被点也越少;反之,则说明判断为剩余植被点的条件越不苛刻,提取的剩余植被点也越多。

1.4基于种子填充法从子步骤1.3获得的剩余点云中提取电力线和电力塔点云。

考虑到电力线和电力塔点云呈线状分布,通常贯穿整个区域,具有较强的延伸性,而噪声点通常是小范围分布,数量较少,具有不规则性。因此,基于种子填充法对剩余点云进行空间分割,获得连通区域,再根据各连通区域长度提取电力线和电力塔点云。

本子步骤的具体实施过程如下:

(1)将剩余点云记为数据集B,对数据集B建立k-d树,初始化计数label为0;

(2)扫描数据集B,若当前点B(i)未被处理,把B(i)加入栈Q,然后,执行如下步骤:

(2a)栈Q中弹出一未处理点作为种子点q,标记为已处理并存入点集C;在数据集B中搜索种子点q的n’邻近点,得n’邻近点集N,将n’邻近点集N中未被处理、且与种子点q的距离和高差均小于对应的第二距离阈值Td2和第二高差阈值H2的点加入栈Q。

第二距离阈值Td2和第二高差阈值H2为与点云质量和点云密度有关的经验值,并且在设值时要尽量保证提取的电力线和电力塔点云的完整性和连续性。当点云质量较好,点云密度较大时,可适当减小第二距离阈值Td2和第二高差阈值H2;当点云质量较差,点云密度较小时,可适当增大第二距离阈值Td2和第二高差阈值H2。本实施例中,第二距离阈值Td2取5,第二高差阈值H2取3。

(2b)重复子步骤(2a)直至栈Q为空,此时,点集C中点即构成一个连通区域,该连通区域内点标记为label。

(2c)获取该连通区域的矩形范围,计算对角线长度,若对角线长度满足长度阈值,则认为该连通区域是电力线和电力塔点云。长度阈值根据所处理电力线的长度进行取值,其取值一般应大于100。

(3)令label+1,重复步骤(2),直至数据集B中所有点均被处理,这样就可以得到剩余点云中的电力线点云和电力塔点云。

二、电力塔点云的提取与定位

基于上述流程可以从原始点云中提取出整个线路上的电力线和电力塔点云,而从中提取出每相邻两基塔间的导线段是电力悬线识别和重建的基础。本发明根据电力塔坐标对电力线点云进行分段,而电力塔的坐标是从点云中计算得到。本方法流程见图3。

2.1对电力线和电力塔点云进行张量分析,提取潜在电力塔点云。

根据电力线点在局部区域内的高程基本相同,而电力塔点在局部区域内高程变化大的特点,采用张量分析的方法从步骤一提取的电力线和电力塔点云中进一步提取潜在电力塔点云。

对步骤一提取的电力线和电力塔点云建立k-d树,对k-d树中各数据点p,分别计算以p为中心的m个邻近点的对称张量矩阵T,如下:

式(1)中:

上述,(xj,yj,zj)为m个邻近点中第j个邻近点的三维坐标,(xp,yp,zp)为数据点p的三维坐标。

m根据需求进行取值,一般m取值大,精度高,但计算量也大。一般情况,m取20即可基本反映邻域特征。

对称张量矩阵T通过拟合以p为中心的m邻域点的误差椭球对p点的几何特性进行描述,其中,mxx描述邻域点在X方向的拟合误差,myy描述邻域点在Y方向的拟合误差,mzz描述邻域点在Z方向的拟合误差。根据电力线在Z方向拟合误差较小、在X和Y方向上拟合误差大,而电力塔在Z方向拟合误差较大、在X和Y方向上拟合误差小的特点,即可提取潜在电力塔点。提取潜在电力塔点时,只需计算各数据点的mxx、myy和mzz,mzz大于mxx和myy的数据点即潜在电力塔点。

2.2对潜在电力塔点云进行空间分割获取电力塔点云,剔除电力塔点云得电力线点云。

采用种子填充法以第三距离阈值Td3对潜在电力塔点云进行空间分割,得到多个可能的电力塔连通区域,高差大于第三高差阈值H3的连通区域则认为是电力塔连通区域。空间分割即从潜在电力塔点云中分割出高度长度与电力塔匹配的连通区域,所以第三距离阈值Td3根据电力塔长度进行取值,第三高差阈值H3根据电力塔高度进行取值。本实施例中,第三距离阈值Td3取为10米,第三高差阈值H3取为5米。

空间分割所提取的电力塔点云并不完整,还需要根据电力塔平面中心坐标提取完整的电力塔点云。

根据电力塔平面中心坐标提取完整的电力塔点云具体为:

获取各电力塔连通区域的数据范围参数,数据范围参数包括(Xmin,Ymin)和(Xmax,Ymax),Xmin和Xmax分别表示电力塔连通区域中点横坐标的最小值和最大值,Ymin和Ymax分别表示电力塔连通区域中点纵坐标的最小值和最大值,计算电力塔平面中心的坐标(X,Y):

在步骤一所得的电力线点电力塔点云中,以电力塔平面中心坐标(X,Y)为圆心,预设半径的圆形范围内点云即电力塔点云。预设半径一般设为电力塔最大长度的一半,本实施例中,预设半径取8米。

对各电力塔连通区域分别执行上述步骤,即可获得完整的电力塔点云。

三、单相电力线的提取与拟合。

采用种子填充算法对电力线点云进行空间分割,得到各相电力线点云。本步骤在进行空间分割时,所设置的第四距离阈值Td4和第四高差阈值H4要确保能区分电力线点云中任意两条电力线,一般情况下,第四距离阈值Td4取5,第四高差阈值H4取3。当电力线附近存在连续噪声时,采用空间分割的方法会将噪声错分为电力线点,因此需要将噪声点剔除以减小拟合误差,采用二维RANSAC对各相电力线点云进行直线检测以剔除噪声点,流程见图4。

随机一致性检测(RANSAC)是用迭代的方法从包含大量异常点的数据中确定数学模型的最优参数,所以也可以用来检测异常点的存在。本步骤采用随机一致性检测来检测噪声点。

本步骤具体为:

对各相电力线点云分别进行如下步骤:

(1)从该相电力线点云中随机产生两个点;

(2)计算两点在XOY水平面上的直线方程,统计到该直线距离小于第四距离阈值Td4的点数;

(3)重复步骤(1)~(2),直至达到预设的迭代次数iter_num;

(4)取统计点数最多的直线方程为该相电力线点云的平面直线方程,到该平面直线方程距离小于第四距离阈值Td4的所有点均为该相电力线点,其他点为异常点,即噪声点,剔除噪声点。

四、多分裂导线的自动提取和精细建模

多分裂导线的精细建模是基于单相导线进行三维重建,是本发明核心内容,流程见图5。

本步骤具体流程如下:

对步骤三提取的各相电力线点云分别进行如下步骤:

(1)对该相电力线点云进行拟合,计算拟合残差,若残差小于残差阈值Tm,则该相电力线点云为单相地线点云,执行步骤(2);否则为多分裂的单相导线点云,执行步骤(3)。残差计算过程中,将误差大于三倍残差阈值Tm的点作为粗差点剔除。本发明中,残差阈值Tm的设值与点云质量有关,点云质量较好时,残差阈值设置为较小值;点云质量较差时,残差阈值设置为较大值。一般情况下残差阈值Tm取0.2。

(2)采用RANSAC法对单相地线点云进行悬链线拟合。

(3)基于二分法从单相导线点云中提取分裂子导线点云,并采用RANSAC法对分裂子导线点云进行悬链线拟合。

将单相导线点云L投影到XOZ竖直面上,将投影等分为若干导线段。当导线段长度较小时,可认为直线。用最小二乘法对每段导线段分别进行线性拟合,得到直线方程Ax+Bz+C=0,对该导线段上所有点(xk,zk)分别计算:

式(4)中,(xk,zk)表示导线段上第k个点。

预设距离阈值DIS,若d>DIS,则认为该点位于该单相导线点云L的上部分,记为上分裂子导线点;若d<-DIS,则认为该点位于该单相导线点云L的下部分,记为下分裂子导线点;若d介于-DIS~DIS间,则认为该点位于该单相导线点云L的中部分,记为中分裂子导线点。

DIS≥0,其值根据多分裂导线规格设置,一般为两根最邻近分裂子导线距离的一半,特别的,对二、四分裂电力线直接设置DIS=0即可。

对各导线段全部划分完后,得到单相导线点云L的上分裂子导线点云S(0)、中分裂子导线点云S(1)和下分裂导线点云S(2)。分别判断各分裂子导线点云S(s)(s=0、1、2)中点数是否大于点数阈值NUM,点数阈值NUM用来判断粗差点。一般来说,一根分裂子导线上点数至少大于200,所以本实施例中,NUM设置为200。若不大于NUM,则将该分裂子导线点云S(s)作为粗差点剔除;否则对该分裂子导线点云S(s)进行基于RANSAC的悬链线拟合并计算残差;若残差小于残差阈值Tm,则认为是该分裂子导线点云S(s)中仅包含一根分裂子导线,直接输出悬链线拟合结果;若不小于残差阈值,则认为该分裂子导线点云S(s)中任然包含多根分裂子导线,对该分裂子导线点云S(s)继续执行如下步骤(a)~(b):

(a)将分裂子导线点云S(s)整体投影到XOY水平面上,根据最小二乘原理对投影进行线性拟合,得到直线方程ax+by+c=0,对该分裂子导线点云S(s)中所有点(x'k,y'k)分别计算:

若t>DIS,则认为该点位于该分裂子导线点云S(s)的左部分,记为左分裂子导线点;若t<-DIS,则认为该点位于该分裂子导线点云S(s)的右部分,记为右分裂子导线点;若t介于-DIS~DIS之间,则认为该点位于该分裂子导线点云S(s)的中间部分,记为中分裂子导线点。

将该分裂导线点云S(s)全部划分完后,得到S(s)的左分裂子导线点云T(0)、中分裂子导线点云T(1)和右分裂子导线点云T(2)。分别判断各分裂子导线点云T(u)(u=0、1、2)上点数是否大于点数阈值NUM,若大于,对该分裂子导线点云T(u)进行基于RANSAC的悬链线拟合并计算残差,输出悬链线拟合结果,将残差作为精度评定结果;否则,将该分裂子导线点云T(u)作为粗差点剔除。

(b)重复步骤(a),直至单相导线点云L中各分裂子导线点云S(s)均被处理。

(4)重复步骤(1)~(3),直至各相电力线点云均被处理。

下面将结合实施例进一步说明本步骤。

对二分裂导线,将各相导线点云投影XOZ竖直面,并分段进行线性拟合。将每相导线点云分为上分裂子导线点云S(0)和下分裂子导线点云S(2),对S(0)和S(2)分别进行基于RANSAC的悬链线拟合并进行误差分析,由于两个分裂子导线点云均仅包含一根分裂子导线,进行误差分析时残差小于残差阈值,因此直接输出S(0)和S(2)的悬链线拟合结果。

对四分裂导线,将各相导线点云投影到XOZ竖直面,并分段进行线性拟合,将每相导线点云分为上分裂子导线点云S(0)和下分裂子导线点云S(2)。对S(0)和S(2)分别进行悬链线拟合并进行误差分析,由于两个分裂子导线点云中均包含两根分裂子导线,进行误差分析时残差大于残差阈值,需继续将两个分裂子导线点云整体投影到XOY水平面上并整体拟合,将两个分裂子导线点云均划分为左分裂子导线点云和右分裂子导线点云,此时,所划分的左分裂子导线点云和右分裂子导线点云中均仅包含一根分裂子导线,对左分裂子导线点云和右分裂子导线点云分别进行基于RANSAC的悬链线拟合并计算残差。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1