




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、2014高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白,在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、网上咨询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的问题。我们知道,抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料(包括网上查到的资料),必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺,严格遵守竞赛规则,以保证竞赛的公正、公平性。如有违反竞赛规则的行为,我们将受到严肃处理。我们参赛选择的题号是(从A/B/C/D中选择一项填写): A 我们的参赛报名号为(如果赛区设置
2、报名号的话): 01059029 所属学校(请填写完整的全名): 北京联合大学 参赛队员 (打印并签名) :1. 靳宝 2. 吴佳怡 3. 魏佳 指导教师或指导教师组负责人 (打印并签名): 张晓晞 日期: 2014 年 9 月 15 日赛区评阅编号(由赛区组委会评阅前进行编号):2014高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号(由赛区组委会评阅前进行编号):赛区评阅记录(可供赛区评阅时使用):评阅人评分备注全国统一编号(由赛区组委会送交全国前编号):全国评阅编号(由全国组委会评阅前进行编号)嫦娥三号软着陆轨道设计与控制策略摘要嫦娥三号任务令中华民族千年登月梦想成真,也标志
3、着中国航天技术再上一层楼,向将来开展载人探月工程又前进了一步。嫦娥三号在高速飞行的情况下,要保证准确地在月球预定区域内实现软着陆,关键问题是着陆轨道与控制策略的设计。为了解决此情况,因此正确计算着陆轨道近月点和远月点的位置以及航天器登月时关键的六阶段的速度和角度调整对嫦娥三号软着陆的轨道设计有着至关重要的意义。在处理问题一时,根据文章中对嫦娥三号从近月点到着陆点的运动状态和附件一中近月轨道示意图,可以将此运动看作近圆运动,利用万有引力定律和开普勒定律可以计算出嫦娥三号在近月点和远月点的速度。从文章中得知嫦娥三号在近月点15公里处以抛物线下降,相对速度从每秒1.7公里逐渐降为零,由此可见,嫦娥三
4、号的下降运动为类平抛运动,因此利用运动学公式可以解决这个问题。一、问题的提出与重述嫦娥三号是中国国家航天局嫦娥工程第二阶段的登月探测器,包括着陆器和月球车。它携带中国的第一艘月球车,并实现中国首次月面软着陆。嫦娥三号由着陆器和巡视探测器组成,进行首次月球软着陆和自动巡视勘察,获取月球内部的物质成分并进行分析,将一期工程的“表面探测”引申至内部探测。其中着陆器定点守候,月球车在月球表面巡游90天,范围可达到5平方公里,并抓取月壤在车内进行分析,得到的数据将直接传回地球。它将是中国发射的第一个地外软着陆探测器和巡视器(月球车),也是阿波罗计划结束后重返月球的第一个软着陆探测器,是探月工程二期(落)
5、的关键任务,起承上启下的作用。叶培建介绍,嫦娥三号探测器将突破月球软着陆、月面巡视勘察、月面生存、深空探测通信与遥控操作、运载火箭直接进入地月转移轨道等关键技术 嫦娥三号在高速飞行的情况下,要保证准确性地在月球预定区域内实现软着陆。其着陆轨道设计的基本要求:着陆准备轨道为近月点15Km,远月点100Km的椭圆形轨道先按椭圆轨道运行有两方面原因:1.节约燃料。嫦娥三号飞向月球的绝对速度是很大的,如果直接减速,需要很大的反推动力和燃料来使其减速,同时很危险使工作人员没有足够的时间来应对突发状况。椭圆形轨道可以起到缓冲的作用。2.确保成功,在椭圆形轨道运行,可以留出多余的时间来确定着陆点,着陆时间,
6、减速参数,等等。根据嫦娥三号的预定着陆点为19.51W,44.12N,海拔为-2641m(见附件1)。着陆准备轨道为近月点15km,远月点100km的椭圆形轨道;着陆轨道为从近月点至着陆点,其软着陆过程共分为6个阶段(见附件2),要求满足每个阶段在关键点所处的状态;尽量减少软着陆过程的燃料消耗的要求,回答下列问题:(1)确定着陆准备轨道近月点和远月点的位置,以及嫦娥三号相应速度的大小与方向。(2)确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。(3)对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。二、问题的分析和解决1.确定着陆准备轨道近月点和远月点的位置,以及嫦娥三号相应速度
7、的大小与方向月球软着陆,是指月球着陆器经取地月转移到达月球附近后,在制动系统的作用下以很小的速度近乎垂直地降落到月面上,以保证宇航员的安全和实验设备的完好。控制要求就是要求当登月探测器高度接近地面时,竖直方向上的速度和加速度也应基本为零,同时在水平方向上也不能有很大的速度和加速度。登月探测器的着陆控制包含竖直方向的运动控制,水平方向的运动控制,以及姿态控制等等着陆器进入着陆准备轨道的近月点是15KM,远月点是100KM。近月点在月心坐标系的位置和软着陆轨道形态共同决定了着陆点的位置。首先进入主减速段的区间是距离月面15km到3km。该阶段的主要是减速,实现到距离月面3公里处嫦娥三号的速度降到5
8、7m/s。快速调整段的主要是调整探测器姿态,需要从距离月面3km到 2.4km处将水平速度减为0m/s,即使主减速发动机的推力竖直向下,之后进入粗避障阶段。粗避障段的范围是距离月面2.4km到100m区间,其主要是要求避开大的陨石坑,实现在设计着陆点上方100m处悬停,并初步确定落月地点。(1) 开普勒定律开普勒第一定律,也称椭圆定律:每一个行星都沿各自的椭圆轨道环绕太阳,而太阳则处在椭圆的一个焦点中。开普勒第二定律,也称面积定律:在相等时间内,太阳和运动着的行星的连线所扫过的面积是相等的。这一定律实际揭示了行星绕太阳公转的角动量守恒。用公式表示为开普勒第三定律,也称调和定律:各个行星绕太阳公
9、转周期的平方和它们的椭圆轨道的半长轴的地方成正比。由这一定律不难导出:行星与太阳间的引力与半径的平方成反比。这是牛顿的万有引力定律的一个重要基础。用公式表示为。(2) 万有引力任意两个质点有通过连心线方向上的力相互吸引。该引力大小与它们质量的乘积成正比与它们距离的平方成反比,与两物体的化学组成和其间介质种类无关。, 其中(牛顿每平方米二次方千克)根据以上的分析,建立以月球赤道平面为平面,月心为原点、为月心与零度经线和零度纬线交线的交点的连线,为极轴(月球的极轴),与和满足右手标架,建立空间直角坐标系(如图5-1所示)。 图5-1 卫星绕月轨迹及软着陆轨迹 由于着陆点在球面上且近月点与远月点是由
10、月球的经度、纬度及高度唯一确定,在此为了便于计算 将极坐标转化为空间直角坐标,并代数题中相关数据,反解出经度。极坐标转化为空间直角坐标即: (5.1.1) (5.1.2) 距离公式: (5.1.3) 其中:为纬度;为经度;为嫦娥三号距月心的距离;为嫦娥三号距着陆点的距离;根据能量守恒、开普勒第二定律(面积定律),建立以下模型即: (5.1.4) 则近月点的速度,近月点的速度: (5.1.5) 其中:为卫星的质量,为海拔高度,近月点距月球表面的距离;,月球半径,远月点距月球表面的距离, 月球重力加速度, 近月点的速度, 近月点的速度。根据题目所给数据以上分析,可知: 将以上数据代入(5.1.1)
11、式可得,着陆点及近月点的空间直角坐标分别为: (5.1.6) (5.1.7)再将(5.1.6)式和(5.1.7)式代入(5.1.3)式可得关于与(近月点和着陆点距离)的函数,?利用Mathematica 5.0编程求解可得:-139.107近月点与远月点的速度方向,即为相应速度在轴与轴方向上的投影(如图5-2所示)图5-2 近月点与远月点的速度方向示意图2.确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略探月器在月球引力场中的动力学方程用极坐标表示为:第一次制动时,开普勒轨道的轨道方程探月器的径向速度和横向速度为:探月器制动变轨的开始和结束位置为待优化参数,设其极角分别为1和2。当给出 一组1
12、和2时,由式(2)可得探月器的位置和速度在制动变轨的初终值条件。探月器质量的初值m0为已知,终值自由。综合起来则可得到式( 1) 的边界条件。第二次制动时,变轨的初始条件与第一次类似,可由开普勒轨道的性质得到。终端条件为:r=1738km,vr =0,v=0 (3)将方程(1)写为以下形式该式称为系统的状态方程,式中X=x1,x2,x3,x4 T=r,vr,v,m T,称为状态变量。设状态变量的初、终值为X(H1)=X1和X(H2)=X2,由上一节分析可得变轨的边界条件:按照燃料最省的要求,取变轨结束时探月器质量的负值为性能指标,即:引入辅助的Hamilton函数:式中 K = 1,2,3,4
13、T,称为协态变量。根据Pontryagin极大值原理可以得到如下最优条件。协态变量满足协态方程:最优推力方向角:横截条件:存在一常向量i使得i 在端点满足:由于变轨的初始位置和终端位置不定,则有: 最优控制问题可以通过求解式( 4) ,( 5) 和( 7) ( 11) 得到,这是一个复杂的两点边值问题。(2) 、模型求解与验证应用遗传算法求解最优控制问题首先将上述最优控制问题转化成一个最优化问题:寻找最优参数1 ,2,1 ( 1 ) ,2 (1 ) ,3 (1 ) ,4 (1 ) ,使状态变量的终值满足边界条件,并使性能指标取最小值。然后引入“ 伴随 控制” 变换,利用协态变量与控制变量存在的
14、依赖关系, 用控制变量 u ( 1) 及其导数 u( 1) 的初始值得到协态变量的初值 1(1 )4(1), 从而使待优化参数的数目从 6 个减少到 4个,减少了运算量。采用实数编码的遗传算法,把1,2,u(1),u(1)作为个体子串,当给出一个初值时,由“伴随控制”变换和初始条件,得到协态变量和状态变量的初值。然后求解初值问题式(4),(8),同时由式(9)确定 u( ),于是=2时,得到状态变量的终值x1(2),x4(2),选择适应函数为:F(1, 2,u(1),u(1)= -x4(2)+Cû g2(X2)û,(12)式中C为自定义的权系数。通过遗传操作进行迭代,当F取
15、极小值时得到最优的*1,*2,u*(*1),u*(*1),解初值问题式( 4 ), (8)可以得到最优轨线x*i ()和*i (),再由式(9)可求得最优控制规律*()。设探月器质量 m= 150 k g,发动机推力F=500N,比冲V=400 s。探月器从200km高的月球停泊轨道出发,经过第一次变轨制动,减速到近地点为1738km(即月球半径),远地点为1938km的椭圆轨道上。在接近月面附近时进行第二次制动,实现在月面的软着陆。计算结果如表1所示。两次制动,发动机推力方向随时间的变化曲线分别如图2,3所示:以燃料消耗最少为目标,提出了通过两次制动来实现探月器软着陆的方案。 采用 Pont
16、ryagin极大值原理,分析了燃料最省变轨控制问题,得到了发动机推力方向角的最优控制规律。 在数值计算上,首先将上述最优控制问题转化成为一个参数最优化问题,然后应用遗传算法搜索最优参数,把最优控制问题的数值计算转化成对参数最优化问题的寻优计算,避免初值猜测,提高计算效率。这一方法适用于其它变轨控制问题。3.对于设计的着陆轨道和控制策略做相应的误差分析和敏感性分析在开始的时候让探测器自由下落,当下落到某种状态时启动发动机,以最大推力对火箭减速。当速度减小到零时正好到达地面。所谓月球软着陆。是指月球着陆器经地月转移到达月 球附近后。在制动系统的作用F以很小的速度近乎垂直地 降落到月面上,以保证宇航
17、员的安全和试验设备的完好。转化为具体的控制要求就是要求当登月探测器高度接近地 面时,竖直方向上的速度和加速度也应基本为零,同时在水 平方向上也不能有很大的速度和加速度。登月探测器的着 陆控制包含竖直方向的运动控制,水平方向的运动控制,以 及姿态控制等等,将在F面分别进行说明。 理想情况F的最优控制方案流程可简单描述为开始时 让探测器自由下落,当下落到某个状态时启动发动机,以最 大推力对火箭减速。当速度减小到零时正好到达地面。然而 这种控制方案对切换时机的准确程度要求极高,同时不允许 存在任何形式的干扰,而这在实际环境中是不可能实现的。 虽然没有采用理想情况下的最优控制方案,但是借鉴其 中开关控
18、制的思想。如果直接对加速度实行开关控制,势必 要频繁开关发动机,而这样做不仅有很大困难而且会很大程 度上减少发动机的寿命,降低系统稳定性和安全性。为了避 免上述情况的发生,将加速度的导数作为开关控制的变量, 这样不仅可以能快速做出反应,还能够使发动机更加平稳的 工作。这是一种采用非线性变结构控制与状态反馈相结合的采用之前提到的采用非线性变结构控制与状态反馈相结合的控制方法。以竖直方向的控制为例,假设登月探测器(包括燃料在内)的质量为m,偏二甲肼燃料燃烧后喷出气体 相对于探测器的速度为k,则喷气发动机产生的推力为棚, 其是燃料消耗的变化率,受不等式 则由它产生的加速度则为又设x为登月探测器离月球
19、表面的高度;则登月探测器垂直运动速度为,又已知月球重力加速度常数g;根据牛顿第二定律,可列出火箭运动方程式,即 根据运动学方程构建状态空间表达式,以竖直向上为正方向,选取状态变量为 (即登月探测器高度(即登月探测器竖直方向速度)、(即 登月探测器竖直方向的加速度) 并设,则系统状态方程为其中初始状态为,火箭终态为,。我们的任务是,根据初步确定的方案选择合适的u(f),使火箭从初始状态,转移到终态控制幅度u(f)需要在系统偏离工作点状态的正负绝对值大小的范围内取值,即,并且使 性能指标极,J、,其中三阶方阵的任意元素 根据性能指标,控制系统对应的哈米尔登函数为 (3) 根据最优控制的极小值原理可
20、知即当 时,达到最优控制由上式可知 又根据:, T 根据最优控制理论的相关知识可知,其中P为矩阵黎卡提方程的解。 将带入控制量u的表达式可得当时,故上式变为 ,即 ,得到了控制量u的表达式 i之后,我们就得到了一个完整控制系统的数学模型,接下来的工作便是通过软件仿真来检验模型的实际可行性。另外在确定e表达式之后,还可以进一步通过公式r计算登月探测器的发动机的实时推力大小, 以及公式三、相关数据的计算和测试1控制过程和仿真结果选择Matlab附带的Simulink作为仿真软件。仿真分为 两步,第一步为暂不各种考虑干扰因素并且假设月球重力加 速度恒定、在较理想的情况下对登月过程的仿真,如图l 所示
21、:图1 为这种条件下基于我们所设计的控制方案的Simulink 仿真程序结构图,其中三个积分器的输出结果从左到右分别 为加速度、速度和高度,上方为状态反馈,下方为变结构控 制的输入部分。 暂不考虑干扰因素并且假设月球重力加速度恒定时的 仿真结果,如图2所示:图2 上图为降落过程中登月探测器高度的仿真曲线,可以看 到下降过程非常的平滑,如图3所示:降落过程速度的仿真曲线显示,速度也没有发生跳变、 并且随高度的减小而减小。降落过程加速度的仿真曲线,如图4所示:可以看出整个控制过程分为3步:首先是开始降落后不 久,探测器处于自由落体状态;然后在50-800秒之间,加 速度竖直向上,探测器减速下降;之
22、后发动机加速度不断切 换直至平稳落地,并在高度小于一定阈值后关闭发动机。从 图像中我们也能看出开关控制的特征。 然而登月探测器在实际降落过程中,由于外太空陨石的 撞击之类的因素,将不可避免的受到各种类型的干扰;另外飞 船的起始降落高度为15千米,在计算重力加速度时需要考虑 探测器高度的影响。考虑上述因素后修改过的Simulink仿 真模型,如图5所示:这里对各个功能模块进行了封装,并且增加了模拟实际 干扰的模块以及重力加速度反馈模块。 下面是我们设计的模拟干扰:假设登陆器在t-400s时 受到一颗陨石撞击,导致高度突然下降100米,向下的速度 大小增大20ms,即对分别对两个物理量的输出在同一
23、时间 增加相应大小的脉冲信号作为干扰,仿真结果,如图6所示:针对目前社会上存在的利用手机短信进行违法活动的 现象,本文提出并设计一个基于贝叶斯分类算法和多模式串 模糊匹配算法的不良短信甄别混合模型。通过对实验数据的 分析,此模型对短信的分类识别具有较高的正确率。 在贝叶斯分类算法模块,本文采用TFIDF方法进行特 征提取,降低了算法的复杂性,使得算法的速度和效率都能 达到不错的效果。贝叶斯分类器在是否不良短信的分类上, 达到了较高的准确率。模型总体效果上,改进的多模式串模 糊匹配算法能够有效地应用于不良短信分类模块。总体实验 证明了该模型的有效性。上图为降落过程中登月探测器高度的仿真曲线,可以
24、看 出与理想情况相差无几。 降落过程速度的仿真曲线,如图7所示:图7 能够看到系统很快地消除了干扰对速度带来的影响,如 图8所示:上图为降落过程加速度的仿真曲线,也能够看出飞船在 在遭受撞击后立刻加大了发动机的推力从而抵消了干扰的 影响。2.总结根据以上分别在理想情况下和存在干扰的情况下 Simulink软件的仿真结果,可以看出这套控制方案不仅能够 控制登月探测器实现软着陆的既定目标,而且对于各种形式 干扰也能够很有效的进行抑制并且快速恢复到正常工作状 态,在理论上能够胜任在外太空环境卜的各类任务,系统的 快速性、稳定性、鲁棒性等都能够满足航天工程所需的要求。 总的来说本方案有一定的实际价值和
25、应用前景,但由于项目 类型和时间的关系,仅仅是在理论上初步进行了分析和验证 工作。亟需进一步的研究和实践。3.月球最优软着陆方案软着陆方案14J是,首先将探月器射入一个约100 km高度的环月停泊圆轨道:满足一定条件后,向飞行器施加一 制动脉冲,使飞行器进入100 kmxl5 km的椭圆轨道:当下降 到约15 km高度的近月点时,制动发动机点火,开始软着陆。探月器质心运动方程组在许多文献中都有介绍ff-6q,为便于比较,这里采用文献【2】的模型。忽略月球自转、月球引力非球 项、13月引力摄动等影响探月器的质心方程组为式中,r为着陆器距月心矢径;秽为着陆器在,方向上的速度; 0为着陆器环绕月球表
26、面的航程角:是航程角的角速度;m为着陆器质量;弘为月球引力常数;F为制动推力器的推力, 为制动推力器的比冲;砂为推力方向角,即推力方 向与当地水平面的夹角,取锐角,向上为正,向下为负。最优软着陆轨道设计的目的是寻找最优控制, ,使性能指标取最小值并且满足软着陆条件式中,t。为软着陆的初始时刻。定义o=0,该时刻的状态参数由椭圆轨道的近月点确定;t为着陆时刻,tf未知;R。为月球半径,Rm=1738 km。根据Pontryagin最大值原理川。无奇异情况下,推力应为开关控制:或者以最大推力工作,或者为零。理论上,需分析 切换点,但为了简化问题。采用常值推力假设闻即认为制动发动机一直以最大推力工作
27、。这样一方面有利于优化,另一 方面呵以降低发动机复杂性。另外,部分间接法的算例也验 证r该假设的合理性嘲。凶此,本文采用常值推力假设,即取4.月球运动方程的积分变换月球最优软着陆问题是一类终端时间自由型且受终端约束的最优控制问题。对于这类问题,一种方法是将终端时 问作为设计变量进行优化但这样做会加大计算量,甚至难以收敛,不利于快速优化;另一种方法是通过引入能量变量替换积分变量将其转化为终端积分变量固定型最优控制问题。但在月球最优软着陆时,能量并不是均匀变化的。在着陆点附近能量变化率趋于零以能量为积分变量的数 值误差较大而微小的能量误差将会带来很大的状态误差。并且,优化迭代过程是未知的。当沿着非
28、最优轨迹时,能量不一定是单调变化的这将会导致基于能量的积分无法进行。本文选择状态量作为积分变量。这样只要推力方向角在(一900,900)范围内的单调性就有保证,且其变化较为均 匀。为了使得积分变量单调递增,引入变量则则将方程组(1)两边同时除以,则转换为对积分, 转化后的方程组为式中为了得到各状态量随时间的变化,需增加微分方程如下:以上两式组成新的状态方程,即,所以实际计算时可以转变目标函数为。变换后的约束条件为:式中,;下标0和f表示初值和终端值, 下同。 至此,原终端积分变量不确定型最优控制问题转化为终端积分变量固定型最优控制问题。转化后的问题一方面更适合于优化数值算法求解;另一方面终端约
29、束条件减少了一 个,终端约束更容易满足,收敛速度更快。并且,转换过程也较为简单。5.模型评论建立的模型方简单易行且易于应用与现实生活,模型具有坚实可靠的数学基础。在建立模型时,考虑的影响因素较少,在处理实际问题时可能存在一些误差,数据具有一定的局限性,考虑的情况比较简单。我们会吸取这次建模的经验与教训以便在以后的学习及生活中有更好地发挥。参考文献:【l】 尹怀勤我国的嫦娥工程和探月卫星天津科技J】 2007,(2)5石 【2】 王鹏基,张熵,曲广吉月球软着陆飞行动力学和制导 控制建模与仿真中国科学E辑:技术科学叨 2009,39(3):521527 【3】 凤凰网嫦娥三号详细资料【EBOL】(2009-0929) 【201 1-09-20 http:apptechifengcombaikeintrophvl?nameAE5A BA6E5A8A5E4B889E58FB7 【4】 新浪网探月工程首席科学家称
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 消费贷款购车合同(3篇)
- 2025年棉花加工成套设备项目合作计划书
- 理财顾问实习报告范文
- 2025年饲料营养型添加剂项目发展计划
- 2025年特种丝制品项目合作计划书
- 教育技术终身学习的助推器
- 2025年浙江省杭州市杭州二中物理高二下期末质量检测试题含解析
- 智慧城市管理与服务的数字化转型之路
- 国际合作在提升教育国际化水平中的贡献
- 专题04 读后续写精彩结尾及主题升华仿写(测试)原卷版-2025年高考英语二轮复习
- 中石化夏季八防培训课件
- 超星尔雅学习通《红色经典影片与近现代中国发展(首都师范大学)》2025章节测试附答案
- 2024届高三生物学科高考备考经验交流与反思
- 2025年河北轨道运输职业技术学院单招职业技能考试题库及答案1套
- 腰椎间盘突出的诊治课件
- 煤矿工作申请书
- 医疗护理医学培训 简易呼吸气囊的介绍及使用课件
- 加油站的运营数据分析
- 《典型生物质颗粒的安全性能分析综述》2200字
- IATF 16949 质量管理手册
- 燃气安全培训课件
评论
0/150
提交评论