CN113361163B - A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction - Google Patents

A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction Download PDF

Info

Publication number
CN113361163B
CN113361163B CN202110619180.1A CN202110619180A CN113361163B CN 113361163 B CN113361163 B CN 113361163B CN 202110619180 A CN202110619180 A CN 202110619180A CN 113361163 B CN113361163 B CN 113361163B
Authority
CN
China
Prior art keywords
irradiance
satellite
target satellite
earth
surface element
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110619180.1A
Other languages
Chinese (zh)
Other versions
CN113361163A (en
Inventor
汪夏
黎吉顺
张雅声
徐灿
李智
方宇强
贺俊
李鹏举
程文华
刁华飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority to CN202110619180.1A priority Critical patent/CN113361163B/en
Publication of CN113361163A publication Critical patent/CN113361163A/en
Application granted granted Critical
Publication of CN113361163B publication Critical patent/CN113361163B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/36—Guiding or controlling apparatus, e.g. for attitude control using sensors, e.g. sun-sensors, horizon sensors
    • B64G1/363—Guiding or controlling apparatus, e.g. for attitude control using sensors, e.g. sun-sensors, horizon sensors using sun sensors

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computer Hardware Design (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Navigation (AREA)
  • Photometry And Measurement Of Optical Pulse Characteristics (AREA)

Abstract

The invention provides a satellite attitude estimation method for correcting terrestrial reflected light, which comprises the following steps: dividing the earth surface into a plurality of surface elements according to longitude and latitude grids, and screening effective surface elements; according to the inverse proportion law of the radiation illumination and the square of the distance, calculating the radiation illumination of the sunlight visible light wave band at each effective surface element, and according to the Lambert reflection principle, calculating the radiation illumination of the effective surface elements at the target satellite; according to the radiation illumination generated by the effective surface elements at the target satellite, calculating the radiation illumination of the target satellite to the detector by taking each effective surface element as a light source according to the Lambert reflection principle; adding the radiation illumination of the target satellite to the detector to obtain the radiation illumination of the target satellite to the detector when the earth is used as a light source; subtracting the real-time radiation illumination of the target satellite to the detector when the earth is used as a light source from the real-time radiation illumination of the target satellite to the detector obtained from the target satellite, and obtaining corrected observation data of the radiation illumination of the satellite to the detector when only the sun is used as the light source; and performing attitude estimation on the satellite by using a tasteless Kalman filtering method to obtain the correction data of the earth reflected light. The invention accurately calculates the illumination of the earth reflected light to the satellite and the illumination of the satellite under the condition, and the satellite attitude estimation based on the illumination can greatly improve the estimation precision.

Description

一种地球反射光校正的卫星姿态估计方法A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction

技术领域technical field

本发明属于光学特性识别技术领域，特别涉及地球反射光对中高轨卫星光学观测影响的计算。The invention belongs to the technical field of optical characteristic identification, and particularly relates to the calculation of the influence of earth reflected light on the optical observation of medium and high orbit satellites.

背景技术Background technique

地基光学观测是获取卫星特征信息的一种重要手段。太阳作为一个恒定光源，发出的辐射光线经过卫星表面的反射后进入探测器入瞳。根据入射光强度不变，结合光线传递的原理，卫星形状、姿态和表面材质的不同均会造成探测器接收到的光学信息发生变化。因此，基于这一原理，通过获取的卫星光度信息可以识别出卫星的形状、大小、工作状态等信息。然而现实观测中存在太阳以外的其它杂散光源，导致探测器接收到的光信号是由卫星对多个角度的入射光反射得到的混合光。Ground-based optical observation is an important means to obtain satellite characteristic information. The sun acts as a constant light source, and the radiated light emitted by the satellite enters the entrance pupil of the detector after being reflected on the surface of the satellite. According to the constant incident light intensity, combined with the principle of light transmission, the difference in the shape, attitude and surface material of the satellite will cause the optical information received by the detector to change. Therefore, based on this principle, the shape, size, working status and other information of the satellite can be identified through the obtained satellite luminosity information. However, there are other stray light sources other than the sun in actual observations, so that the light signal received by the detector is a mixed light reflected by the satellite to incident light from multiple angles.

目前已知的主要杂散光源为地球，地球是比较特殊的光源，地球不具备自身发光的能力，而是对太阳发出的光进行反射以此来照亮卫星。由于太阳在同一时刻内只能照亮地球的一半，因此，卫星接收到的间接太阳辐射光强度与太阳、地球、卫星的几何位置关系有着密不可分的关系。The main stray light source known at present is the earth. The earth is a special light source. The earth does not have the ability to emit light by itself, but reflects the light emitted by the sun to illuminate the satellite. Since the sun can only illuminate half of the earth at the same time, the indirect solar radiation intensity received by the satellite is inseparable from the geometrical relationship between the sun, the earth and the satellite.

研究初期，基于地球较低的反照率，研究者们普遍认为它们对卫星光学观测造成的影响可以忽略不计，因而只关注太阳的直接照射。但是随着观测技术的发展及研究的深入，一些不合理现象出现在了日常卫星观测中，在对不同的卫星进行多相位的光学观测后，发现它们的光度曲线在大相位角处的变化趋势均趋于平缓，这说明在大相位角处，地球反射光对卫星光度的贡献要大于太阳。另外，大多卫星形状、材料不规则，导致探测器接收到的光信号对观测几何位置关系(光源、目标卫星、探测器的相对位置关系)极为敏感，往往需要同时考虑太阳、地球、卫星、探测器四者的相对位置关系。At the beginning of the study, based on the low albedo of the Earth, researchers generally believed that their impact on satellite optical observations was negligible, so they only focused on direct sunlight. However, with the development of observation technology and in-depth research, some unreasonable phenomena appear in daily satellite observations. After multi-phase optical observation of different satellites, it is found that their luminosity curves change at large phase angles. All tend to be flat, which means that at large phase angles, the contribution of the Earth's reflected light to the satellite's luminosity is greater than that of the sun. In addition, most of the satellites have irregular shapes and materials, which makes the optical signals received by the detectors extremely sensitive to the geometrical position relationship of the observation (the relative positional relationship between the light source, the target satellite, and the detector). The relative positional relationship of the four devices.

近年来，人们已经逐渐意识到间接太阳辐射给卫星光学观测带来的不可忽视的影响。基于现有的天体反射模型，分析了地球反射光对卫星观测的影响。但大多研究者们探究的目的在于验证白天利用地球反射光进行卫星观测的可行性。他们着重研究的是复合光线的谱形，仅对单一时刻或某几个时刻下的目标反射光信号进行了分析。并且没有详细比较在不同时刻下地球反射光相对于直接太阳辐射照度的影响程度，对于可否忽略间接辐射没有给出一个定量的标准。因此，在实际的应用中，未经校正的卫星观测数据将导致卫星状态估计精度降低。In recent years, people have gradually realized the non-negligible impact of indirect solar radiation on satellite optical observations. Based on the existing celestial reflection model, the influence of the Earth's reflected light on satellite observations is analyzed. However, the purpose of most researchers' exploration is to verify the feasibility of satellite observations using the reflected light of the earth during the day. They focused on the spectral shape of the composite light, and only analyzed the reflected light signal of the target at a single moment or at certain moments. And there is no detailed comparison of the influence of the Earth's reflected light relative to the direct solar irradiance at different times, and there is no quantitative standard for whether the indirect radiation can be ignored. Therefore, in practical applications, uncorrected satellite observation data will lead to a decrease in the accuracy of satellite state estimation.

发明内容SUMMARY OF THE INVENTION

本发明提供了一种地球反射光校正的卫星姿态估计方法，解决了地球反射光的定量计算的技术问题以及对于传统的基于光度数据的卫星姿态估计准确率较低的技术问题，本发明的目的在于对获得的卫星观测数据进行校正，消除地球反射光对卫星观测数据的影响，从而提高卫星的姿态估计精度。在研究卫星的光学特性过程中，地球反射光的存在会改变探测器探测到的卫星的亮度，会对卫星特征分析过造成较大影响。本发明精确计算了地球反射光在不同时段影响的大小，基于计算结果对卫星观测数据进行校正，通过无味卡尔曼滤波的方法进行姿态估计后，提高了估计精度。The invention provides a satellite attitude estimation method for earth reflection light correction, which solves the technical problem of quantitative calculation of earth reflection light and the technical problem that the accuracy of traditional satellite attitude estimation based on photometric data is low. The purpose is to correct the obtained satellite observation data to eliminate the influence of the earth's reflected light on the satellite observation data, thereby improving the accuracy of satellite attitude estimation. In the process of studying the optical characteristics of satellites, the existence of reflected light from the earth will change the brightness of the satellites detected by the detector, which will have a great impact on the analysis of satellite characteristics. The invention accurately calculates the influence of the earth reflected light in different time periods, corrects the satellite observation data based on the calculation results, and improves the estimation accuracy after the attitude estimation is performed by the tasteless Kalman filtering method.

为达到上述目的，本发明提供的一种地球反射光校正的卫星姿态估计方法，包括以下步骤：In order to achieve the above purpose, a satellite attitude estimation method for earth reflection light correction provided by the present invention comprises the following steps:

步骤一：将地球表面按照经纬网格划分为多个面元，计算面元面积；基于太阳、地球和目标卫星之间的遮挡关系筛选有效面元；Step 1: Divide the earth's surface into multiple surfels according to the latitude and longitude grid, and calculate the surfel area; screen the effective surfels based on the occlusion relationship between the sun, the earth and the target satellite;

步骤二：根据辐射照度与距离的平方成反比定律，计算太阳可见光波段在各有效面元处的辐射照度，并根据朗伯反射原理计算有效面元在目标卫星处产生的辐射照度；Step 2: According to the law of irradiance being inversely proportional to the square of the distance, calculate the irradiance of the solar visible light band at each effective surface element, and calculate the irradiance generated by the effective surface element at the target satellite according to the Lambertian reflection principle;

步骤三：根据有效面元在目标卫星处产生的辐射照度，通过朗伯反射原理计算出以各个有效面元为光源时，目标卫星对探测器的辐射照度；将目标卫星对探测器的辐射照度相加得到以地球为光源时，目标卫星对探测器的辐射照度；Step 3: According to the irradiance generated by the effective surface element at the target satellite, the irradiance of the target satellite to the detector is calculated by the Lambertian reflection principle when each effective surface element is used as the light source; the irradiance of the target satellite to the detector is calculated. Add up to get the irradiance of the target satellite to the detector when the earth is the light source;

步骤四：从目标卫星处获得实时的目标卫星对探测器的辐射照度减去以地球为光源时，目标卫星对探测器的辐射照度，得到了仅以太阳为光源时的目标卫星对探测器的辐射照度的校正观测数据；Step 4: Obtain the real-time irradiance of the target satellite to the detector from the target satellite minus the irradiance of the target satellite to the detector when the earth is used as the light source, and obtain the irradiance of the target satellite to the detector when only the sun is used as the light source. Corrected observation data of irradiance;

步骤五：基于校正观测数据，使用无味卡尔曼滤波的方法对目标卫星进行姿态估计，得到地球反射光校正数据。Step 5: Based on the corrected observation data, use the method of tasteless Kalman filtering to estimate the attitude of the target satellite, and obtain the correction data of the earth's reflected light.

所述步骤一中，将地球表面按照经纬网格划分为多个面元的方法包括：基于地球表面的经纬网格，将面元中心位于经纬点，长和宽分别为一个纬度距离和一个经度距离。In the step 1, the method of dividing the surface of the earth into a plurality of bins according to the latitude and longitude grid includes: based on the latitude and longitude grid of the earth surface, the center of the bin is located at the latitude and longitude point, and the length and width are respectively a latitude distance and a longitude. distance.

计算面元面积的方法包括：每个地表面元的长均为lj＝2πRE/360；纬度i处的地表面元，宽为di＝2πREcosi/360；地表面元面积为sij＝ljdi；其中，lj为经度j处地表面元的长，di为纬度i处的地表面元的宽，sij为面元面积，RE为地球半径，i、j分别表示纬度和经度。The method for calculating the area of the surface element includes: the length of each surface element is l j = 2πRE /360; the width of the surface element at latitude i is d i = 2πRE cosi /360; the area of the surface element is s ij =l j d i ; where l j is the length of the ground surface element at longitude j, d i is the width of the ground surface element at latitude i, s ij is the area of the surface element, RE is the radius of the earth, i , j represent latitude and longitude, respectively.

筛选有效面元的方法包括：Methods for filtering valid bins include:

Figure BDA0003099037250000031
Figure BDA0003099037250000031

其中，θin为太阳光对地表面元的入射天顶角，θout为反射天顶角。Among them, θ in is the incident zenith angle of sunlight to the ground surface element, and θ out is the reflection zenith angle.

所述步骤二中，计算太阳可见光波段在各有效面元处的辐射照度的方法包括：根据辐射照度与距离的平方成反比定律，在太阳光不被遮挡的情况下，得到太阳可见光波段在各有效面元处的辐射照度：In the second step, the method for calculating the irradiance of the solar visible light band at each effective surface element includes: according to the law that the irradiance is inversely proportional to the square of the distance, under the condition that the sunlight is not blocked, obtain the solar visible light band in each effective area. The irradiance at the effective bin:

Figure BDA0003099037250000032
其中，R0为一年内日地的平均距离；Q0是在距离太阳R0处，可见光波段(0.4～0.7μm)部分的辐射照度，R为地表面元到太阳的距离，Q是在距离太阳R的有效面元处，对可见光波段(0.4～0.7μm)的辐射照度。
Figure BDA0003099037250000032
Among them, R 0 is the average distance between the sun and the earth in one year; Q 0 is the irradiance in the visible light band (0.4-0.7 μm) at the distance from the sun at R 0 , R is the distance from the surface element to the sun, and Q is the distance between the sun and the sun. The irradiance of the visible light band (0.4-0.7μm) at the effective surface element of the sun R.

根据朗伯反射原理计算有效面元在目标卫星处产生的辐射照度，得出目标卫星所受单个地表面元的辐射照度的方法包括：地球表面太阳辐射照度为：

Figure BDA0003099037250000041
其中，R0为一年内日地的平均距离，Rse为太阳到地球的距离，Qse为太阳在地球表面的辐射照度，Q0为在距离太阳R0处，对可见光波段的辐射照度；According to the principle of Lambertian reflection, the irradiance generated by the effective surface element at the target satellite is calculated, and the method to obtain the irradiance of a single ground surface element received by the target satellite includes: the solar irradiance on the earth's surface is:
Figure BDA0003099037250000041
Among them, R 0 is the average distance between the sun and the earth in one year, R se is the distance from the sun to the earth, Q se is the irradiance of the sun on the surface of the earth, and Q 0 is the irradiance of the visible light band at the distance R 0 from the sun;

地表面元对各方向产生的辐射亮度为：

Figure BDA0003099037250000042
其中，θin为太阳光对地表面元的入射天顶角，Qsecosθin为地表面元所接收的太阳辐射照度，地球表面及地表大气总的平均反照率值为
Figure BDA0003099037250000043
Lij为地表面元对各方向产生的辐射亮度；The radiance generated by the surface element in each direction is:
Figure BDA0003099037250000042
Among them, θ in is the incident zenith angle of sunlight to the surface element, Q se cos θ in is the solar irradiance received by the surface element, and the total average albedo value of the earth's surface and the surface atmosphere is
Figure BDA0003099037250000043
L ij is the radiance generated by the ground surface element in all directions;

单个地表面元(i,j)在

Figure BDA0003099037250000044
方向目标卫星处产生的辐射照度为：
Figure BDA0003099037250000045
其中，
Figure BDA0003099037250000046
为面元指向目标卫星的方向向量，
Figure BDA0003099037250000047
为地表面元相对目标卫星所呈立体角，sij为经纬度位于(i,j)的地表面元面积，
Figure BDA0003099037250000048
为地表面元至目标卫星的距离，Qij为单个地表面元(i,j)在
Figure BDA0003099037250000049
方向目标卫星处产生的辐射照度。A single ground element (i, j) is
Figure BDA0003099037250000044
The irradiance generated at the target satellite is:
Figure BDA0003099037250000045
in,
Figure BDA0003099037250000046
is the direction vector of the surfel pointing to the target satellite,
Figure BDA0003099037250000047
is the solid angle of the surface element relative to the target satellite, s ij is the area of the surface element located at (i, j) in latitude and longitude,
Figure BDA0003099037250000048
is the distance from the surface element to the target satellite, and Q ij is the single surface element (i, j) in
Figure BDA0003099037250000049
The irradiance generated at the target satellite.

所述步骤三中，根据有效面元在目标卫星处产生的辐射照度，通过朗伯反射原理计算出以各个有效面元为光源时，目标卫星对探测器的辐射照度的方法包括：

Figure BDA00030990372500000410
其中，太阳指向目标的向量即入射光向量为λS；目标指向探测器的向量即观测向量为λF；目标卫星某面元的法向量为n；入射向量与面元法向量的夹角为到入射天顶角，记作θS，观测向量与面元法向量的夹角为观测天顶角，记作θF；根据入射天顶角θS、观测天顶角θF的大小关系判断探测器是否能收到目标卫星反射光；In the third step, according to the irradiance generated by the effective surface element at the target satellite, the method of calculating the irradiance of the target satellite to the detector when each effective surface element is used as the light source by the Lambertian reflection principle includes:
Figure BDA00030990372500000410
Among them, the vector that the sun points to the target, that is, the incident light vector is λ S ; the vector that the target points to the detector, that is, the observation vector is λ F ; the normal vector of a certain panel of the target satellite is n; the angle between the incident vector and the normal vector of the panel is To the incident zenith angle, denoted as θ S , the angle between the observation vector and the normal vector of the surface element is the observation zenith angle, denoted as θ F ; Judging according to the relationship between the incident zenith angle θ S and the observation zenith angle θ F Whether the detector can receive the reflected light from the target satellite;

目标卫星上第i个面元所接收的单个地表面元的辐射照度为：Qls(i)＝Qls·cos(θS(i))；其中，Qls表示光源在目标卫星处的辐射照度为，θS(i)表示第i个地表面元的光线入射天顶角，Qls(i)为目标卫星上第i个面元所接收的单个地表面元的辐射照度；The irradiance of a single ground surface element received by the i-th surface element on the target satellite is: Q ls (i)=Q ls ·cos(θ S (i)); where Q ls represents the radiation of the light source at the target satellite The illuminance is, θ S (i) represents the incident zenith angle of the i-th surface element, and Q ls (i) is the irradiance of a single surface element received by the i-th surface element on the target satellite;

根据辐射照度传递公式，目标卫星第i面元对探测器的辐射照度为：

Figure BDA0003099037250000051
其中，αi为第i面元的反射率，S(i)为第i面元的面积，Rsf为目标卫星距探测器的距离，Qsf(i)为目标卫星第i面元对探测器的辐射照度，θF(i)为观测天顶角；目标卫星面元遮挡判别系数η为：According to the irradiance transfer formula, the irradiance of the i-th surface element of the target satellite to the detector is:
Figure BDA0003099037250000051
Among them, α i is the reflectivity of the i-th surface element, S(i) is the area of the i-th surface element, R sf is the distance between the target satellite and the detector, and Q sf (i) is the detection of the i-th surface element of the target satellite. The irradiance of the sensor, θ F (i) is the observation zenith angle; the target satellite panel occlusion discriminant coefficient η is:

Figure BDA0003099037250000052
其中，θS为入射天顶角、θF为观测天顶角。
Figure BDA0003099037250000052
Among them, θ S is the incident zenith angle, and θ F is the observation zenith angle.

以单个地球有效面元为唯一光源时目标卫星对探测器的辐射照度为：

Figure BDA0003099037250000053
其中，Qsf(i)为目标卫星第i面元对探测器的辐射照度，η(i)为目标卫星第i面元遮挡判别系数，n为目标卫星总面元数，Qsf为目标卫星对探测器的总辐射照度。When a single effective surface element of the earth is the only light source, the irradiance of the target satellite to the detector is:
Figure BDA0003099037250000053
Among them, Q sf (i) is the irradiance of the ith panel of the target satellite to the detector, η(i) is the occlusion discrimination coefficient of the ith panel of the target satellite, n is the total number of panels of the target satellite, and Q sf is the target satellite The total irradiance to the detector.

所述步骤三中，将目标卫星对探测器的辐射照度相加得到以地球为光源时，目标卫星对探测器的辐射照度，包括：根据以单个有效面元为唯一光源时目标卫星对探测器的辐射照度，针对所有有效面元进行叠加可以得到基于整个地球反射光的目标卫星辐射照度为

Figure BDA0003099037250000054
S有效为地表有效面元区域。In the third step, adding the irradiance of the target satellite to the detector to obtain the irradiance of the target satellite to the detector when the earth is the light source, including: according to the single effective surface element as the only light source, the target satellite to the detector. The irradiance of the target satellite based on the reflected light of the whole earth can be obtained by superimposing all the effective surface elements as
Figure BDA0003099037250000054
S is effectively the surface effective bin area.

所述步骤四中，校正观测数据的校正方法包括：

Figure BDA0003099037250000055
其中，从目标卫星处获得的实时的目标卫星对探测器的辐射照度为QA，校正观测数据为QJ，以地球为光源时，目标卫星对探测器的辐射照度为
Figure BDA0003099037250000061
In the step 4, the correction method for correcting the observed data includes:
Figure BDA0003099037250000055
Among them, the real-time irradiance of the target satellite to the detector obtained from the target satellite is Q A , the corrected observation data is Q J , and when the earth is used as the light source, the irradiance of the target satellite to the detector is
Figure BDA0003099037250000061

所述步骤五中，目标卫星的姿态估计方法包括：状态向量

Figure BDA0003099037250000062
其中，δp表示误差罗德里格参数；ω为目标卫星翻滚角速度，目标姿态指目标本体系相对地球惯性坐标系的姿态；状态向量经过无味变换得到2n+1个sigma点，其表达式为：In the step 5, the attitude estimation method of the target satellite includes: a state vector
Figure BDA0003099037250000062
Among them, δp represents the error Rodrigue parameter; ω is the roll angular velocity of the target satellite, and the target attitude refers to the attitude of the target system relative to the earth's inertial coordinate system; the state vector is tastelessly transformed to obtain 2n+1 sigma points, and its expression is:

Figure BDA0003099037250000063
Figure BDA0003099037250000063

Figure BDA0003099037250000064
Figure BDA0003099037250000064

其中，λ＝α2(n+κ)-n，为缩放比例系数，通常有κ＝3-n；

Figure BDA0003099037250000065
为当前状态量的第i个sigma点，σ(i)为第i个sigma点；将第k时刻的滤波值
Figure BDA0003099037250000066
作为当前的均值sigma点，记为
Figure BDA0003099037250000067
Figure BDA0003099037250000068
为k时刻状态协方差矩阵的滤波值；用Pi表示矩阵P的第i列；chol(P)表示对矩阵P进行cholesky分解，其对应权值为：Among them, λ=α 2 (n+κ)-n, is the scaling factor, usually κ=3-n;
Figure BDA0003099037250000065
is the i-th sigma point of the current state quantity, and σ(i) is the i-th sigma point;
Figure BDA0003099037250000066
As the current mean sigma point, denoted as
Figure BDA0003099037250000067
Figure BDA0003099037250000068
is the filter value of the state covariance matrix at time k; P i represents the i-th column of the matrix P; chol(P) represents the cholesky decomposition of the matrix P, and its corresponding weight is:

Figure BDA0003099037250000069
Figure BDA0003099037250000069

其中，n为状态向量的维数，Wmean表示均值，Wcov表示协方差，α为控制sigma点扩散程度的系数，取值范围通常为[10-4,1]；β为权值调整系数，通常取值为2；λ＝α2(n+κ)-n，为缩放比例系数，通常有κ＝3-n；获取误差四元数的sigma，Among them, n is the dimension of the state vector, W mean is the mean value, W cov is the covariance, α is the coefficient that controls the spread of the sigma point, the value range is usually [10 -4 ,1]; β is the weight adjustment coefficient , usually the value is 2; λ=α 2 (n+κ)-n, is the scaling factor, usually κ=3-n; to obtain the sigma of the error quaternion,

Figure BDA0003099037250000071
Figure BDA0003099037250000071

Figure BDA0003099037250000072
Figure BDA0003099037250000072

其中，a＝1,f＝2(a+1)，

Figure BDA0003099037250000073
表示四元数的前三个分量
Figure BDA0003099037250000074
表示四元数的第四个分量
Figure BDA0003099037250000075
表示四元数误差分量；再进行误差四元数向四元数的转换，
Figure BDA0003099037250000076
Among them, a=1, f=2(a+1),
Figure BDA0003099037250000073
represents the first three components of a quaternion
Figure BDA0003099037250000074
represents the fourth component of the quaternion
Figure BDA0003099037250000075
represents the quaternion error component; then convert the error quaternion to quaternion,
Figure BDA0003099037250000076

其中，

Figure BDA0003099037250000077
表示当前四元数sigma点的均值，
Figure BDA0003099037250000078
表示第i个sigma点的误差四元数根据四元数乘法转换得到的四元数值，
Figure BDA0003099037250000079
表示第i个sigma点的误差四元数；四元数乘法定义如下，in,
Figure BDA0003099037250000077
represents the mean of the current quaternion sigma point,
Figure BDA0003099037250000078
Represents the quaternion value obtained by converting the error quaternion of the i-th sigma point according to quaternion multiplication,
Figure BDA0003099037250000079
represents the error quaternion of the i-th sigma point; the quaternion multiplication is defined as follows,

Figure BDA00030990372500000710
其中，
Figure BDA00030990372500000711
对于三维向量a，[a×]表示向量叉积矩阵，
Figure BDA00030990372500000712
qa和qb表示任意两个四元数；四元数sigma点和角速度sigma点传递，
Figure BDA00030990372500000713
其中，
Figure BDA00030990372500000710
in,
Figure BDA00030990372500000711
For a three-dimensional vector a, [a×] represents the vector cross product matrix,
Figure BDA00030990372500000712
q a and q b represent any two quaternions; the quaternion sigma point and the angular velocity sigma point transfer,
Figure BDA00030990372500000713
in,

Figure BDA00030990372500000714
Figure BDA00030990372500000714

Δt为探测器采样时间间隔；Γ(ωk)由动力学方程在时域离散化得到；基于以上传递后的四元数sigma点，再进行误差罗德里格参数点的转换，

Figure BDA00030990372500000715
其中，共轭四元数q-1定义为
Figure BDA0003099037250000081
继而，
Figure BDA0003099037250000082
其中，a＝1,f＝2(a+1)。令其均值sigma点
Figure BDA0003099037250000083
得到传递后状态向量的sigma点为
Figure BDA0003099037250000084
最后通过加权求和的方式进行状态更新和协方差更新，
Figure BDA0003099037250000085
式中，Qk+1为离散时间的过程噪声协方差；由观测方程可计算得到每个sigma点的观测预测值，
Figure BDA0003099037250000086
h表示观测方程；
Figure BDA0003099037250000087
为sigma点对应的观测预测值；再通过加权求和得到观测预测均值以及量测协方差和互协方差：
Figure BDA0003099037250000088
Δt is the sampling time interval of the detector; Γ(ω k ) is obtained by discretizing the dynamic equation in the time domain; based on the quaternion sigma point after the above transfer, the error Rodrigue parameter point is converted,
Figure BDA00030990372500000715
where the conjugate quaternion q -1 is defined as
Figure BDA0003099037250000081
Then,
Figure BDA0003099037250000082
Among them, a=1, f=2(a+1). Let its mean sigma point
Figure BDA0003099037250000083
The sigma point of the state vector after the transfer is obtained as
Figure BDA0003099037250000084
Finally, the state update and covariance update are performed by means of weighted summation,
Figure BDA0003099037250000085
In the formula, Q k+1 is the process noise covariance in discrete time; the observed predicted value of each sigma point can be calculated from the observation equation,
Figure BDA0003099037250000086
h represents the observation equation;
Figure BDA0003099037250000087
is the observed predicted value corresponding to the sigma point; then the mean of the observed predicted and the measured covariance and cross-covariance are obtained by weighted summation:
Figure BDA0003099037250000088

Figure BDA0003099037250000089
Rk+1表示
Figure BDA00030990372500000810
表示量测协方差，
Figure BDA00030990372500000811
表示互协方差，
Figure BDA00030990372500000812
表示第k+1步的观测预测值；卡尔曼增益计算，
Figure BDA00030990372500000813
Kk+1表示第k+1次滤波的卡尔曼增益；状态更新和协方差更新，
Figure BDA00030990372500000814
式中，
Figure BDA00030990372500000815
为第k+1时刻的设备观测值，即校正后的目标卫星观测数据QJ；将状态滤波值中的误差罗德里格参数转化为四元数，并在下一次滤波开始前进行误差罗德里格参数置零，即
Figure BDA00030990372500000816
Figure BDA00030990372500000817
表示第k+1次滤波时的状态向量，03×1表示3行1列的零矩阵。
Figure BDA0003099037250000089
R k+1 means
Figure BDA00030990372500000810
represents the measurement covariance,
Figure BDA00030990372500000811
represents the cross-covariance,
Figure BDA00030990372500000812
Indicates the observed predicted value of the k+1 step; Kalman gain calculation,
Figure BDA00030990372500000813
K k+1 represents the Kalman gain of the k+1th filtering; state update and covariance update,
Figure BDA00030990372500000814
In the formula,
Figure BDA00030990372500000815
is the equipment observation value at the k+1th moment, that is, the corrected target satellite observation data Q J ; the error Rodrigue parameter in the state filtering value is converted into a quaternion, and the error Rodrigue is carried out before the next filtering starts. The parameter is set to zero, i.e.
Figure BDA00030990372500000816
Figure BDA00030990372500000817
Represents the state vector during the k+1th filtering, and 0 3×1 represents a zero matrix with 3 rows and 1 column.

本发明具有如下有益效果：The present invention has the following beneficial effects:

1、本发明通过对地球模型进行网格化处理，精确计算了地球反射光对卫星的照度，以及探测器接收到的卫星总照度，模型详细考虑了地球的体积及卫星的轨道高度因素对仿真结果的影响。1. The present invention accurately calculates the illuminance of the reflected light from the earth to the satellite and the total illuminance of the satellite received by the detector by performing grid processing on the earth model. The model considers the volume of the earth and the orbital height of the satellite in detail. impact on results.

2、仿真了一年内由地球反射光引起的卫星亮度变化情况，并做了以下两种情况的考虑：①因月球遮挡导致太阳光无法到达地球的情况。②卫星处于地球阴影区，无法受到地球反射光照射的情况。③卫星自身面元之间的遮挡关系。2. Simulate the change of satellite brightness caused by the reflected light of the earth in one year, and consider the following two situations: ① The situation that sunlight cannot reach the earth due to the occlusion of the moon. ② The satellite is in the shadow area of the earth and cannot be irradiated by the reflected light of the earth. ③ The occlusion relationship between the satellites themselves.

3、基于所建立的地球反射光模型，可以对获得的卫星观测数据进行校正，将校正数据输入到无味卡尔曼滤波器中，可以大大提高卫星姿态的估计精度。3. Based on the established earth reflected light model, the obtained satellite observation data can be corrected, and the correction data can be input into the tasteless Kalman filter, which can greatly improve the estimation accuracy of satellite attitude.

4、另一方面，给出了不同时刻地球反射光对卫星观测影响的占比，如图为了消减影响可以选择在影响较弱的时刻进行观测，对卫星观测者们的研究具有一定的指导意义。4. On the other hand, the proportion of the impact of the Earth's reflected light on satellite observations at different times is given, as shown in the figure, in order to reduce the impact, you can choose to observe at a time when the impact is weak, which has certain guiding significance for the research of satellite observers. .

附图说明Description of drawings

图1是本发明的一种地球反射光校正的卫星姿态估计方法流程示意图。FIG. 1 is a schematic flowchart of a method for estimating satellite attitude based on earth reflected light correction according to the present invention.

图2是地球表面按照经纬网格划分为多个面元的示意图。FIG. 2 is a schematic diagram showing that the earth surface is divided into a plurality of surface elements according to the latitude and longitude grid.

图3是地表面元(i,j)长和宽计算示意图。Figure 3 is a schematic diagram of the calculation of the length and width of the ground surface element (i, j).

图4是地表面元反射示意图。Figure 4 is a schematic diagram of ground surface element reflection.

图5是探测器接收卫星反射光示意图。Figure 5 is a schematic diagram of the detector receiving the reflected light from the satellite.

图6-1是地表有效面元正视图。Figure 6-1 is a front view of the effective surface element on the surface.

图6-2是地表有效面元俯视图。Figure 6-2 is a top view of the effective surface element on the surface.

图7是一年内地球反射光照射的GEO卫星对探测器的照度。Figure 7 shows the illuminance of the GEO satellite to the detector illuminated by the Earth's reflected light in one year.

图8是一年内不同时刻地球反射光对GEO卫星观测结果的影响占比。Figure 8 shows the proportion of the Earth's reflected light at different times of the year on the GEO satellite observations.

图9-1是未校正的卫星欧拉角θ1滤波估计图。Figure 9-1 is a filtered estimate of the uncorrected satellite Euler angle θ 1 .

图9-2是未校正的卫星欧拉角θ2滤波估计图。Figure 9-2 is a plot of the uncorrected satellite Euler angle θ 2 filter estimate.

图9-3是未校正的卫星欧拉角θ3滤波估计图。Figure 9-3 is a plot of the uncorrected satellite Euler angle θ 3 filter estimate.

图10-1是校正后的卫星欧拉角θ1滤波估计图。Figure 10-1 is the corrected satellite Euler angle θ 1 filter estimate.

图10-2是校正后的卫星欧拉角θ2滤波估计图。Figure 10-2 is the corrected satellite Euler angle θ 2 filter estimate.

图10-3是校正后的卫星欧拉角θ3滤波估计图。Figure 10-3 is the corrected satellite Euler angle θ 3 filter estimate.

具体实施方式Detailed ways

本发明实施例中，假设探测器位于地表(75.5966°E,30.0386°N,0m)，以一年为观测时间(1 Jan 2021 00:00:00.000UTCG～1 Jan 2022 00:00:00.000UTCG)，60s为观测间隔，对GEO卫星进行观测。设置卫星尺寸为5×5×5m3，姿态为失稳翻滚，初始角速度[0 -0.10]°/s，初始状态四元数设为[-0.646924 0.401929 -0.333014 0.555917]。观测噪声设为标准差为0.05m2的零均值白噪声。卫星各面的反射率均为0.3。如图1所示，本发明公开了一种地球反射光校正的卫星姿态估计方法，包括以下步骤：In this embodiment of the present invention, it is assumed that the detector is located on the surface (75.5966°E, 30.0386°N, 0m), and one year is used as the observation time (1 Jan 2021 00:00:00.000UTCG～1 Jan 2022 00:00:00.000UTCG) , 60s is the observation interval, and the GEO satellite is observed. The satellite size is set to 5×5×5m 3 , the attitude is unstable roll, the initial angular velocity is [0 -0.10]°/s, and the initial state quaternion is set to [-0.646924 0.401929 -0.333014 0.555917]. The observation noise was set as zero-mean white noise with a standard deviation of 0.05 m 2 . The reflectivity of all sides of the satellite is 0.3. As shown in FIG. 1, the present invention discloses a satellite attitude estimation method for earth reflection light correction, which includes the following steps:

步骤一：将地球表面按照经纬网格划分为多个面元，计算面元面积；基于太阳、地球和目标卫星之间的遮挡关系筛选有效面元；Step 1: Divide the earth's surface into multiple surfels according to the latitude and longitude grid, and calculate the surfel area; screen the effective surfels based on the occlusion relationship between the sun, the earth and the target satellite;

如图2所示，将地球表面按照经纬网格划分为网格面元，面元中心位于经纬点(i,j)(i,j∈N)，长和宽分别为一个纬度距离和一个经度距离。As shown in Figure 2, the earth's surface is divided into grid panels according to the latitude and longitude grid, the center of the panel is located at the latitude and longitude point (i, j) (i, j∈N), and the length and width are a latitude distance and a longitude respectively. distance.

如图3所示，地球的半径为RE，每一个面元可以视为中心位于经纬点(i,j)(i,j∈N)，长和宽分别为一个纬度距离和一个经度距离的长方形平面，每个地表面元的长均为：lj＝2πRE/360；纬度i处的地表面元，宽为：di＝2πREcosi/360；地表面元面积为：sij＝ljdi；其中，lj为经度j处地表面元的长，di为纬度i处的地表面元的宽，sij为面元面积，RE为地球半径，i、j分别表示纬度和经度。As shown in Figure 3, the radius of the earth is RE, and each surface element can be regarded as the center is located at the latitude and longitude point ( i ,j) (i,j∈N), and the length and width are respectively a latitude distance and a longitude distance. For a rectangular plane, the length of each ground element is: l j =2πR E /360; the width of the ground element at latitude i is: d i =2πR E cosi/360; the area of the ground element is: s ij = l j d i ; where, l j is the length of the ground surface element at longitude j, d i is the width of the ground surface element at latitude i, s ij is the area of the surface element, RE is the radius of the earth, i and j represent the Latitude and longitude.

如图4所示，筛选有效面元筛选方法包括：太阳光对地表面元的入射天顶角为θin，反射天顶角为θout，筛选有效面元的方法包括：As shown in Figure 4, the screening method for screening effective bins includes: the incident zenith angle of sunlight to the ground surface element is θ in , the reflection zenith angle is θ out , and the method for screening effective bins includes:

Figure BDA0003099037250000111
Figure BDA0003099037250000111

其中，θin为太阳光对地表面元的入射天顶角，θout为反射天顶角。Among them, θ in is the incident zenith angle of sunlight to the ground surface element, and θ out is the reflection zenith angle.

步骤二：根据辐射照度传递规律，即辐射照度与距离的平方成反比定律，计算太阳可见光波段在各有效面元处的辐射照度，并根据朗伯反射原理计算有效面元在目标卫星处产生的辐射照度，包括：根据辐射照度与距离的平方成反比定律，在太阳光不被遮挡的情况下，得到太阳可见光波段在各有效面元处的辐射照度：

Figure BDA0003099037250000112
其中，R0为一年内日地的平均距离；Q0是在距离太阳R0处，可见光波段(0.4～0.7μm)部分的辐射照度，R为地表面元到太阳的距离，Q是在距离太阳R的有效面元处，对可见光波段(0.4～0.7μm)的辐射照度。Step 2: According to the irradiance transfer law, that is, the irradiance is inversely proportional to the square of the distance, calculate the irradiance of the solar visible light band at each effective surface element, and calculate the effective surface element generated at the target satellite according to the Lambertian reflection principle. The irradiance includes: according to the law of inverse proportionality between the irradiance and the square of the distance, when the sunlight is not blocked, the irradiance of the visible light band of the sun at each effective surface element is obtained:
Figure BDA0003099037250000112
Among them, R 0 is the average distance between the sun and the earth in one year; Q 0 is the irradiance in the visible light band (0.4-0.7 μm) at the distance from the sun at R 0 , R is the distance from the surface element to the sun, and Q is the distance between the sun and the sun. The irradiance of the visible light band (0.4-0.7μm) at the effective surface element of the sun R.

例如：假设一年内日地的平均距离R0为1.495×108km，则在距离R0处，对可见光波段(0.4～0.7μm)积分的辐射照度Q0为438.894W/m2。根据太阳辐射照度与距离的平方成反比，在太阳光不被遮挡的情况下，任意位置处太阳可见光波段的辐射照度可以表示为

Figure BDA0003099037250000113
For example: Assuming that the average distance R 0 between the sun and the earth in one year is 1.495×10 8 km, then at the distance R 0 , the integrated irradiance Q 0 in the visible light band (0.4-0.7μm) is 438.894W/m 2 . According to the inverse ratio of the solar irradiance to the square of the distance, in the case where the sunlight is not blocked, the irradiance in the visible light band of the sun at any position can be expressed as
Figure BDA0003099037250000113

步骤二中，计算有效面元在目标卫星处产生的辐射照度，得出卫星所受单个地表面元的辐射照度的方法包括：地球表面太阳辐射照度为：

Figure BDA0003099037250000114
其中，R0为一年内日地的平均距离，Rse为太阳到地球的距离，Qse为太阳在地球表面的辐射照度，Q0为在距离太阳R0处，对可见光波段的辐射照度；假设地球表面及地表大气总的平均反照率值为
Figure BDA0003099037250000115
可以计算出地表面元对各方向产生的辐射亮度为：
Figure BDA0003099037250000116
其中，θin为太阳光对地表面元的入射天顶角，Qsecosθin为地表面元所接收的太阳辐射照度；单个地表面元(i,j)在
Figure BDA0003099037250000121
方向目标卫星处产生的辐射照度为：
Figure BDA0003099037250000122
其中，
Figure BDA0003099037250000123
为面元指向目标卫星的方向向量，
Figure BDA0003099037250000124
为地表面元相对目标卫星所呈立体角，sij为经纬度位于(i,j)的地表面元面积，
Figure BDA0003099037250000125
为地表面元至目标卫星的距离，Qij为单个地表面元(i,j)在
Figure BDA0003099037250000126
方向目标卫星处产生的辐射照度。In step 2, the irradiance generated by the effective surface element at the target satellite is calculated, and the method for obtaining the irradiance of a single ground surface element received by the satellite includes: the solar irradiance on the surface of the earth is:
Figure BDA0003099037250000114
Among them, R 0 is the average distance between the sun and the earth in one year, R se is the distance from the sun to the earth, Q se is the irradiance of the sun on the surface of the earth, and Q 0 is the irradiance of the visible light band at the distance R 0 from the sun; Assume that the total average albedo of the earth's surface and the surface atmosphere is
Figure BDA0003099037250000115
It can be calculated that the radiance generated by the surface element in each direction is:
Figure BDA0003099037250000116
Among them, θ in is the incident zenith angle of sunlight to the ground surface element, Q se cosθ in is the solar radiation illuminance received by the ground surface element; a single ground surface element (i, j) is in
Figure BDA0003099037250000121
The irradiance generated at the target satellite is:
Figure BDA0003099037250000122
in,
Figure BDA0003099037250000123
is the direction vector of the surfel pointing to the target satellite,
Figure BDA0003099037250000124
is the solid angle of the surface element relative to the target satellite, s ij is the area of the surface element located at (i, j) in latitude and longitude,
Figure BDA0003099037250000125
is the distance from the surface element to the target satellite, and Q ij is the single surface element (i, j) in
Figure BDA0003099037250000126
The irradiance generated at the target satellite.

步骤三：根据有效面元在目标卫星处产生的辐射照度，通过朗伯反射原理计算出以各个有效面元为光源时，目标卫星对探测器的辐射照度；将目标卫星对探测器的辐射照度相加得到以地球为光源时，目标卫星对探测器的辐射照度；其中，根据有效面元在目标卫星处产生的辐射照度，通过朗伯反射原理计算出以各个有效面元为光源时，卫星对探测器的辐射照度的方法包括：卫星反射模型如图5所示，以边长5m的立方星为例，设太阳指向目标的向量即入射光向量为λS；设目标指向探测器的向量即观测向量为λF；卫星某面元的法向量为n；定义入射向量与面元法向量的夹角为到入射天顶角，记作θS，定义观测向量与面元法向量的夹角为观测天顶角，记作θF。表达式为：Step 3: According to the irradiance generated by the effective surface element at the target satellite, the irradiance of the target satellite to the detector is calculated by the Lambertian reflection principle when each effective surface element is used as the light source; the irradiance of the target satellite to the detector is calculated. Add up to get the irradiance of the target satellite to the detector when the earth is used as the light source; among them, according to the irradiance generated by the effective surface element at the target satellite, the Lambertian reflection principle is used to calculate when each effective surface element is used as the light source, the satellite The method for the irradiance of the detector includes: the satellite reflection model is shown in Figure 5, taking a cube star with a side length of 5m as an example, set the vector of the sun pointing to the target, that is, the incident light vector as λ S ; set the vector of the target pointing to the detector That is, the observation vector is λ F ; the normal vector of a certain panel of the satellite is n; the angle between the incident vector and the normal vector of the panel is defined as the angle to the incident zenith, denoted as θ S , and the angle between the observation vector and the normal vector of the panel is defined. The angle is the observed zenith angle, denoted as θ F . The expression is:

Figure BDA0003099037250000127
根据入射天顶角θS、观测天顶角θF的大小关系判断探测器是否能收到卫星反射光。
Figure BDA0003099037250000127
According to the relationship between the incident zenith angle θ S and the observation zenith angle θ F , it is judged whether the detector can receive the reflected light from the satellite.

目标卫星上第i个面元所接收的单个地表面元的辐射照度为：The irradiance of a single ground surface element received by the i-th surface element on the target satellite is:

Qls(i)＝Qls·cos(θS(i))；其中，Qls表示光源在卫星处的辐射照度为，θS(i)表示第i个地表面元的光线入射天顶角。根据辐射照度传递公式，卫星第i面元对探测器的辐射照度为：

Figure BDA0003099037250000128
Q ls (i)=Q ls ·cos(θ S (i)); wherein, Q ls represents the irradiance of the light source at the satellite, and θ S (i) represents the incident zenith angle of the i-th surface element . According to the irradiance transfer formula, the irradiance of the i-th surface element of the satellite to the detector is:
Figure BDA0003099037250000128

其中，αi为第i面元的反射率，S(i)为第i面元的面积，Rsf为卫星距探测器的距离，Qsf(i)为卫星第i面元对探测器的辐射照度，θF(i)为观测天顶角。卫星面元遮挡判别系数η为：

Figure BDA0003099037250000131
Among them, α i is the reflectivity of the i-th surface element, S(i) is the area of the i-th surface element, R sf is the distance from the satellite to the detector, and Q sf (i) is the ith surface element of the satellite to the detector. Irradiance, θ F (i) is the observed zenith angle. The satellite occlusion discrimination coefficient η is:
Figure BDA0003099037250000131

其中，θS为入射天顶角、θF为观测天顶角。Among them, θ S is the incident zenith angle, and θ F is the observation zenith angle.

以单个地球有效面元为唯一光源时卫星对探测器的辐射照度为：When a single effective surface element of the earth is the only light source, the irradiance of the satellite to the detector is:

Figure BDA0003099037250000132
其中，Qsf(i)为卫星第i面元对探测器的辐射照度，η(i)为卫星第i面元遮挡判别系数，n为卫星总面元数，Qsf为卫星对探测器的总辐射照度。
Figure BDA0003099037250000132
Among them, Q sf (i) is the irradiance of the i-th panel of the satellite to the detector, η(i) is the occlusion discrimination coefficient of the i-th panel of the satellite, n is the total number of satellite panels, and Q sf is the irradiance of the satellite to the detector. total irradiance.

步骤三中，将所求得的各辐射照度进行叠加得到以整个地球为光源时，卫星对探测器的辐射照度，包括：根据以单个有效面元为唯一光源时卫星对探测器的辐射照度，针对所有有效面元进行叠加可以得到基于整个地球反射光的卫星辐射照度为

Figure BDA0003099037250000133
S有效为地表有效面元区域，有效区域如图6-1和图6-2所示。In step 3, the obtained irradiances are superimposed to obtain the irradiance of the satellite to the detector when the whole earth is used as the light source, including: according to the irradiance of the satellite to the detector when a single effective surface element is the only light source, By superimposing all effective bins, the satellite irradiance based on the reflected light of the whole earth can be obtained as
Figure BDA0003099037250000133
S is effectively the effective surface area of the surface, and the effective area is shown in Figure 6-1 and Figure 6-2.

步骤四：从目标卫星处获得实时的目标卫星对探测器的辐射照度减去以地球为光源时，目标卫星对探测器的辐射照度，得到了仅以太阳为光源时的目标卫星对探测器的辐射照度的校正观测数据；Step 4: Obtain the real-time irradiance of the target satellite to the detector from the target satellite minus the irradiance of the target satellite to the detector when the earth is used as the light source, and obtain the irradiance of the target satellite to the detector when only the sun is used as the light source. Corrected observation data of irradiance;

根据步骤一、二、三的过程，以一年为仿真时间，计算出以地球为唯一光源时的GEO轨道卫星对探测器的辐射照度。以GEO轨道卫星为研究对象时的仿真结果如图7所示，不同地球反射光对卫星观测结果的影响占比如图8所示。According to the process of steps 1, 2 and 3, and taking one year as the simulation time, calculate the irradiance of the GEO orbiting satellite to the detector when the earth is the only light source. The simulation results when GEO orbiting satellites are used as the research object are shown in Figure 7, and the proportion of the influence of different Earth reflected light on the satellite observation results is shown in Figure 8.

步骤四中，仅以太阳为光源时的目标卫星对探测器的辐射照度的校正观测数据的校正方法包括：

Figure BDA0003099037250000134
其中，从目标卫星处获得的实时的目标卫星对探测器的辐射照度为QA，校正观测数据为QJ，以地球为光源时，目标卫星对探测器的辐射照度为
Figure BDA0003099037250000135
In step 4, the correction method for correcting the observation data of the irradiance of the detector by the target satellite only when the sun is used as the light source includes:
Figure BDA0003099037250000134
Among them, the real-time irradiance of the target satellite to the detector obtained from the target satellite is Q A , the corrected observation data is Q J , and when the earth is used as the light source, the irradiance of the target satellite to the detector is
Figure BDA0003099037250000135

步骤五：基于校正观测数据，使用无味卡尔曼滤波的方法对目标卫星进行姿态估计，得到地球反射光校正数据。例如：选取观测站当地时间1月1日傍晚16时40分至18时40分的观测数据用于姿态估计。此时计算得到的地球反射光影响占比为0.1～0.2。为了避免“三分量”参数计算过程中可能出现的奇异值问题，本发明选用运算相对线性的四元数作为描述参量。四元数定义为q＝[εT q4]T，εT为三维矢量，q4为标量，并满足约束qTq＝1。四元数表示的姿态运动方程与角速度动力学方程为：Step 5: Based on the corrected observation data, use the method of tasteless Kalman filtering to estimate the attitude of the target satellite, and obtain the correction data of the earth's reflected light. For example, the observation data of the observatory from 16:40 to 18:40 in the evening of January 1, local time, is selected for attitude estimation. At this time, the calculated influence ratio of the reflected light from the earth is 0.1 to 0.2. In order to avoid the singular value problem that may occur in the "three-component" parameter calculation process, the present invention selects the quaternion with relatively linear operation as the description parameter. A quaternion is defined as q=[ε T q 4 ] T , where ε T is a three-dimensional vector, q 4 is a scalar, and satisfies the constraint q T q=1. The attitude motion equation and angular velocity dynamic equation represented by quaternion are:

Figure BDA0003099037250000141
式中，
Figure BDA0003099037250000142
E3×3表示3×3的单位矩阵，ω(t)＝[ωx ωy ωz]T表示本体系下目标的角速度，Π(t)为目标所受外部力矩和自身控制力拒之和，J为目标的转动惯量。w(t)和v(t)分别为过程噪声和观测噪声，以均值为0的高斯白噪声代替。对于三维向量a，[a×]表示向量叉积矩阵，
Figure BDA0003099037250000143
Figure BDA0003099037250000141
In the formula,
Figure BDA0003099037250000142
E 3×3 represents the unit matrix of 3×3, ω(t)=[ω x ω y ω z ] T represents the angular velocity of the target under this system, Π(t) is the external torque and self-control force rejected by the target and, J is the moment of inertia of the target. w(t) and v(t) are process noise and observation noise, respectively, and are replaced by Gaussian white noise with mean 0. For a three-dimensional vector a, [a×] represents the vector cross product matrix,
Figure BDA0003099037250000143

四元数姿态转移矩阵为：A(q)＝Ξ(q)TΨ(q)；式中，

Figure BDA0003099037250000144
The quaternion attitude transfer matrix is: A(q)=Ξ(q) T Ψ(q); in the formula,
Figure BDA0003099037250000144

观测方程建立为：The observation equation is established as:

Figure BDA0003099037250000145
式中，k为观测时刻，Q0为太阳常数，Rsf为卫星离探测器的距离，θi(k)为卫星第i个面元在第k个时刻的入射天顶角，θr(k)为第i个面元在第k个时刻的观测天顶角，η(k)为k时刻的面元遮挡判别系数，S(i)为卫星第i个面元的面积，n为卫星面元总数，这里取6。状态向量
Figure BDA0003099037250000146
其中，δp表示误差罗德里格参数；ω为卫星翻滚角速度，目标姿态指目标本体系相对地球惯性坐标系的姿态；状态向量经过无味变换得到2n+1个sigma点，其表达式为：
Figure BDA0003099037250000145
In the formula, k is the observation time, Q 0 is the solar constant, R sf is the distance from the satellite to the detector, θ i (k) is the incident zenith angle of the i-th panel of the satellite at the k-th moment, θ r ( k) is the observation zenith angle of the i-th surface element at the k-th time, η(k) is the surface occlusion discrimination coefficient at the k time, S(i) is the area of the i-th surface element of the satellite, and n is the satellite The total number of bins, here is 6. state vector
Figure BDA0003099037250000146
Among them, δp represents the error Rodrigue parameter; ω is the satellite roll angular velocity, and the target attitude refers to the attitude of the target system relative to the earth's inertial coordinate system; the state vector is tastelessly transformed to obtain 2n+1 sigma points, and its expression is:

Figure BDA0003099037250000151
Figure BDA0003099037250000151

其中，λ＝α2(n+κ)-n，为缩放比例系数，通常有κ＝3-n；

Figure BDA0003099037250000152
为当前状态量的第i个sigma点，σ(i)为第i个sigma点；将第k时刻的滤波值
Figure BDA0003099037250000153
作为当前的均值sigma点，记为
Figure BDA0003099037250000154
Figure BDA0003099037250000155
为k时刻状态协方差矩阵的滤波值；用Pi表示矩阵P的第i列；chol(P)表示对矩阵P进行cholesky分解，其对应权值为：Among them, λ=α 2 (n+κ)-n, is the scaling factor, usually κ=3-n;
Figure BDA0003099037250000152
is the i-th sigma point of the current state quantity, and σ(i) is the i-th sigma point;
Figure BDA0003099037250000153
As the current mean sigma point, denoted as
Figure BDA0003099037250000154
Figure BDA0003099037250000155
is the filter value of the state covariance matrix at time k; P i represents the i-th column of the matrix P; chol(P) represents the cholesky decomposition of the matrix P, and its corresponding weight is:

Figure BDA0003099037250000156
Figure BDA0003099037250000156

其中，n为状态向量的维数，Wmean表示均值，Wcov表示协方差，α为控制sigma点扩散程度的系数，取值范围通常为[10-4,1]；β为权值调整系数，通常取值为2；λ＝α2(n+κ)-n，为缩放比例系数，通常有κ＝3-n。获取误差四元数的sigma，Among them, n is the dimension of the state vector, W mean is the mean value, W cov is the covariance, α is the coefficient that controls the spread of the sigma point, the value range is usually [10 -4 ,1]; β is the weight adjustment coefficient , usually takes a value of 2; λ=α 2 (n+κ)-n, is a scaling factor, usually κ=3-n. Get the sigma of the error quaternion,

Figure BDA0003099037250000157
Figure BDA0003099037250000157

其中，a＝1,f＝2(a+1)，

Figure BDA0003099037250000158
表示四元数的前三个分量，
Figure BDA0003099037250000159
表示四元数的第四个分量
Figure BDA00030990372500001510
表示四元数误差分量；再进行误差四元数向四元数的转换，
Figure BDA00030990372500001511
其中，
Figure BDA00030990372500001512
表示当前四元数sigma点的均值，
Figure BDA00030990372500001513
表示第i个sigma点的误差四元数根据四元数乘法转换得到的四元数值，
Figure BDA00030990372500001514
表示第i个sigma点的误差四元数；Among them, a=1, f=2(a+1),
Figure BDA0003099037250000158
represents the first three components of the quaternion,
Figure BDA0003099037250000159
represents the fourth component of the quaternion
Figure BDA00030990372500001510
represents the quaternion error component; then convert the error quaternion to quaternion,
Figure BDA00030990372500001511
in,
Figure BDA00030990372500001512
represents the mean of the current quaternion sigma point,
Figure BDA00030990372500001513
Represents the quaternion value obtained by converting the error quaternion of the i-th sigma point according to quaternion multiplication,
Figure BDA00030990372500001514
represents the error quaternion of the i-th sigma point;

四元数乘法可定义如下，

Figure BDA0003099037250000161
其中，
Figure BDA0003099037250000162
对于三维向量a，[a×]表示向量叉积矩阵，
Figure BDA0003099037250000163
qa和qb表示任意两个四元数。Quaternion multiplication can be defined as follows,
Figure BDA0003099037250000161
in,
Figure BDA0003099037250000162
For a three-dimensional vector a, [a×] represents the vector cross product matrix,
Figure BDA0003099037250000163
q a and q b represent any two quaternions.

四元数sigma点和角速度sigma点传递，Quaternion sigma point and angular velocity sigma point transfer,

Figure BDA0003099037250000164
其中，
Figure BDA0003099037250000164
in,

Figure BDA0003099037250000165
Figure BDA0003099037250000165

Δt为探测器采样时间间隔；Γ(ωk)由动力学方程在时域离散化得到；Δt is the sampling time interval of the detector; Γ(ω k ) is obtained by discretizing the dynamic equation in the time domain;

基于以上传递后的四元数sigma点，再进行误差罗德里格参数点的转换，

Figure BDA0003099037250000166
其中，共轭四元数q-1定义为
Figure BDA0003099037250000167
Based on the quaternion sigma points passed above, the conversion of the error Rodrigue parameter points is performed,
Figure BDA0003099037250000166
where the conjugate quaternion q -1 is defined as
Figure BDA0003099037250000167

继而，

Figure BDA0003099037250000168
其中，a＝1,f＝2(a+1)。令其均值sigma点
Figure BDA0003099037250000169
得到传递后状态向量的sigma点为
Figure BDA00030990372500001610
Then,
Figure BDA0003099037250000168
Among them, a=1, f=2(a+1). Let its mean sigma point
Figure BDA0003099037250000169
The sigma point of the state vector after the transfer is obtained as
Figure BDA00030990372500001610

最后通过加权求和的方式进行状态更新和协方差更新，

Figure BDA00030990372500001611
式中，Qk+1为离散时间的过程噪声协方差。由观测方程可计算得到每个sigma点的观测预测值，
Figure BDA0003099037250000171
h表示观测方程。Finally, the state update and covariance update are performed by means of weighted summation,
Figure BDA00030990372500001611
where Q k+1 is the discrete-time process noise covariance. The observed predicted value of each sigma point can be calculated from the observation equation,
Figure BDA0003099037250000171
h represents the observation equation.

Figure BDA0003099037250000172
为sigma点对应的观测预测值；再通过加权求和得到观测预测均值以及量测协方差和互协方差：
Figure BDA0003099037250000173
Figure BDA0003099037250000172
is the observed predicted value corresponding to the sigma point; then the mean of the observed predicted and the measured covariance and cross-covariance are obtained by weighted summation:
Figure BDA0003099037250000173

Figure BDA0003099037250000174
Figure BDA0003099037250000174

Rk+1表示

Figure BDA0003099037250000175
表示量测协方差，
Figure BDA0003099037250000176
表示互协方差，
Figure BDA0003099037250000177
表示第k+1步的观测预测值。R k+1 means
Figure BDA0003099037250000175
represents the measurement covariance,
Figure BDA0003099037250000176
represents the cross-covariance,
Figure BDA0003099037250000177
Represents the observed predicted value at step k+1.

卡尔曼增益计算，

Figure BDA0003099037250000178
Kk+1表示第k+1次滤波的卡尔曼增益。状态更新和协方差更新，
Figure BDA0003099037250000179
Kalman gain calculation,
Figure BDA0003099037250000178
K k+1 represents the Kalman gain of the k+1th filtering. state updates and covariance updates,
Figure BDA0003099037250000179

式中，

Figure BDA00030990372500001710
为第k+1时刻的设备观测值，即校正后的卫星观测数据QJ。最后，将状态滤波值中的误差罗德里格参数转化为四元数，并在下一次滤波开始前进行误差罗德里格参数置零，即
Figure BDA00030990372500001711
最终得到的姿态估计结果如图9-1、图9-2和图9-3所示。In the formula,
Figure BDA00030990372500001710
is the equipment observation value at the k+1th time, that is, the corrected satellite observation data Q J . Finally, convert the error Rodrigue parameter in the state filtering value into a quaternion, and set the error Rodrigue parameter to zero before the next filtering starts, that is,
Figure BDA00030990372500001711
The final attitude estimation results are shown in Figure 9-1, Figure 9-2 and Figure 9-3.

作为对比实验，将未校正的卫星观测数据QA作为无味滤波器中的观测数据，即步骤五中的状态更新中的式中，

Figure BDA00030990372500001712
估计结果如图10-1、图10-2和图10-3所示。从结果图中可知，基于修正后的卫星观测数据进行姿态估计，能够获得更快的估计准确度以及收敛速度。因此本发明具有较大的应用价值。As a comparative experiment, the uncorrected satellite observation data QA is used as the observation data in the odorless filter, that is, in the formula in the state update in step 5,
Figure BDA00030990372500001712
The estimation results are shown in Figure 10-1, Figure 10-2 and Figure 10-3. It can be seen from the result graph that the attitude estimation based on the corrected satellite observation data can obtain faster estimation accuracy and convergence speed. Therefore, the present invention has great application value.

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。The above are only specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art can easily think of changes or substitutions within the technical scope disclosed by the present invention. should be included within the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims (10)

1.一种地球反射光校正的卫星姿态估计方法，其特征在于，包括以下步骤：1. a satellite attitude estimation method of earth reflected light correction, is characterized in that, comprises the following steps: 步骤一：将地球表面按照经纬网格划分为多个面元，计算面元面积；基于太阳、地球和目标卫星之间的遮挡关系筛选有效面元；Step 1: Divide the earth's surface into multiple surfels according to the latitude and longitude grid, and calculate the surfel area; screen the effective surfels based on the occlusion relationship between the sun, the earth and the target satellite; 步骤二：根据辐射照度与距离的平方成反比定律，计算太阳可见光波段在各有效面元处的辐射照度，并根据朗伯反射原理计算有效面元在目标卫星处产生的辐射照度；Step 2: According to the law of irradiance being inversely proportional to the square of the distance, calculate the irradiance of the solar visible light band at each effective surface element, and calculate the irradiance generated by the effective surface element at the target satellite according to the Lambertian reflection principle; 步骤三：根据有效面元在目标卫星处产生的辐射照度，通过朗伯反射原理计算出以各个有效面元为光源时，目标卫星对探测器的辐射照度；将目标卫星对探测器的辐射照度相加得到以地球为光源时，目标卫星对探测器的辐射照度；Step 3: According to the irradiance generated by the effective surface element at the target satellite, the irradiance of the target satellite to the detector is calculated by the Lambertian reflection principle when each effective surface element is used as the light source; the irradiance of the target satellite to the detector is calculated. Add up to get the irradiance of the target satellite to the detector when the earth is the light source; 步骤四：从目标卫星处获得实时的目标卫星对探测器的辐射照度减去以地球为光源时，目标卫星对探测器的辐射照度，得到了仅以太阳为光源时的目标卫星对探测器的辐射照度的校正观测数据；Step 4: Obtain the real-time irradiance of the target satellite to the detector from the target satellite minus the irradiance of the target satellite to the detector when the earth is used as the light source, and obtain the irradiance of the target satellite to the detector when only the sun is used as the light source. Corrected observation data of irradiance; 步骤五：基于校正观测数据，使用无味卡尔曼滤波的方法对目标卫星进行姿态估计，得到地球反射光校正数据。Step 5: Based on the corrected observation data, use the method of tasteless Kalman filtering to estimate the attitude of the target satellite, and obtain the correction data of the earth's reflected light. 2.如权利要求1所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤一中，将地球表面按照经纬网格划分为多个面元的方法包括：基于地球表面的经纬网格，将面元中心位于经纬点，长和宽分别为一个纬度距离和一个经度距离。2. a kind of satellite attitude estimation method of earth reflected light correction as claimed in claim 1 is characterized in that, in described step 1, the method for dividing the earth surface into a plurality of surface elements according to the latitude and longitude grid comprises: based on the earth The latitude and longitude grid of the surface, the center of the bin is located at the latitude and longitude point, and the length and width are a latitude distance and a longitude distance respectively. 3.如权利要求1或2所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤一中，计算面元面积的方法包括：3. The method for estimating the satellite attitude of a kind of earth reflected light correction as claimed in claim 1 or 2, characterized in that, in the step 1, the method for calculating the surface area comprises: 每个地表面元的长均为：lj＝2πRE/360；The length of each surface element is: l j =2πR E /360; 纬度i处的地表面元，宽为：di＝2πREcosi/360；The ground surface element at latitude i, the width is: d i =2πR E cosi/360; 地表面元面积为：sij＝ljdi；The surface element area is: s ij =l j d i ; 其中，lj为经度j处地表面元的长，di为纬度i处的地表面元的宽，sij为面元面积，RE为地球半径，i、j分别表示纬度和经度。Among them, l j is the length of the surface element at longitude j, d i is the width of the surface element at latitude i, s ij is the area of the surface element, RE is the radius of the earth, and i and j represent the latitude and longitude, respectively. 4.如权利要求1-3之一所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤一中，筛选有效面元的方法包括：4. The method for estimating satellite attitude for earth reflected light correction according to one of claims 1 to 3, wherein in the step 1, the method for screening effective surface elements comprises:
Figure FDA0003099037240000021
Figure FDA0003099037240000021
 其中，θin为太阳光对地表面元的入射天顶角，θout为反射天顶角。Among them, θ in is the incident zenith angle of sunlight to the ground surface element, and θ out is the reflection zenith angle. 5.如权利要求1所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤二中，计算太阳可见光波段在各有效面元处的辐射照度的方法包括：5. The method for estimating satellite attitude corrected by earth reflected light as claimed in claim 1, wherein in the step 2, the method for calculating the irradiance of the solar visible light band at each effective surface element comprises: 根据辐射照度与距离的平方成反比定律，在太阳光不被遮挡的情况下，得到太阳可见光波段在各有效面元处的辐射照度：According to the law that the irradiance is inversely proportional to the square of the distance, when the sunlight is not blocked, the irradiance of the solar visible light band at each effective surface element is obtained:
Figure FDA0003099037240000022
Figure FDA0003099037240000022
 其中，R0为一年内日地的平均距离；Q0是在距离太阳R0处，可见光波段(0.4～0.7μm)部分的辐射照度，R为地表面元到太阳的距离，Q是在距离太阳R的有效面元处，对可见光波段(0.4～0.7μm)的辐射照度。Among them, R 0 is the average distance between the sun and the earth in one year; Q 0 is the irradiance in the visible light band (0.4-0.7 μm) at the distance from the sun at R 0 , R is the distance from the surface element to the sun, and Q is the distance between the sun and the sun. The irradiance of the visible light band (0.4-0.7μm) at the effective surface element of the sun R. 6.如权利要求1或5所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤二中，根据朗伯反射原理计算有效面元在目标卫星处产生的辐射照度，得出目标卫星所受单个地表面元的辐射照度的方法包括：6. The satellite attitude estimation method of a kind of earth reflected light correction as claimed in claim 1 or 5, it is characterized in that, in described step 2, according to Lambertian reflection principle, calculate the irradiance that effective surface element produces at the target satellite , the methods for deriving the irradiance of a single surface element on the target satellite include: 地球表面太阳辐射照度为：
Figure FDA0003099037240000031
The solar irradiance on the Earth's surface is:
Figure FDA0003099037240000031
 其中，R0为一年内日地的平均距离，Rse为太阳到地球的距离，Qse为太阳在地球表面的辐射照度，Q0为在距离太阳R0处，对可见光波段的辐射照度；Among them, R 0 is the average distance between the sun and the earth in one year, R se is the distance from the sun to the earth, Q se is the irradiance of the sun on the surface of the earth, and Q 0 is the irradiance of the visible light band at the distance R 0 from the sun; 地表面元对各方向产生的辐射亮度为：
Figure FDA0003099037240000032
The radiance generated by the surface element in each direction is:
Figure FDA0003099037240000032
 其中，θin为太阳光对地表面元的入射天顶角，Qsecosθin为地表面元所接收的太阳辐射照度，地球表面及地表大气总的平均反照率值为
Figure FDA0003099037240000033
Lij为地表面元对各方向产生的辐射亮度；Among them, θ in is the incident zenith angle of sunlight to the surface element, Q se cos θ in is the solar irradiance received by the surface element, and the total average albedo value of the earth's surface and the surface atmosphere is
Figure FDA0003099037240000033
L ij is the radiance generated by the ground surface element in all directions; 单个地表面元(i,j)在
Figure FDA0003099037240000034
方向目标卫星处产生的辐射照度为：A single ground element (i, j) is
Figure FDA0003099037240000034
The irradiance generated at the target satellite is:
Figure FDA0003099037240000035
Figure FDA0003099037240000035
 其中，
Figure FDA0003099037240000036
为面元指向目标卫星的方向向量，
Figure FDA0003099037240000037
为地表面元相对目标卫星所呈立体角，sij为经纬度位于(i,j)的地表面元面积，
Figure FDA0003099037240000038
为地表面元至目标卫星的距离，Qij为单个地表面元(i,j)在
Figure FDA0003099037240000039
方向目标卫星处产生的辐射照度。in,
Figure FDA0003099037240000036
is the direction vector of the surfel pointing to the target satellite,
Figure FDA0003099037240000037
is the solid angle of the surface element relative to the target satellite, s ij is the area of the surface element located at (i, j) in latitude and longitude,
Figure FDA0003099037240000038
is the distance from the surface element to the target satellite, and Q ij is the single surface element (i, j) in
Figure FDA0003099037240000039
The irradiance generated at the target satellite. 7.如权利要求1所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤三中，根据有效面元在目标卫星处产生的辐射照度，通过朗伯反射原理计算出以各个有效面元为光源时，目标卫星对探测器的辐射照度的方法包括：7. The satellite attitude estimation method of a kind of earth reflected light correction as claimed in claim 1, is characterized in that, in described step 3, according to the irradiance that effective surface element produces at target satellite, calculate by Lambertian reflection principle When each effective surface element is used as the light source, the method for determining the irradiance of the target satellite to the detector includes:
Figure FDA00030990372400000310
Figure FDA00030990372400000310
 其中，太阳指向目标的向量即入射光向量为λS；目标指向探测器的向量即观测向量为λF；目标卫星某面元的法向量为n；入射向量与面元法向量的夹角为到入射天顶角，记作θS，观测向量与面元法向量的夹角为观测天顶角，记作θF；根据入射天顶角θS、观测天顶角θF的大小关系判断探测器是否能收到目标卫星反射光；Among them, the vector that the sun points to the target, that is, the incident light vector is λ S ; the vector that the target points to the detector, that is, the observation vector is λ F ; the normal vector of a certain panel of the target satellite is n; the angle between the incident vector and the normal vector of the panel is To the incident zenith angle, denoted as θ S , the angle between the observation vector and the normal vector of the surface element is the observation zenith angle, denoted as θ F ; Judging according to the relationship between the incident zenith angle θ S and the observation zenith angle θ F Whether the detector can receive the reflected light from the target satellite; 目标卫星上第i个面元所接收的单个地表面元的辐射照度为：The irradiance of a single ground surface element received by the i-th surface element on the target satellite is: Qls(i)＝Qls·cos(θS(i))；Q ls (i)=Q ls ·cos(θ S (i)); 其中，Qls表示光源在目标卫星处的辐射照度为，θS(i)表示第i个地表面元的光线入射天顶角，Qls(i)为目标卫星上第i个面元所接收的单个地表面元的辐射照度；根据辐射照度传递公式，目标卫星第i面元对探测器的辐射照度为：Among them, Q ls represents the irradiance of the light source at the target satellite, θ S (i) represents the light incident zenith angle of the i-th surface element, and Q ls (i) is received by the i-th surface element on the target satellite The irradiance of a single ground surface element; according to the irradiance transfer formula, the irradiance of the i-th surface element of the target satellite to the detector is:
Figure FDA0003099037240000041
Figure FDA0003099037240000041
 其中，αi为第i面元的反射率，S(i)为第i面元的面积，Rsf为目标卫星距探测器的距离，Qsf(i)为目标卫星第i面元对探测器的辐射照度，θF(i)为观测天顶角；目标卫星面元遮挡判别系数η为：Among them, α i is the reflectivity of the i-th surface element, S(i) is the area of the i-th surface element, R sf is the distance between the target satellite and the detector, and Q sf (i) is the detection of the i-th surface element of the target satellite. The irradiance of the sensor, θ F (i) is the observation zenith angle; the target satellite panel occlusion discriminant coefficient η is:
Figure FDA0003099037240000042
Figure FDA0003099037240000042
 其中，θS为入射天顶角、θF为观测天顶角；Among them, θ S is the incident zenith angle, and θ F is the observation zenith angle; 以单个地球有效面元为唯一光源时目标卫星对探测器的辐射照度为：
Figure FDA0003099037240000043
其中，Qsf(i)为目标卫星第i面元对探测器的辐射照度，η(i)为目标卫星第i面元遮挡判别系数，n为目标卫星总面元数，Qsf为目标卫星对探测器的总辐射照度。When a single effective surface element of the earth is the only light source, the irradiance of the target satellite to the detector is:
Figure FDA0003099037240000043
Among them, Q sf (i) is the irradiance of the ith panel of the target satellite to the detector, η(i) is the occlusion discrimination coefficient of the ith panel of the target satellite, n is the total number of panels of the target satellite, and Q sf is the target satellite The total irradiance to the detector. 8.如权利要求1或7所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤三中，将目标卫星对探测器的辐射照度相加得到以地球为光源时，目标卫星对探测器的辐射照度，包括：8. the satellite attitude estimation method of a kind of earth reflected light correction as claimed in claim 1 or 7, it is characterized in that, in described step 3, the irradiance illuminance of target satellite to detector is added to obtain when taking earth as light source , the irradiance of the target satellite to the detector, including: 根据以单个有效面元为唯一光源时目标卫星对探测器的辐射照度，针对所有有效面元进行叠加可以得到基于整个地球反射光的目标卫星辐射照度为
Figure FDA0003099037240000051
S有效为地表有效面元区域。According to the irradiance of the target satellite to the detector when a single effective surface element is the only light source, the irradiance of the target satellite based on the reflected light of the whole earth can be obtained by superimposing all the effective surface elements as:
Figure FDA0003099037240000051
S is effectively the surface effective bin area. 9.如权利要求1所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤四中，校正观测数据的校正方法包括：9. The method for estimating satellite attitude of earth reflected light correction according to claim 1, wherein in the step 4, the correction method for correcting the observation data comprises:
Figure FDA0003099037240000052
Figure FDA0003099037240000052
 其中，从目标卫星处获得的实时的目标卫星对探测器的辐射照度为QA，校正观测数据为QJ，以地球为光源时，目标卫星对探测器的辐射照度为
Figure FDA0003099037240000053
Among them, the real-time irradiance of the target satellite to the detector obtained from the target satellite is Q A , the corrected observation data is Q J , and when the earth is used as the light source, the irradiance of the target satellite to the detector is
Figure FDA0003099037240000053
 10.如权利要求1所述的一种地球反射光校正的卫星姿态估计方法，其特征在于，所述步骤五中，目标卫星的姿态估计方法包括：10. The method for estimating the attitude of a satellite based on earth reflected light correction according to claim 1, wherein in the step 5, the method for estimating the attitude of the target satellite comprises: 状态向量
Figure FDA0003099037240000054
其中，δp表示误差罗德里格参数；ω为目标卫星翻滚角速度，目标姿态指目标本体系相对地球惯性坐标系的姿态；state vector
Figure FDA0003099037240000054
Among them, δp represents the error Rodrigue parameter; ω is the roll angular velocity of the target satellite, and the target attitude refers to the attitude of the target system relative to the earth's inertial coordinate system; 状态向量经过无味变换得到2n+1个sigma点，其表达式为：The state vector is tastelessly transformed to obtain 2n+1 sigma points, and its expression is:
Figure FDA0003099037240000055
Figure FDA0003099037240000055
 其中，λ＝α2(n+κ)-n，为缩放比例系数，通常有κ＝3-n；
Figure FDA0003099037240000061
为当前状态量的第i个sigma点，σ(i)为第i个sigma点；将第k时刻的滤波值
Figure FDA0003099037240000062
作为当前的均值sigma点，记为
Figure FDA0003099037240000063
为k时刻状态协方差矩阵的滤波值；用Pi表示矩阵P的第i列；chol(P)表示对矩阵P进行cholesky分解，其对应权值为：Among them, λ=α 2 (n+κ)-n, is the scaling factor, usually κ=3-n;
Figure FDA0003099037240000061
is the i-th sigma point of the current state quantity, and σ(i) is the i-th sigma point;
Figure FDA0003099037240000062
As the current mean sigma point, denoted as
Figure FDA0003099037240000063
is the filter value of the state covariance matrix at time k; P i represents the i-th column of the matrix P; chol(P) represents the cholesky decomposition of the matrix P, and its corresponding weight is:
Figure FDA0003099037240000064
Figure FDA0003099037240000064
 其中，n为状态向量的维数，Wmean表示均值，Wcov表示协方差，α为控制sigma点扩散程度的系数，取值范围通常为[10-4,1]；β为权值调整系数，通常取值为2；λ＝α2(n+κ)-n，为缩放比例系数，通常有κ＝3-n；获取误差四元数的sigma，Among them, n is the dimension of the state vector, W mean is the mean value, W cov is the covariance, α is the coefficient that controls the spread of the sigma point, the value range is usually [10 -4 ,1]; β is the weight adjustment coefficient , usually the value is 2; λ=α 2 (n+κ)-n, is the scaling factor, usually κ=3-n; to obtain the sigma of the error quaternion,
Figure FDA0003099037240000065
Figure FDA0003099037240000065
Figure FDA0003099037240000066
Figure FDA0003099037240000066
 其中，a＝1,f＝2(a+1)，
Figure FDA0003099037240000067
表示四元数的前三个分量
Figure FDA0003099037240000068
表示四元数的第四个分量
Figure FDA0003099037240000069
表示四元数误差分量；再进行误差四元数向四元数的转换，
Figure FDA00030990372400000610
Among them, a=1, f=2(a+1),
Figure FDA0003099037240000067
represents the first three components of a quaternion
Figure FDA0003099037240000068
represents the fourth component of the quaternion
Figure FDA0003099037240000069
represents the quaternion error component; then convert the error quaternion to quaternion,
Figure FDA00030990372400000610
 其中，
Figure FDA00030990372400000611
表示当前四元数sigma点的均值，
Figure FDA00030990372400000612
表示第i个sigma点的误差四元数根据四元数乘法转换得到的四元数值，
Figure FDA00030990372400000613
表示第i个sigma点的误差四元数；四元数乘法定义如下，in,
Figure FDA00030990372400000611
represents the mean of the current quaternion sigma point,
Figure FDA00030990372400000612
Represents the quaternion value obtained by converting the error quaternion of the i-th sigma point according to quaternion multiplication,
Figure FDA00030990372400000613
represents the error quaternion of the i-th sigma point; the quaternion multiplication is defined as follows,
Figure FDA00030990372400000614
Figure FDA00030990372400000614
 其中，
Figure FDA0003099037240000071
对于三维向量a，[a×]表示向量叉积矩阵，
Figure FDA0003099037240000072
qa和qb表示任意两个四元数；in,
Figure FDA0003099037240000071
For a three-dimensional vector a, [a×] represents the vector cross product matrix,
Figure FDA0003099037240000072
q a and q b represent any two quaternions; 四元数sigma点和角速度sigma点传递，
Figure FDA0003099037240000073
其中，Quaternion sigma point and angular velocity sigma point transfer,
Figure FDA0003099037240000073
in,
Figure FDA0003099037240000074
Figure FDA0003099037240000074
 Δt为探测器采样时间间隔；Γ(ωk)由动力学方程在时域离散化得到；基于以上传递后的四元数sigma点，再进行误差罗德里格参数点的转换，
Figure FDA0003099037240000075
Δt is the sampling time interval of the detector; Γ(ω k ) is obtained by discretizing the dynamic equation in the time domain; based on the quaternion sigma point after the above transfer, the error Rodrigue parameter point is converted,
Figure FDA0003099037240000075
 其中，共轭四元数q-1定义为
Figure FDA0003099037240000076
where the conjugate quaternion q -1 is defined as
Figure FDA0003099037240000076
 继而，
Figure FDA0003099037240000077
Then,
Figure FDA0003099037240000077
 其中，a＝1,f＝2(a+1)。令其均值sigma点
Figure FDA0003099037240000078
得到传递后状态向量的sigma点为
Figure FDA0003099037240000079
Among them, a=1, f=2(a+1). Let its mean sigma point
Figure FDA0003099037240000078
The sigma point of the state vector after the transfer is obtained as
Figure FDA0003099037240000079
 最后通过加权求和的方式进行状态更新和协方差更新，
Figure FDA00030990372400000710
Finally, the state update and covariance update are performed by means of weighted summation,
Figure FDA00030990372400000710
 式中，Qk+1为离散时间的过程噪声协方差；由观测方程可计算得到每个sigma点的观测预测值，
Figure FDA0003099037240000081
h表示观测方程；
Figure FDA0003099037240000082
为sigma点对应的观测预测值；再通过加权求和得到观测预测均值以及量测协方差和互协方差：
Figure FDA0003099037240000083
In the formula, Q k+1 is the process noise covariance in discrete time; the observed predicted value of each sigma point can be calculated from the observation equation,
Figure FDA0003099037240000081
h represents the observation equation;
Figure FDA0003099037240000082
is the observed predicted value corresponding to the sigma point; then the mean of the observed predicted and the measured covariance and cross-covariance are obtained by weighted summation:
Figure FDA0003099037240000083
Figure FDA0003099037240000084
Figure FDA0003099037240000084
 Rk+1表示
Figure FDA0003099037240000085
表示量测协方差，
Figure FDA0003099037240000086
表示互协方差，
Figure FDA0003099037240000087
表示第k+1步的观测预测值；卡尔曼增益计算，
Figure FDA0003099037240000088
Kk+1表示第k+1次滤波的卡尔曼增益；R k+1 means
Figure FDA0003099037240000085
represents the measurement covariance,
Figure FDA0003099037240000086
represents the cross-covariance,
Figure FDA0003099037240000087
Indicates the observed predicted value of the k+1 step; Kalman gain calculation,
Figure FDA0003099037240000088
K k+1 represents the Kalman gain of the k+1th filtering; 状态更新和协方差更新，
Figure FDA0003099037240000089
state updates and covariance updates,
Figure FDA0003099037240000089
 式中，
Figure FDA00030990372400000810
为第k+1时刻的设备观测值，即校正后的目标卫星观测数据QJ；将状态滤波值中的误差罗德里格参数转化为四元数，并在下一次滤波开始前进行误差罗德里格参数置零，即
Figure FDA00030990372400000811
表示第k+1次滤波时的状态向量，03×1表示3行1列的零矩阵。In the formula,
Figure FDA00030990372400000810
is the equipment observation value at the k+1th moment, that is, the corrected target satellite observation data Q J ; the error Rodrigue parameter in the state filtering value is converted into a quaternion, and the error Rodrigue is carried out before the next filtering starts. The parameter is set to zero, i.e.
Figure FDA00030990372400000811
Represents the state vector during the k+1th filtering, and 0 3×1 represents a zero matrix with 3 rows and 1 column.
CN202110619180.1A 2021-06-03 2021-06-03 A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction Active CN113361163B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110619180.1A CN113361163B (en) 2021-06-03 2021-06-03 A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110619180.1A CN113361163B (en) 2021-06-03 2021-06-03 A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction

Publications (2)

Publication Number Publication Date
CN113361163A CN113361163A (en) 2021-09-07
CN113361163B true CN113361163B (en) 2022-09-30

Family

ID=77531766

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110619180.1A Active CN113361163B (en) 2021-06-03 2021-06-03 A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction

Country Status (1)

Country Link
CN (1) CN113361163B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115077500B (en) * 2022-05-07 2024-03-29 中国人民解放军国防科技大学 A method for determining ground sunlight reflection points and its related components

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102564574B (en) * 2011-12-29 2014-04-02 北京控制工程研究所 Method for measuring radiant illumination of earth albedo
CN105352609B (en) * 2015-11-13 2018-06-01 北京空间飞行器总体设计部 A kind of Optical remote satellite absolute radiation calibration method based on space lambert's sphere

Also Published As

Publication number Publication date
CN113361163A (en) 2021-09-07

Similar Documents

Publication Publication Date Title
CN111563962B (en) Remote sensing image simulation method based on geometric radiation integrated sampling
De et al. Palomar gattini-ir: Survey overview, data processing system, on-sky performance and first results
Calcabrini et al. A simplified skyline-based method for estimating the annual solar energy potential in urban environments
Rahmati et al. MAVEN measured oxygen and hydrogen pickup ions: Probing the Martian exosphere and neutral escape
CN104573251B (en) A kind of star-loaded optical remote sensing device full filed Apparent spectral radiance determines method
CN101598797B (en) Method for realizing rugged topography remote sensing scene simulation
Andrews et al. The effect of spectral albedo on amorphous silicon and crystalline silicon solar photovoltaic device performance
Brunger et al. Anisotropic sky radiance model based on narrow field of view measurements of shortwave radiance
Frieß et al. Intercomparison of aerosol extinction profiles retrieved from MAX-DOAS measurements
Notton et al. Neural network approach to estimate 10-min solar global irradiation values on tilted planes
Luque et al. Diffusing reflectors for bifacial photovoltaic panels
CN103675794A (en) Spaceflight optical remote sensing imaging simulation method based on space-time unified feature
CN104101297B (en) Space object dimension acquisition method based on photoelectric observation
Salo et al. Photometric modeling of Saturn's rings: I. Monte Carlo method and the effect of nonzero volume filling factor
Machado et al. Mapping zonal winds at Venus’s cloud tops from ground-based Doppler velocimetry
CN101450716A (en) Fault photo-detection method for earth synchronous transfer orbit satellite in orbit
Redweik et al. 3D local scale solar radiation model based on urban LiDAR data
Andres et al. Time‐varying, ray tracing irradiance simulation approach for photovoltaic systems in complex scenarios with decoupled geometry, optical properties and illumination conditions
CN113361163B (en) A Satellite Attitude Estimation Method Based on Earth Reflected Light Correction
CN104880701A (en) Satellite-borne sensor imaging simulation method and device
Mangano et al. Dynamical evolution of sodium anisotropies in the exosphere of Mercury
Santos-Martin et al. SoL–A PV generation model for grid integration analysis in distribution networks
Stern et al. SPATIALLY RESOLVING THE KINEMATICS OF THE $\lesssim 100\;\mu {\rm as} $ QUASAR BROAD-LINE REGION USING SPECTROASTROMETRY
CN101876700A (en) A method for simulating radiative transfer in complex terrain areas based on radiosity
CN102564574B (en) Method for measuring radiant illumination of earth albedo

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant