2. 河北工业大学土木与交通学院,天津市西平道5340号,300401
三维激光扫描仪在获取桥梁几何空间信息方面具有高效率、高精度、非接触等独特优势[1]。由于采集到的桥梁点云数据量较大,不利于开展后续桥梁几何参数的测量分析和模型重建,因此需要将桥梁点云分割为具有语义信息的离散点集。
点云是包含三维坐标信息的离散点,体素包含若干离散点,超体素面片包含若干体素(简称“面片”),融合区域包含若干面片。点云分割可应用于植物叶片分割[2]、地面分割[3]、屋顶分割[4]、建筑物立面分割[5]和室内语义分割[6]等领域,主要方法包括随机采样一致性算法[7]、聚类法[8]、区域生长法[9]和机器学习法[10]。随机采样一致性算法对点云的噪声点和异常点表现出更好的鲁棒性,但过于依赖模型系数;聚类法能根据点云之间的距离、颜色等属性分组,空间聚集能力强,但对物体边界区分能力较弱;区域生长法生长结果的好坏取决于种子点的选取和相邻点区域的生长条件,对于桥梁的细分割很容易越过边界,造成错误分割;机器学习的分割算法能得到很好的语义分割效果,但训练时间较长且前期需要大量的数据集。对于建筑物点云的分割,Xu等[11]提出一种结合体素结构和区域增长的策略,利用平滑度、连续性和凹凸性作为几何线索,通过区域生长法对屋顶体素进行融合,但边界的检测仍存在困难;董震等[12]提出一种融合颜色、反射强度和空间距离的多尺度超体素方法,利用图像分割法分割点云,但需要点云本身带有多种特征,不适用于只包含空间特征的点云分割。
综上所述,本文利用桥面在桥梁结构中的位置特点和超体素面片的局部特性,基于桥梁的法线特征、曲率、空间位置和点云数量占比,提出一种邻接区域平面元融合的桥面分割方法。实验证明,该方法的分割效果较好,可实现桥梁点云中曲型桥面区域的分割。
1 桥面分割方法本文方法主要流程见图 1,图中每个处理过程的基本单位都包括点云、体素、超体素面片、融合区域,其中超体素面片和融合区域的最小组成单元都是体素。
超体素过分割的最小单元是体素,因此需要先对点云进行体素化处理。体素化是指给所有点云建立一个轴对齐包围盒,以体素分辨率Rvoxel为单位轴向划分为多个小包围盒,用小包围盒代表其内部的点云,该包围盒即为体素,后续计算过程中用体素的质点坐标表示体素的空间位置。在3×3×3的空间内一个体素周围最多有26个邻接体素,因此使用26邻接关系建立体素间的邻接图。
体素可有效减少后续算法的计算量,但难以代表一块区域的局部特征。超体素过分割能将具有相似属性的体素聚类成一块面片,相较于体素分割,其表示局部特征的能力更强。超体素过分割的步骤如下:
1) 以种子分辨率Rseed为单位建立若干初始种子。
2) 以Rsearch为搜索半径,计算种子范围内体素的数量,过滤掉小于阈值的种子,去除孤立点。
3) 以种子为中心逐层向外搜索,使用局部K均值聚类迭代生长法更新聚类中心,计算邻接体素与聚类中心所有体素之间相似性距离的平均值,将平均距离最近的邻接体素加入聚类。
4) 重复步骤3)直到无其他邻接体素或达到最大搜索范围,结束增长,获得超体素面片。面片中心点由靠近超体素质点的体素表示。
生长过程中的相似性距离D由空间位置和法线特征的加权方程决定:
$ D=\sqrt{w_{s} \frac{D_{s}^{2}}{3 R_{\text {seed }}^{2}}+w_{n} D_{n}^{2}} $ | (1) |
$ D_{s}=\sqrt{\left(x_{i}-x_{j}\right)^{2}+\left(y_{i}-y_{j}\right)^{2}+\left(z_{i}-z_{j}\right)^{2}} $ | (2) |
$ D_{n}=1-\cos \left(\boldsymbol{n}_{i}, \boldsymbol{n}_{j}\right) $ | (3) |
式中,ws为空间权值,Ds为体素间的空间距离,wn为法线权值,Dn为体素间的法线偏差,(xi, yi, zi)为第i个体素的位置,ni为第i个体素的法线,(xj, yj, zj)为第j个体素的位置,nj为第j个体素的法线,Rseed用于归一化处理。
1.2 面片过滤面片的属性特征是后续融合条件的重要依据,因此本文将面片的法线和曲率作为2个基本属性。首先将面片内所有体素的法线重定向,对所有体素法线加权融合获得面片法线,然后计算面片内所有体素的曲率平均值获得面片曲率。同时为防止融合过程越过区域边界,还需对面片进行过滤,过滤指标为面片法线的方向和模长,最终获得包含桥面的候选面片。
将面片中每个体素的法向量以y轴正方向为参考方向进行重定向:
$ \boldsymbol{n}_{i}^{*}=\left\{\begin{array}{l} \boldsymbol{n}_{i}, \boldsymbol{n}_{i} \cdot \boldsymbol{e}_{y} \geqslant 0 \\ -\boldsymbol{n}_{i}, \boldsymbol{n}_{i} \cdot \boldsymbol{e}_{y}<0 \end{array} i=1, 2, 3, \cdots, M\right. $ | (4) |
式中,M为面片的体素数量,ni为第i个体素的单位法向量,ey=(0, 1, 0)为y轴正方向的单位法向量,ni*为第i个体素重定向后的法向量。
重定向后利用分层加权融合的方法计算面片的法向量。将面片分为中心点和外围的L个面片层:
$ \left\{\begin{array}{l} \boldsymbol{n}_{\text {面片 }}^{*}=w_{0} \boldsymbol{n}_{0}^{*}+\sum\limits_{j=1}^{L} w_{j} \frac{1}{N_{j}} \sum\limits_{i=1}^{N_{j}} \boldsymbol{n}_{(j, i)}^{*} \\ w_{0}+\sum\limits_{j=1}^{L} w_{j}=1 \end{array}\right. $ | (5) |
式中,n面片*为面片中心点处法向量,w0为中心点权重,wj为第j层的权重,n0*为中心点法线,n(j, i)*为第j层第i个体素的法向量,Nj为第j层的体素数量。
桥梁不同部位的面片法向量在模长和方向上具有明显差异,图 2为面片过滤的2种指标。
第一个过滤指标是面片法线模长。由图 2(a)可见,虽然桥面部分是曲面,但超体素过分割后的桥面部分面片曲率几乎为0(Ⅰ类面片),设置的法线模长阈值εm用于过滤非桥面面片(Ⅱ类面片)。
第二个过滤指标是面片法线方向。由图 2(b)可见,桥面(区域1)融合遇到边界时法线变化不明显,符合融合条件。越过区域1和区域2的边界后,2个区域合二为一,会影响分割效果。由于桥面的面片法线和y轴正方向的夹角不超过桥面的最大坡度,因此设置法线夹角阈值εθ保留桥面区域、过滤边界,达到面片过滤分层的效果。最终获得的候选面片需满足公式:
$ \left\{\begin{array}{l} \left\|\boldsymbol{n}_{\text {面片 }}^{*}\right\|>\varepsilon_{m} \\ R\left(\boldsymbol{n}_{\text {面片 }}^{*}, \boldsymbol{e}_{y}\right)<\varepsilon_{\theta} \end{array}\right. $ | (6) |
平面元融合是指将候选面片融合为多个区域,是桥面分割的关键步骤。虽然桥面同时具有凸曲面和凹曲面,但经过过分割后每个面片对于桥面而言相当于一块小平面,称为“面元”,使用多个面元表示曲面,对候选部分进行平面元融合。
邻接区域平面元融合可以获取多块具有几何意义的桥梁区域,包括桥面区域以及桥梁的其他曲面区域。设置种子面片Sci和邻接面片Bcj,输入为初始面片区域A、面片法向量n和曲率γ、面片之间的邻接关系索引表,输出为多组融合区域B。面片融合的具体步骤如下:
1) 为提高面片的增长效率,选取所有曲率最小的面片作为初始种子。
2) 从区域A中取出初始种子,将初始种子加入到一组区域Rc中,并加入种子队列Sc。
3) 开始区域增长:
① 从种子队列Sc中取出队头的种子面片Sci;
② 根据面片邻接关系索引表获得种子面片Sci的所有邻接面片,存于邻接区域Bc中;
③ 计算种子面片Sci与其邻接面片Bcj的相似性度量Sij,若Sij大于相似性度量阈值Sth,则加入该组区域Rc中;若面片曲率γ大于曲率阈值γth,则将该邻接面片Bcj加入种子队列Sc的队尾;
④ 清空邻接区域Bc,返回步骤①,直到种子队列Sc为空。
4) 将增长后的区域Rc添加至区域B中,清空区域Rc,返回步骤1),直到区域A为空。至此,平面元融合结束,获得多组融合区域B。
上述方法中的相似性度量Sij由面片间的法线方向和中心点判决结果决定。法线方向用于确定面片之间的法线相似性Sn,在满足法线相似性的情况下,面片夹角可能会出现锐角或钝角的情况。由于桥面区域面片夹角为钝角,因此中心点判决可用于排除面片夹角为锐角的情况。
根据余弦定理,面片法线之间的余弦值公式为:
$ S_{n}=|\cos \alpha|=\left|\frac{\boldsymbol{n}_{i} \cdot \boldsymbol{n}_{j}}{\left|\boldsymbol{n}_{i}\right|\left|\boldsymbol{n}_{j}\right|}\right| $ | (7) |
式中,α为2个面片法线的夹角,ni为面片Sci的法线向量,nj为面片Bcj的法线向量。由于进行了法线重定向,因此法线余弦值取值范围为[0, 1]。
中心点判决过程如图 3所示。首先计算种子面片Sci和邻接面片Bcj之间的交线ℓ;然后从种子面片中心点M对交线ℓ作垂线,交点为M0,从邻接面片中心点P对交线ℓ作垂线,交点为P0;最后根据向量
$ k= \begin{cases}\text { true } & R\left(\overrightarrow{M_{0} M}, \overrightarrow{P_{0} P}\right)>\frac{\pi}{2} \\ \text { false } & R\left(\overrightarrow{M_{0} M}, \overrightarrow{P_{0} P}\right) \leqslant \frac{\pi}{2}\end{cases} $ | (8) |
受到体素分辨率的限制,超体素过分割后的面片间夹角不会出现极小锐角的情况。因此,当2个面片法线方向几乎相同时,无需进行中心点判决。设置面片法线间的夹角阈值αth=5°,计算公式为:
$ S_{i j}=\left\{\begin{array}{l} S_{n}, \alpha \leqslant \alpha_{\mathrm{th}} \\ k \cdot S_{n}, \alpha>\alpha_{\mathrm{th}} \end{array}\right. $ | (9) |
式中,k为邻接面片中心点判决结果,Sn为面片间法向量夹角余弦值,α为2个面片法线的夹角。
1.4 桥面点云为分割出最终的桥面点云,需要将以体素为最小单位的融合区域映射到以点云为最小单位的融合区域上,并计算每组融合区域点云数占原始桥梁点云数的比例,剔除占比不足1%的微小区域,剩余的融合区域包括桥面区域和非桥面区域。由于点云具有离散性,因此需要使用统计分布方法从区域中筛选出桥面点云。分析桥面在桥梁中的结构和位置特点可知,桥面位于桥梁较高处,向地面的投影类似矩形,且投影面积较大。最终使用矩形相似性、相对区域高度和投影面积指数3个指标来判断融合区域是桥面点云的可能性。第n组区域的3个指标计算步骤如下:
1) 桥梁点云y轴垂直于地面,y坐标最小值为Ymin,统计融合区域点云的y坐标,计算平均值Ymean_n。
2) 融合区域的点云向地面xoz投影,计算x轴方向的最大值Xmax_n和最小值Xmin_n,计算z轴方向的最大值Zmax_n和最小值Zmin_n。
3) 对投影区域进行正方形网格划分,设置每个网格边长为Gsize,统计每个网格中投影点的数量,生成点云分布矩阵。
4) 将点云分布矩阵转化为灰度图像,为更直观地展示投影形状,使用大津法[13]设置阈值,进行二值化处理,即图像中大于阈值的像素设为1,最后生成二值图。
5) 统计二值图中值为1的像素数,记为投影面积Sp_n。使用旋转主轴法[14]计算二值图的最小外接矩形,统计最小外接矩形所占用的像素数,记为外接矩形面积Srect_n。
6) 3个指标的计算公式为:
$ D_{\text {矩 }\_n}=\frac{S_{p_{-} n}}{S_{\text {rect_ } n}} $ | (10) |
$ H_{\text {mean_ } n}=\frac{Y_{\text {mean }\_n}-Y_{\text {min }}}{Y_{\text {mean_max }}-Y_{\text {min }}} $ | (11) |
$ R_{n}^{\prime}=1.0-0.02\left(R_{n}-1\right) $ | (12) |
式中,D矩_n为矩形相似性,Hmean_n为相对区域高度,Ymean_max为所有融合区域Ymean_n的最大值,R′n为投影面积指数,Rn为第n组区域投影面积在所有融合区域中的名次。由于投影面积不是决定性因素,为减弱该指标对结果的影响,将其与其他指标处于同一度量下,进行式(12)的计算。以上3个指标的取值范围都是[0, 1],指标越接近于1,是桥面点云的可能性越大。最终评价每一组融合区域的公式见式(13),通过比较每一组区域的Pn,筛选出最大值对应区域,即桥面点云:
$ P_{n}=D_{\text {矩_ } n} \cdot H_{\text {mean_ } n} \cdot R_{n}^{\prime} $ | (13) |
实验平台为Intel Core i5-9400F CPU@2.90GHz、VMware Workstation 15 Pro、Ubuntu18.04,开源点云库为PCL 1.11.1。为得到较为真实的实验数据,对拱形桥、钢筋桥和斜拉桥的3D模型进行采样,获得模型点云。由于模型内部也会被采样,因此需要使用点云软件Cloud Compare的Hidden Point Removal功能,从桥梁的7个视点模拟扫描出模型点云数据(图 4(a)),最后融合7组结果并删除重复点,获得实验数据。以拱形桥点云作为实验对象,结果如图 4(b)所示。
首先对原始点云数据进行体素化处理,既能保持桥梁外型轮廓不变,又能降低点云数据量、提高后续处理效率。体素分辨率为0.02 m,得到771 732个体素。然后进行超体素过分割,种子分辨率为0.2 m,获得18 627块面片。接着以ey=(0, 1, 0)为基准向量对面片的法向量作重定向处理,重定向后法向量与y轴的夹角范围为[0, π/2]。考虑到曲型桥面的坡度,为面片法线偏离y轴预留π/6的阈值。同时,为使桥面在候选面片中有更高的占比,设置法线模长阈值为0.8,过滤法线模长小于阈值的面片。图 5为面片过滤前后的效果对比,过滤后的候选区域剩余14 667块面片,法线方向垂直于y轴的面片被过滤,部分桥墩和护栏也被过滤,分割后获得较好的分层效果。
最后按照§1.3的步骤进行邻接面片的平面元融合,设置生长条件的相似性度量阈值为0.8,加入种子队列条件的曲率阈值为0.08,融合后共计4 546组区域,部分融合区域具有一定的语义信息,便于后续桥面的获取。
图 6为融合区域的筛选过程,借助PCL点云库中的getLabeledCloud函数将以体素为单位的融合区域映射到以点云为单位的融合区域上,计算每组区域点云数占桥梁点云数的比例,剔除不足1%的区域。图 6(b)中大部分微小区域被剔除,剩下5组融合区域,分别计算对应的评价指标,得到Pn。详细数据如表 1所示,虽然1号区域在投影面积方面不占优势,但另外2项指标得分较高,因此1号区域的Pn最高,被认为是桥面点云。结合图 6(c)可知,1号区域确实是桥面点云。
为进一步验证本文方法的可行性和稳定性,对钢筋桥和斜拉桥进行桥面分割实验(图 7),记录下3种桥梁(1号为拱形桥,2号为钢筋桥,3号为斜拉桥)的调试参数和结果(表 2、表 3)。由图 7可见,在面片过滤阶段,大部分钢筋、绳索和桥墩等非桥面区域被过滤掉,使得后续融合区域更加准确,最后统计各区域的点云分布,筛选出桥面点云。
由表 2可见,体素分辨率Rvoxel和种子分辨率Rseed的设置主要与桥梁点云的密度和结构尺寸有关,二者之间的关系可设置为Rseed≈(6~10)Rvoxel,Rvoxel和Rseed是影响分割时间的主要参数,其几何意义具有参考性,其他参数则相差不大。表 3为桥面点云的数量,为直观地展示本文方法的分割效果,使用Cloud Compare软件的Closest Point Set功能和Remove Duplicate Points功能计算算法分割桥面的结果以及人工分割桥面结果的交并比。由表可见,交并比均在95%以上。因此,本文方法对以上3种类型的桥面分割均具有较高的准确性。
2.3 与其他方法的比较由于在融合桥面区域的过程中会出现法线变化不明显的情况,从而导致融合的区域越界。为验证本文方法在特殊边界融合区域上的优势,采用区域生长法、超体素区域生长法[9]和本文方法进行比较,并截取S3DIS(stanford large-scale 3D indoor spaces)数据集的墙面制作点云数据,包括3块区域(R1、R2、R3)和2处边界(B1、B2),如图 8(a)所示。实验对象的数据无颜色信息,在MATLAB中使用渐变色展示其结构特征。T为处理数据所用的时间。
分析边界B2的越界情况可知,在点或面片过滤前,3种方法都会将区域R1、边界B2、区域R3融合到一起,这是由于边界B2点与点、面与面之间的法线方向变化不明显所致。在点或面片过滤后,点法线或面片法线过滤掉部分边界B2,边界B2被截断,区域R1和区域R3融合为单独区域。
分析边界B1的越界情况可知,在点或面片过滤后,仍然不能截断边界B1,这是因为该边界的点法线或面片法线与y轴的夹角在阈值εθ范围内,无法被过滤,最终导致区域R1和区域R2融合在一起。区域生长法依靠点法线变化判断生长条件,在融合区域时无法分辨边界B1的情况(图 8(b))。从图 8(a)的邻接关系图可以看出,区域R1和边界B1交界处的面片具有邻接关系,超体素区域生长法的生长条件只依靠相邻面片中心点法线的变化,因此面片融合时同样无法分辨边界B1的情况(图 8(c))。本文方法在法线变化的基础上加入中心点判决准则,防止相邻面片间夹角为锐角的情况出现(图 8(d)),区域R1和区域R2融合为单独区域。
综上所述,相较于2种传统方法,本文方法使3块区域单独融合的同时不会越过边界,在跨边界融合方面具有较好的抑制作用。在速度方面,以面片为单位的处理速度要优于以点为单位的处理速度。
3 结语本文提出一种邻接区域平面元融合的桥面分割方法,主要用于分割建筑物点云中的大型曲面,可极大减少人工干预。相较于传统分割方法,本文方法在融合过程中加入面片过滤和面片判决准则,可进一步避免融合区域越过边界的情况出现。
虽然本文方法的阈值参数具有较好的通用性,但超体素过分割种子点的选择仍然会影响过分割的效果,进而导致后续融合区域的边界较为粗糙。如何得到融合区域的清晰边界,将是后续的研究重点。
[1] |
徐进军, 郭鑫伟, 廖骅, 等. 基于地面三维激光扫描的桥梁挠度变形测量[J]. 大地测量与地球动力学, 2017, 37(6): 609-613 (Xu Jinjun, Guo Xinwei, Liao Hua, et al. The Test on Bridge Deflection Deformation Monitoring by Terrestrial Laser Scanning[J]. Journal of Geodesy and Geodynamics, 2017, 37(6): 609-613)
(0) |
[2] |
赖亦斌, 陆声链, 钱婷婷, 等. 植物三维点云分割[J]. 应用科学学报, 2021, 39(4): 660-671 (Lai Yibin, Lu Shenglian, Qian Tingting, et al. Three-Dimensional Point Cloud Segmentation for Plants[J]. Journal of Applied Sciences, 2021, 39(4): 660-671 DOI:10.3969/j.issn.0255-8297.2021.04.013)
(0) |
[3] |
张凯, 于春磊, 赵亚丽, 等. 基于自适应阈值的三维激光点云地面分割算法研究[J]. 汽车工程, 2021, 43(7): 1 005-1 012 (Zhang Kai, Yu Chunlei, Zhao Yali, et al. Research on Ground Segmentation Algorithm Based on Adaptive Thresholds for 3D Laser Point Clouds[J]. Automotive Engineering, 2021, 43(7): 1 005-1 012)
(0) |
[4] |
周钦坤, 岳建平, 李乐乐, 等. 一种分步式建筑物屋顶面点云高精度提取方法[J]. 大地测量与地球动力学, 2021, 41(6): 633-638 (Zhou Qinkun, Yue Jianping, Li Lele, et al. A Stepwise High-Precision Extraction for Point Cloud of Building Roof[J]. Journal of Geodesy and Geodynamics, 2021, 41(6): 633-638)
(0) |
[5] |
Hamid L F. Structural-Based Point Cloud Segmentation of Highly Ornate Building Facades for Computational Modelling[J]. Automation in Construction, 2019, 108
(0) |
[6] |
熊汉江, 郑先伟, 丁友丽, 等. 基于2D-3D语义传递的室内三维点云模型语义分割[J]. 武汉大学学报: 信息科学版, 2018, 43(12): 2 303-2 309 (Xiong Hanjiang, Zheng Xianwei, Ding Youli, et al. Semantic Segmentation of Indoor 3D Point Cloud Model Based on 2D-3D Semantic Transfer[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 2 303-2 309)
(0) |
[7] |
Poz A P D, Ywata M S Y. Adaptive Random Sample Consensus Approach for Segmentation of Building Roof in Airborne Laser Scanning Point Cloud[J]. International Journal of Remote Sensing, 2020, 41(6): 2 047-2 061 DOI:10.1080/01431161.2019.1683644
(0) |
[8] |
Zhao X K, Li H, Zhu Q B, et al. Automatic Sweet Pepper Detection Based on Point Cloud Images Using Subtractive Clustering[J]. International Journal of Agricultural and Biological Engineering, 2020, 13(3): 154-160 DOI:10.25165/j.ijabe.20201303.5460
(0) |
[9] |
杨琳, 翟瑞芳, 阳旭, 等. 结合超体素和区域增长的植物器官点云分割[J]. 计算机工程与应用, 2019, 55(16): 197-203 (Yang Lin, Zhai Ruifang, Yang Xu, et al. Segmentation of Plant Organs Point Clouds through Super Voxel-Based Region Growing Methodology[J]. Computer Engineering and Applications, 2019, 55(16): 197-203 DOI:10.3778/j.issn.1002-8331.1805-0221)
(0) |
[10] |
Boulch A, Guerry J, Saux B, et al. SnapNet: 3D Point Cloud Semantic Labeling with 2D Deep Segmentation Networks[J]. Computers and Graphics, 2018, 71: 189-198 DOI:10.1016/j.cag.2017.11.010
(0) |
[11] |
Xu Y S, Yao W, Hoegner L, et al. Segmentation of Building Roofs from Airborne LiDAR Point Clouds Using Robust Voxel-Based Region Growing[J]. Remote Sensing Letters, 2017, 8(11): 1 062-1 071 DOI:10.1080/2150704X.2017.1349961
(0) |
[12] |
董震, 杨必胜. 车载激光扫描数据中多类目标的层次化提取方法[J]. 测绘学报, 2015, 44(9): 980-987 (Dong Zhen, Yang Bisheng. Hierarchical Extraction of Multiple Objects from Mobile Laser Scanning Data[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(9): 980-987)
(0) |
[13] |
卜文斌, 游福成, 李泉, 等. 一种基于大津法改进的图像分割方法[J]. 北京印刷学院学报, 2015, 23(4): 76-78 (Bu Wenbin, You Fucheng, Li Quan, et al. An Improved Image Segmentation Method Based on Otsu[J]. Journal of Beijing Institute of Graphic Communication, 2015, 23(4): 76-78 DOI:10.3969/j.issn.1004-8626.2015.04.018)
(0) |
[14] |
张法全, 王国富, 曾庆宁, 等. 利用重心原理的图像目标最小外接矩形快速算法[J]. 红外与激光工程, 2013, 42(5): 1 382-1 387 (Zhang Faquan, Wang Guofu, Zeng Qingning, et al. New Algorithm for Minimum Enclosing Rectangle of the Object in the Image Region Based on Center-of-Gravity Principle[J]. Infrared and Laser Engineering, 2013, 42(5): 1 382-1 387)
(0) |
2. School of Civil and Transportation Engineering, Hebei University of Technology, 5340 Xiping Road, Tianjin 300401, China