赞
踩
LSM6DSV16X包含三轴陀螺仪与三轴加速度计。
姿态有多种数学表示方式,常见的是四元数,欧拉角,矩阵和轴角。他们各自有其自身的优点,在不同的领域使用不同的表示方式。在四轴飞行器中使用到了四元数和欧拉角。
姿态解算选用的旋转顺序为ZYX,即IMU坐标系初始时刻与大地坐标系重合,然后依次绕自己的Z、Y、X轴进行旋转:
绕IMU的Z轴旋转:航向角yaw
绕IMU的Y轴旋转:俯仰角pitch
绕IMU的X轴旋转:横滚角row
最近在弄ST和瑞萨RA的课程,需要样片的可以加群申请:615061293 。
横滚roll,俯仰pitch,偏航yaw的实际含义如下图:
由于需要解析姿态角,故将陀螺仪速度修改快一点。
https://www.bilibili.com/video/BV18M411d7Wg/
陀螺仪LSM6DSV16X与AI集成(2)----姿态解算
https://www.wjx.top/vm/OhcKxJk.aspx#
https://download.csdn.net/download/qq_24312945/88598263
横滚角φ:机体绕OBXB转动,轴Y’B与平面OBXBYB构成的夹角。
俯仰角θ:机体绕OBYB转动,轴Z’B与平面OBYBZB构成的夹角。
偏航角ψ:机体绕OBZB转动,轴X’B与平面OBXBZB构成的夹角。
将姿态角从机体坐标系转换到惯性坐标系中是为了便于分析无人机状态,反映无人机在惯性坐标系下的姿态运动状态,利用齐次线性变换可实现坐标系的转换,旋转矩阵就是在线性变化中产生的,用REB表示惯性坐标系{E}到机体坐标系{B}的变换。
例如,绕OBXB旋转必角,此时两个坐标系存在必的角度差,不再重合。点(x, y, z)的转换方程为:
可提取转换矩阵:
同理,绕口OBYB旋转θ角得:
而绕OBZB旋转ψ角得:
不同旋转顺序有不同的旋转矩阵,按照偏航,俯仰,横滚的顺序,即分别绕X-Y-Z旋转,就可计算出旋转矩阵REB,REB等于依次旋转所得的矩阵连乘,且顺序为从右向左排列。
当俯仰角θ=±Π/2时,横滚运动与偏航运动的旋转轴重合,出现万向节死锁现象,在空间失去了一个自由度。如式所示,φ或ψ的变化具有相同的效果,因此不再具有唯一性啊。
本文选择的是四元数法进行姿态解算。无人机姿态解算方法主要有四种,它们各自的优缺点如下图所示。欧拉角法不能用于计算飞行器的全姿态角,且难以实时计算而不易于工程应用。方向余弦法不会出现“奇点”现象,但计算量大,效率低。四元数法避免了复杂的三角函数运算,变为求解线性微分方程,算法简单易操作,且不存在角度奇异性问题,可以更好的线性化系统,是一种更实用的工程方法。
四元数的概念诞生在1843年的爱尔兰,是数学家哈密顿研究空间几何时提出。在如今的导航技术领域,四元数的优势逐渐被发现,得到了研究者们的广泛关注,并逐渐应用在姿态解算领域。
四元数是由四个元构成的数Q(q0,q1,q2,q3) = q0 + q1i + q2j + q3k;其中,q0,q1,q2,q3是实数,i,j,k既是互相正交的单位向量,又是虚单位根号-1。四元数即可看作四维空间中的一个向量,又可以看做一个超复数。对于后续有一个重要的变化需要记住:
Q=q0 + q1i + q2j + q3k
可视为一个超复数,Q 的共轭复数记为
Q’=q0 - q1i - q2j - q3k
Q°称为Q的共轭四元数。
同时,有
ij=k,jk=i,ki=j,ji=-k,kj=-i,ik=-j
i2 = j2 = k2 =ijk=-1
其中,i、j、k是相互正交的单位向量,其几何意义可理解为分别绕三个坐标轴的旋转,q0、q1、q2、q3为常数,有
通过四元数进行欧拉角求解,可以减少芯片运算负担,提高运算速度。
一个矢量V相对于坐标系OXYZ固定:V = xi + yj + zk;坐标系OXYZ转动了Q得到一个新坐标系OX’Y’Z’:V = x’i’ + y’j‘ + z’k’;设四元数Ve、Ve‘
Ve = xi + yj + zk;
Ve’ = x’i + y’j + z’k;
则Ve’ = Q* Ve * Q’;
设Q = q0 + q1i + q2j + q3k;则Q’ = q0 - q1i - q2j - q3k;
则Ve’ = Q* Ve * Q’=(q0 + q1i + q2j + q3k) * (0+xi + yj + zk) + (q0 - q1i - q2j - q3k)
可以算出
x’=(q0 ^2+q1 ^2-q2 ^2-q3 ^2)x+2(q1q2+ q1q3)y+2(q1q3-q0q2)z
y’ = 2(q1q2-q0q3)x+(q0 ^2-q1 ^2+q2 ^2-q3 ^2)y+2(q2q3+q0q1)z
z’ = 2(q1q3+q0q2)x+2(q2q3-q0q1)y+(q0 ^2-q1 ^2-q2 ^2+q3 ^2)z
结合
可以反推
Pitch = asin(2 * q2 * q3 + 2 * q0* q1)* 57.3; // pitch ,转换为度数
Roll = atan2(-2 * q1 * q3 + 2 * q0 * q2, q0*q0-q1*q1-q2*q2+q3*q3)* 57.3; // rollv
Yaw = atan2(2*(q1*q2 - q0*q3),q0*q0-q1*q1+q2*q2-q3*q3) * 57.3; //偏移太大,
将加速度的三维向量转为单位向量
// 测量正常化
norm = sqrt(ax*ax + ay*ay + az*az);
ax = ax / norm; //单位化
ay = ay / norm;
az = az / norm;
世界坐标系重力分向量是通过方向旋转矩阵的最后一列的三个元素乘上加速度就可以算出机体坐标系中的重力向量。
// 估计方向的重力
vx = 2*(q1*q3 - q0*q2);//由下向上方向的加速度在加速度计X分量
vy = 2*(q0*q1 + q2*q3);//由下向上方向的加速度在加速度计X分量
vz = q0*q0 - q1*q1 - q2*q2 + q3*q3;//由下向上方向的加速度在加速度计Z分量
陀螺仪能够迅速响应设备的旋转,在短时间内误差较小且可靠。然而,因为温度漂移、零漂移和积分误差会随时间累积,陀螺仪的长时间精度受到影响。在静止状态下,加速度计的漂移很小,其倾角求解过程中不存在积分误差,但在飞行过程中,加速度计受到发动机和机架振动以及转动和运动加速度的干扰。磁罗盘测量的地磁向量在特定地理范围内可视为不变,但磁罗盘易受硬磁场和软磁场干扰。
因此,若系统外环采用九轴姿态传感器(包括三轴加速度计、三轴磁罗盘和三轴陀螺仪)进行数据融合,磁罗盘易受干扰可能导致融合后的数据仍有较大误差。为此,在内环使用六轴姿态传感器(包括三轴加速度计和三轴陀螺仪)进行数据融合,对融合后的传感器姿态偏差进行二次修正,以提高整体精度。
外环九轴姿态传感器数据融合,记在飞行器机体坐标系下an=[ax ay az]T和mn=[mx my mz]T分别为加速度计和磁罗盘实际测量得到的重力向量和地磁向量。
记vn=[vx vy vz]T和wn=[mx my mz]T是将地理坐标系下重力向量kb=[0 0 1g]T和地磁向量nb=[nx 0 nz]T(不考虑地理磁偏角因素,将机头固定向北)通过四元数坐标换算成机体坐标系下的重力向量和地磁向量。向量之间的误差为坐标轴的旋转误差,可以用向量的叉积en=[ex ey ez]T表示,如下所示。
由于我的LSM6DS3TR-C为六轴,不带三轴陀螺仪,故代码如下。
//这个叉积向量仍旧是位于机体坐标系上的,而陀螺积分误差也是在机体坐标系,而且叉积的大小与陀螺积分误差成正比,正好拿来纠正陀螺。
//(你可以自己拿东西想象一下)由于陀螺是对机体直接积分,所以对陀螺的纠正量会直接体现在对机体坐标系的纠正。
ex = (ay*vz - az*vy);
ey = (az*vx - ax*vz);
ez = (ax*vy - ay*vx);
由于陀螺仪是对机体直接积分,所以,陀螺仪的误差可以体现为机体坐标的误差。因此修正坐标轴的误差可以达到修正陀螺仪误差的目的,从而将加速度计和磁罗盘进行修正陀螺仪,实现了九轴的数据融合。即如果陀螺仪按照叉积误差的轴,转动叉积误差的角度,就可以消除机体坐标上实际测量的重力向量和地磁向量和坐标换算后的重力向量和地磁向量之间的误差。
PI调节器的比例部分用于迅速纠正陀螺仪误差,积分部分用于消除稳态偏差。PI调节器的比例系数和积分系数自己去修正。陀螺仪经过外环PI控制器修正姿态误差后输出值为了gn =[gx gy gz]T
// 积分误差比例积分增益,计算陀螺仪测量的重力向量与估计方向的重力向量之间的误差。
exInt = exInt + ex*Ki;
eyInt = eyInt + ey*Ki;
ezInt = ezInt + ez*Ki;
// 调整后的陀螺仪测量,使用叉积误差来进行比例-积分(PI)修正陀螺仪的零偏。将修正量乘以比例增益Kp,并加上之前计算的积分误差exInt、eyInt和ezInt。
gx = gx + Kp*ex + exInt;
gy = gy + Kp*ey + eyInt;
gz = gz + Kp*ez + ezInt;
内环的六轴姿态传感器数据融合是将地理坐标系下的重力场向量与加速度计在机体坐标系下采集到的重力向量进行叉乘,求出两者向量误差。并通过PI控制器修正向量误差,从而达到修正外环九轴数据融合后的陀螺仪的偏差的目的。在每个姿态解算周期读取出机体坐标系下双环PI控制后的陀螺仪的角速率
整合四元数率和正常化,根据陀螺仪的测量值和比例-积分修正值,对四元数进行更新。
// 整合四元数率和正常化,根据陀螺仪的测量值和比例-积分修正值,对四元数进行更新。根据微分方程的离散化形式,将四元数的每个分量加上相应的微分项乘以采样周期的一半(halfT)。
q0 = q0 + (-q1*gx - q2*gy - q3*gz)*halfT;
q1 = q1 + (q0*gx + q2*gz - q3*gy)*halfT;
q2 = q2 + (q0*gy - q1*gz + q3*gx)*halfT;
q3 = q3 + (q0*gz + q1*gy - q2*gx)*halfT;
// 正常化四元数
norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
q0 = q0 / norm;
q1 = q1 / norm;
q2 = q2 / norm;
q3 = q3 / norm;
六轴传感器(包括三轴加速度计和三轴陀螺仪)可以用于估算设备在空间中的姿态,包括俯仰角(Pitch)、横滚角(Roll)和偏航角(Yaw)。然而,六轴传感器仅依赖陀螺仪和加速度计数据,可能无法准确测量偏航角(Yaw),原因如下:
无磁场参考:六轴传感器缺少磁罗盘,没有固定的参考方向。因此,在长时间内,陀螺仪的积分误差可能导致偏航角估计漂移。
陀螺仪误差累积:陀螺仪测量的是角速度,要得到偏航角,需要将角速度积分。由于陀螺仪存在零漂、噪声和温度漂移等误差,这些误差在积分过程中会累积,使得偏航角估计产生较大的漂移。
虽然六轴传感器可能无法准确测量偏航角,但可以通过将其与磁罗盘(三轴磁场传感器)结合,形成九轴传感器(包括三轴加速度计、三轴磁罗盘和三轴陀螺仪),以提高偏航角估计的准确性。九轴传感器融合了磁场信息,为偏航角提供了一个稳定的参考方向,有助于减小陀螺仪误差对偏航角估计的影响。
//加速度单位g,陀螺仪rad/s void IMUupdate(float gx, float gy, float gz, float ax, float ay, float az) { float norm; float vx, vy, vz; float ex, ey, ez; // 测量正常化,把加计的三维向量转成单位向量。 norm = sqrt(ax*ax + ay*ay + az*az); ax = ax / norm; //单位化 ay = ay / norm; az = az / norm; // 估计方向的重力,世界坐标系重力分向量是通过方向旋转矩阵的最后一列的三个元素乘上加速度就可以算出机体坐标系中的重力向量。 vx = 2*(q1*q3 - q0*q2);//由下向上方向的加速度在加速度计X分量 vy = 2*(q0*q1 + q2*q3);//由下向上方向的加速度在加速度计X分量 vz = q0*q0 - q1*q1 - q2*q2 + q3*q3;//由下向上方向的加速度在加速度计Z分量 //这个叉积向量仍旧是位于机体坐标系上的,而陀螺积分误差也是在机体坐标系,而且叉积的大小与陀螺积分误差成正比,正好拿来纠正陀螺。 //(你可以自己拿东西想象一下)由于陀螺是对机体直接积分,所以对陀螺的纠正量会直接体现在对机体坐标系的纠正。 ex = (ay*vz - az*vy); ey = (az*vx - ax*vz); ez = (ax*vy - ay*vx); // 积分误差比例积分增益,计算陀螺仪测量的重力向量与估计方向的重力向量之间的误差。 exInt = exInt + ex*Ki; eyInt = eyInt + ey*Ki; ezInt = ezInt + ez*Ki; // 调整后的陀螺仪测量,使用叉积误差来进行比例-积分(PI)修正陀螺仪的零偏。将修正量乘以比例增益Kp,并加上之前计算的积分误差exInt、eyInt和ezInt。 gx = gx + Kp*ex + exInt; gy = gy + Kp*ey + eyInt; gz = gz + Kp*ez + ezInt; // 整合四元数率和正常化,根据陀螺仪的测量值和比例-积分修正值,对四元数进行更新。根据微分方程的离散化形式,将四元数的每个分量加上相应的微分项乘以采样周期的一半(halfT)。 q0 = q0 + (-q1*gx - q2*gy - q3*gz)*halfT; q1 = q1 + (q0*gx + q2*gz - q3*gy)*halfT; q2 = q2 + (q0*gy - q1*gz + q3*gx)*halfT; q3 = q3 + (q0*gz + q1*gy - q2*gx)*halfT; // 正常化四元数 norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); q0 = q0 / norm; q1 = q1 / norm; q2 = q2 / norm; q3 = q3 / norm; Pitch = asin(2 * q2 * q3 + 2 * q0* q1)* 57.3; // pitch ,转换为度数 Roll = atan2(-2 * q1 * q3 + 2 * q0 * q2, q0*q0-q1*q1-q2*q2+q3*q3)* 57.3; // rollv Yaw = atan2(2*(q1*q2 - q0*q3),q0*q0-q1*q1+q2*q2-q3*q3) * 57.3; //偏移太大,等我找一个好用的 }
这里使用的是匿名助手的上位机
https://gitee.com/anotc/AnoAssistant
有专门的通讯协议
串口通讯协议格式如下所示,需要注意传输为小端模式传输。
对应的源地址和目标地址分别为0xFD和0xFE。
我们只需要上报加速度和陀螺仪数据,所以功能码为0x01,数据长度为0x0D,需要主要为小端模式传输。
实测移动模块分别为X、Y、Z轴向下,可以看见数值基本上为1000mg。
加速度计测量线性加速度,而陀螺仪测量角旋转。为此,他们测量了科里奥利效应产生的力。
陀螺仪是一种运动传感器,能够感测物体在一轴或多轴上的旋转角速度。它能够精确地感测自由空间中复杂的移动动作,因此成为追踪物体移动方位和旋转动作的必要设备。与加速计和电子罗盘不同,陀螺仪不需要依赖外部力量(如重力或磁场),可以自主地发挥其功能。因此,从理论上讲,只使用陀螺仪就可以完成姿态导航的任务。
陀螺仪的每个通道检测一个轴的旋转。也就是说陀螺仪通过测量自身的旋转状态,判断出设备当前运动状态,是向前、向后、向上、向下、向左还是向右呢,是加速(角速度)还是减速(角速度)呢,都可以实现,但是要判断出设备的方位(东西南北),陀螺仪就没有办法。
MEMS陀螺仪主要利用科里奥利力(旋转物体在有径向运动时所受到的切向力)原理,公开的微机械陀螺仪均采用振动物体传感角速度的概念,利用振动来诱导和探测科里奥利力。
MEMS陀螺仪的核心是一个微加工机械单元,在设计上按照一个音叉机制共振运动,通过科里奥利力原理把角速率转换成一个特定感测结构的位移。
两个相同的质量块以方向相反的做水平震荡。当外部施加一个角速率,就会出现一个科氏力,力的方向垂直于质量运动方向,如垂直方向箭头所示。产生的科氏力使感测质量发生位移,位移大小与所施加的角速率大小成正比,科氏力引起的电容变化即可计算出角速率大小。
科里奥利效应指出,当质量 (m) 以速度 (v) 沿特定方向移动并施加外部角速率 (Ω)(红色箭头)时,科里奥利效应会产生一个力(黄色箭头),导致质量垂直移动。该位移的值与应用的角速率直接相关。
/* USER CODE BEGIN 2 */ printf("123123123"); lsm6dsv16x_reset_t rst; stmdev_ctx_t dev_ctx; /* Initialize mems driver interface */ dev_ctx.write_reg = platform_write; dev_ctx.read_reg = platform_read; dev_ctx.handle = &SENSOR_BUS; HAL_GPIO_WritePin(CS_GPIO_Port, CS_Pin, GPIO_PIN_SET); HAL_GPIO_WritePin(SA0_GPIO_Port, SA0_Pin, GPIO_PIN_RESET); /* Wait sensor boot time */ platform_delay(BOOT_TIME); /* Check device ID */ lsm6dsv16x_device_id_get(&dev_ctx, &whoamI); printf("LSM6DSV16X_ID=0x%x,whoamI=0x%x",LSM6DSV16X_ID,whoamI); if (whoamI != LSM6DSV16X_ID) while (1); /* Restore default configuration */ lsm6dsv16x_reset_set(&dev_ctx, LSM6DSV16X_RESTORE_CTRL_REGS); do { lsm6dsv16x_reset_get(&dev_ctx, &rst); } while (rst != LSM6DSV16X_READY); /* Enable Block Data Update */ lsm6dsv16x_block_data_update_set(&dev_ctx, PROPERTY_ENABLE); /* Set Output Data Rate. * Selected data rate have to be equal or greater with respect * with MLC data rate. */ lsm6dsv16x_xl_data_rate_set(&dev_ctx, LSM6DSV16X_ODR_AT_960Hz); lsm6dsv16x_gy_data_rate_set(&dev_ctx, LSM6DSV16X_ODR_AT_960Hz); /* Set full scale */ lsm6dsv16x_xl_full_scale_set(&dev_ctx, LSM6DSV16X_2g); lsm6dsv16x_gy_full_scale_set(&dev_ctx, LSM6DSV16X_2000dps); /* Configure filtering chain */ filt_settling_mask.drdy = PROPERTY_ENABLE; filt_settling_mask.irq_xl = PROPERTY_ENABLE; filt_settling_mask.irq_g = PROPERTY_ENABLE; lsm6dsv16x_filt_settling_mask_set(&dev_ctx, filt_settling_mask); lsm6dsv16x_filt_gy_lp1_set(&dev_ctx, PROPERTY_ENABLE); lsm6dsv16x_filt_gy_lp1_bandwidth_set(&dev_ctx, LSM6DSV16X_GY_ULTRA_LIGHT); lsm6dsv16x_filt_xl_lp2_set(&dev_ctx, PROPERTY_ENABLE); lsm6dsv16x_filt_xl_lp2_bandwidth_set(&dev_ctx, LSM6DSV16X_XL_STRONG); int16_t acc_int16[3] ={0,0,0}; int16_t gyr_int16[3] ={0,0,0}; float acc[3] = {0}; float gyr[3] = {0}; uint8_t data[21]={0}; data[0]=0xAB;//帧头 data[1]=0xFD;//源地址 data[2]=0xFE;//目标地址 data[3]=0x01;//功能码ID data[4]=0x0D;//数据长度LEN data[5]=0x00;//数据长度LEN 13 uint8_t sumcheck = 0; uint8_t addcheck = 0; int16_t angular_rate_raw[3]={0,0,0}; //pitch,roll,yaw uint8_t data_angular_rate_raw[16]={0}; data_angular_rate_raw[0]=0xAB;//帧头 data_angular_rate_raw[1]=0xFD;//源地址 data_angular_rate_raw[2]=0xFE;//目标地址 data_angular_rate_raw[3]=0x03;//功能码ID data_angular_rate_raw[4]=0x08;//数据长度LEN data_angular_rate_raw[5]=0x00;//数据长度LEN 8 data_angular_rate_raw[6]=0x01;//mode = 1 data_angular_rate_raw[13]=0x00;//FUSION _STA:融合状态 /* USER CODE END 2 */
/* Infinite loop */ /* USER CODE BEGIN WHILE */ while (1) { lsm6dsv16x_data_ready_t drdy; /* Read output only if new xl value is available */ lsm6dsv16x_flag_data_ready_get(&dev_ctx, &drdy); if (drdy.drdy_xl) { /* Read acceleration field data */ memset(data_raw_acceleration, 0x00, 3 * sizeof(int16_t)); lsm6dsv16x_acceleration_raw_get(&dev_ctx, data_raw_acceleration); acceleration_mg[0] = lsm6dsv16x_from_fs2_to_mg(data_raw_acceleration[0]); acceleration_mg[1] = lsm6dsv16x_from_fs2_to_mg(data_raw_acceleration[1]); acceleration_mg[2] = lsm6dsv16x_from_fs2_to_mg(data_raw_acceleration[2]); // printf("Acceleration [mg]:%4.2f\t%4.2f\t%4.2f\r\n",acceleration_mg[0], acceleration_mg[1], acceleration_mg[2]); } /* Read output only if new xl value is available */ if (drdy.drdy_gy) { /* Read angular rate field data */ memset(data_raw_angular_rate, 0x00, 3 * sizeof(int16_t)); lsm6dsv16x_angular_rate_raw_get(&dev_ctx, data_raw_angular_rate); angular_rate_mdps[0] = lsm6dsv16x_from_fs2000_to_mdps(data_raw_angular_rate[0]); angular_rate_mdps[1] = lsm6dsv16x_from_fs2000_to_mdps(data_raw_angular_rate[1]); angular_rate_mdps[2] = lsm6dsv16x_from_fs2000_to_mdps(data_raw_angular_rate[2]); // printf("Angular rate [mdps]:%4.2f\t%4.2f\t%4.2f\r\n",angular_rate_mdps[0], angular_rate_mdps[1], angular_rate_mdps[2]); } if (drdy.drdy_temp) { /* Read temperature data */ memset(&data_raw_temperature, 0x00, sizeof(int16_t)); lsm6dsv16x_temperature_raw_get(&dev_ctx, &data_raw_temperature); temperature_degC = lsm6dsv16x_from_lsb_to_celsius( data_raw_temperature); // printf("Temperature [degC]:%6.2f\r\n", temperature_degC); } IMUupdate(angular_rate_mdps[0]/1000,angular_rate_mdps[1]/1000,angular_rate_mdps[2]/1000,acceleration_mg[0]/1000,acceleration_mg[1]/1000,acceleration_mg[2]/1000); Roll=Roll*100; Pitch=Pitch*100; Yaw=Yaw*100; data_angular_rate_raw[8]=(int16_t)Roll>>8;//roll data_angular_rate_raw[7]=(int16_t)Roll; data_angular_rate_raw[10]=(int16_t)Pitch>>8;//pitch data_angular_rate_raw[9]=(int16_t)Pitch; data_angular_rate_raw[12]=(int16_t)Yaw>>8;//yaw data_angular_rate_raw[11]=(int16_t)Yaw; sumcheck = 0; addcheck = 0; for(uint16_t i=0; i < 14; i++) { sumcheck += data_angular_rate_raw[i]; //从帧头开始,对每一字节进行求和,直到 DATA 区结束 addcheck += sumcheck; //每一字节的求和操作,进行一次 sumcheck 的累加 } data_angular_rate_raw[14]=sumcheck; data_angular_rate_raw[15]=addcheck; HAL_UART_Transmit(&huart1 , (uint8_t *)&data_angular_rate_raw, 16, 0xFFFF); //printf("Roll=%.2f,Pitch=%.2f,Yaw=%.2f\n",Roll,Pitch,Yaw); HAL_Delay(100); /* USER CODE END WHILE */ /* USER CODE BEGIN 3 */ } /* USER CODE END 3 */
上报匿名助手能正常进行解析。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。