Defect generation of small sample laser welding based on generative adversarial network
-
摘要: 为了提升深度学习模型在不均衡小样本激光焊接表面缺陷数据集上的性能，优化了小样本数据量输入下的对抗生成网络(Generative Adversarial Network, GAN)模型. 通过对比激光焊接缺陷和其他用于测试对抗生成网络的公共数据集在特征复杂度上的区别，设计了一种全新的OCM (one class mixup)模块，并将其引入至针对有限样本的stylegan2-ada中，以提升GAN的性能，加快其收敛. 试验结果表明，在分类模型上，通过OCM-stylegan2-ada生成的数据集，比原始数据集性能提升40%，比使用mixup和stylegan2-ada增强后的数据集性能上提升20%，同时生成的焊接缺陷图片质量大幅提升.Abstract: To improve the performance of deep learning models on an unbalanced dataset of small samples of laser welded surface defects, an adversarial generative network (GAN) model using small datasets as input is optimized. By comparing the difference in feature complexity between laser welding defects and other public datasets used to test adversarial generative networks, a new OCM (one class mixup) module is designed and introduced into the stylegan2-ada for a limited number of samples to improve the performance of the adversarial generative network and accelerate its convergence. The results show that the dataset generated by OCM-stylegan2-ada improves the performance of the classification model by 40% over the original dataset and by 20% over the dataset enhanced with mixup and stylegan2-ada. Also the quality of the visually generated images of weld defects is greatly improved.
-
Keywords:
- welding defects /
- generative adversarial network /
- small sample /
- unbalanced dataset
-
0. 序言
焊接结构的主S-N曲线法已经在焊接领域开展了广泛应用，取得了很好的应用效果. 该方法采用等效结构应力将不同接头形式、厚度和加载模式的疲劳S-N数据有效地统一起来，解决了S-N曲线选取的困扰，同时由于采用与外载荷平衡的结构应力，所以也具有网格不敏感性等优点[1-3].
ASME标准中给出的主S-N曲线法通常先采用典型工况计算单位外载荷下的等效结构应力，再乘以典型工况的载荷谱实现应力谱计算，最后结合主S-N曲线法采用损伤叠加的方式完成寿命评估，该流程实际是一种准静态计算流程[4]. 准静态流程实施时，要将构件从系统中隔离出来，在构件的载荷传递接口部位施加单位载荷和约束，当结构承受动态载荷时，这种静态平衡的处理方法与实际承载是存在一定差异的. 另外，准静态法要获得典型工况载荷谱就需要对实测载荷进行工况分解，处理过程中容易导致载荷成分的丢失. 以铁路货车车体疲劳评估为例，由于我国载荷谱的积累较少，通常采用美国AAR标准中的载荷谱结合主S-N曲线法开展车体设计阶段的疲劳可靠性分析. 由于载荷不一致性导致疲劳评估精度大打折扣，这与近年来发生的多起疲劳失效故障有很大关系，给制造企业和运营造成了较大的经济损失[5-6].
2012年铁路货车全尺寸车体疲劳试验台建成并投入使用后，为了获得国内线路载荷，先后对中国通用干线、神华线、鹰夏线等线路开展了线路动态响应测试工作，取得了大量试验数据，采用加速试验方法和时域迭代的线路模拟方法形成了国内铁路线路的台架载荷数据[7]. 当开展台架载荷数据应用研究时发现，由于载荷数据是施加到试验台架作动器上的位移和力的时域曲线，是动态加载过程，传统将车体隔离出来计算的准静态分析流程无法使用该载荷数据. 为了解决载荷数据加载问题，考虑将台架模型作为车体边界条件引入到车体仿真分析中，建立刚柔耦合分析系统实现复合工况的动应力计算，计算的时域应力能够与试验测试结果较好吻合[8-10]. 然而如何将刚柔耦合计算获得的动应力计算过程转化为动态结构应力计算过程呢，再采用主S-N曲线法开展疲劳评估，这是值得思考的问题.
鉴于此，文中从模态叠加法的动力学理论入手思考该问题，基于振型向量的线性叠加计算位移思想，提出模态结构应力的叠加实现动态结构应力计算方法，并采用通用软件进行程序模块开发，再通过实际工程应用验证其合理性，模态结构应力法有效地拓宽主S-N曲线疲劳预测方法的应用领域，为焊接结构的研究提供理论基础.
1. 模态结构应力法
1.1 模态叠加法的动力学计算理论
大型结构的多自由的系统，外载荷作用下的动力学方程为
$$ {\boldsymbol{M}\ddot u}(t) + {\boldsymbol{C}\dot u}(t) + {\boldsymbol{K}u}(t) = {{f}}(t) $$ (1) 式中：
$ {\boldsymbol{M}} $ 为质量矩阵；$ {\boldsymbol{C}} $ 为阻尼矩阵；$ {\boldsymbol{K}} $ 为刚度矩阵；$ {{u}}(t) $ 为随时间变化的位移；$ {{f}}(t) $ 为随时间变化的外力. 当采用模态叠加法计算时，首先采用模态集$ {{\varPhi }} $ 对动力学方程的位移第一次坐标变换，即$$ {{u}}(t) = \sum\limits_{i = 1}^n {{\xi _i}} (t)\left\{ {{\phi _i}} \right\},\begin{array}{*{20}{c}} {}&{{{\varPhi }} = \left\{ {{\phi _1},{\phi _2}, \cdots ,{\phi _n}} \right\}} \end{array} $$ (2) 式中：
$ \left\{ \phi \right\} $ 为模态集${{\varPhi }}$ 各阶振型；$ {\xi _i}(t) $ 为随时间变化的模态坐标.将式(2)带入式(1)中，得到解耦的动力学方程，即
$$ {\boldsymbol{\tilde {M} }\ddot \xi}(t) + {\boldsymbol{\tilde {C} }\dot \xi}(t) + {\boldsymbol{\tilde {K} }\xi}(t) = {{{\varPhi }}^{\text{T}}}{{f}}(t) $$ (3) 式中：
${\boldsymbol{\tilde{ M}}}$ ，${\boldsymbol{\tilde{ K}}}$ 和${\boldsymbol {\tilde {C}}}$ 分别为主坐标刚度质量矩阵、主坐标刚度矩阵和比例阻尼矩阵.对解耦后的振动方程式(3)进行求解，可获得模态坐标
$ \xi $ 的时域历程，将模态坐标带入式(2)可以获得随时间变化的位移值$ {{u}}(t) $ .1.2 模态结构应力计算理论
主S-N曲线法中的结构应力与外载荷平衡，与之对应的是有限元分析中的节点力. 目前，多数刚柔耦合的动力学分析软件是无法直接获得动态节点力的，根据有限元的基本理论，可以将动力学计算后的
$ u(t) $ 与单元刚度矩阵$ {{\boldsymbol{K}}^{\bf{e}}} $ 相乘获得动态节点力，再通过动态节点力计算动态结构应力. 操作时发现，这个流程操作复杂，计算量大，程序也不易实现.模态振型向量的叠加既然能够获得节点位移，那么模态结构应力叠加是否能获得动态结构应力. 由于模态集
$ {{\varPhi }} $ 中的各阶振型是结构的固有属性，是不随时间变化的. 在$ {{\varPhi }} $ 中的第$ i $ 振型$ \left\{ {{\phi _i}} \right\} $ 由于节点位移相对变化会产生节点力$ \left\{ {{F_i}} \right\} $ 和$ \left\{ {{M_i}} \right\} $ ，既然振型是结构的固有属性，很显然模态节点力也是常数矩阵，其不随时间变化的，那么全局坐标系下的模态节点力与模态振型向量之间就存在线性关系，所以线性叠加原理能够在结构应力计算时使用[11].模态计算完成后即可获得模态节点力，需要将结构总体直角坐标系下求解的焊趾处模态节点力向量
$ \left\{ {{F_j},{M_j}} \right\} $ 变换到焊线局部坐标系下，坐标变换矩阵$ {\boldsymbol{T}} $ 可以通过有限元网格模型的节点坐标算出，得到变换后的焊线局部坐标系下的力向量$ \left\{ {{{F'}_j},{{M'}_j}} \right\} $ ，局部焊线坐标系$ y' $ 始终与焊线的切向方向垂直，$ x' $ 方向与焊线相切.基于节点力向量及相邻节点距离，利用长度等效矩阵
$ {\boldsymbol{L}} $ ，将每一阶模态焊线局部坐标轴$ y' $ 方向的节点力转化为该方向单元边上的线载荷$ {f_{iy}} $ ，将每一阶模态焊线局部坐标轴$ x' $ 方向的节点力转化为该方向单元边上的线载荷$ {m_{ix}} $ .$$ \begin{gathered} {\left\{ {{f_{iy}}(t)} \right\}^T} = {\left\{ {{{F'}_{iy}}(t)} \right\}^T}{{\boldsymbol{L}}^{ - 1}} \\ {\left\{ {{m_{ix}}(t)} \right\}^T} = {\left\{ {{{M'}_{ix}}(t)} \right\}^T}{{\boldsymbol{L}}^{ - 1}} \\ \end{gathered} $$ (4) 在已知焊趾处各节点线力和线弯矩情况下，采用结构应力计算公式，可以得第
$ i $ 阶模态对应焊线上得模态结构应力向量$ \left\{ {\sigma _{\text{s}}^i} \right\} $ 为$$ \left\{ {\sigma _{{\rm{s}}}^i} \right\} = \left\{ {\sigma _{{\rm{m}}}^i} \right\} + \left\{ {\sigma _{{\rm{b}}}^i} \right\} = \frac{{\left\{ {{f_i}} \right\}}}{d} + \frac{{6\left\{ {{m_i}} \right\}}}{{{d^2}}} $$ (5) 式中：
$ d $ 为板厚；$ \left\{ {\sigma _{\text{m}}^j} \right\} $ 为模态膜应力；$ \left\{ {\sigma _{\text{b}}^j} \right\} $ 为模态弯曲应力.1.3 动态结构应力计算及疲劳评估
将动力学计算的模态坐标时间历程
$ {\xi _i}(t) $ 与式(5)中模态结构应力$ q_{0} $ 进行叠加，即可获得动态结构应力，即$$ \left\{ {\sigma _s^{}(t)} \right\} = \sum\limits_{i = 1}^n {\left\{ {\sigma _{\rm{s}}^i} \right\}} {\xi _{\rm{t}}}(t) $$ (6) 对复合工况加载下获得的动态结构应力进行统计，单轴应力响应下可采用雨流技术，多轴载荷时可采用PDMR方法，获得不同等级的结构应力变化范围及其循环次数，分别将不同等级的结构应力变化范围带入，计算动态等效结构应力
$ \Delta {S_{\rm{s}}} $ 变化范围为$$ \Delta S_{\rm{s}}=\frac{\Delta \sigma_{\rm{s}}}{d^{(2-m) / 2 m} \cdot I(r)^{1 / m}} $$ (7) 式中：
$ I(r) $ 是弯曲比$ r $ 的无量纲函数$r = {{\Delta {\sigma _{\rm{b}}}} {\Delta {\sigma _{\rm{s}}}}}$ ；$ m = 3.6 $ ；$ d $ 为板厚.根据线性Miner线性损伤累计法则可以得出最终寿命为
$$ N = \frac{1}{{ \displaystyle\sum\limits_{i = 1}^k {{{\left( {\frac{{{C_{\rm{d}}}}}{{\Delta {S_{\rm{s}}}(i)}}} \right)}^{\dfrac{1}{h}}} \times {n_i}} }} $$ (8) 式中：
$ k $ 为应力幅划分的等级.2. 软件开发
通过上述原理及软件的需求分析，设计了程序流程，如图1所示，具体为：①选择截止频率，计算柔性体的模态计算；②定义的焊线组件和焊线坐标，实现焊线路径排序；③计算焊线的模态节点力，包含3个力和3个弯矩分量；④根据焊线的节点坐标计算长度等效矩阵
$ {\boldsymbol{L}} $ ，再通过$ {\boldsymbol{L}} $ 矩阵获得线载荷，根据线载荷计算模态结构应力；⑤根据模态结构应力获得模态等效结构应力；⑥导入动力学计算的模态坐标，通过模态坐标与模态结构应力叠加获得时域动态结构应力；⑦结合主S-N曲线法完成寿命评估.在图1的模态结构应力实施流程中，焊线组件定义、焊线排序、模态节点力计算、等效长度变换L矩阵的求逆、模态线载荷计算和模态等效结构应力计算均与准静态计算方法相同，只是由于模态阶数的选择导致计算量不同. 2种方法的不同点在于，模态结构应力法分析时输入的是时域载荷，是动力学计算过程中，同时兼顾模态参与效应，考虑系统的频率和阻尼影响.
根据流程和计算公式，可在通用软件中对模态结构应力进行程序开发，以某应用软件平台进行二次开发的程序如图2所示.
3. 应用实例
为了说明模态结构应力法如何在载荷输入和疲劳评估2个方面解决焊接结构复合工况下的动态寿命预测，以铁路货车车体的疲劳评估为例，从试验等效载荷、仿真系统建模、疲劳寿命预测验证和结构改进等方面进行应用介绍.
3.1 试验等效载荷
准静态结构应力法是将实测的载荷变成载荷谱，通过载荷谱换算出应力谱进行寿命评估. 模态结构应力法将实测的时域载荷直接动态加载，获得结构应力的时域曲线，再统计出应力谱实现寿命评估. 在载荷边界处理上是二者的主要不同.
模态结构应力法实施时仍需考虑如下2个问题：①当时域载荷曲线长度较长时，如何提高计算效率；②当时域载荷较难直接获得时，如何间接获得等效载荷.
以铁路货车车体分析为例，一方面，由于不间断的线路测试，数据量大，小载荷循环较多，为了提高效率，通过识别车体所有测试数据中疲劳损伤较小的区段，设定阈值同步删除疲劳损伤影响较小的小载荷事件，进行试验数据压缩[7]；另一方面，铁路货车车体与转向架通过心盘、旁承和车钩进行载荷传递，如图3所示. 上、下心盘通过刚性圆面接触，旁承通过弹性橡胶接触，较难直接通过心盘、旁承载荷进行动态加载. 为了能够获得车体准确模拟线路状态的等效载荷，用于开展车体的疲劳评估，以压缩试验数据中的车体枕梁加速度为目标，依托车体疲劳试验台架，采用时域迭代的方法，先识别系统的频响函数，在通过反复激励试验台架使车体响应与压缩后数据的响应的一致，从而获得与实际线路相当的等效载荷.
车体疲劳试验台架如图4所示，获得的等效载荷包含垂、横向作动器的位移数据和纵向车钩力数据.
值得注意的是，压缩后的测试数据保留的是大载荷的低频响应，通过迭代方法获得的等效载荷也是低频输入. 转向架和车体组成系统中由于转向架中包含斜楔阻尼减震器，在线路载荷下的高频载荷会发生衰减，而试验台架与车体组成的系统，由于没有阻尼减震，在低频输入时也会存在一定共振现象，故而在等效载荷获取前通常先进行台架边界下的模态试验，为线路模拟时的消除动力放大效应提供参数依据，敞车的台架边界下的车体频率和阻尼见表1[12].
表 1 台架边界下的车体频率及阻尼比Table 1. Body frequency and damping ratio under the boundary of car body振型 C70E C80 BHP矿石车 频率 f/Hz 阻尼比 $\xi $(%) 频率 f/Hz 阻尼比 $\xi $(%) 频率 f/Hz 阻尼比 $\xi $(%) 侧滚 3.58 4.7 1.37 0.98 1.25 1.97 扭转 9.03 1.78 3.73 1.52 4.52 1.69 垂向一阶弯曲 9.94 0.94 8.80 0.82 — — 横向一阶弯曲 13.45 3.26 9.24 0.92 — — 3.2 仿真系统建模
台架试验获得的等效载荷实际是液压油缸的驱动，传统的准静态方法将车体隔离出来无法加载，需要台架的机械结构来传递载荷，因此建立了刚性台架机械结构模型，机械结构包含4个垂向、2个横向、4个纵向及1个车钩作动器，如图5所示.
焊接结构研究对象选用的是矿石车车体，将装载的散体煤按质量单元附着在地板上[9]，采用固定交界面法获得柔性体结构，再将台架和车体模型连接起来，形成刚柔耦合多体模型，如图6所示，通过调整连接部位刚度，使仿真系统的频率与实际台架基本一致，其中低频模态的阻尼比按表1设置，高频模态的阻尼比按1%设置.
![]()
图 6 刚柔耦合多体模型和台架试验Figure 6. Rigid-flexible coupling simulation model and physical test rig. (a) rigid-flexible coupling model; (b) physical test将物理台架试验获得的多通道等效载荷作为刚柔耦合多体模型的输入，该载荷为澳大利亚的纽曼到黑德兰港之间的主线，单程为430 km，往返860 km，压缩后的载荷长度为280.4 s. 结果如图7所示，车体中梁中部焊缝动应力经过3 s后瞬态振动发生衰减，进入稳态振动，计算结果表明实测应力和仿真应力的波形基本一致.
刚柔耦合动力学计算后，导出柔性车体模态坐标时域曲线，为计算动态结构应力提供权值.
3.3 疲劳寿命预测验证
通过开发的软件计算了矿石车车体的模态等效结构应力，再将模态等效结构应力与动力学计算的模态坐标进行叠加，获得等效结构应力时间历程，采用钢材的中值主S-N曲线数据，对图8中10条车体关键焊缝进行了寿命评估. 结果发现10号的压车梁焊缝是所有焊缝中的相对薄弱部位.
10号焊缝路径从1节点开始到另一侧端部节点，节点数量为30，如图9所示. 在第5节点、10节点及对称的21节点、26节点为薄弱部位.
![]()
图 9 压车梁U形焊缝的寿命Figure 9. Life of U-shaped welds in the roller beam. (a) weld path; (b) weld life选取A部位节点的动态等效结构应力时域曲线的峰值点，将各阶模态的模态坐标和模态等效结构应力绘制到图10中，可见该焊缝主要是41阶模态产生的作用，该模态对应于车体的纵向加载过程，车体结构应该根据纵向载荷进行优化.
![]()
图 10 焊缝模态坐标和模态等效结构应力Figure 10. Modal coordinates and modal equivalent structural stresses of weld驱动文件总长为280.4 s，按25年全寿命考核，共5814次循环，用时长约为453 h.
等效15 ~ 20年试验时，在车体压车梁上的U形焊缝的A，B，C 3个部位均发现了疲劳裂纹，如图11所示，裂纹的出现位置与采用模态结构应力法计算结果吻合.
![]()
图 11 车体疲劳裂纹Figure 11. Fatigue crack of car body. (a) cracks at position A; (b) enlarged view of position A; (c) cracks at position B; (d) cracks at position C可见，采用试验台架作为边界条件，实现了以线路等效载荷开展动力学计算，再使用模态结构应力法与主S-N曲线法评估结合，能够实现复合工况下车体的关键焊缝进行疲劳可靠性预测，可视为物理台架的数字孪生. 通过该方法能够在一定程度上部分替代实物台架，开展结构的改进验证，以矿石车为例，改进结果如图12所示，通过再次疲劳评估表明改进后U形焊缝的损伤明显变小.
![]()
图 12 车体结构改进前后比较Figure 12. Structural improvement of car body. (a) before improvement; (b) after improvement4. 结论
(1)加载的准确性是保证焊接结构疲劳寿命预测的关键，准静态的简化处理容易使载荷成分丢失或无法考虑模态参与效应，可考虑建立引入边界条件模型，使仿真模型与实际结构承载具有可比性.
(2)以模态坐标和模态结构应力叠加实现动态结构应力计算的疲劳评估方法，能够兼顾复合工况的加载和模态参与效应，更能有效的识别出动态加载下车体的疲劳薄弱部位.
(3)模态结构应力法有效地拓宽了主S-N曲线疲劳预测方法的应用领域，为研究焊接结构动态载荷下疲劳可靠性提供技术基础.
-
![]()
图 1 Stylegan2-ada与OCM-stylegan2-ada网络结构对比
Figure 1. Comparison of Stylegan2-ada and OCM-stylegan2-ada network structures. (a) Stylegan2-ada; (b) OCM-stylegan2-ada
![]()
图 2 激光焊接缺陷
Figure 2. Laser welding defects. (a) normal; (b) porosity; (c) lack of penetration; (d) black
![]()
图 3 不同算法生成图片对比
Figure 3. Different algorithms generated image comparisons. (a) Original image; (b) Mixup; (c) Stylegan2-ada; (d) OCM-stylegan2-ada; (e) OCM-stylegan2-ada; (f) OCM-stylegan2-ada
![]()
图 7 OCM0和OCM1生成图片对比
Figure 7. OCM0 and OCM1 generated image comparisons (a) OCM0-stylegan2-ada; (b) OCM1-stylegan2-ada
表 1 不同数据集的测试准确率
Table 1 Test accuracy of different datasets
数据集种类 正常 发黑 孔洞 未焊透 整体 原始 0.854 0.446 0.438 0.868 0.661 Mixup 0.794 0.973 0.889 0.600 0.838 Stylegan2-ada 0.955 0.951 0.291 0.676 0.755 OCM0-Stylegan2-ada 1.000 0.996 0.994 0.976 0.991 OCM1-Stylegan2-ada 0.992 1.000 1.000 0.987 0.996 
下载: 导出CSV
表 2 不同数据集的AUC值
Table 2 AUC values of different datasets
数据集种类 正常 发黑 孔洞 未焊透 整体 原始 0.839 0.693 0.638 0.863 0.758 Mixup 0.881 0.954 0.895 0.783 0.878 Stylegan2-ada 0.877 0.932 0.633 0.832 0.819 OCM-Stylegan2-ada 1.000 0.996 0.994 0.976 0.991 OCM-Stylegan2-ada 0.999 0.997 0.999 0.996 0.998
下载: 导出CSV
-
[1] Mery D, Riffo V, Zscherpel U, et al. GDXray: The database of X-ray images for nondestructive testing[J]. Journal of Nondestructive Evaluation, 2015, 34(4): 1 − 12.
[2] Zhang Hongyi, Moustapha Cisse, Yann N. , et al. Mixup: beyond empirical risk minimization[C].//International Conference on Learning Representations. 2018.
[3] Creswell A, White T, Dumoulin V, et al. Generative adversarial networks: An overview[J]. IEEE Signal Processing Magazine, 2018, 35(1): 53 − 65. doi: 10.1109/MSP.2017.2765202
[4] 黄旭丰. 基于深度迁移学习的焊接质量在线监测方法研究[D]. 南宁: 广西大学, 2019. Huang Xufeng. Deep transfer learning for online welding quality monitoring[D]. Nanning: Guangxi University, 2019.
[5] Gurumurthy S, Kiran Sarvadevabhatla R, Venkatesh Babu R. Deligan: Generative adversarial networks for diverse and limited data[C]//Proceedings of the IEEE conference on computer vision and pattern recognition. 2017: 166 − 174.
[6] 谷静, 张可帅, 朱漪曼. 基于卷积神经网络的焊缝缺陷图像分类研究[J]. 应用光学, 2020, 41(3): 531 − 537. GU Jing, ZHANG Keshuai, ZHU Yima. Research on weld defect image classification based on convolutional neural network[J]. Journal of Applied Optics, 2020, 41(3): 531 − 537.
[7] Radford A , Metz L , Chintala S . Unsupervised Representation Learning with Deep Convolutional Generative Adversarial Networks[C]//International Conference on Computer Vision. 2016.
[8] Karras T, Aittala M, Hellsten J, et al. Training generative adversarial networks with limited data[J]. Advances in Neural Information Processing Systems, 2020, 33: 12104 − 12114.
[9] Karras T, Laine S, Aila T. A style-based generator a-rchitecture for generative adversarial networks[C]//Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2019: 4401 − 4410.
[10] Karras T, Laine S, Aittala M, et al. Analyzing and improving the image quality of stylegan[C]//Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2020: 8110 − 8119.
[11] Heusel M, Ramsauer H, Unterthiner T, et al. GANS trained by a two time-scale update rule converge to a local nash equilibrium[J]. Advances in Neural Information Processing Systems, 2017, 30: 1 − 38.
[12] Howard A, Sandler M, Chu G, et al. Searching for mobilenetv3[C]//Proceedings of the IEEE/CVF Interna-tional Conference on Computer Vision. 2019: 1314 − 1324.
[13] Howard A G, Zhu M, Chen B, et al. Mobilenets: Efficient convolutional neural networks for mobile vision applications[C]//Proceedings of the IEEE conference on computer vision and pattern recognition. 2017: 432 − 445.
[14] Sandler M, Howard A, Zhu M, et al. Mobilenetv2: Inverted residuals and linear bottlenecks[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2018: 4510 − 4520.
计量
- 文章访问数: 503
- HTML全文浏览量: 119
- PDF下载量: 97
