|
关于我们
 新书资讯 新书推荐 |
空间飞行器GNC系统简化算法 
定 价：198 元
9 7 9 8 2 7 4 5 9 1 3 5 9  本书主要介绍了飞行器 GNC系统研制工程应用算法、建模和仿真验证,包括组合定姿、组合导航、姿态控制、轨道机动和深空探测 GNC等相关的算法、建模和仿真.空间飞行器轨道机动大都基于数学仿真,特别是超实时数学仿真,需要动力学快速迭代推算,基于挠性和晃动动力学模型解耦及局部迭代算法,可以大幅缩短动力学迭代计算时间.空间飞行器 GNC系统工程应用算法简单可靠,物理概念清晰,适合工程应用,可以作为飞行器导航、制导与控制专业研究生教材,也可以供从事飞行器 GNC 系统设计与研制的技术人员及高等院校相关专业的师生参考。 考虑到卡尔曼滤波算法相对复杂、物理概念不够清晰，提出基于PI滤波的陀螺与磁强计组合定姿，陀螺与星敏感器组合定姿，GNSS兼容机与加表组合导航，伪距及伪距率与加表组合导航，基于跟瞄和加表近程相对导航等组合定姿和组合导航算法，并进行了仿真和在轨验证。
王献忠，工学博士，专注于高精度指向控制、远近程相对导航、空间自主交会、在轨操控与服务等控制技术研究。从事深空探测飞行器研制，主要负责国家首次火星探测工程中的环绕器研制。 第1章概述 1.1空间飞行器GNC系统概述 1.1.1GNC系统含义 近地轨道、月球等载人探索和无人科学任务，火星等深空探测任务，都涉及GNC技术，包括自主交会对接（ARD）、自主精确着陆（APL）、进入下降着陆（EDL）等。GNC系统技术是空间探测关键需求和未来技术发展方向。 空间飞行器GNC系统涵盖制导、导航与控制，典型的GNC系统框图如图11所示，GNC是一个内部相互作用和耦合的系统，也是一个具有很强外部耦合作用的集成系统，耦合存在于飞行机械、电子设备、软件、推进、飞行科学、载荷等相关分系统，也与载人/无人自主飞行操控交互作用。 图11典型的GNC系统框图 制导：一个运动的飞行器从当前的位置/速度/姿态到一个期望的位置/速度/姿态状态的路径决策，同时满足规定的制约因素，如推进剂消耗、安全、动态热负荷和时间等约束。 导航：确定一个运动飞行器在指定参考系下的当前运动状态，包括姿态、位置、速度等。 控制：作用在一个运动飞行器上，能够稳定和调节飞行器的运动，使飞行器状态朝向期望导引状态的力和力矩控制指令，通常采用闭环控制方式。 1.1.2载人空间探索领域GNC GNC系统是确保飞行器在上升段、逃逸段、在轨和再入任务阶段可靠性和鲁棒性的关键分系统，以支持航天员安全和任务成功为目标，GNC系统相关技术包括：
1.1.3空间科学探测领域GNC 把人类的探索范围延伸到太空，通过地球轨道和深空放置天文台，飞行器探访月球和行星，探测器着陆、环绕和采样返回等，飞行器可以位于低地球轨道、地球静止轨道、高地球轨道、地月转移轨道、绕月轨道、地球行星转移轨道、绕行星轨道、拉格朗日点等，需要GNC具备高容错性和高可靠性，相关的GNC系统技术包括： 1）天文级精密望远镜捕获、跟踪和指向控制。哈勃望远镜指向精度要求为0.007，挠性振荡对韦伯望远镜影响更大，但韦伯望远镜指向精度约为哈勃望远镜的1/18，韦伯望远镜的指向精度比哈勃望远镜更高。 2）行星大气层减速，大气捕获，进入、下降和着陆控制，包括行星探测的定点着陆和自主障碍规避。 3）行星取样返回自主交会对接。 4）多飞行器精确编队飞行。涉及编队飞行敏感器和执行机构，编队飞行任务设计工具，地面验证测试平台和在轨技术演示验证。 5）推进剂在轨贮存和补加。 6）用于测量行星结构和自然资源的超高精度引力敏感器。 7）优于10 cm精度的激光深空测距系统。光学跟踪技术将用于深空探测网，以实现更高精度的三维测量。 8）先进的控制系统技术。研制安全可靠的新一代无人飞行控制系统，集成了先进的通信和指控技术，扩展了在极具挑战性的太空环境中的操作能力。 9）高度集成的自主故障管理多功能模块化GNC，可用于深空和近地探测。 10）小型化、低功耗GNC敏感器，用于提高载人和无人飞行器的安全性和可靠性，并可用于微小卫星/纳米卫星。微技术和纳米技术应用于超低功耗电子产品，有助于减小飞行器的重量和体积。 11）可靠的长寿命执行机构，包括高灵敏度的控制力矩陀螺、先进的电推等。 1.1.4天文观测和激光通信领域GNC 随着对望远镜指向控制和抖动抑制要求越来越高，高速激光通信需要光束捕获、指向和稳定跟踪控制，指向传感与控制相关的GNC系统技术如下： 1）高分辨率、低抖动导航敏感器，包括高精度相对导航敏感器和惯组。 2）高分辨率和高动态执行机构。 3）影像稳定和防抖控制，包括惯性伪星参考稳像技术，图像匹配识别技术等。 4）波前传感与控制，包括先进的光学波前传感器、快速转向反射镜、高带宽驱动、光路长度延迟线控制等。 5）振动/抖动主动传感与控制。 6）多变量自适应控制和自主重构。开发自适应飞行控制系统，根据操作环境的变化，对飞行器进行指令重新配置。 7）多学科优化设计。 8）端到端集成的动态建模技术。研究光学、结构、散热、电源等相关系统之间的相互作用。 9）由传感器、算法和执行机构组成的超精密抗干扰系统。 1.1.5GNC系统发展趋势 随着微电子和微机械应用，星敏感器、跟瞄相机、惯组、太阳敏感器、反作用飞轮、动量轮、减振平台、快反镜、力矩陀螺、焦平面阵列、制冷机、推力器等空间飞行器GNC系统敏感器和执行机构不断小型化。随着集成芯片技术的发展，控制器计算速度不断提高，未来GNC系统可能的发展方向如下： 1）编队飞行和自主交会对接敏感器及执行机构小型化； 2）无人空间飞行器自适应飞行控制； 3）抗干扰高低精度敏感器组合定姿和组合导航； 4）低功耗、轻质、小型化、高可靠敏感器； 5）力矩陀螺、飞轮等高可靠、长寿命执行机构； 6）微振动抑制高精度指向、捕获和跟踪； 7）空间组网运行星座自主相位保持及快速机动； 8）深空探测自主导航和自主轨控。 1.2主要内容 1.2.1数学基础 GNC系统数学建模需要一定的数学基础，主要包括坐标系转换、四元数乘法、不同转序姿态解算、球面三角形公式、卡尔曼滤波和龙格库塔积分等，这些数学公式是后续建模和仿真的基础。 1.2.2仿真基础 GNC系统数学建模也需要一定的工程基础，主要包括姿态和轨道简易动力学、章动角和章动角速率、格林尼治恒星时计算、微波和光学基本公式、单/双目视觉、大气密度、轨道参数转换、日月轨道及日月指向等算法。考虑到超实时仿真对动力学快速迭代推算的需求，结合前期发表的论文，给出挠性动力学模型解耦及局部迭代方法、挠性和晃动动力学模型解耦及局部迭代等算法。姿态和轨道动力学及运动学是姿态和轨道确定与控制的基础，也是数学仿真的基础。 1.2.3组合定姿 组合定姿包括常规的陀螺漂移与姿态四元数偏差关系、地磁场模型、基于卡尔曼滤波磁强计与陀螺组合定姿、双矢量和多矢量定姿等姿态确定算法。考虑到卡尔曼滤波算法相对复杂、物理概念不够清晰，结合前期发表的论文，给出基于PI滤波陀螺与磁强计组合定姿，陀螺与星敏感器组合定姿及陀螺漂移估计，地平仪与陀螺轨道罗盘组合定姿，以及星敏外场观星姿态及极性测试算法，并经仿真和在轨验证。组合定姿是GNC系统的关键，没有姿态基准，姿态控制、导航和轨控都将无法进行。 1.2.4组合导航 组合导航包括常规的轨道参数转换、轨道平根与瞬根相互转换、WGS84地固系与J2000惯性系相互转换、WGS84位置/速度与经纬高及地平系速度相互转换、轨道根数与J2000惯性系位置/速度相互转换、基于最小二乘地面静态对准、惯组建模及标定、纯惯导、GNSS兼容机与惯组卡尔曼滤波组合导航、基于跟瞄近程相对导航卡尔曼滤波算法等组合导航算法。考虑到卡尔曼滤波算法相对复杂、物理概念不够清晰，结合前期发表的论文，给出基于PI滤波估计加表漂移的GNSS兼容机与加表组合导航、伪距及伪距率与加表组合导航及漂移估计、基于跟瞄和加表PI滤波近程相对导航、基于光电测角相对导航等算法，以及基于星敏与地平仪的惯性天文组合导航，星光折射敏感器谱宽选择及探测能力分析，基于星光折射敏感器和惯组组合导航，基于准北东地系地面静态对准等算法，并经仿真和在轨验证。组合导航是轨道控制的基础，地面测定轨为主逐步向星上自主定轨转变。 1.2.5姿态控制 姿态控制也是GNC系统设计和研制的重要内容之一，主要包括挠性动力学及模态增益、常规的干扰力矩分析、飞轮加磁稳定控制、磁消旋及磁速率阻尼控制、力矩陀螺控制、六支腿Stewart减振平台动力学、基于姿态偏差四元数姿态机动等姿控算法。考虑到干扰力矩补偿和部件故障时姿态控制的可靠性和安全性，结合前期发表的论文，给出低轨偏置动量卫星气动干扰力矩补偿控制、近地零动量卫星干扰力矩飞轮补偿控制、非偏置动量单飞轮加磁控制、异轨交会过程高精度指向跟踪控制、微振动抑制和高精度指向控制等算法，经仿真和在轨验证的稳定平台是轨控的基础。 1.2.6轨道控制 轨道控制研究包括常规的地球非球形摄动分析、太阳同步等典型轨道分析、改进春分点根数转换、轨道面内轨控、轨道面外轨控、兰伯特轨道转移等轨控相关算法，考虑到载荷对轨道机动要求不断提高，结合前期发表的论文，给出远程快速轨道机动、再入返回离轨时机及制动策略、基于再入误差预测精确返回控制、近程自主伴飞及防撞控制等算法，并经仿真和在轨验证。随着飞行器组网和编队需求的增强，轨道机动控制越来越重要。 1.2.7深空探测GNC系统 深空探测越来越受到重视，结合火星探测研制经验，给出深空探测不同坐标系及轨道转换、行星轨道推算、行星探测器对地指向、火星时、火影、对火中继、引力影响球、近火捕获和再入制动策略、火星表面静态对准，以及小行星探测轨道转移策略等GNC系统相关算法及仿真验证结果。 1.3GNC系统建模及仿真 GNC系统工程研制涉及大量的模型及仿真，模型包括敏感器模型、姿态动力学模型、轨道动力学模型、执行机构模型等，仿真包括超实时数学仿真、实时仿真、半实物仿真和全物理仿真等。 航天早期研制过程由于缺乏经验和计算机技术限制，较多采用半实物仿真和全物理仿真，随着计算机技术的发展，数学仿真和实时仿真的重要性越来越突出，基于数字单机的虚拟仿真在研制过程中的作用正在逐步增强。 GNC系统数学建模及仿真一方面有助于熟悉和了解GNC系统研制理论和特点，另一方面为后续进入工程研制打下建模和仿真基础。 从工程应用角度，在满足研制需求的前提下GNC系统模型越简单越好，精简的模型不但有利于提高仿真效率，还便于分析和查找问题，模型的精简需要一定的工程经验，特别是参数选取。 仿真模型的验证至关重要，因为很多航天事故是由于模型不正确导致的，本书是在前期型号研制的基础上编写的，在符合工程研制需求的前提下尽量采用精简的模型，这些模型在地面经过了大量的复核，大部分模型已经过在轨考核验证。 本书提供了适合于GNC系统设计与研制的工程应用算法，并经数学仿真和在轨验证。早期GNC系统设计与研制主要依靠物理仿真验证，随着计算机技术的发展和飞行器动力学及运动学建模准确性的提高，GNC系统设计与研制越来越依靠数学仿真验证，特别是超实时仿真，不但能加快研制速度，而且也能够大幅降低研制成本。 1.4GNC系统研制关注重点 GNC系统研制关注重点如下： （1）理解公式物理概念 姿态方面包括坐标转换、欧拉轴角与四元数关系、四元数非唯一性、偏差四元数与陀螺角速率关系、双矢量定姿、刚体动力学、动不平衡力矩、静不平衡力矩、磁场与磁矩、挠性动力学、模态增益、基频与刚度关系等物理概念。 轨道方面包括质点轨道动力学、力学时与世界时、跳秒、单/双目视觉、跟瞄视线距和视线角、相对位置和相对姿态确定、星光折射、星等、天线、轨道六要素、太阳同步轨道、回归轨道、冻结轨道等物理概念。 （2）建立基本模型库 姿态模块包括四元数乘法、龙格库塔积分、四元数求312转序欧拉角、四元数求321转序欧拉角、基于四元数求转换矩阵、基于312转序欧拉角求转换矩阵、基于321转序欧拉角求转换矩阵、惯组、飞轮、地平仪、磁场表等单机建模等。 轨道模块包括地球非球形J2至J4项摄动、太阳章动和岁差参数、儒略日和儒略世纪数计算、格林尼治真恒星时计算、日月指向简单计算、平流层和轨道空间大气密度计算、轨道根数不同表示形式相互转换、近点角MEf相互转换、轨道根数与J2000惯性系位置/速度相互转换、WGS84坐标系与地平坐标系相互转换、轨道平根与瞬根相互转换、轨道交点周期、再入点纵程与横程计算等。 建模的关键是模型要正确，基于标准数据校验模型正确性，如WGS84坐标系与J2000惯性系位置/速度转换、轨道平根与瞬根相互转换等。 （3）学会基于仿真分析问题 姿态方面包括外场观星数据正确性判别、基于地平仪轨道罗盘组合定姿仿真、基于磁强计与陀螺卡尔曼滤波组合定姿仿真、多矢量定姿、干扰力矩建模及对姿态影响仿真分析、零动量加磁卸载控制、偏置动量加磁章进动控制、双自由度偏置动量加磁卸载控制、非偏置动量单飞轮加磁控制、基于偏差四元数姿态机动控制等。 轨道方面包括摄动力建模及对轨道影响仿真分析、静态粗对准和精对准、兼容机输出位置/速度与惯性组合导航、基于跟瞄输出近程相对导航、小偏心率轨道推算、大偏心率轨道推算、轨道面内控制、轨道面外控制、升交点赤经漂移控制、再入点迹向与法向误差分析等。 （4）提高解决建模和仿真问题能力 复杂的建模与仿真包括高阶地球非球形摄动模型与STK误差比对、基于地固系惯导与基于J2000惯性系惯导比对、地面重力加速度与引力和离心力关系、兼容机输出伪距与惯性组合导航等，通过这些复杂建模和仿真提高解决问题的能力。 具有飞行器姿态和轨道控制基础，且后续还将从事飞行器GNC系统设计与研制的读者，可以结合习题仿真验证模型，为后续从事该专业工作奠定基础。
［1］ ［2］王献忠,张肖.挠性和晃动动力学模型解耦及局部迭代算法［J］.航天控制，2021，38（6）：913，33. ［3］王献忠.陀螺与磁强计组合定姿及工程应用分析［J］.航天控制，2009，26（5）：4347. ［4］王献忠,张肖,刘艳.陀螺与磁强计组合定姿及陀螺漂移估计［J］.航天控制，2017，34（4）：1519. ［5］王献忠,张肖,刘艳.陀螺与磁强计组合定姿及漂移估计算法优化［J］.航天控制，2020，37（4）：913. ［6］王献忠,张肖.陀螺与星敏感器组合定姿及陀螺漂移估计［J］.航天控制，2018，35（4）：711. ［7］王献忠,张肖,张丽敏，等.基于PI滤波估计加表漂移的GNSS兼容机与加表组合导航［J］.航天控制，2018，35（6）：3641. ［8］王献忠,张肖,张丽敏，等.基于跟瞄和加表PI滤波近程相对导航方法［J］.上海航天，2019，36（1）：4347. ［9］王献忠,张肖,张丽敏.基于地平仪的惯性天文组合导航［J］.导航定位学报，2016，22（3）：2630. ［10］王献忠,张肖,张丽敏，等.伪距及伪距率与加表组合导航及漂移估计［J］.航天控制，2020，37（6）：38. ［11］王献忠,刘禹,汤敏兰，等.星光折射敏感器谱宽选择及探测能力分析［J］.系统仿真学报，2018，30（2）：722730. ［12］王献忠,张肖.基于准北东地系地面静态对准算法［J］.航天控制，2022年，39（1）：3741. ［13］王献忠,张肖,张丽敏.低轨偏置动量卫星气动干扰力矩补偿控制研究［J］.上海航天，2017，34（2）：7478. ［14］王献忠,张肖,张丽敏.近地零动量卫星干扰力矩飞轮补偿控制［J］.航天控制，2016，33（2）：6065. ［15］王献忠,张肖.一种非偏置动量单飞轮加磁控制算法［J］.空间控制技术与应用，2014，40（6）：4852. ［16］王献忠,张肖,张国柱.挠性振荡在轨自主识别算法［J］.控制工程，2017，S1：5658. ［17］王献忠,汤敏兰,张丽敏.再入返回离轨时机及制动策略［J］.航天控制，2016，33（4）：5358.
第2章数学基础 2.1坐标系旋转及欧拉轴角 2.1.1坐标系旋转 分别绕三轴旋转角度的转换矩阵如下： Rx（）=1000cossin0-sincos（21） Ry（）=cos0-sin010sin0cos（22） Rz（）=cossin0-sincos0001（23） 设坐标系绕X轴旋转角度，原坐标系中的矢量r在新坐标系中的矢量rn为 rn=Rx（）·r=1000cossin0-sincos·r（24） 新坐标系中的矢量rn在原坐标系中的矢量r为 r=Rx（）T·rn=1000cos-sin0sincos·rn（25） 坐标系绕Y轴和Z轴旋转角度，计算过程类似。 2.1.2欧拉轴角 以欧拉轴角定义的姿态四元数如下： q=q0q1q2q3=cos2exsin2eysin2ezsin2（26） 式中 绕旋转轴旋转角度，以欧拉轴角定义的姿态四元数具有清晰的姿态角含义。 绕e旋转角度和绕-e旋转角度2-得到的姿态是相同的，如图21所示。 图21姿态不同旋转方向 ·· ·· 从数学上计算如下： q=cos2-2exsin2-2eysin2-2ezsin2-2=-cos2exsin2eysin2ezsin2=-q0q1q2q3=-q（27） 说明四元数非唯一性，一般对q0＜0进行正则化处理如下： q=-q0q1q2q3（28） 正则化处理后的姿态四元数用于姿态机动转动角度q0＜0情况，有利于快速姿态机动。 2.2四元数乘法 2.2.1四元数乘法形式 姿态从惯性坐标系转到轨道坐标系，再从轨道坐标系转到星体坐标系，姿态转移矩阵采用左乘形式 Abi=Abo·Aoi（29） 姿态四元数乘法采用右乘形式，如图22所示，初始姿态四元数p经四元数q转动后为姿态四元数r，姿态四元数采用右乘形式 r=pq（210） 图22姿态旋转 姿态转换矩阵采用左乘形式比较直观，姿态四元数乘法采用右乘形式比较直观，应用时姿态转换矩阵一般采用左乘形式，姿态四元数乘法一般采用右乘形式。 2.2.2四元数乘法推导 基于右乘形式推导四元数乘法公式，四元数乘法可以看作将3维叉乘扩展到4维空间，令 p=p0Pi，P=p1p2p3（211） q=q0Qi，Q=q1q2q3（212） 求得右乘形式四元数乘法公式为 pq=p0Piq0Qi=p0·q0-P·Q（p0·Q q0·P PQ）i（213） 四元数书写时一般省去矢量标识i，得到四元数乘法如下 pq=p0Pq0Q=p0·q0-P·Qp0·Q q0·P PQ =p0-PTP［P］ p0·I33·q0Q =q0-QTQ-［Q］ q0·I33·p0P（214） 其中 ［Q］=0-q3q2q30-q1-q2q10 ［P］=0-p3p2p30-p1-p2p10 同理求得 qp=q0Qp0P=q0·p0-Q·Pq0·P p0·Q QP =q0-QTQ［Q］ q0·I33·p0P =p0-PTP-［P］ p0·I33·q0Q（215） 3维叉乘QP=-PQ，由式（214）和式（215）得 qp=pq 20013031［Q］·p0P（216） 2.2.3姿态四元数运动学推导 设姿态角速率为，姿态推算时间为dt，则姿态角变化A为 A=·dt 姿态变化四元数qe为 qe=cos（/2）sin（/2）A（217） 其中 =A 推算时间短，较小，sin（/2）/2，姿态四元数qe近似为 qe1·dt/2（218） 令初始姿态四元数q为 q=q0Q（219） 求得新姿态四元数qn为 qn=qqe=q0Q1·dt/2（220） 即 qn=q0-Q··dt/2Q q0·dt/2 Q·dt/2（221） 求得dq为 dq=qn-q=-Q··dt/2q0·dt/2 Q·dt/2（222） 求得q·为 q·=dqdt=-Q·/2q0/2 Q/2（223） 应用四元数乘法形式把式（223）写成 q·=12q0Q0（224） 其中，0不是真正四元数，但满足四元数乘法规则。 得到姿态四元数运动学方程如下： q·=12q0（225） 2.2.4姿态四元数应用处理 实际应用中要对姿态四元数进行归一化处理，即 p=p0p1p2p3=1p20 p21 p22 p23·p0p1p2p3（226） 姿态四元数写成欧拉轴角形式为 p=cos2sin2·e（227） 其中 |e|=1 绕欧拉轴e旋转角与绕矢量轴-e旋转2-结果是一样的，姿态四元数不具有唯一性，即 p=p0P与p=-p0P=-p0-P 是等效的。 从姿态机动快捷性出发，姿态机动四元数选择p00情况，应用中需要对姿态四元数进行正则化处理，当p0＜0时 p=|p0|sign（p0）·P（228） 姿态四元数的共轭四元数为 p*=p0-P（229）
2.3.1姿态转换矩阵 姿态四元数表示的姿态转换矩阵如下： A=q20 q21-q22-q232q0q3 2q1q2-2q0q2 2q1q3-2q0q3 2q1q2q20-q21 q22-q232q0q1 2q2q32q0q2 2q1q3-2q0q1 2q2q3q20-q21-q22 q23（230） 令 A=a11a12a13a21a22a23a31a32a33（231） 2.3.2姿态四元数按312转序转为姿态角 姿态按312转序表示的姿态转换矩阵如下： A=cos0-sin010sin0cos·1000cossin0-sincos·cossin0-sincos0001（232） A=cos0-sin010sin0cos·cossin0-cossincoscossinsinsin-sincoscos（233） A=coscos-sinsinsincossin sinsincos-cossin-cossincoscossinsincos sincossinsinsin-sincoscoscoscos（234） 姿态四元数按312转序求三轴姿态角，如果|a23|0.999 9，则 sin=a23，=arcsin（a23）（235） tan=-a13/cosa33/cos，=arctan2-a13/cosa33/cos（236） tan=-a21/cosa22/cos，=arctan2-a21/cosa22/cos（237） 否则 =arcsin（a23）（238） =0（239） =arctan2a12a11（240） 2.3.3姿态四元数按313转序转为姿态角 姿态按313转序表示的姿态转换矩阵如下： A=cossin0-sincos0001·1000cossin0-sincos·cossin0-sincos0001（241） A=cossin0-sincos0001·cossin0-cossincoscossinsinsin-sincoscos（242） A=coscos-sincossincossin sincoscossinsin-sincos-coscossin-sinsin coscoscoscossinsinsin-sincoscos（243） 姿态四元数按313转序求三轴姿态角，如果|a23|0.999 9，则 cos=a33，=arccos（a33）（244） tan=a31/sin-a32/sin，=arctan2a31/sin-a32/sin（245） tan=a13/sina23/sin，=arctan2a13/sina23/sin（246） 否则 =arctan2a12a11（247） =arccos（a33）（248） =0（249） 2.3.4姿态四元数按321转序转为姿态角 姿态按321转序表示的姿态转换矩阵如下： A=1000cossin0-sincos·cos0-sin010sin0cos·cossin0-sincos0001（250） A=coscoscossin-sinsinsincos-cossincoscos sinsinsinsincoscossincos sinsincossinsin-sincoscoscos（251）
你还可能感兴趣
我要评论
|
