Two-dimensional coordinate image measurement system based on angle intersection
-
摘要:
针对基准平面障碍物检测定位问题，设计了一种采用角度交会法测量二维坐标的嵌入式图像检测系统。二维坐标测量系统采用双线法标定线阵CCD图像传感器，并采用角度交会法对被测对象进行坐标计算，从而测得未修正的二维坐标测量结果。采用控制变量法分别测量X轴和Y轴坐标，使用Matlab软件对数据进行处理，并对X轴和Y轴测量误差分别进行多项式线性拟合，进而得到坐标修正公式，修正的二维坐标误差明显变小。实验结果表明：基于角度交会的二维坐标图像测量系统能够实时、准确、快速和可靠地测量二维坐标，为基准平面障碍物二维坐标测量定位提供了一种可行的方案，具有一定的应用价值和意义。
Abstract:Aiming at the problem of detection and positioning of reference plane obstacles, an embedded image detection system using angle intersection method to measure two-dimensional coordinates was designed. The two-line method was adopted to calibrate the linear array CCD image sensor by two-dimensional coordinate measurement system, and the angle intersection method was used to calculate the coordinates of the measured object, so that the uncorrected two-dimensional coordinate measurement results were obtained. The control variable method was used to measure the coordinates of X axis and Y axis, respectively, the Matlab software was used to process the data, and the polynomial linear fitting of the measurement errors of X axis and Y axis was carried out respectively to obtain the coordinate correction formula, in which the corrected two-dimensional coordinate errors became significantly smaller. The experimental results show that the two-dimensional coordinate image measurement system based on angle intersection can measure two-dimensional coordinates in real time, accurately, quickly and reliably, and provides a feasible scheme for two-dimensional coordinate measurement of reference plane obstacles, which has certain application values and significance.
-
引言
飞行目标的偏航角、俯仰角等姿态特征参数是火箭发射装备、发射火药及国防工程设计、定型、验收和故障诊断等过程的主要参数，准确地测量这些参数对提高武器效能和国防工程的防御能力意义重大[1-2]。目前，靶场姿态测量以两站以上交会测量为主，但在设备故障及捕获目标失误的情况下，单站测姿必须发挥重要作用，故研究单站测姿非常有必要。
综合现有文献的单站测姿方法，基本可划分为需要建立精确目标模型前提下的单站模型匹配法; 不需要相机内参数的DLT方法[3]，但该方法至少需要5对以上的物像点对应; 仅适用于投影轴对称目标的类似比例方法包括长宽比、椭圆长短轴比例法[4]等，该方法将透视投影做了平行投影的近似; 为突破该方法仅适合投影轴对称目标的局限，文献[5]利用了主轴以外的侧翼等信息，提出利用弱透视投影理论求解飞机等非投影轴对称目标的姿态角，但对模型缺乏系统、完整性分析; 文献[6]中使用的像点坐标差值比例法，该方法保持了透视模型的准确度，且利用了经纬仪较为准确的测角信息，但至少需要4个特征点，并对特征点的选取有一定要求，另外用到了站点的位置信息; 需要利用相机内外参数的PNP[7-9]以及光雷一体化方法[10]，其中PNP也需要3对以上的控制点对应，且大多基于迭代的方法。
本文测姿方法立足靶场实际，不需要建立目标模型，在只有单站图像序列，仅仅依靠某些先验特征长度、摄像系统相机内外参数以及成像信息的情况下，不需要利用测站位置信息。在保持透视投影模型准确度的前提下，对于投影轴对称目标，可以精确定位弹首、弹尾2个特征点作为输入条件，较准确地获取目标的姿态角。
1 基本原理
本文以中轴线为透视投影对象，将已知长度的目标向已知内、外相机参数的姿态测量系统像面投影，获得模拟投影像长; 同时提取实际图像序列中的目标实际像长，两者匹配相等，即可获取包括参考物距、目标深度等关键中间参数，然后利用关键参数，可计算出目标姿态角。大致流程图如图 1所示。
2 算法基础—透视投影2种等效方法
透视投影2种等效方式是本文获取关键中间参数的重要环节，可等效为2种方式。
2.1 目标跨度内的物面等效为垂直于光轴的若干连续物面组
经过目标垂直于主光轴的平面作为近似目标平面组，这样的平面连续且有无限个，故本来连续的目标物可以离散化为无数个连续的近似目标平面，透视投影则为无数个离散目标微元在像面上投影的累加。
2.2 以任意的其中一个等效物面为基准
可以近似目标平面组的任意一个平面为等效目标平面，按照透视投影的关系，连接光心及目标点的直线与等效目标平面相交，交点则为该目标点在等效目标平面的等效物点，如此获得等效物面内的等效物长。
3 目标姿态求解
3.1 方法分析及关键过程设置
本文方法涉及的主要坐标系，除了基准坐标系、摄像机坐标系外，还涉及目标体坐标系，与测量坐标系平行的过渡坐标系(零姿态坐标系)，O0-XpYpZp所有坐标系均符合右手坐标系。其中目标体坐标系O0-X0Y0Z0以体现目标姿态状态的中轴线为主轴建立; 设置初始姿态下的体坐标系与摄像机坐标系平行。
3.2 算法推导
已知条件：
物方 LAB=L;
像方 AB在像面的成像情况：可利用的包括A点成像AA(x1, y1)，B点成像BB(x2, y2);
相机内参数 像主点(x0, y0)以及焦距f已知，畸变暂不考虑;
相机外参数 光轴指向α、β角，即为摄像机坐标系相对于基准坐标系的旋转角度。
求解：目标中轴线相对于基准坐标系的姿态ψ和φ。
透视投影如图 2所示，始终设B所在物面为参考物面，A向B所在的参考物面投影点为D，∠ABD=θ，由上面的直线透视投影原理，计算过程为
1) 计算中间变量u及sinθ
$$ \left\{ \begin{align} & \sqrt{{{({{x}_{1}}-{{x}_{2}})}^{2}}+{{({{y}_{1}}-{{y}_{2}})}^{2}}}\times \mu \times {{10}^{-6}}=1 \\ & \frac{u}{{{10}^{3}}f}=\frac{\sqrt{{{(\frac{u}{\text{cos}{{\alpha }_{1}}})}^{2}}+{{(\frac{u}{\text{cos}{{\alpha }_{2}}})}^{2}}-\frac{u}{u+L~\text{sin}\theta }\left[ \frac{{{(u+L\text{sin}\theta )}^{2}}}{{{(\text{cos}{{\alpha }_{1}})}^{2}}}+\frac{{{u}^{2}}}{(\text{cos}{{\alpha }_{2}})}-{{L}^{2}} \right]}}{\mu \sqrt{{{({{x}_{1}}-{{x}_{2}})}^{2}}+{{({{y}_{1}}-{{y}_{2}})}^{2}}}} \\ \end{align} \right. $$ (1) 2) 计算关键特征点在摄像机坐标系下的坐标;
3) 计算目标相对于摄像机坐标系的姿态角;
4) 计算目标中轴线矢量在过渡坐标系下的坐标，即计算过渡坐标系下的目标中轴线矢量;
5) 经进一步坐标系旋转，计算单位中轴线矢量在基准坐标系下的坐标(X, Y,Z), 姿态角为
$$ \left\{ \begin{array}{l} \psi = {\tan ^{ - 1}}\frac{{{Z^\prime }}}{{{X^\prime }}} + k\pi \\ \phi = {\tan ^{ - 1}}\frac{{{Y^\prime }}}{{\sqrt {{X^\prime }^2 + {Z^\prime }^2} }} \end{array} \right.(k以实际象限取值) $$ (2) 4 试验验证与误差分析
4.1 试验验证
为了验证姿态模型的精度，我们在靶场距离姿态测量系统测试点1.7 km左右处放置了一门120火炮，以炮管为目标进行测量和验证，炮管的的方位和俯仰角度值由全站仪测得，姿态测量系统对炮管进行静态或动态扫描拍摄，获得了系列不同角度值的测量图像，以其中2组为例，如图 3所示。
其中姿态测量系统内外参数为：拍摄频率为100 fps，可见光高速电视分辨率1 536×1024，焦距值为1 523.1 mm，光轴指向如表 1所示，被研究的关键段长度约为2.655 m，用本文的算法进行姿态处理结果如表 1所示。
表 1 算法结果单位（°） Table 1. Results of algorithmα β 偏航绝对误差 俯仰角误差 1 323.063 3 -0.913 1 -1.975 6 1.5 2 323.053 2 -0.940 2 1.575 6 -0.8 3 323.102 1 -0.878 9 1.224 4 0.8 4 323.047 3 -0.861 6 -1.075 6 0.6 5 323.018 3 -0.870 3 1.775 6 0.3 6 323.019 9 -0.824 8 1.424 4 0.8 7 323.029 -0.867 1.924 4 1.2 8 323.052 -0.844 2.024 4 1.1 RMSE 1.7 1.0 从表 1中可得，姿态角度误差值偏航角约为1.7°，俯仰角误差约1.0°，主要原因在于偏航角主要靠对目标深度的鉴别力，而实际物距相对目标的比值大，故系统对偏航角的分辨力差。下面按照输入参数对模型进行误差分析与仿真。
4.2 误差分析
从模型(2)分析
$$ \left\{ \begin{align} & \psi =f(l, {{x}_{1}}, {{y}_{1}}, {{x}_{2}}, {{y}_{2}}, f, \alpha , \beta ) \\ & \phi =g(L, {{x}_{1}}, {{y}_{1}}, {{x}_{2}}, {{y}_{2}}, f, \alpha , \beta ) \\ \end{align} \right. $$ (3) 下面逐一分析关键模型参数对姿态角的影响。
1) 分析经纬仪指向误差对姿态角误差的影响。以表 1中经纬仪的测角数据为基准，按照目前设备的指向精度为指标参考，左右放缩0.3°为极限标准，分别以光轴指向的方位角和俯仰角2种情况进行分析，经纬仪指向误差对姿态角的测量误差趋势如图 4所示。
![]()
图 4 经纬仪指向误差对姿态角度误差的影响Figure 4. Influence of pointing error of theodoliteon pose error从图 4可知，测角精度对整个姿态结果有一定的影响，以目前的指向精度指标不超过0.3°为依据，对姿态角的影响不超过0.3°，光轴指向的方位角测角误差对偏航角影响大于俯仰角; 光轴指向的俯仰角测角误差对俯仰角的影响则大于偏航角。
2) 分析作用距离对姿态角误差的影响。从模型上分析，对作用距离敏感性较强的是偏航角，故作用距离对姿态角的影响主要以偏航角展开，作用距离对姿态角误差的影响如图 5所示。
从图 5中可以看出，作用距离对姿态角误差影响比较大，当横轴比值为4 800时，则姿态角绝对误差可达1°，当比值为7 200时，绝对误差可达2°。
3) 分析像点判读误差对姿态角的影响。以文中投影轴对称目标试验为基础，物距大约为1 700 m左右，在此物距及其它系统参数条件下，以表 1中的某一组数据为基础，2个端点像点反向移动下图中横轴的像素数，像素判读误差对姿态角误差影响如图 6所示。
![]()
图 6 特征像点提取误差对姿态角的影响Figure 6. Influence of feature points extractionerror on pose error从图 6可知，像素判读误差对姿态角误差影响占绝对比重，故在测量中，尽量利用目前比较成熟的亚像素定位精度技术，提高特征像点定位精度。
5 结论
文章提出了直线透视投影的2种形式，以长度已知的特征直线为研究对象，通过直线透视投影，可直接获取像长，或在某等效物面透视投影得到等效物长。以此为基础，通过像机内参数和外参数建立了投影轴对称的姿态测量模型; 并以炮管试验为依托，获得序列图像和相对真值，对透视算法进行了数据验证，结果为在试验设置的条件下投影轴对称姿态精度偏航角约为1.7°，俯仰角精度约为1.0°; 最后对算法进行了误差分析及实验仿真，分析了关键因素对姿态角误差的影响。
对于靶场姿态测量，文章的算法简单、方便，不需要3D建模及测试站点信息，测试精度满足靶场单站测姿需求，且该算法同样适用于非投影轴对称目标，为靶场姿态测量提供了一种新的思路和方法。
-
![]()
图 1 二维坐标测量系统总体框图
Figure 1. Overall block diagram of two-dimensional coordinate measurement system
![]()
图 2 角度交会法求解坐标示意图
Figure 2. Schematic diagram of solving coordinates by angle intersection method
![]()
图 5 TSL1401线阵CCD驱动时序图
Figure 5. Drive time-sequence diagram of TSL1401 linear array CCD
![]()
图 10 $ m_{i} $和$ n_{i} $采样点展示图
Figure 10. Display diagram of sampling points $ m_{i} $ and $ n_{i} $
![]()
图 11 $ \beta_{i} $与像元值的最小二乘法线性回归方程拟合图
Figure 11. $ \beta_{i} $ and least square linear regression equation for pixel values
![]()
图 13 Y轴坐标的测量误差多项式拟合曲线图
Figure 13. Polynomial fitting curve of measurement error for Y-axis coordinates
![]()
图 14 X轴坐标的测量误差多项式拟合曲线图
Figure 14. Polynomial fitting curve of measurement error for X-axis coordinates
表 1 $ m_{i} $采样点与对应像元值测试数据
Table 1 Sampling points $ m_{i} $ and corresponding pixel value test data
${\boldsymbol{m} }_{ {\boldsymbol{i} } }\_{ {\boldsymbol{X} } }_{ {\boldsymbol{i} } }$ ${C_{{x_i}}}$ ${\boldsymbol{m} }_{ {\boldsymbol{i} } }\_{ {\boldsymbol{X} } }_{ {\boldsymbol{i} } }$ ${C_{{x_i}}}$ ${\boldsymbol{m} }_{ {\boldsymbol{i} } }\_{ {\boldsymbol{X} } }_{ {\boldsymbol{i} } }$ ${C_{{x_i}}}$ 1 114 8 72 15 26 2 109 9 66 16 20 3 103 10 62 17 16 4 98 11 54 18 9 5 91 12 45 19 5 6 84 13 40 20 0 7 78 14 32 
下载: 导出CSV
表 2 $ n_{i} $采样点与对应像元值测试数据
Table 2 Sampling points $ n_{i} $ and corresponding pixel value test data
${\boldsymbol{n} }_{ {\boldsymbol{i} } }\_{ {\boldsymbol{X} } }_{ {\boldsymbol{i} } }$ ${C_{{x_i}}}$ ${\boldsymbol{n} }_{ {\boldsymbol{i} } }\_{ {\boldsymbol{X} } }_{ {\boldsymbol{i} } }$ ${C_{{x_i}}}$ ${\boldsymbol{n} }_{ {\boldsymbol{i} } }\_{ {\boldsymbol{X} } }_{ {\boldsymbol{i} } }$ ${C_{{x_i}}}$ 5.0 0 8.5 78 12.0 33 5.5 114 9.0 72 12.5 26 6.0 108 9.5 66 13.0 20 6.5 103 10.0 61 13.5 15 7.0 98 10.5 54 14.0 9 7.5 90 11.0 45 14.5 5 8.0 84 11.5 40 15.0 1
下载: 导出CSV
表 3 像元值与对应$ \beta_{i} $数据
Table 3 Image value and the corresponding $ \beta_{i} $ data
$ {C_{{x_i}}} $ $ {{\boldsymbol{\beta}} _{\boldsymbol{i}}} $ $ {C_{{x_i}}} $ $ {{\boldsymbol{\beta}} _{\boldsymbol{i}}} $ $ {C_{{x_i}}} $ $ {{\boldsymbol{\beta}} _{\boldsymbol{i}}} $ 0 45.0 78 73.3 32 111.8 114 48.0 72 78.8 26 116.6 108 51.3 66 84.3 20 121.0 103 55.0 61 90.0 15 125.0 98 59.0 54 95.7 9 128.7 90 63.4 45 101.3 5 132.0 84 68.2 40 106.7 1 135.0
下载: 导出CSV
表 4 Y=10时坐标测量数据
Table 4 Coordinate test data when Y=10
$ {\boldsymbol{Y}} $ $ {\boldsymbol{X}} $ ${\boldsymbol{Y} }_{ {{\rm{test}} } }$ $ {\boldsymbol{X}} $ ${\boldsymbol{Y} }_{ {{\rm{test}} } }$ $ {\boldsymbol{X}} $ ${\boldsymbol{Y} }_{ {{\rm{test}} } }$ 10 1 10.22 8 9.86 15 10.20 10 2 10.07 9 9.87 16 10.17 10 3 10.00 10 9.91 17 10.10 10 4 9.96 11 9.97 18 10.09 10 5 9.92 12 10.04 19 9.94 10 6 9.89 13 10.12 10 7 9.87 14 10.17
下载: 导出CSV
表 5 不同Y轴坐标测量数据
Table 5 Different Y-axis coordinate measurement data
$ {\boldsymbol{Y}} $ $\bar{\boldsymbol{Y} }_{ {{\rm{test}} } }$ $ \Delta {\boldsymbol{Y}} $ $ {\boldsymbol{Y}} $ $\bar{\boldsymbol{Y} }_{ {{\rm{test}} } }$ $ \Delta {\boldsymbol{Y}} $ 2.0 2.045 0.045 7.0 6.994 −0.006 2.5 2.515 0.015 7.5 7.508 0.008 3.0 3.000 0.000 8.0 8.024 0.024 3.5 3.491 −0.009 8.5 8.531 0.031 4.0 3.983 −0.017 9.0 9.038 0.038 4.5 4.481 −0.020 9.5 9.534 0.034 5.0 4.975 −0.025 10.0 10.020 0.020 5.5 5.474 −0.026 10.5 10.506 0.006 6.0 5.975 −0.025 11.0 10.988 −0.012 6.5 6.483 −0.017
下载: 导出CSV
表 6 不同X轴坐标测量数据
Table 6 Different X-axis coordinate measurement data
$ {{\boldsymbol{X}}} $ $\overline{ { {\boldsymbol{X} } } }_{{\rm{test}}}$ $ \Delta {{\boldsymbol{X}}} $ $ {{\boldsymbol{X}}} $ $\overline{ { {\boldsymbol{X} } } }_{{\rm{test}}}$ $ \Delta {{\boldsymbol{X}}} $ 2 1.982 −0.018 11 10.987 −0.013 3 2.984 −0.016 12 11.977 −0.023 4 3.999 −0.001 13 12.974 −0.026 5 5.013 0.013 14 13.976 −0.024 6 6.022 0.022 15 14.988 −0.012 7 7.026 0.026 16 16.001 0.001 8 8.021 0.021 17 17.013 0.013 9 9.012 0.012 18 18.026 0.026 10 10.00 0.000 19 19.034 0.034
下载: 导出CSV
-
[1] 卢佳宝, 韩学辉, 王彩霞. 基于线阵CCD的精密尺寸测量系统[J]. 光电子·激光,2019,30(9):935-940. LU Jiabao, HAN Xuehui, WANG Caixia. Precise dimension measurement system based on linear CCD[J]. Journal of Optoelectronics·Laser,2019,30(9):935-940.
[2] 宫二敏, 黄强先, 余夫领. 共平面二维高精度工作台误差修正与实验研究[J]. 中国测试,2014,40(4):10-14. doi: 10.11857/j.issn.1674-5124.2014.04.003 GONG Ermin, HUANG Qiangxian, YU Fuling. Research of error correction and experiment of co-planar 2D high precision stage[J]. China Measurement & Test,2014,40(4):10-14. doi: 10.11857/j.issn.1674-5124.2014.04.003
[3] HOCKEN R J. Coordinate measuring machines and systems[D]. Boca Raton: Crc Press, 2016.
[4] CHEN Z Y, WANG Y C, YANG Y, et al. PSD: principled synthetic-to-real dehazing guided by physical priors[C]//2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). Nashville: IEEE, 2021: 7180-7189.
[5] 朱博源, 任永杰, 段显龙, 等. 旋转激光扫描无源多靶点定位方法[J]. 光学学报,2022,42(5):66-76. ZHU Boyuan, REN Yongjie, DUAN Xianlong, et al. Passive multi-target location method based on rotating laser scanning[J]. Acta Optica Sinica,2022,42(5):66-76.
[6] LI J, QIN H, WANG J Z, et al. Open street map-based autonomous navigation for the four wheel-legged robot via 3D-lidar and CCD camera[J]. IEEE Transactions on Industrial Electronics,2022,69(3):2708-2717. doi: 10.1109/TIE.2021.3070508
[7] 白玮本, 曹慧峰, 刘欣宜. 基于结构光的嵌入式图像采集系统[J]. 仪表技术与传感器,2021(8):55-58. doi: 10.3969/j.issn.1002-1841.2021.08.012 BAI Weiben, CAO Huifeng, LIU Xinyi. Embedded image acquisition system based on structured light[J]. Instrument Technique and Sensor,2021(8):55-58. doi: 10.3969/j.issn.1002-1841.2021.08.012
[8] PAN Y, DONG Y Q, WANG D L, et al. Three-dimensional reconstruction of structural surface model of heritage bridges using UAV-based photogrammetric point clouds[J]. Remote Sensing,2019,11(10):1204. doi: 10.3390/rs11101204
[9] SHAH V, KENIYA R, SHRIDHARANI A, et al. Diagnosis of COVID-19 using CT scan images and deep learning techniques[J]. Emergency Radiology,2021,28(3):497-505. doi: 10.1007/s10140-020-01886-y
[10] RUMBOLDT Z, HUDA W, ALL J W. Review of portable CT with assessment of a dedicated head CT scanner[J]. American Journal of Neuroradiology,2009,30(9):1630-1636. doi: 10.3174/ajnr.A1603
[11] 张同双, 傅敏辉, 钟德安, 等. 基于MLE算法的海上角度交会测量方法及其精度分析[J]. 电讯技术,2013,53(8):1033-1038. doi: 10.3969/j.issn.1001-893x.2013.08.012 ZHANG Tongshuang, FU Minhui, ZHONG Dean, et al. Marine angle intersection method based on MLE algorithm and its precision analysis[J]. Telecommunication Engineering,2013,53(8):1033-1038. doi: 10.3969/j.issn.1001-893x.2013.08.012
[12] 田又源, 程瑶, 贾宁, 等. 高速线阵CCD驱动与数据采集系统设计[J]. 仪表技术与传感器,2022(3):84-87. doi: 10.3969/j.issn.1002-1841.2022.03.016 TIAN Youyuan, CHENG Yao, JIA Ning, et al. Design of high speed linear CCD driving and data acquisition system[J]. Instrument Technique and Sensor,2022(3):84-87. doi: 10.3969/j.issn.1002-1841.2022.03.016
[13] 高毅, 于瀛, 杨絮, 等. 基于多波段光源的智能痕迹分色方法[J]. 应用光学,2023,44(1):46-54. GAO Yi, YU Ying, YANG Xu, et al. Intelligent trace color separation method based on multi-band light source[J]. Journal of Applied Optics,2023,44(1):46-54.
[14] 伍川辉, 尹纪磊, 郭辉, 等. 基于二维激光位移传感器的通过式轮对测量系统设计[J]. 仪表技术与传感器,2020(11):50-53. doi: 10.3969/j.issn.1002-1841.2020.11.011 WU Chuanhui, YIN Jilei, GUO Hui, et al. Design of wheel-set measurement system based on two-dimensional laser displacement sensor[J]. Instrument Technique and Sensor,2020(11):50-53. doi: 10.3969/j.issn.1002-1841.2020.11.011
[15] 顾天奇, 林述温, 陈剑雄, 等. 坐标测量数据稳健拟合算法与实验研究[J]. 机械设计与制造,2019(7):165-167. doi: 10.3969/j.issn.1001-3997.2019.07.041 GU Tianqi, LIN Shuwen, CHEN Jianxiong, et al. The robust fitting algorithm and experimental study of measurement data of coordinate measurement machine[J]. Machinery Design & Manufacture,2019(7):165-167. doi: 10.3969/j.issn.1001-3997.2019.07.041
-
期刊类型引用(4)
1. 钟志，王琛，刘磊，孙岩，刘彬，单明广. 数字全息声场成像测量技术综述（特邀）. 激光与光电子学进展. 2024(24): 11-24 . 
百度学术
2. 席锋，秦华锋，胡莉，林睿. 可视化聚焦纹影系统实现流场中的三维扰动成像. 光学学报. 2021(09): 90-95 . 百度学术
3. 陈水桥，程诗卓，郭鹍，陆子毅，许熠辉. 非均匀介质驻波声场纹影法成像及其声速测量的改善. 实验室研究与探索. 2020(12): 46-50 . 百度学术
4. 罗劲明，薛凌峰，郑景辉，曾锦富，李子文，纪桂彬. 空气折射率梯度分布的可视化实验. 物理实验. 2018(01): 19-21+28 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 80
- HTML全文浏览量: 15
- PDF下载量: 33
- 被引次数: 6
陕公网安备 61011302001501号 