星敏感器姿态确定仿真综合实验.doc

。北京航空航天大学专业综合实验报告 学 院 宇航学院 班 级 111514 学 号 11151146 姓 名 高荣荣 指导老师 王海涌 2015年1月3日星敏感器姿态确定仿真综合实验摘要：通过对电子星图模拟器和星敏感器PC仿真平台的操作，实现星敏感器姿态确定，以及借助星象天文馆，来实现电子星图模拟器对星的标定。完成星敏感器系统仿真。关键词：星敏感器定姿 星图模拟 星图姿态矩阵一 实验目的通过电子星图模拟器（ESS）和星敏感器PC仿真平台的操作，熟悉星图模拟和星敏感器姿态基本流程及各模块功能，完成星敏感器系统仿真。二 实验原理1.星图模拟原理（1）星图模拟系统是一种近似模拟星空的仿真系统。为星敏感器算法调试、星敏感器产品测试及天文导航半物理仿真系统运行提供标准的星图输入，并提供已知参考星光矢量及星像中心的理想映射坐标。分光学物理星模、电子星模和计算机软件星模三种。（2）星图模拟系统的实现旋转关系：OZ是光轴指向： 星图模拟是根据光轴指向及旋角(0,0,)所确定的既定视场，将其范围内的星空目标映射到星敏感器CCD面阵上并模拟出目标图像的过程。其中，需要确定第二赤道坐标系、航天器本体坐标系、星敏感器坐标系的转换矩阵。星敏感器固联（安装矩阵为常数阵），那么只考虑第二赤道坐标系和星敏感器坐标系之间的转换关系。令O-UVW为第二赤道坐标系，令O-XYZ为星敏感器坐标系。那么星光矢量在两个坐标系下的分量列阵的关系可以表示为：X,Y,ZT = TsiU,V,WT，其中Tsi为转换矩阵。2.星敏感器定姿基本原理OsXsYsZs 星敏感器坐标系Ouv CCD 成像面坐标系OsO之间距离 f 为光学透镜的焦距由图中的几何关系可得：第 n 颗星的单位矢量在星敏感器坐标系中的分量列阵：或，根据星像点质心坐标直接计算得单位矢量：式中VS为星敏感器测量误差矢量。星敏视场内 n 颗恒星在星敏感器坐标系Sb中的单位矢量坐标分别为 (X1, Y1, Z1), (X2, Y2, Z2), , (Xn, Yn, Zn)。经过星图匹配，获知这 n 颗恒星在惯性空间 i 系中单位矢量坐标：U1, V1, W1, U2, V2, W2, , Un, Vn, Wn，是由星表中赤经赤纬经过直角坐标转换求得。则存在以下矩阵转换关系式：上式各矩阵赋予命名，对应关系式：S = CTis当观测星数 n 等于3颗时，星历矩阵C33为非奇异矩阵，那么转移矩阵Tis C-1S。n 3时，可以采用最小二乘法得到计算式：Tis = (CTC)-1CTS 姿态矩阵Tsi = TsbTblTli，则载体相对于发射点惯性系l系的姿态矩阵 Tbl = (Tsb)-1 (Tis)-1(Tli)-1 =(Tli Tis Tsb)-1，这里 Tsb 为星敏感器固联常量安装矩阵，发射点惯性坐标系 l系相对于赤道惯性系i系的转移矩阵：其中：A 为轴xl的方位角；S 为发射时刻的格林尼治恒星时，或发射时刻的春分点的格林时角GHA；(, ) 为发射点经纬度。如果安装矩阵与火箭本体坐标系方向重合，则安装矩阵Tbs = 单位阵 I那么：Tbl = (Tli Tis)-1欧拉角法表示火箭姿态发惯系l，原点在发射点L，轴yl在当地铅垂线向上，轴 xl 在当地水平，且在名义射击平面内。发射点惯性坐标系 Lxlylzl转动 3 次转到本体坐标系 oxbybzb 。转动顺序：那么，从Sl到Sb的坐标变换矩阵：矩阵元素对照得主值：姿态角和的取值范围都在-90,90；的取值范围都在-180,180然后基于主值再获得真值。三 实验步骤1. 设定观星地址和观星时间，采用静基座天顶观测方式，解算出相机光轴的赤经、赤纬和旋角。（1）设定观星地点为北京市延庆县马匹营村，北纬403014.571，东经1160420.815（2）设定北京时间为2011年8月20日22h22min222s（3）2011年8月20日当天世界时0时对应的恒星时为12h16m40.6407s; GHAr=S0+(1+)M=26.3156由此得RA=300.2372 Dec=40.5041指定 =02开启ESS，进行参数设置，选择“星图模拟” “单帧重复”菜单命令，将（ ， ， ）计算值填入ESS的人机界面，观察“发送的星图”，此时ESS持续不断地将模拟星图数字视频传输给星敏感器PC仿真平台。设置：指定仿真星等阈值：6峰值饱和灰度对应星等：8得到模拟星图：其中，星模软件使用方法见附录；3.在PC仿真平台，存储接收的模拟星图，并选定六颗颗亮星。结果如图所示。得到亮星的坐标，星等，赤经和赤纬：310.358041845.280358371.3243.0664099873.2034322305.557124840.256699942.3514.3430786718.3319065311.552858233.970255722.6810.0597353993.0520476296.243895845.130739033.0272.523317369.3284183303.407916546.741358354.0190.4634936621.7853499299.076604335.083419324.0783.1009519466.73437074.运行虚拟天文馆软件stellarium.exe，选定（1）中地址和观星时间，目视对比模拟星图和stellarium界面，找出这六颗选定亮星的匹配星，记录其星号、赤经、赤纬和星等属性信息，计算在春分点赤道惯性系i系下的恒星方位矢量。运行虚拟天文馆结果如图所示：Cyg-18O1 Cyg-31天津四Cyg-21天津九天津一根据经纬度关系，并借助主观肉眼判断，匹配情况如下：星名/星等赤经/（J2000）赤纬天津四/1.25（脉动变星）310.25/20h41m25s45.16天津一/2.20（恒星）305.52/20h22m13.7s40.17天津九/2.45（恒星）311.50/20h16m41s34.01Cyg-18 Cyg-HIP 97165 A/2.9（双星）296.25/19h45m20s45.07O1 Cyg-31 Cyg-V0695/3.8（食双星系统）303.404/20h13m37.9s46.46Cyg-21 Cyg-HIP 98110/3.85(恒星)299.011/19h56m45s35.06注意事项：（1）在运用虚拟天文馆软件时，要注意调节观察视角，一般为20度左右事宜。本次实验，调节观察视角读数为22.8。（2）根据纬度关系显示，将虚拟天文馆的视场调节到指定的范围内，以便容易找到星星。本次实验，大约根据纬度关系调节到45的范围内。（3）在记录赤经、赤纬的信息时，有两种表示的方法：一种是J2000，另一种是J2011.2。两者之间有非常细微的差距。因为地轴的进动，使得春分点不断西退，而赤经赤纬又是依靠春分点的位置来确定，所以随着春分点位置的不断变化，赤经赤纬也会不断变化。J2000版指的是以2000年春分点作为基准天体的坐标，日期版就是以现在的春分点作为基准的坐标。差距是非常细微的，对一般的爱好者来说，忽略这个微小的差距也无妨，但是对于专业的观测来说，这点误差是要命的。在本实验中，记录的是J2000下的赤经赤纬信息。记录星等信息时，有目视星等和绝对星等两种星等信息。绝对星等是假定把恒星放在距地球10秒差距（32.6光年）的地方测得的恒星的亮度，用以区别于视星等。它反映天体的真实发光本领。由于本实验中记录的是目视星等。5.给定该相机焦距f的标定值，查找星敏感器PC平台下的日志文件获取6颗选定亮星的星象质心，是在星图预处理环节由灰度重心法计算获得的，计算星敏感器本体系s系下的星光矢量。 (1).焦距f = 2889.2643，单位：pixel，主点O的值(508.3, 513.3) pixel。 (2).阵列平面坐标系和星敏感器本体坐标系的关系如图所示： 其中星敏感器本体系s系下的星光矢量由下式计算得到 6.基于i系下的恒星方位矢量和s系下的星光矢量，计算出观星时刻星敏感器相对于i系的姿态矩阵，并换算出星敏感器的光轴赤经赤纬和旋角，对比星敏感器PC仿真平台的计算结果。Matlab程序见附录；得到赤经赤纬及旋角为：RA=300.2372Dec=40.5041四 实验结论及误差分析：1.以观测者所在位置，运用查表法得到的赤经赤纬为Ra=300.2372，Dec=40.5041，Gamma=0；再运用模拟星图ESS和虚拟天文馆得到的赤经赤纬为RaDeter=300.2385，DecDeter=40.5180 ，GammaDeter=-0.3402;所以，经运算得到绝对误差和相对误差如下：绝对误差为：|Ra-RaDeter|=0.0013 |Dec-DecDeter|=0.0139 |Gamma-GammaDeter|=0.3402相对误差为： (Ra)=0.004%(Dec)=0.031%(Gamma)=0.340%2.在星图模拟和虚拟天文馆的对比找星过程中，也存在一定的误差，主要是，星等误差，赤经和赤纬对比误差等等。误差定量计算如下：（1）赤经绝对误差：绝对误差310.25/20h41m25s310.35804180.100305.52/20h22m13.7s305.55712480.030311.50/20h16m41s311.55285820.050296.25/19h45m20s296.24389580.010303.404/20h13m37.9s303.40791650.003299.011/19h56m45s299.07660430.065分析：最大误差为0.1，因此，在合理误差范围之内，说明星星确认正确；（2）赤纬绝对误差：绝对误差45.1645.280358370.1240.1740.256699940.0834.0133.970255720.0445.0745.130739030.0646.4646.741358350.2835.0635.083419320.02分析：最大误差为0.28，计算相对误差为=0.6%，误差总体较小，因此也在合理范围之内。说明确认星星正确；（3）星等绝对误差绝对误差1.251.30.052.202.30.102.452.60.052.93.00.103.84.00.203.854.00.15分析：最大误差为0.20，计算相对误差一定比较大，但是不影响星星判定，星等的误差，大多还是由于虚拟天文馆的制作造成的，因此，从整体绝对误差上看，星星的误差还是在合理范围之内的。因此，认为确定星星正确。综上所述：经过赤经，赤纬，星等的误差对比，确认模拟星图所找的星星在虚拟天文馆中得到一一对应。误差可视为合理误差。五 实验总结 为期一周的综合实验到这里结束了，首先，要感谢王海涌老师的悉心教导，在王老师的帮助指导下，将天文导航这门课程的精髓在实验中得到了应用，收获颇丰。 本实验，主要就图模拟原理和星敏感器定姿原理展开了实验上的运用。我们先利用观星地址和观星时间解算出光轴的赤经、赤纬和旋角。再计算出由春分点确定的子午线和格林尼治子午线的时角。第一步是根据观星时间查紫金山天文台出版的天文年历确定当日格林尼治零时的春分点时角，第二步由观星时间确定观星时刻春分点所在子午线和格林尼治子午线的时角。进而计算由春分点第二赤道坐标系i系和地球坐标系e系间的基元旋转矩阵，由假定测着所在地理位置计算地理坐标系t系向地球坐标系e系的转移矩阵，从而确定地理坐标系t系向第二赤道惯性坐标系i系的转移矩阵，即可算出光轴的赤经、赤纬和旋角。通过计算观星时刻的世界时计算儒略日从而计算地球转角、TT儒略世纪数t ，最后得到格林尼治平恒星时GMST，再用上测者的位置信息，即可求得t系向i系的转移矩阵。这种算法的好处是避免了查天文年历，比较方便。 在运行虚拟天文馆时，一定要注意光轴是和地面垂直这一条件。按住键盘的向上的按键一直到头，则界面显示的就是正上方的星空。同时也要注意调整视角的大小，因为星图模拟器的视角是2020，要调整虚拟天文馆的FOV至合适值，大约是22.8。在使用虚拟天文馆时，要注意用已知信息和实际信息进行匹配，如，运用在模拟星图中的赤经和赤纬的信息在虚拟天文馆中进行信息对比，如果位置正确，注意调节视角等再进行匹配。直至匹配完成。 附录：Matlab程序：1电算法% 假定在北京古观象台（北纬403014.571，东经1160420.815）进行观星Lat=(39+54/60)*pi/180;Lon=(116+18/60)*pi/180;% 电算法解算格林尼治平恒星时%时间为北京时2011年8月20日22点22分22.2秒format long;year = 2011; % 按 UTC和UT1近似相等处理month = 8;day = 22;hour = 22;minute = 22;second = 22.2;%估计预测DUT1a = floor(14-month)/12);y = year+4800-a;m = month+12*a-3;JDN = day+floor(153*m+2)/5)+floor(365+97/400)*y)-32045.5;%求地球自转角UTC_TAI = -34; % 需查2011年的bulletina-xxiv-017.txt,文档中已给JdUT1 = JDN+(hour+minute/60+second/3600)/24; Tu = JdUT1 - 2451545;Theta = 2*pi*(0.7790572732640+1.00273781191135448*Tu)*180/pi*3600; %计算力学时TT，求GMSTsecTT = second - UTC_TAI + 32.184; % 力学时 0hTTJdTT = JDN+(hour+minute/60+secTT/3600)/24; t = (JdTT-2451545);t = t/36525;GMST = 0.014506+Theta+4612.15739966*t+1.39667721*t.*t-0.00009344*t.3+0.00001882*t.4; % 注意：单位是角度秒，arcsecondGMSTarc = deg2rad(mod(GMST/3600,360) % 得到格林尼治平恒星时temp = rem(GMST/15,24*3600); % GMST/15则为时间秒H = floor(temp/3600); % 恒星时“h”M = floor(mod(temp,3600)/60); % 恒星时“m”S = mod(temp,3600)-M*60; % 恒星时“s”% 计算天顶方向在惯性系i系下的经纬度%旋转矩阵Fai=pi/180;Pfai=Fai;Thet=Fai;L1=1 0 0; 0 cos(Fai) sin(Fai); 0 -sin(Fai) cos(Fai);L2=cos(Pfai) 0 -sin(Pfai); 0 1 0; sin(Pfai) 0 cos(Pfai);L3=cos(Thet) sin(Thet) 0; -sin(Thet) cos(Thet) 0; 0 0 1; %转移矩阵 Cts=L1*L2*L3; Cet=-sin(Lon) -sin(Lat)*cos(Lon) cos(Lat)*cos(Lon); cos(Lon) -sin(Lat)*sin(Lon) cos(Lat)*sin(Lon); 0 cos(Lat) sin(Lat); Cie=cos(GMSTarc) -sin(GMSTarc) 0; sin(GMSTarc) cos(GMSTarc) 0; 0 0 1; R=Cie*Cet*Cts*0;0;1; %解算赤经赤纬 Alpha=degrees2dms(rad2deg(atan(R(2,:)/R(1,:), Delta=degrees2dms(rad2deg(asin(R(3,:),2.matlab下打开raw文件；clc;clear all;fid = fopen(0.raw, r);mk=1024;nk=1024;z=uint16(zeros(mk,nk);for i=1:mk for j=1:nk z(i,j)=fread(fid, 1,uint16); end end g=double(max(max(z); temp=double(255/g);for i=1:mk for j=1:nk A(i,j)=z(i,j)*temp; endendA=uint8(A);figure;imshow(A)（3）计算赤经赤纬：u0=508.3;v0=513;f=2889.2643;star1 = 113.6499751 31.88848347;star2 = 134.8018832 48.04160012;star3 = 135.9063711 47.15658602;star4 = 112.2777956 31.78433636;star5 = 135.1599624 41.78273636;star6 = 125.7089622 43.18811131;Vstar1 = cosd(star1(1)*cosd(star1(2),sind(star1(1)*cosd(star1(2),sind(star1(2);Vstar2 = cosd(star2(1)*cosd(star2(2),sind(star2(1)*cosd(star2(2),sind(star2(2);Vstar3 = cosd(star3(1)*cosd(star3(2),sind(star3(1)*cosd(star3(2),sind(star3(2);Vstar4 = cosd(star4(1)*cosd(star4(2),sind(star4(1)*cosd(star4(2),sind(star4(2);Vstar5 = cosd(star5(1)*cosd(star5(2),sind(star5(1)*cosd(star5(2),sind(star5(2);Vstar6 = cosd(star6(1)*cosd(star6(2),sind(star6(1)*cosd(star6(2),sind(star6(2); uv1 = 895.2076186 74.50844256; uv2 = 66.40130886 885.4437202;uv3 = 105.8599846 929.8961347;uv4 = 894.6814767 13.83111223; uv5 = 382.7642809 940.7513796;uv6 = 341.0580967 581.6978552;VsensorBody1 = uv1(2)-v0,-(uv1(1)-u0),f/sqrt(uv1(1)-u0)2+(uv1(2)-v0)2+f*f);VsensorBody2 = uv2(2)-v0,-(uv2(1)-u0),f/sqrt(uv2(1)-u0)2+(uv2(2)-v0)2+f*f);VsensorBody3 = uv3(2)-v0,-(uv3(1)-u0),f/sqrt(uv3(1)-u0)2+(uv3(2)-v0)2+f*f);VsensorBody4 = uv4(2)-v0,-(uv4(1)-u0),f/sqrt(uv4(1)-u0)2+(uv4(2)-v0)2+f*f);VsensorBody5 = uv5(2)-v0,-(uv5(1)-u0),f/sqrt(uv5(1)-u0)2+(uv5(2)-v0)2+f*f);VsensorBody6 = uv6(2)-v0,-(uv6(1)-u0),f/sqrt(uv6(1)-u0)2+(uv6(2)-v0)2+f*f); S = VsensorBody1;VsensorBody2;VsensorBody3;VsensorBody4;VsensorBody5;VsensorBody6;C = Vstar1;Vstar2;Vstar3;Vstar4;Vstar5;Vstar6;Tsi = (S*S)S*C; % 采用最小二乘法 % * 浮点计算条件下，解算星敏光轴的赤经赤纬和旋角(0,0,)*DecDeter = asin(Tsi(3,3); % 毫无疑问!真值就是主值，单位：弧度ref = sind(89.99);if abs(Tsi(3,3)ref % 赤纬接近90 RaDeter = 0; % 赤纬位于天南极或天北极，赤经可以是任意值，直接令其指向春分点吧！else if Tsi(3,2)/cos(DecDeter)ref % 判断sin(0)，0接近90 RaDeter = 90; elseif Tsi(3,2)/cos(DecDeter)-ref % 0接近270 RaDeter = 270; else RaDeter = atand(Tsi(3,2)/Tsi(3,1