两种典型异构星座摄动轨道偏置与保持控制*

贺波勇, 姜宇, 李恒年

航天控制 ›› 2022, Vol. 40 ›› Issue (2) : 41-46.

PDF(919 KB)
2025年4月22日 星期二 Email Alert  RSS
PDF(919 KB)
航天控制 ›› 2022, Vol. 40 ›› Issue (2) : 41-46.
制导、导航与控制技术

两种典型异构星座摄动轨道偏置与保持控制*

作者信息 +

Perturbed Orbit Offset Deployment and Configuration Maintenance of Two Typical Heterogeneous Constellations

Author information +
文章历史 +

摘要

分析了低中高不同轨道高度受摄演化规律,针对2种典型的异构混合星座,提出了星座部署偏置量计算方法及星座构型保持控制策略。算例表明,合理利用入轨轨道参数偏置和摄动运动保持控制,可以实现异构星座卫星最优协同工作。

Abstract

The perturbation evolution laws of low, medium and high orbital heights are analyzed. Regarding two typical heterogeneous constellations, the optimal offset calculation method of orbit deployment and the constellation configuration keeping control strategy are designed. Examples show that the optimal coordination of heterogeneous constellation satellites can be realized using orbit development bias and perturbation motion keeping control.

关键词

异构星座 / 摄动分析 / 偏置部署 / 保持控制

Key words

Heterogeneous constellations / Perturbation analysis / Orbit bias / Keep control

引用本文

导出引用
贺波勇 , 姜宇 , 李恒年. 两种典型异构星座摄动轨道偏置与保持控制*[J]. 航天控制, 2022, 40(2): 41-46
He Boyong , Jiang Yu , Li Hengnian. Perturbed Orbit Offset Deployment and Configuration Maintenance of Two Typical Heterogeneous Constellations[J]. Aerospace Control, 2022, 40(2): 41-46
中图分类号: V412.4 (航天器飞行力学)   

0 引言

异构混合星座是指携带功能相似或功能不同的有效载荷,为了协同完成某项任务而组合在一起的2个或2个以上子星座构成的复合星座[1]。例如,我国环境-1A和环境-1B两颗光学卫星组网能够以50m的对地分辨率在2天内覆盖全球[2],若与环境-1C电子侦察卫星组网,则星座能够以更高的时间分辨率对地协同观测[3]。张雅声等[4]利用椭圆冻结轨道和赤道中轨道特殊几何关系组合,设计了一种只需6颗卫星的高性价比异构预警星座。Zhao等[5]给出了一般性异构星座构型的重构思路。王茂才等[6]设计了一种双层异构协同对地观测星座。加拿大发射部署的“雷达卫星星座任务(Radar constellation mission)”由3颗完全相同的卫星组成,雷达卫星-1/2轨道高为798km,RCM卫星轨道高为600km,均为太阳同步晨昏轨道[7]。我国2020年6月建成的北斗三代导航星座采用异构星座增强亚太区域导航[8]。美国Space-X公司正在建设的互联网星座Star-link预计2025年有12000颗卫星组网,之后可扩充至42000颗,Star-link大致可以分为3个不同轨道高度的Walker星座,轨道高度和倾角分别为340km(倾角:42~53°)、550km(倾角:53°)和1150km(倾角:53~81°)[9]。可见,异构星座因其独特的优势逐渐被重视。
航天器所处轨道高度不同,受到的摄动力种类和大小不同[10]。传统习惯以1000km和20000km作为低、中、高轨道高度分界线,则低轨卫星除考虑地球非球形摄动外,还须考虑大气阻力引起的半长轴、偏心率和相位变化。中轨卫星和高轨道受大气摄动可以忽略,但三体摄动和光压摄动效果增强。李恒年等[11-12]量化分析了中高轨道摄动运动规律及星座构型发散的主要原因,并分析了绝对控制与相对控制的优缺点。姜宇等[13]研究了Walker星座摄动及构型保持控制策略,表明星座入轨偏置量数值修正求解法略优于摄动补偿法。钱山等[14]提出了一种解耦的中轨道星座位置保持控制策略,并分析了星座部署时机对位置保持控制量的影响。陈雨等[15]结合我国某首例低轨Walker星座轨道实测数据,利用遍历寻优基准星的相对控制策略,实现了卫星间的相对相位保持控制。陈长春等[16]提出采用最小二乘法求解Walker星座整体偏移方法,获得燃耗最小控制策略。相比较于异构星座,同构Walker星座由于其全球覆盖均匀性特点,研究较为广泛。
本文针对2种典型的异构星座,分析其受摄运动规律,并设计入轨偏置部署与保持控制策略。

1 摄动轨道

地球轨道卫星受到的主要摄动力为地球扁率J2项、大气阻力、日月中心引力以及太阳光压等[17]

1.1 地球扁率摄动

由Langrage方程知,将地球扁率J2项势函数代入,可得卫星平均轨道根数的长期变化率为
a˙=e˙=i˙=0Ω·=-32J2Rep2n-cosiω˙=34J2Rep2n-4-5sin2iM˙=n-=n1+34J2Rep22-3sin2i1-e2
(1)
式中: n-为考虑地球扁率摄动的平运动角速率,n= μa-3为不考虑摄动的平均运动角速率,a为半长轴,e为偏心率,i为倾角,Ω为升交点赤经,ω为近地点幅角,M为平近点角,J2=1.08263×10-3,地球赤道半径Re=6378.137 km,地球引力常数μ=3.98600436×105m3·s-2,p=a(1-e2)为半通径。
可见,考虑地球扁率摄动作用的轨道升交点赤经、近地点幅角及平近点角变化率均有长期影响,变化量与卫星的轨道半长轴、偏心率和倾角有关。

1.2 大气阻力摄动分析

面质比为 sm的低轨航天器受到大气的阻力加速度为
aD=- s2mCDρvv
(2)
式中:CD为阻力系数,ρ为航天器当前位置的大气密度,有很多大气模型可供选择。v为航天器相对大气速度矢量,v为其大小。大气阻力主要引起航天器轨道能量衰减,造成半长轴和偏心率摄动变化,其长期变化率为
a˙=-CDsmρna21+e2+2ecosf1-e232e˙=-CDsmρna(e+cosf)1+e2+2ecosf1-e2
(3)
式中:f为航天器轨道真近点角。大气阻力摄动作用取决于航天器面质比、轨道高度和大气密度,大气密度模型受时间、季节、太阳活动和地磁活动等影响。

1.3 日月三体摄动分析

太阳中心引力摄动引起的卫星轨道倾角长期变化率为
i·= 3ns22n(cos β scos Ω+sin β scos issin Ω)·[cos β ssin isin Ω-sin β scos issin icos Ω+sin β ssin iscos i]
(4)
式中:βs为太阳视运动的黄经,is为黄道倾角,ns为地球绕太阳公转的角速率。如果不计βs的周期项,则上式变为
i·= 3ns28n[sin 2Ω sin i+sin 2issin Ω cos i-sin 2Ω cos 2issin i]
(5)
对于太阳同步轨道特例
i˙=- 3ns216nsin i(1+cos is)2sin(2βs-2Ω)
(6)
月球中心引力摄动引起的航天器轨道倾角长期变化率为
i·= 38mmme+mmnm2nsin(Ω-Ωm)·[2cos(Ω-Ωm)sini+sin2imcosi-2cos2imsin icos(Ω-Ωm)]
(7)
式中:mm为月球质量,me为地球质量,im为月球轨道倾角,即当时白赤夹角,Ωm为月球轨道的升交点赤经,nm为月球公转的角速率。
可见,日月中心引力引起航天器轨道倾角摄动量与航天器半长轴、倾角和升交点赤经有关,轨道越高,轨道倾角受摄长期变化率幅值越大。

1.4 太阳光压摄动分析

利用偏心率矢量 ξ=ecosω,η=-esinω,则太阳光压引起的星座卫星轨道偏心率矢量变化率为
ξ˙=-3CPspssin(βs+α)2nam·(cosisinΩ)2+sinisinis+cosicosiscosΩ2η˙=-3CPspscos(βs+α)2nam·cos2Ω+(cosissinΩ)2
(8)
式中:
cosα= sinisinis+cosicosiscosΩ[cosisinΩ,sinisinis+cosicosiscosΩ]
sinα= -cosisinΩ[cosisinΩ,sinisinis+cosicosiscosΩ]
(9)
这里, sm为航天器受晒的面质比,CP=0~1为光压系数,对应全反射和全吸收,一般取均值0.5,地球附近光压强度ps=4.56×10-6N/m2
可见,ξ与轨道倾角相关,而η与轨道倾角无关。太阳光压对中轨道航天器轨道偏心率矢量具有周期为一年的长周期项影响,幅值为 3CPsps2nam。特别地,地球静止轨道卫星倾角为0,设 α˜=α,有
cos α˜= cosΩ[cos2Ω+cos2issin2Ω]
sin α˜ cosissinΩ[cos2Ω+cos2issin2Ω]
(10)

2 典型异构星座偏置部署

以异构预警星座和双层协同星座2个典型异构星座为例,研究异构星座偏置部署方法。

2.1 异构预警星座偏置部署

文献[1]和[4]设计了一种典型的异构预警星座,如图1所示,采用4颗地球倾斜冻结轨道和2颗赤道轨道组合的异构星座可以实现北半球连续一重覆盖,高纬度地区二重覆盖,且具备星间链路条件,卫星轨道六根数如表1所示。
图1 异构预警星座构型

Full size|PPT slide

表1 文献[1]中6颗卫星轨道参数
a/km e i/(°) Ω ω/(°) f/(°)
16763 0.53 63.435 69 270 90
16763 0.53 63.435 69 270 191
16763 0.53 63.435 249 270 90
16763 0.53 63.435 249 270 191
16763 0 0 30 0 0
16763 0 0 210 0 0
按照式(1)~(8)分析I号卫星与V号卫星摄动运动120个整周期(30天)后,轨道参数漂移量如表2所示。
表2 120个整周期后轨道参数增量
Δi/(°) ΔΩ/(°) Δω/(°) Δλ/(°)
0.0083 -8.8565 0 -
0.0064 - - 10.26
表2中:ΔλΩωf为衡量零倾角圆轨道卫星相位变化的赤经增量。该异构星座中,Ⅴ号卫星和Ⅵ卫星的作用主要是弥补Ⅰ~Ⅳ号卫星分别处于远地点和近地点时,2个轨道平面所夹的中间部分出现了的2个覆盖空隙,如文献[1]中图5.9 所示。只需控制Ⅴ号星和Ⅵ号星赤经漂移率与Ⅰ号星升交点赤经漂移率一致即可,又由于Ⅰ号卫星为冻结轨道,无法采用倾角偏置实现升交点赤经漂移率控制。所以,应对Ⅴ号和Ⅵ号卫星采用平半长轴偏置策略
Δ a-= 2a523tμ(Δλ-ΔΩ)
(11)
使用式(11)迭代计算3~4次可得,在Ⅴ号和Ⅵ号卫星入轨部署时,使之平半长轴比设计值增大Δ a-=+4.8993km,则可实现在Ⅰ~Ⅳ号卫星分别处于远地点和近地点时,Ⅴ号和Ⅵ号卫星对赤道附近2个轨道平面所夹的中间部分覆盖的功能。该异构星座中6颗卫星均不为太阳同步轨道,日月中心引力摄动引起的轨道倾角呈现回归年和朔望月两种复合周期性,且摄动运动幅值较小。所以,只需考虑更高阶地球带谐项作用,在卫星入轨部署时,精确部署轨道参数[18]

2.2 双层协同星座偏置部署

以文献[6]双层协同对地观测星座为例,如图2所示,研究偏置部署策略。该星座中10颗卫星轨道六根数如表3所示。
图2 双层协同星座构型

Full size|PPT slide

表3 文献[6]中10颗卫星轨道参数
a/km e i/(°) Ω/(°) ω/(°) f/(°)
A1 8714.027 0 107.1545 47.7809 0 0
A2 8714.027 0 107.1545 47.7809 0 120
A3 8714.027 0 107.1545 47.7809 0 240
B1 8714.027 0 107.1545 87.1559 0 60
B2 8714.027 0 107.1545 87.1559 0 180
B3 8714.027 0 107.1545 87.1559 0 330
C1 7012.537 0 97.9284 67.468 0 0
C2 7012.537 0 97.9284 67.468 0 90
C3 7012.537 0 97.9284 67.468 0 180
C4 7012.537 0 97.9284 67.468 0 270
假设卫星面质比 sm=0.005,按照式(1)~(8)分析A1号、B1号和C1号卫星摄动运动360天轨道参数漂移量如表4 所示。
表4 360天后轨道参数增量
Δi/(°) ΔΩ/(°) Δu/(°)
A1 0.0072 355.18 322.25
B1 0.0215 355.18 322.25
C1 0.0020 355.29 241.71
表4中Δuωf为纬度幅角增量。可见,日月三体摄动引起轨道升交点赤经摄动差异较小,360天后A1号卫星和C1号卫星升交点赤经差异0.11°,验证了文献[6]设计的星座具备较为稳定构型特性。尽管相位差距较大,但任意时刻C1~C4号卫星均在A1~A3和B1~B3号卫星构成的中轨星间链路框架中。该星座构型稳定运行的条件是精确部署入轨,也可以采用倾角偏置方式进一步提升星座升交点赤经漂移一致性
Δi=- 72Δaatani
(12)
使用式(12)迭代计算3~4次可得,在C1~C4号卫星入轨部署时,使之平倾角比设计值减小Δi=-0.0024°,则可实现升交点赤经与A1~A4号和B1~B2号卫星升交点赤经漂移率一致。

3 典型异构星座保持控制

仍以上述2个典型异构星座为例,研究各自构型保持控制策略

3.1 异构预警星座保持控制

如对文献[1]中异构预警星座倾角有更高精度要求,需进一步分别对Ⅰ~Ⅳ号卫星和Ⅴ~Ⅵ号卫星进行构型保持控制。按照式(1)~(8)分析Ⅰ号卫星与Ⅴ号卫星摄动运动1440个整周期(约1年)后,Ⅰ号卫星倾角增量为0.1°,Ⅴ号卫星倾角增量为0.08°,则会影响卫星星座性能。由于Ⅴ号和Ⅵ号卫星采用半长轴偏置方式使相位漂移率配合Ⅰ~Ⅳ号卫星,则只需依据Gauss摄动方程,在升/降交点处采用法向速度增量解耦调整轨道倾角
Δvz= μprcosuΔi
(13)
将控制周期设为1个月、3个月、6个月、和1年的倾角保持控制速度增量如表5所示。
表5 异构预警星座倾角保持控制频次与速度增量
1个月 3个月 6个月 1年
I Δi/(°) 0.0083 0.0249 0.0499 0.1001
Δvz/(m/s) 0.833 2.499 5.008 10.046
V Δi/(°) 0.0064 0.0195 0.0393 0.0800
Δvz/(m/s) 0.545 1.659 3.345 6.808
可见,异构预警星座倾角摄动主要是由日月三体摄动引起的长周期项,每1个月、3个月、6个月和1年的控制速度增量具有较好的线性叠加性,星座保持控制频次可以根据实际对倾角精度要求决定。

3.2 双层协同星座保持控制

如对文献[6]中双层协同星座倾角有更高精度要求,则需进一步控制星座中卫星轨道面协同变化。例如:以1 Jul 2020 00:00:00 (UTCG) 时刻为初始时刻,A1、B1和C1号(倾角偏置0.0024°后)卫星受摄运动1个月、3个月、6个月和1年轨道倾角漂移量如表6 所示。
表6 双层协同星座轨道面保持控制频次与速度增量
1个月 3个月 6个月 1年
A1 Δi/(°) 0.0016 0.0046 0.0091 0.0173
Δvz/(m/s) 0.1889 0.5430 1.0742 2.0421
B1 Δi/(°) 0.0031 0.0089 0.0175 0.0331
Δvz/(m/s) 0.3660 1.0506 2.0657 3.9072
C1 Δi/(°) 0.0027 0.0078 0.0155 0.0300
Δvz/(m/s) 0.3553 1.0264 2.0396 3.9476
同样依据Gauss摄动方程,在升/降交点处采用法向速度增量解耦调整轨道倾角,利用式(13)计算的倾角控制量如表6所示。
可见,由于双层协同星座轨道半长轴小于异构预警星座中椭圆冻结轨道半长轴,由日月三体摄动引起倾角摄动运动较小,每1个月、3个月、6个月和1年的控制速度增量具有较好的线性叠加性,单颗卫星年倾角保持控制速度增量小于4m/s,倾角保持控制频次可以根据实际对倾角精度要求决定。

4 结论与展望

通过对异构星座摄动因素和摄动运动分析,研究了对应的入轨参数偏置策略和保持控制方法。在实际应用中,应以异构星座中不同种类卫星相互协同配合的目的为切入点,研究对应的入轨轨道参数偏置策略和摄动运动保持控制方法,例如:异构预警星座利用零倾角中轨道半长轴偏置产生的相位滞后匹配椭圆冻结轨道升交点赤经西漂,只需定期对轨道倾角保持控制,以减小燃料消耗;双层协同星座中所有卫星自身具备太阳同步轨道特性,只需低层C1~C4号共4颗卫星入轨时倾角微量偏置即可达到升交点赤经漂移率与A1~A4号卫星和B1~B4号卫星升交点赤经漂移率一致,但需定期对所有卫星进行轨道倾角保持控制。
由于异构星座构型形式多样,不同异构星座中不同卫星轨道协同配合的目的和方式不同,应具体问题具体对待,但入轨轨道参数偏置方式和保持控制策略都应以不同种类卫星轨道协同配合的目的为切入点,借助自然摄动力,尽可能以较小的轨道调整代价控制少数卫星来匹配多数卫星轨道自然摄动漂移,协同配合共同完成工作。
由于精密轨道确定误差和入轨参数部署偏差导致的异构星座中,卫星轨道摄动运动相对漂移量一般大于同构Walker星座相应入轨参数部署偏差导致的摄动运动漂移量,故应尽快消除入轨偏差,避免可能引起的星座构型偏差非线性增长,甚至构型破坏,这是后续研究的一个重点内容。

参考文献

[1]
张雅声, 冯飞. 卫星星座轨道设计方法[M]. 北京: 国防工业出版社, 2019.
(Zhang Yasheng, Feng Fei. Orbital design method of satellite constellations[M]. Beijing: National Defense Industry Press, 2019.)
[2]
谭田, 杨芳. 环境减灾-1A、1B卫星星座轨道设计[J]. 航天器工程, 2009, 18(6): 27-30.
(Tan Tian, Yang Fang. HJ-1A/1B constellation orbit design[J]. Spacecraft Engineering, 2009, 18(6): 27-30.)
[3]
张润宁, 姜秀鹏. 环境一号C卫星系统总体设计及其在轨验证[J]. 雷达学报, 2014, 3(3): 249-255.
(Zhang Running, Jiang Xiupeng. System design and in-orbit verification of the HJ-1C SAR satellite[J]. Journal of Radars, 2014, 3(3): 249-255.)
[4]
张雅声, 姚勇. 异构预警卫星星座设计与分析[J]. 装备指挥技术学院学报, 2009, 20(3): 47-51.
(Zhang Yasheng, Yao Yong. Design and analysis of non-isomorphic early satellites constellation[J]. Journal of the Academy of Equipment Command & Technology, 2009, 20(3): 47-51.)
[5]
Zhao S, Xu Y L, Dai H Y. Research on the configuration design method of heterogeneous constellation reconstruction under the multiple objective and multiple constraint[C]// AIP Conference Proceedings, 1839,(1), 2017.
[6]
王茂才, 罗鑫, 宋志明, 等. 双层协同对地观测卫星星座设计[J]. 华中科技大学学报(自然科学版), 2018, 46(2): 100-105.
(Wang Maocai, Luo Xin, Song Zhiming, et al. Double-layered satellite constellation design for earth cooperative observation[J]. Journal of Huazhong University of Science & Technology(Natural Science Edition), 2018, 46(2): 100-105.)
[7]
李意, 徐冰. 加拿大“雷达卫星星座任务”及应用领域[J]. 国际太空, 2019, 7: 30-34.
[8]
蒋虎, 邓雷, 余金培. 北斗导航星座现状仿真分析与定量评估[J]. 天文研究与技术, 2020, 17(2): 171-177.
(Jiang Hu, Deng Lei, Yu Jinpei. Simulation of current Bei-dou navigation constellation and its quantitative assessment[J]. Astronomical Research and Technology, 2020, 17(2): 171-177.)
[9]
孙俞, 沈红新. 基于TLE的低轨巨星座控制研究[J]. 力学与实践, 2020, 42(2): 156-162.
(Sun Yu, Shen Hongxin. The control of mega-constellation at low Earth orbit based on TLE[J]. Mechanics in Engineering, 2020, 42(2): 156-162.)
[10]
何丽娜. 不同摄动力对低中高轨航天器轨道的影响分析[J]. 大地测量与地球动力学, 2017, 37(11): 1156-1160.
(He Lina. Perturbation forces analysis for spacecraft of different orbit altitudes[J]. Journal of Geodesy and Geodynamics, 2017, 37(11): 1156-1160.)
[11]
李恒年, 李济生, 焦文海. 全球星摄动运动及摄动补偿运控策略[J]. 宇航学报, 2010, 31(7): 1756-1761.
(Li Hengnian, Li Jisheng, Jiao Wenhai. Analyzing perturbation motion and studying configuration maintenance strategy for compass-M navigation constellation[J]. Journal of Astronautics, 2010, 31(7): 1756-1761.)
[12]
李恒年. 地球静止卫星轨道与共位控制技术[M]. 北京: 国防工业出版社,2010年10月第1版.
(Li Hengnian. Geostationary satellite orbital analysis and collocation strategies[M]. Beijing: National Defense Industry Press, 1st Print in Oct 2010.)
[13]
姜宇, 李恒年, 宝音贺西. Walker星座摄动分析与保持控制策略[J]. 空间控制技术与应用, 2013, 39(2): 36-41.
(Jiang Yu, Li Hengnian, Baoyin Hexi. On perturbation and orbital maintenance control strategy for Walker constellation[J]. Aerospace Control and Application, 2013, 39(2): 36-41.)
[14]
钱山, 李恒年, 伍升钢. MEO非共振轨道导航星座摄动补偿控制[J]. 国防科技大学学报, 2014, 36(2): 53-60.
(Qian Shan, Li Hengnian, Wu Shenggang. Perturbation compensation strategy for MEO non-resonant navigation constellation maintenance[J]. Journal of National University of Defense Technology, 2014, 36(2): 53-60.)
[15]
陈雨, 赵灵峰, 刘会节, 等. 低轨Walker星座构型演化及维持控制[J]. 宇航学报, 2019, 40(11): 1296-1303.
(Chen Yu, Zhao Lingfeng, Liu Huijie, et al. Analysis of configuration and maintenance strategy of LEO Walker constellation[J]. Journal of Astronautics, 2019, 40(11): 1296-1303.)
[16]
陈长春, 林滢, 沈鸣, 等. 一种考虑摄动影响的星座构型稳定性设计方法[J]. 上海航天, 2020, 37(1): 33-37.
(Chen Changchun, Lin Ying, Shen Ming, et al. A novel design method for the constellation configuration stability considering the perturbation influence[J]. Aerospace Shanghai, 2020, 37(1): 33-37.)
[17]
McGrath C N, Macdonald M. General perturbation method for satellite constellation reconfiguration using low-thrust maneuvers[J]. Journal of Guidance, Control, and Dynamics, 2019, 42(8): 1676-1692.
A general perturbation solution to a restricted low-thrust Lambert rendezvous problem, considering circular-to-circular in-plane maneuvers using tangential thrust and including a coast arc, is developed. This provides a fully analytical solution to the satellite reconnaissance problem. The solution requires no iteration. Its speed and simplicity allow problems involving numerous spacecraft and maneuvers to be studied; this is demonstrated through two case studies. In the first, a range of maneuvers providing a rapid flyover of Los Angeles is generated, giving an insight to the trade space and allowing the maneuver that best fulfills the mission to be selected. A reduction in flyover time from 13.8 to 1.6 days is possible using a less than 17 m/s velocity change. A comparison with a numerical propagator including atmospheric friction and an 18th-order tesseral model shows 4 s of difference in the time of flyover. A second study considers a constellation of 24 satellites that can maneuver into repeating ground track orbits to provide persistent coverage of a region. A set of maneuvers for all satellites is generated for four sequential targets, allowing the most suitable maneuver strategy to be selected. Improvements in coverage of greater than 10 times are possible as compared to a static constellation using 35% of the propellant available across the constellation.
[18]
Cutting E, Born G H, Frautnick J G. Orbit analysis for Seat-A[J]. Journal of the Astronautical Sciences, 1978, 26(4): 315-342.

基金

*国家自然科学基金(11902362)
中国博士后科学基金第68批面上资助(2020M683764)
PDF(919 KB)

743

Accesses

0

Citation

Detail

段落导航
相关文章

/