一种改进的中途修正变结构控制方法.docx

上传人:rrsccc 文档编号:9305590 上传时间:2021-02-17 格式:DOCX 页数:8 大小:248.59KB
返回 下载 相关 举报
一种改进的中途修正变结构控制方法.docx_第1页
第1页 / 共8页
一种改进的中途修正变结构控制方法.docx_第2页
第2页 / 共8页
一种改进的中途修正变结构控制方法.docx_第3页
第3页 / 共8页
一种改进的中途修正变结构控制方法.docx_第4页
第4页 / 共8页
一种改进的中途修正变结构控制方法.docx_第5页
第5页 / 共8页
点击查看更多>>
资源描述

《一种改进的中途修正变结构控制方法.docx》由会员分享,可在线阅读,更多相关《一种改进的中途修正变结构控制方法.docx(8页珍藏版)》请在三一文库上搜索。

1、一种改进的中途修正变结构控制方法晁宁 1,罗晓英 2CHAO Ning1, LUO Xiaoying21.中国空间技术研究院西安分院,西安 7101002.西安电力电子技术研究所,西安 7100611.China Academy of Space Technology Xian Branch, Xian 710100, China 2.Xian Power Electronic Research Institute, Xian 710061, ChinaCHAO Ning, LUO Xiaoying. Improved variable structure control method in

2、midcourse correction. Computer Engi- neering and Applications, 2014, 50(5):265-270.Abstract:Combining with characteristic of small midcourse correction level of transfer orbit, approximate error linear model of restricted 3-body dynamic problem in the earth-moon system is established. On this founda

3、tion, variable struc- ture adjuster theory and model reference variable structure control theory of linear system are used to deduce midcourse control law separately. The optimal choosing method for control sub matrix b2 under different orbit injection error level and calculating method for input ve

4、ctor of reference model under different state are proposed. Simulation experiments compare the performance of standard orbit and re-designed correcting orbit under two kinds of variable structure theories, and analyze affection on control effect brought by r corresponding to orbit state under differ

5、ent characteristic points. The experimental results show validity of proposed methods.Key words:libration point; transfer orbit; midcourse correction; variable structure control摘要:结合转移轨道中途修正量级较小的特点,建立了地月系限制性三体动力学问题的近似误差线性模型。在此 基础上,分别利用线性系统的变结构调节器理论与模型参考变结构控制器理论推导了对应的中途修正控制律。提 出了控制子矩阵 b2 在不同入轨误差等级下的最

6、优选择方法和不同状态下参考模型输入向量 r 的计算方法。仿真实 验对比了分别利用两种变结构理论进行控制时,重新规划出的修正轨道与理论轨道的性能指标,并分析了选择不同 特征点处轨道状态对应的 r 对控制效果的影响情况。实验结果验证了所提方法的有效性。 关键词:平动点;转移轨道;中途修正;变结构控制文献标志码:A中图分类号:V412.4doi:10.3778/j.issn.1002-8331.1205-01801引言从 20 世纪 60 年代美苏两国开展月球及深空探测开 始,国内外许多学者开始转移轨道中途修正方面的深入 研究。郗晓宁1 等以中途制导算例为基础,分析了摄动 法和显示制导法的修正策略,

7、并进行了统计学分析;周 文艳2等研究了修正速度增量同入轨误差和修正时刻间 的关系。Serban3等应用动力系统理论和最有控制方法 给 出 了 地 月 平 动 点 的 中 途 修 正 策 略 。 Breakwell4 提 出 的用于探测器进行行星际飞行的间距比策略。其基本思想是尽快进行首次修正,以后每次修正的时间间隔与 前一时间间隔为常数。利用变结构控制实现中途修正 通常得到的是连续小推力控制律,由于其算法繁琐、硬 件实现复杂及对执行机构要求高等特点,目前各国对于 这方面的工程应用尚不多见。而将滑模变结构控制应 用于飞控系统和伺服控制系统有许多研究成果。针对 变结构系统的能耗与振颤,洪剑波5 提

8、出了一种新的负 指数分布变参数滑模变结构控制,以自适应调节边界层 系 数 、滑 动 模 系 数 和 控 制 增 益 系 数 等 变 结 构 控 制 器 参基金项目:国家自然基金(No.61174204)。 作者简介:晁宁(1982),男,工学博士,工程师,研究方向为星载天线控制系统设计及轨道预报;罗晓英(1984),女,工学硕士。E-mail:收稿日期:2012-05-21修回日期:2012-06-26文章编号:1002-8331(2014)05-0265-06CNKI 网络优先出版:2012-08-16, http:/ 采用变结构控制方法设计了控制器,并 利用干扰观测器估计了摩擦力矩并自适应

9、地抑制了其 对控制系统的影响。本文以地月系限制性三体动力学问题的近似误差 线性模型为基础,分别推导了基于线性系统变结构调节 器理论与模型参考变结构控制器理论的中途修正控制 律。选择转移轨道初始点为中途修正时刻,在不同入轨性系统可近似化为线性系统8:Z = AZ + bu(2)方便起见,下面均用 X 代替 Z 。4变结构控制调节器(Variable Structure Con- trol Adjuster,VSCA)从前面近似误差线性化方法中知道,限制性三体动误 差 等 级 条 件 下 ,重 新 规 划 出 满 足 终 端 条 件 的 修 正 轨道,并对比了修正轨道与标称轨道的性能指标。针对参力

10、学模型在地月 L1点至月球的转移轨道入轨点处进行考输入和控制矩阵取值对控制效果的影响,提出了控制 子矩阵 b2 在不同入轨误差等级下的最优选择方法和不一阶泰勒展开,并取状态误差 X = X - X ,标准状态为 X。令 X = X。原系统可以写作如下线性系统 X = AX + Bu ,或 同状态下参考模型输入向量 r 的计算方法,对比了选取 不同终端条件下的轨道控制情况。xyz 00010 0 x 00001 0 yI 00000 1 0X = + 2y = 02 0 z + 3u(3)xXXXYXZx 3y - 2xYXYYYZ -2 0 0 y2圆限制性三体模型(Circular Rest

11、ricted Three-Body Problem,CRTBP) z ZXZYZZ 00 0 z 对于圆限制性三体问题,可描述为一个质量忽略不 计 的 小 天 体 P 在 两 个 相 互 环 绕 运 动 的 大 天 体 P1、P2这是线性系统式的简约形式,可根据控制量数目写作分块形式:A11A120(P1 P2) 引力作用下的运动状态,通常利用以两大天体A = AA ,B = b (4) 2122 2质 心 为 圆 心 旋 转 的 会 合 坐 标 系 或 称 旋 转 坐 标 系 来 研究。假定 P 不影响 P1、P2 的运动,两大天体共同绕其质b2 为 m m 非 奇 异 矩 阵 ,A11为

12、(n - m) (n - m) 矩 阵 ,X xx 1心做角速度为 的圆周运动。 P1 指向 P2 的方向为会合坐标系的 x 轴, 方向为 z 轴,y 轴与 x、z 轴成右手X = X ,X21 = y ,X z2 = y 。这是变结构系统的简约型。 z 系。这种假设下的动力学模型不考虑摄动因素,即无摄 圆限制性三体模型,小天体的无量纲化运动对应一个二 阶常微初值问题7,其分量形式为:广义上讲,变结构系统主要有两类:一类是具有滑动模态的变结构系统;另一类是不具有滑动模态的变结 构系统。通常说的变结构系统均为指前者。这是由于x - 2y = x = x - (1 - ) x + r31 y- x

13、 + - 1r32 y变结构系统对外界干扰和参数摄动的鲁棒性完全通过滑动模态来实现,变结构控制系统中滑动模态的设计很 关键,它决定了最终的控制效果。设计时保证滑动模态rry + 2x = y = y - (1 - ) 3 -312(1)运动稳定;且具有良好的运动品质。通常利用闭环系统 z z 特征根配置要求进行设计,即滑动模态的极点配置法。z = z = - (1 - ) r 3 - r 312其中,(x y z) = (x2 + y2)/2 + U (r r )引力势函数U (r r ) =本文利用理论轨道初始入轨点来构造近似线性系统闭环极点集,并通过闭环极点配置方法构造滑动超平面。1 21

14、 22(1 - )/r1 + /r2 = m2 /(m1 + m2) r1 = (x + )2+ y22+ z2 是航2对上式来说,选择滑动超平面为:S(X t) = CX(5)天 器 到 大 天 体 的 距 离 ;r2 = (x + - 1) + y + z器到小天体的距离。是 航 天式中 C = C1 C2 Rm nC RC Rm (n - m) 21m m且满秩。3线性化用于深空探测的无摄圆限制性三体模型表现出很 强的非线性,在考虑摄动因素后这种非线性特性就更加 明显。而在一定时间范围内,当实际轨道与标准轨道的 状 态 误 差 较 小 时 ,其 误 差 系 统 模 型 具 有 较 强 的

15、 线 性 特 性。这种线性化方法将误差状态作为新的状态变量,利 用多元向量函数的雅克比矩阵作为线性系统的系数矩令 S = 0 ,则系统滑动模态方程为:X1 = (A11 - A12 K )X1(6)极点配置增益阵 K = C -1C ,K Rm (n - m) 。21考虑到线性系统模型的近似性、未建模动态及天体 摄动力等因素影响,近似线性系统写作不确定性多变量 系统8:X (t) = A + DA(t)X (t) + B + DB(t)u(t) + Df (t)(7)式中状态变量 X Rn 控制量 U Rm 外界干扰 f Rl。阵 A。设状态变量 X = x y z x y z Tu = ux

16、 uyu T另外 A Rn n、B Rn m、D Rn l、DA、DB 分别为标称z晕轨道标准状态为 X ,状态误差 X = Z = X - X ,则非线系统矩阵、标称控制矩阵、扰动分配矩阵以及 A 和 B 的摄动矩阵。假设该系统满足如下条件:(1)(A B) 为完全可控对;系统具有简约标准型,即BT = 0BT ,且 BT 为非奇异。(1)Am Bm 均有界;(2)(Am Bm ) 为可控对,且 rank(Bm ) = l m 。 要使得被控对象完全跟踪参考模型,应满足:lim e(t) = 0(14)22t (2)DA DB 对于任意 t 均连续且有界,DA DB在 上一直有界。(3)参数

17、摄动及外界扰动上界均已知,即 DA a DB b Df f 其 中 a b f 为 已 知 正 常 数 , 为矩阵或向量的无穷范数。于是,取变结构控制律 u 如下形式:u = -g(t)(CB)-1 sgn(S)(8) 其中:式中 e(t) = Xm (t) - Xp (t) 为系统误差状态变量。 选择滑动模态切换平面为:S = Ce(t)(15)式中 C = C1 C2 为滑动模型参数矩阵,并令 S = S = 0 同 样可以使用式(5)中 C 的极点配置方法来确定。计算方 法详见作者晕轨道极点配置控制方法一文中。CRTBP 近似线性系统方程的 A、B 具有简约标准形式,因此不 需要进行初等

18、行变换即有:-1Am 11Am 12g(t) = (1 - a3)a1 X + a2 + (9)mA = ,B= 0 (16)式中 为以小的正常数,其他系数为:Am 21Am 22pBp2-1a1 = CA + aC ,a2 = CD f ,a3 = bC (CB) 式中 Am 11 R(n - m) (n - m),Am 22 Rm m,Bp Rm m 。a ,b ,f 为已知正常数,且为系统不确定性的界,即根据完全模型跟踪匹配条件,模型匹配控制律为:-1-1 DA a , DB a , f f(10) 这里的矩阵或向量范数 均采用诱导范数,即 范数。5模型参考变结构控制方法(Model R

19、eferrence Variable Structure Control,MRVSC)从 前 面 的 分 析 中 可 以 看 到 ,将 受 控 系 统 的 初 始 误差、参数变化及非线性等因素作为扰动处理时,可以通u m = Bp20Im(A m - A p)X p + Bp20ImB mr (17)而 变 结 构 控 制 律 u v 的 选 择 应 通 过 使 得 Lyapunov 函数变化率 v = ST QS 0 为滑动模态存在的充分条件。参数系数问题。设计思想是通过匹配模型与变结构两部分输入 为 已知正 常数 ,如 前所述 ,均 为不确 定系 统的量来实现控制,即abfu = um +

20、 uv(11)其中 um 为模型参考闭环控制系统的匹配控制律,uv 为界,即 DA a, DB b, f f(21)变结构控制律。以 CRTBP 中 基 于 小 偏 差 的 误 差 线 性 系 统 为 例 ,参 照式(7)不确定性多变量系统,被控对象取为:X通过后面的实验验证可知,滑模变结构控制方法与模型参考变结构控制方法均可实现对处于转移轨道段 航天器的中途修正控制,使其沿着一条新的轨迹趋近目 标点。而其中模型参考变结构控制方法更加适用于对p (t) = Ap + DAp (t)Xp (t) + Bp + DBp (t)U (t) + Dp fp (t)(12)式中诸变量定义同式(7),下标

21、 p 表示受控系统,且系统 式满足第 4 章中系统的假设(1)(3)。而参考模型的数学描述可表示为一个确定型多变 量系统:X中途修正的问题描述,这种方法中需要对可控子系统的 控制矩阵取值进行确定。6输入子矩阵 b2 的确定由于近似线性三体动力学模型中控制矩阵的选择m (t) = Am Xm (t) + Bm r(t)(13)式 中 带 下 标 m 的 变 量 均 表 示 模 型 参 考 量 。 r(t) Rl 表示 了 参 考 模 型 一 致 有 界 的 外 部 输 入 量 ,且 B Rn l 。m模型参考式满足6:是自主的,不同元素取值会产生不同的控制效果。控制矩阵 B 符合变结构控制的简约

22、标准型,因此其中子矩阵 b2 元素就成为了主要确定对象。本文通过对入轨点起 始误差的分析,根据简约线性系统 A12 子块三轴比例关系来确定控制子矩阵中三通道元素的比例关系,达到了 更具针对性的控制效果。结 合 前 面 误 差 状 态 的 定 义 ,取 初 始 状 态 偏 差 为DX = X - XP ,其 中 X 为 理 论 状 态 向 量 ,Xp 为 实 际 状 态 向 量 ,并 用 X 代 替 DX 。 控 制 子 矩 阵 b2 元 素 为diag (b b b ) 根 据 式(4),A 子 矩 阵 为 R3 3 对 称21 2223 12阵,其中元素为:a11 a12a13以不加区分地处

23、理为具有一定界限的随机变量。8.1参数选择以地月系 L1 平动点晕轨道与月球间转移轨道为研 究对象。标准初始状态为 X (0) 实际状态向量为 Xr (0) 。 若在入轨点处对原系统进行线性近似,则将 X (0)、Xr (0) 代入式(3)可以分别得到参考模型和实际系统的系统矩 阵,取:X (0) = 0.861 863 27 0.004 395 03 0.006 098 44A = a12 21a22a 23(22)-0.154 212 51 0.001 423 41 -0.000 039 31Ta31 a32 a33Xr (0) = X (0) + DX (0)(26)由 于 三 轴 控

24、制 量 是 针 对 A12 子 矩 阵 对 应 数 据 进 行 DX (0) X为 随 机 向 量 且 有 界 X,仿 真 中 分 别 取 值a的,所以若取 b21 = 1 则 B2 可变换为 B2 = kb diag(1 b22 /b2110-5 、10-6 和 10-7 进 行 计 算 。 参 数 摄 动 为 DA ,b23 /b21) 。对应的比例关系为: DB , f 。333bfb21:b22:b23 = a1i:a2i:a3i i = 1 2 3(23)各参数及扰动取值如下:i = 1i = 1i = 1(1)a = 2.0736 ,b = 0.3 ,f = 0.06 。式中比例系

25、数 kb 用来决定控制的力度。通 过 这 种 方 法 ,可 以 实 现 MRVSC 参 数 的 优 化 配 置,提高控制效率。(2)开关函数中小的正常数取为 i = 10 。-4(3)扰动分配矩阵为:D = 037参考输入 r 的选择从前面的控制律中看到,参考模型中的输入向量 r也是影响控制效果的重要因素 。对于 CRTBP 来说,理diag(d1 d2 d3)-3其中扰动矩阵元素 d1 d2 d3 = 10-58.2性能对比1 2 1 。论模型中的输入必然是呈现非线性特性的时变量。而选择 X 分别为 10和 10-7对两种变结构控制方法MRVSC 系统的匹配控制律需要知道 r 的诱导范数,为

26、简化计算通常取 r 为常向量。这样一来,选择哪个特征 点处的 r 就成为了必需解决的问题。本文选择理论轨 道上的若干特征点为对象,计算出对应 r ,然后在仿真 中比较了不同 r 对应的不同的控制效果。转移轨道的起点、中点与终点具有一定代表性,因 此本文选择这三处作为特征点进行分析。不同特征点 对应了不同的终端条件,因此需要对以上三点为终端条 件的控制情况进行分析。在圆限制性三体问题框架下, 探测器偏离理论轨道可以近似认为两个主天体对第三进行比较。图 1(a)、(b)与(c)、(d)分别为 10-5 和 10-7 等 级下两种方法的轨道控制效果。从图中可见,入轨误差 级小的受控轨道与理论十分接近

27、,并且具有耗能小、飞行 时间相近的优点。且直观来看,MRVSC 方法与前者相 比,其受控轨道收敛速度较快,距离理论轨道距离较近。图 2 为模型参考变结构方法的控制加速度,可见在 经过一段时间的振荡后,三轴加速度均趋于收敛。下面 将两种相同误差等级下(10-5 和 10-7 )两种控制器 VSCA 与 MRVSC 的性能指标进行比较,见表 1。表 1 不同误差级下的指标XTVSTbtcAccdAcca10-50.082 41.290 91.276 10.168 12.482 310-70.111 11.208 51.161 00.223 72.233 210-50.673 61.290 91.3

28、06 90.241 75.492 610-70.036 11.208 51.220 90.041 00.361 9体的引力失衡所致。因此本文以探测器在空间受到的 天体引力加速度差为输入量,并取当前特征点对应特征向量为:跟踪器TX (k) = X1(k)X2 (k) X6 (k) =调节器x(k)y(k)z(k)x (k)y (k)z (k)T(24)式中 k 为序号,则该点输入向量 r(k) 按下式计算:224466r(k) = X (k) - X (k - 1) X (k) - X (k - 1) X (k) - X (k - 1)T (25)8仿真实验对于不确定性,包括参数摄动与入轨偏差。

29、参数摄表中均为无量纲数据,tb、tc、Ac cd、Acca 分别为标准 飞行时间、受控飞行时间、最大加速度和总加速度。跟 踪方差和(TVS)表示三轴分量方差在控制区间 T 内的总 和。TVS 按照式(28)计算:nTVS = var(X动的原因主要由环境条件变化引起,如大气密度、气压、i = 1i - Xconi)(28)温度、高度等;而带来入轨偏差的原因很多,如发动机推 力偏差、导航信息误差、系统随机噪声等。两者通常可式中 X 为标准轨道状态,Xcon 为受控轨道状态,n 为 T内 的 累 加 次 数 。 表 中 数 据 除 误 差 级 外 均 为 无 量 纲 量 。3z/103AU2101

30、23 0.050L1入轨点标准轨道 受控轨道z/103AU实际轨道 20月球20.050y/AU标准轨道 受控轨道 实际轨道L1入轨点月球1.000.05y/AU0.750.100.800.85x/AU0.900.951.000.050.10 0.750.800.850.950.90x/AU(a)VSCA(105)(b)MRVSC(105)2z/103AU1012 0.05入轨点L10月球0.85标准轨道 受控轨道 实际轨道1.000.950.9032z/103AU101230.050入轨点L1标准轨道 受控轨道 实际轨道月球0.95 1.00y/AU0.050.800.10 0.75x/AU

31、y/AU0.050.10 0.75 0.800.85 0.90x/AU(c)VSCA(107)(d)MRVSC(107)图 1 不同误差等级下两种变结构控制情况三轴偏差加速度(/ 2.730 75e-6kms2)0.50.40.30.20.100.10.20.30.4x 轴加速度 y 轴加速度 z 轴加速度三轴偏差加速度(/ 2.730 75e-6kms2)0.50.40.30.20.100.10.20.30.4x 轴加速度 y 轴加速度 z 轴加速度0.5 050100150时间 t/0.189 7d0.5 050100150时间 t/0.189 7d(a)误差等级为 105(b)误差等级为

32、 107图 2 MRVSC 加速度假 设 用 下 标“00”表 示 无 量 纲 量 ,d、V、a、t 分 别 表 示 距 离、速度、加速度和时间,则其与实际物理量对应关系为:5VSCA 的优势。8.3输入 r 选择在不同点时的影响1d00 = 3.844 01 10km 1V00 = 1.024 55 km/s模型匹配时采用不同点固定输入 r 的不同控制效-61a00 = 2.730 75 10km/s2 1t = 3.7516 107 s果 ,以 10-7为 对 象 ,分 别 选 择 为 :(1)转 移 轨 道 入 轨 点 ;00可以看到当入轨误差级别较大时,仅靠 VSCA 对多方面 扰动以

33、及入轨偏差进行调节的效果不很理想,能耗也大 幅增加。而在加入模型参数匹配控制律后再进行变结 构 控 制 的 跟 踪 精 度 有 了 很 大 提 高 ,控 制 能 耗 也 大 幅 减 少 。 由 于 匹 配 加 速 度 的 存 在 ,因 此 其 加 速 度 值 比 较 平 稳。可见,MRVSC 系统适用于对系统参数进行大范围 调 整 的 情 况 。 而 对 于 入 轨 偏 差 较 小 的 情 况(以 10-7 为 例),MRVSC 的 控 制 效 果 与 能 耗 方 面 并 未 体 现 出 对 于(2)转移轨道中点;(3)近月点,并对这三点对应的输入进行控制效果比较,结果见表 2。表 2 107

34、 等级下不同 R 的控制性能TVStcAccdAcca10.191 01.166 20.240 32.545 620.116 11.161 10.227 32.310 1 3 0.111 1 1.161 1 0.223 7 2.233 3 非理想 CRTBP 为一个非线性系统。对参考模型来说,需要一个时变引力加速度差来作为控制输入,选择 一个定常输入向量进行跟踪必然存在固有的缺陷。从 表 2 中可以看出,r 从起点至终点的选择过程中,4 个性 能指标均逐级递减,并且这些指标都是越小表示控制效 果就越好。于是可以推断,参考模型变结构控制中的 r 如果选择定常向量,那么选择轨道控制目标点将是控制

35、输入的最佳选择。9 结束语对 无 摄 CRTBP 在 L1 点 附 近 晕 轨 道 至 月 球 的 理 论 转移轨道入轨点进行了一阶泰勒展开,得到误差近似线 性系统。在此系统基础上,分别利用线性模型的变结构 调节器与模型参考变结构控制跟踪器理论推导了对应控制律,提出了控制子矩阵 b2 在不同入轨误差下的最 优选择方法和参考模型输入向量 r 的计算方法,并分析 了选择不同特征点终端条件对应的 r 对修正轨道的影 响情况。根据研究分析,由于控制律均是在小偏差假设 下近似线性模型基础上推导的线性控制律,在较大轨道 偏差、非线性及参数摄动等不利因素的影响下具有一定 局限。针对这个问题,对精确线性化方法

36、和非线性滑模变结构控制方法的分析将是课题的下一步研究方向。(上接 242 页)参考文献:1 Luck R,Ray A.Experimental verification of a delay com- pensation algorithm for integrated communication and control systemsJ.International Journal of Control,1994, 59(6):1357-1372.2 Walsh G C,Ye H,Bushnell G.Stability analysis of net-worked control syste

37、msJ.IEEE Transactions on Control Systems Technology,2002,10(3):438-442.3 Zhang W,Branicky M S,Phillips S M.Stability of net- worked control systemsJ.IEEE Control Systems Maga- zine,2001,21(1):84-99.参考文献:1 郗晓宁,曾国强,任萱,等.月球探测器轨道设计M.北京: 国防工业出版社,2001.2 周文艳,杨维廉.月球探测器转移轨道的中途修正J.宇航 学报,2004,25(1):89-92.3 Ser

38、ban R,Koon W S,Lo M W,et al.Halo orbit mission correction maneuvers using optimal controlJ.Automatica, 2002,38:571-583.4 Breakwell J V,Kamel A A,Ratner M J.Station-keeping of a translunar communication stationJ.Celestial Mechanics, 1974,10(3):357-373.5 胡剑波,王坚浩.一种变参数滑模变结构控制及飞行控制 应用研究J.飞行力学,2010,28(4):46-49.6 刘金琨,孙富春.滑模变结构控制理论及其算法研究与进 展J.控制理论与应用,2007,24(3):407-418.7 俞辉,宝音贺西,李俊峰.双三体系统不变流形拼接成的低 成本探月轨道J.宇航学报,2007,28 (3):637-642.8 晁宁,李言俊.基于分段连续推力的晕轨

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 社会民生


经营许可证编号:宁ICP备18001539号-1