首页> 中国专利> 基于可见光地球敏感器和太阳敏感器的卫星姿态确定方法

基于可见光地球敏感器和太阳敏感器的卫星姿态确定方法

摘要

基于可见光地球敏感器和太阳敏感器的卫星姿态确定方法,利用可见光地球敏感器和太阳敏感器联合确定卫星对地的三轴姿态。首先通过对CMOS成像器件所成的地球图像处理,提取到地球图像的有效边缘,然后应用最小二乘法拟合提取到的有效边缘的中心。其次在卫星对地姿态角为小角度的假设下,近似计算卫星的对地滚动角和俯仰角。最后应用太阳敏感器和近似计算得到的滚动角和俯仰角姿态信息计算得到卫星对地的偏航角。相对于传统的利用地球红外敏感器和太阳敏感器组合确定卫星对地的三轴姿态角,本发明方法提高了卫星姿态角测量精度,并且该方法适用于小型化的小卫星姿态测量,顺应当前小卫星的发展潮流。

著录项

  • 公开/公告号CN105928527A

    专利类型发明专利

  • 公开/公告日2016-09-07

    原文格式PDF

  • 申请/专利权人 航天东方红卫星有限公司;

    申请/专利号CN201610262687.5

  • 发明设计人 万小波;张晓敏;

    申请日2016-04-25

  • 分类号G01C21/24(20060101);

  • 代理机构11009 中国航天科技专利中心;

  • 代理人陈鹏

  • 地址 100094 北京市海淀区北京市5616信箱

  • 入库时间 2023-06-19 00:28:54

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2019-04-30

    授权

    授权

  • 2016-10-05

    实质审查的生效 IPC(主分类):G01C21/24 申请日:20160425

    实质审查的生效

  • 2016-09-07

    公开

    公开

说明书

技术领域

本发明属于卫星控制领域,涉及一种确定卫星姿态的方法。

背景技术

目前,在以地球作为参考矢量确定卫星姿态角的敏感器中,红外地球敏感器得到了广泛的应用,并且日趋成熟,因此利用红外地球敏感器和太阳敏感器联合确定卫星对地的三轴姿态成为了常见方法。

现有的利用红外地球敏感器和太阳敏感器联合确定卫星对地三轴姿态的方法是:一方面,利用红外地球敏感器扫描机构进出地球二氧化碳红外辐射带以及红外扫描机构到达基准点的时间,计算出红外地球敏感器扫描地球的弦宽和地心矢量E在红外地球敏感器的测量坐标系xSEoySE平面中的投影与xSE轴的夹角,由上述计算所得的夹角以及球面三角的知识可以计算出地心矢量E在红外地球敏感器测量坐标系中的坐标表示,又由红外地球敏感器的安装矩阵可得地心矢量E在星体坐标系中的坐标表示;另一方面,地心矢量E在卫星轨道坐标系的坐标已知,由未知卫星姿态角计算的旋转矩阵可计算得地心矢量E在星体坐标系的坐标表示,由上述两种途径转换成的地心矢量在星体坐标系下的不同表达式,即可求出卫星的滚动角和俯仰角,然后,太阳敏感器利用已计算得到的卫星滚动角和俯仰角求得卫星的偏航角。

由上述红外地球敏感器和太阳敏感器联合确定卫星对地三轴姿态的原理可知,卫星姿态角的测量精度受限于时间测量精度和扫描机构扫描角速度的测量精度,以及卫星轨道参数测量精度(视不同的安装矩阵)。由于上述测量物理量均具有不确定性,因此造成卫星姿态的测量精度波动较大。

另外,当前卫星正朝着小型化、智能化、高精度、高稳定性、低功耗和长寿命方向发展,而主要的红外地球敏感器有圆锥扫描和摆动扫描等,这类红外 地球敏感器由于存在扫描机构,其质量大、功耗多,也不适合当前卫星的发展潮流。

发明内容

本发明解决的技术问题是:克服现有技术的不足,提供了一种基于可见光地球敏感器和太阳敏感器的卫星姿态确定方法,解决了现有红外地球敏感器和太阳敏感器联合确定卫星姿态中时间测量精度、扫描机构角速度测量精度、卫星轨道参数测量精度导致卫星姿态角测量精度下降的问题,并且该方法适用小型化的小卫星姿态测量,顺应当前小卫星的发展潮流。

本发明的技术解决方案是:基于可见光地球敏感器和太阳敏感器的卫星姿态确定方法,包括如下步骤:

(1)利用可见光地球敏感器中的CMOS成像器件对地球进行成像,获取输入图像;

(2)对所述的输入图像进行预处理,去除噪声,平滑图像;

(3)对步骤(2)获取的图像进行边缘检测,获取有效边缘;

(4)对所述有效边缘进行拟合,获取有效边缘中心和半径;

(5)计算获得卫星对地的滚动角φ和俯仰角θ,其中

φ=arctan(-Cy*cos2()f)θ=arctan(Cx*cos2()f)

(Cx,Cy)为有效边缘中心坐标,R为有效边缘的半径,f为CMOS成像器件的焦距;

(6)计算获得卫星对地的偏航角其中DssB为太阳敏感器相对 卫星本体坐标系的安装矩阵,式中SAX、SAY、SAZ分别表示太阳敏感器测量坐标系xss轴与卫星本体坐标系xb、yb、zb三轴的夹角余弦值,SBX、SBY、SBZ分别表示太阳敏感器测量坐标系yss轴与卫星本体坐标系xb、yb、zb三轴的夹角余弦值,SCX、SCY、SCZ分别表示太阳敏感器测量坐标系zss轴与卫星本体坐标系xb、yb、zb三轴的夹角余弦值;So=[SoX>oY>oZ]T为太阳矢量在轨道坐标系下的三轴分量形式,α是太阳方向矢量在太阳敏感器测量坐标系yssozss平面投影与ozss轴的夹角。

所述步骤(1)中的输入图像应满足:CMOS成像器件所成地球图像地平轮廓有效弧长占到总弧长的50%。所述步骤(2)中平滑图像的方法为3×3卷积模板的中值滤波。所述步骤(3)中边缘检测的方法为:首先采用Sobel算子检测所有的边缘,然后采用Hough变换检测出有效边缘。所述步骤(4)中对有效边缘进行拟合的方法为最小二乘法。

本发明与现有技术相比的优点在于:

(1)本发明方法通过可见光地球敏感器和太阳敏感器联合确定卫星对地的三轴姿态角,相比于传统的红外地球敏感器和太阳敏感器确定卫星的对地姿态角,该方法消除了时间测量精度、扫描机构角速度测量精度、轨道参数测量精度对红外地球敏感器测量卫星姿态角的影响,从而提高卫星姿态测量的精度,而且该方法实现容易,算法简单,对卫星平台的控制要求不高,是一套完整的确定卫星对地三轴姿态方法;

(2)本发明方法的重点在于有效边缘的提取,边缘提取的好坏直接决定方法的精度。本发明采取的提取边缘方法不仅能精确提取有效边缘,满足计算量的要求,而且算法的鲁棒性很好;

(3)本发明方法可集成到控制器中,与其他部组件共同形成轻质量、精度满足要求的可见光敏感器,为确定卫星姿态提供另外一种解决方案。

附图说明

图1为本发明方法的流程框图;

图2为本发明CMOS相机成像示意图;

图3为本发明卫星姿态坐标系示意图;

图4为本发明太阳敏感器坐标系示意图。

具体实施方式

随着微小型化CMOS成像器件的快速发展,不仅为地球可见光敏感器的真正使用提供了极大可能,而且顺应了卫星的发展趋势。本发明就基于可见光地球敏感器和太阳敏感器如何联合确定卫星对地三轴姿态,提出一套完整的方法。

如图1流程所示,本发明方法主要包括以下几个步骤:第一,CMOS成像器件对地球成像后输入图像;第二,图像预处理;第三,提取图像有效边缘;第四,应用最小二乘拟合图像有效边缘中心;第五,在卫星对地姿态角为小角度假设前提下,根据有效边缘中心近似计算卫星滚动和俯仰姿态角;第六,应用太阳敏感器和近似计算得到的滚动角和俯仰角姿态信息,计算卫星对地的偏航角。以下分别详细介绍上述六个步骤。

第一,可见光地球敏感器中的CMOS成像器件对地球成像后输入图像。

为了保证最终确定的卫星俯仰角和滚动角精度,要求CMOS成像器件所成地球图像地平轮廓有足够长的弧度。一般要求成像平面有效弧长占到总弧长的50%。

CMOS成像器件所拍摄的图像满足弧长要求后输入图像,进行下一步的处理。

第二,图像预处理。

对图像进行预处理的目的在于去除噪声,平滑图像。图像预处理的方法有很多,如中值滤波,均值滤波等,本发明中采用中值滤波。

中值滤波采用3×3卷积模板,具体做法如下:将中心像素及其3×3邻域内8个像素的灰度值按从小到大排序,取第5个元素灰度值为中心像素的灰度值。中值滤波不仅能去除噪声点,为以后边缘提取减少计算量,还能保护图像的边 缘,获得较满意的图像预处理效果。

第三,检测有效边缘。

这一步是本发明方法为解决上述所提出问题的关键步骤,也是核心部分。检测有效边缘方法的好坏直接决定了本发明方法所需计算量的大小和所确定姿态角的精度。

目前检测边缘的算法有很多种,常见的有:微分算子、拉普拉斯高斯算子、Canny算子、拟合法等。本发明采用Sobel微分算子和Hough变换检测有效边缘,Sobel算子是检测所有的边缘,包括有效边缘和无效边缘,而Hough变换则是剔除Sobel算子所检测到的无效边缘,保留其有效边缘。下面分别介绍这两种方法。

Sobel算子首先计算某一像素点处的灰度梯度幅值,然后将此值与所设的阀值进行比较,如果该点处灰度梯度幅值大于阀值,则保留该点,否则剔除该点。该方法的中心思想是图像边缘像素点邻域内灰度值变化比较剧烈,应用此方法确定出邻域内灰度值变化比较大的像素点,从而检测出图像边缘。Sobel算子中所设阀值的大小不变,阀值大小需要根据实际图像确定,一般阀值的大小选取为不大于整个图像灰度值平均值的最大正整数。

下面介绍Sobel算子计算像素点处灰度梯度幅值的方法。

任一点处梯度幅值计算如下:

M=sx2+sy2---(1)

其中sx,sy分别表示在该点处x、y方向导数,而对于数字图像,可以用一阶差分代替一阶微分。通常计算某一像素点处灰度梯度幅值有3×3模板、5×5模板等,本发明采用3×3模板,故sx,sy可写为:

sx=(a2+2a3+a4)-(a0+a7+a6)(2)

sy=(a0+2a1+a2)-(a6+2a5+a4)(3)

其中a0…a7代表该像素点处3×3邻域内元素的灰度值,a0…a7各像素点灰度值在3×3邻域内位置如下:

a0a1a2a7[i,j]a3a6a5a4

注:[i,j]代表所求灰度梯度幅值像素点坐标。

如前所述,Sobel算子检测到的边缘包括有效边缘和无效边缘,其中有效边缘是地球的轮廓边缘,无效边缘包括地球晨昏线边缘、图像背景明亮星体边缘、月球边缘等。为了提高所确定姿态角的精度,需要对Sobel算子检测到的边缘进一步检测,剔除无效边缘,保留有效边缘。下面介绍检测有效边缘方法之Hough变换。

Hough变换有已知图像圆半径Hough变换和未知图像圆半径Hough变换两种。对于本问题,当然属于后者,因为地球成像后所形成图像圆的半径是无法提前精确获取的,一种快速的未知图像圆半径的Hough变换方法如下:

首先,对于圆的检测,利用梯度信息,Ballard用下式检测圆

x0=x±r>

y0=y±r>

式中,(x0,y0)为圆心坐标,h为圆边缘点(x,y)处梯度方向角,r为图像圆半径。

基于上述公式,检测步骤如下:

1.首先根据先验知识或者已知数据估计参数半径r的变化范围r∈[rmin,rmax]以及步长Δr大小。

对于本问题,地球成像后成像面上近似圆形的半径r范围确定如下:

首先,相机成像示意图如图2所示,图中,R为相机可分辨到的大气层高度,变化范围R∈[6371,6371+120];H为卫星高度,已知值;为半张角,大小RE为地球半径,则由此可确定半张角的范围;f为相机焦距,已知值;r为地球成像后成像面上近似圆形的半径,由可确定r的范围r∈[rmin,rmax];步长Δr由敏感器控制器的计算能力以及精度要求等决定,如果控制器的计算能力强,那步长就可以选的小一点,不仅满足姿态精度 要求,同时满足姿态实时性的要求。

从半径r0=rmin开始,每一次循环半径增加Δr,直至r=rmax,对于任一循环,半径ri=rmin+Δr×i,(i=0,1…m),其中m为循环的总次数,表达式将每一边缘坐标点(xj,yj)(j=1,…n)及其梯度方向角hj依次代入上述(4)和(5)Ballard公式中,计算出圆心坐标(x0ij,y0ij),其中j表示有效和无效边缘点总个数,然后将计算所得的圆心坐标四舍五入取整,上述各个边缘点坐标是在以图像左下角为坐标原点,向右为x轴,向上为y轴的坐标系下表示。为防止计算得到的整数圆心坐标为负数,将原有坐标系分别沿x轴,y轴负方向平移1024个像素点,平移后得到新的正的整数圆心坐标值,以新坐标值为二维数组的坐标值,二维数组的值统计累加该二维坐标值出现的次数,遍历所有的边缘点,最后由累加次数峰值确定此次循环r=ri时的圆参数ci={x0i,y0i,ri},每一次循环确定一组参数ci={x0i,y0i,ri}和累加次数Ki。通过此次循环累加次数Ki和上一次循环累加次数Ki-1的比较,记录比较高的累加次数Kmax以及相应的圆参数cmax={X0,Y0,R},当ri=rmax时,最终确定的圆参数cmax={X0,Y0,R}为有效边缘像素级圆心坐标{X0,Y0}以及半径R。

(4)和(5)Ballard公式均取正号即可,灰度梯度方向角h的正弦和余弦计算如下:

sin(h)=sy/s2x+s2y---(6)

cos(h)=sx/s2x+s2y---(7)

其中sx,sy由(2)和(3)计算可得。

2.未知半径Hough变换得到像素级圆心坐标{X0,Y0}以及半径R,如果某边缘点为有效边缘点,即圆周上的边界点,则该点与像素级圆心坐标{X0,Y0}之间距离满足:

R-1(xj-X0)2+(yj-Y0)2R+1---(8)

不符合关系式(8)的干扰边缘点将被删除,从而获得有效边缘。

第四,应用最小二乘拟合有效边缘中心;

得到有效边缘之后,接下来应用最小二乘拟合有效边缘中心。

对给定的一组数据点(xk,yk)(k=1,…l),l为数据点个数,至少为3,设与该组数据点的距离的平方和为最小的圆的方程为:

(x-x0)2+(y-y0)2=r2(9)

式中(x,y)表示上述圆周上的坐标值,(x0,y0)表示上述圆的圆心坐标,r表示上述圆的半径。

又设a1=∑xk,b1=∑yk,a2=∑(xk)2,b2=∑(yk)2,a3=∑(xk)3,b3=∑(yk)3,c11=∑xk*yk,c12=∑xk*yk2,c21=∑(xk)2*yk

设f(x0,y0,R)=∑((xk-x0)2+(yk-y0)2-r2)2,分别令即

-4∑((xk-x0)2+(yk-y0)2)(xk-x0)=0

-4∑((xk-x0)2+(yk-y0)2-R2)(yk-y0)=0

-4∑((xk-x0)2+(yk-y0)2)r=0

将a1,b1,a2,b2,a3,b3,c11,c12,c21代入上面三个方程整理得:

a1((x0)2+(y0)2)-2a2x0-2c11y0+a3+c12-r2a1=0

b1((x0)2+(y0)2)-2c11x0-2b2y0+c21+b3-r2b1=0

n((x0)2+(y0)2)-2a1x0-2a1x0+a2+b2-lr2=0

解得圆心坐标及半径为:

x0=(2a1b1/l-2c11)(b3+c21-b1(a2+b2)/l)-(2(b1)2/l-2b2)(a3+c12-a1(a2+b2)/l)(2(b1)2/l-2b2)(2(b1)2/l-2b2)-(2a1b1/l-2c11)2---(10)

y0=(2a1b1/l-2c11)(a3+c12-a1(a2+b2)/l)-(2(a1)2/l-2a2)(b3+c21-b1(a2+b2)/l)(2(a1)2/l-2a2)(2(b1)2/l-2b2)-(2a1b1/l-2c11)2---(11)

R=Σ(xi-x0)2+Σ(yi-y0)2/l---(12)

(10)、(11)、(12)计算所得结果是坐标数值结果,没有实际物理意义,第五步中需要有物理意义的结果,设单位像素长度为d,给上述(10)、(11)、(12)计算结果都乘以单位像素长度d,所得结果为有长度物理意义的圆心坐标和半径:

Cx=x0×d(13)

Cy=y0×d(14)

R=r×d(15)

第五,计算卫星对地的滚动角和俯仰角。

在计算卫星姿态角之前,首先建立卫星姿态坐标系O-XYZ(如附图3)和敏感器坐标器O-XTYTZT。卫星姿态坐标系O-XYZ定义如下:原点O为卫星质心,OX轴指向卫星飞行方向,OZ轴指向地心,OY轴由右手定则确定。敏感器坐标器OT-XTYTZT各个坐标轴均与卫星姿态坐标系坐标轴平行,唯一不同原点OT位于敏感器的质心。

根据上述坐标系的定义,姿态角有如下说明:卫星绕X轴旋转为滚动角,绕Y轴旋转为俯仰角,绕Z轴旋转为偏航角,均以绕X轴,Y轴,Z轴正向旋转为正的姿态角。

在此基础上,确定卫星的滚动角φ和俯仰角θ:附图3中,为半张角,由附图2可得式中,R由(15)求得,f已知,则半张角有:

=arctan(Rf)---(16)

设地心在成像面中的坐标为(Mx,My),由图像几何中心即最小二乘所得的有效边缘中心坐标为(Cx,Cy),其中地心在成像面中的坐标(Mx,My)有:

Mx=f*tanθ(17)

My=f*tanφ(18)

推导有效边缘中心坐标(Cx,Cy)的计算表达式,为说明问题简单起见,假设卫星滚动角φ=0,俯仰角θ≠0,则成像面上在X轴方向距离有效边缘中心最远点和最近点的坐标分别为则有效边缘中心X坐标为:Cx=Xmax-Xmin,即:

Cx=f2(tan(+θ)-tan(-θ))---(19)

同理可得:

Cy=f2(tan(-φ)-tan(+φ))---(20)

在滚动角φ和俯仰角θ为小角度下,由(17)和(19)、(18)和(20)可得:

MxCx*cos2()---(21)

My-Cy*cos2()---(22)

故由(17)和(21)、(18)和(22)以及半张角计算公式(16)可得卫星滚动角φ和俯仰角θ有:

φ=arctan(-Cy*cos()f)---(23)

θ=arctan(Cx*cos2()f)---(24)

第六,计算卫星对地偏航姿态角。

首先如图4所示,建立数字式太阳敏感器的测量坐标系o-xssysszss,其中,ozss轴沿码盘平面的法线方向,oxss轴沿太阳光入射的狭缝方向,依右手正交定确定oyss轴,则太阳方向矢量在测量坐标系中的坐标为

Sss=[Sssx>ssy>ssz]T(25)

基于太阳敏感器的测量模型,太阳敏感器的测量值满足:

tanα=SssySssz---(26)

α是太阳方向矢量在太阳敏感器的测量坐标系yssozss平面投影与ozss轴的夹角。

根据太阳星历表数据可以计算出测量时刻太阳在黄道上的赤经和赤纬,依此可以确定太阳矢量S在地心赤道惯性坐标系中的坐标为SI,即:

SI=[cosδs>s>s>s>s]T(27)

式中αss分别为太阳矢量方向的赤经和赤纬,地心赤道惯性坐标系Oe-xIyIzI原点在地心,xI轴指向赤道平面与黄道平面相交节线的升交点(指向春分点γ),zI轴与地球的自旋轴平行,yI轴与xI和yI轴组成右手正交坐标系。

然后,根据测量时刻的轨道根数可以计算出轨道坐标系相对于地心赤道惯性坐标系的转换矩阵,计算公式如下:

Coi=01000-1-100cos>usin>u0-sin>ucos>u0001×1000cos>isin>i0-sin>icos>icosΩsinΩ0-sinΩcosΩ0001---(28)

式中i为轨道倾角,Ω为升交点赤经,u为轨道纬度幅角,轨道坐标系O-xoyozo其原点在卫星的质心,卫星的轨道平面是坐标平面,由质心指向地心的坐标轴是zo轴,xo轴在轨道平面上与zo垂直,指向卫星的速度方向,yo轴与xo,zo组成右手正交坐标系。

则可以导出太阳方向矢量S在轨道坐标系中的坐标S0为:

S0=CoiSI(29)

Coi表示地心惯性坐标系到轨道坐标系的旋转矩阵。

把太阳矢量由轨道坐标系转换成太阳敏感器的测量坐标系的表达式为:

Sss=DssBCboS0(30)

式中:S0由(29)求得,Cbo为由轨道坐标系到卫星本体坐标系的旋转姿态矩阵,如果首先按zB旋转,接着按中间坐标系y轴旋转,最后以新坐标系的x轴旋转,则旋转坐标系Cbo表达式如下;

Cbo(Z-Y-X)=1000cosφsinφ0-sinφcosφcosθ0-sinθ010sinθ0cosθcosψsinψ0-sinψcosψ0001

=cosθcosψcosθsinψ-sinθsinφsinθcosψ-cosφsinψsinφsinθsinψ+cosφcosψsinφcosθcosφsinθcosψ+sinφsinψcosφsinθsinψ-sinφcosψcosφcosθ---(31)

式中卫星的滚动角φ和俯仰角θ已由地球可见光敏感器求得,只有偏航角ψ为未知量,又DssB为太阳敏感器相对卫星本体坐标系的安装矩阵,为已知量,可有如下表达式:

DssB=SAXSAYSAZSBXSBYSBZSCXSCYSCZ---(32)

式中SAX、SAY、SAZ分别表示太阳敏感器xss轴与星体坐标系xb、yb、zb三轴的夹角余弦值,SBX、SBY、SBZ分别表示太阳敏感器yss轴与星体坐标系xb、yb、zb三轴的夹角余弦值,SCX、SCY、SCZ分别表示太阳敏感器zss轴与星体坐标系xb、yb、zb三轴的夹角余弦值。

太阳矢量在轨道坐标系下坐标S0表示成为三轴分量形式为:

So=[SoX>oY>oZ]T(33)

又记

Cyx(θ,φ)=1000cosφsinφ0-sinφcosφcosθ0sinθ010-sinθ0cosθ---(34)

将式(31)、(32)带入(30),可得:

SssY=SBXSBYSBZCbo(Z-Y-X)SoXSoYSoZ=SBXSBYSBZCyx(θ,φ)SoXSoY0SoY-SoX000S0Zcosψsinψ1---(35)

SssZ=SCXSCYSCZCbo(Z-Y-X)SoXSoYSoZ

=SCXSCYSCZCyx(θ,φ)SoXSoY0SoY-SoX000S0Zcosψsinψ1---(36)

将式(35)、(36)带入(26),可得:

A1cosψ+A2sinψ+A3=0(37)

求此方程可以求得卫星对地的偏航角为:

ψ=-arcsinA3A12+A22-arctanA1A2---(38)

式中A1,A2,A3可由下式求得:

A1A2A3=SoXSoY0SoY-SoX000S0Z×cosθ0-sinθ010sinθ0cosθ1000cosφsinφ0-sinφcosφ{SBXSBYSBZ-tanαSCXSCYSCZ}---(39)

首先卫星携带的可见光地球敏感器的CMOS成像器件对空间中地球成像,然后输入图像,其次检测所得图像边缘以及剔除干扰边缘,从而得到地球成像图像的最终有效边缘,再者应用最小二乘拟合出有效边缘的中心位置,接着,根据图像中心位置坐标及卫星姿态角为小角度假设下近似计算卫星滚动和俯仰两轴姿态角,最后,利用太阳敏感器的测量原理以及已计算得到的卫星对地滚动和俯仰姿态角信息确定卫星的偏航角。至此,卫星对地的三轴姿态角均已确定,最终可实现卫星对地球的可持续观察。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号