CN105928527B - Satellite attitude determination method based on visible light earth sensor and sun sensor - Google Patents
Satellite attitude determination method based on visible light earth sensor and sun sensor Download PDFInfo
- Publication number
- CN105928527B CN105928527B CN201610262687.5A CN201610262687A CN105928527B CN 105928527 B CN105928527 B CN 105928527B CN 201610262687 A CN201610262687 A CN 201610262687A CN 105928527 B CN105928527 B CN 105928527B
- Authority
- CN
- China
- Prior art keywords
- satellite
- angle
- axis
- circle
- coordinate
- 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
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Astronomy & Astrophysics (AREA)
- Automation & Control Theory (AREA)
- General Physics & Mathematics (AREA)
- Navigation (AREA)
Abstract
Satellite attitude determination method based on visible light earth sensor and sun sensor is combined using visible light earth sensor and sun sensor and determines the three-axis attitude of satellite over the ground.First by extracting the efficient frontier of earth image to the processing of earth image formed by cmos imaging device, the center of the efficient frontier then extracted using least square method fitting.Secondly under the satellite hypothesis that attitude angle is low-angle over the ground, the roll angle over the ground and pitch angle of approximate calculation satellite.The yaw angle of satellite over the ground is calculated in the roll angle and pitch angle posture information finally obtained using sun sensor and approximate calculation.The triaxial attitude angle of satellite over the ground is determined relative to traditional combine using earth infrared senstive device and sun sensor, the method of the present invention improves attitude of satellite angular measurement precision, and this method is suitable for the small satellite attitude measurement of miniaturization, complies with the trend of current moonlet.
Description
Technical field
The invention belongs to satellite control fields, are related to a kind of method of determining attitude of satellite.
Background technique
Currently, infrared earth sensor obtains in the sensor for determining attitude of satellite angle using the earth as reference vector
It is widely applied, and reaches its maturity, therefore combined using infrared earth sensor and sun sensor and determine satellite over the ground
Three-axis attitude become common methods.
The existing method for determining satellite three-axis attitude over the ground of being combined using infrared earth sensor and sun sensor is:
On the one hand, earth carbon dioxide infra-red radiation band is passed in and out using infrared earth sensor sweep mechanism and infrared scan mechanism arrives
Up to the time of datum mark, the chord width for calculating the infrared earth sensor scanning earth and geocentric vector E are in infrared earth sensor
Measurement coordinate system xSEoySEProjection and x in planeSEThe angle of axis, by the resulting angle of above-mentioned calculating and spherical trigonometry
Knowledge can calculate coordinate representation of the geocentric vector E in infrared earth sensor measurement coordinate system, and quick by infrared earth
The installation matrix of sensor can obtain coordinate representation of the geocentric vector E in celestial body coordinate system；On the other hand, geocentric vector E is in satellite
The coordinate of orbital coordinate system it is known that by the spin matrix that unknown attitude of satellite angle calculates can calculate geocentric vector E is sat in celestial body
Mark the coordinate representation of system, different expression formulas of the geocentric vector being converted by above two approach under celestial body coordinate system
The roll angle and pitch angle of satellite are found out, then, sun sensor is asked using the satellite roll angle and pitch angle being calculated
Obtain the yaw angle of satellite.
Combine the principle for determining satellite three-axis attitude over the ground by above-mentioned infrared earth sensor and sun sensor it is found that defending
The measurement accuracy of star attitude angle is limited to the measurement accuracy and satellite rail of time resolution and sweep mechanism angular scanning speed
Road parameter measurement precision (depending on different installation matrixes).Since above-mentioned measurement physical quantity all has uncertainty, cause to defend
The measurement accuracy fluctuation of star posture is larger.
In addition, present satellites are just sent out towards miniaturization, intelligence, high-precision, high stability, low-power consumption and long-life direction
Exhibition, and main infrared earth sensor has conical scanning and oscillatory scanning etc., this kind of infrared earth sensor is swept due to existing
Mechanism is retouched, quality is big, power consumption is more, is also not suitable for the trend of present satellites.
Summary of the invention
Technical problem solved by the present invention is having overcome the deficiencies of the prior art and provide a kind of quick based on the visible light earth
The satellite attitude determination method of sensor and sun sensor, solves existing infrared earth sensor and sun sensor joint is true
Determining time resolution in the attitude of satellite, sweep mechanism angular velocity measurement precision, satellite orbit parameter measurement accuracy leads to satellite
The problem of attitude angle accuracy decline, and this method is applicable in the small satellite attitude measurement of miniaturization, complies with current moonlet
Trend.
The technical solution of the invention is as follows: being determined based on the attitude of satellite of visible light earth sensor and sun sensor
Method includes the following steps:
(1) using the cmos imaging device in visible light earth sensor, ball is imaged over the ground, obtains input picture；
(2) input picture is pre-processed, removes noise, smoothed image；
(3) edge detection is carried out to the image that step (2) obtain, obtains efficient frontier；
(4) efficient frontier is fitted, obtains efficient frontier center and radius；
(5) the roll angle φ obtained satellite over the ground and pitching angle theta are calculated, wherein
(Cx,Cy) it is efficient frontier centre coordinate,R is the radius of efficient frontier, and f is cmos imager
The focal length of part；
(6) it calculates and obtains the yaw angle of satellite over the groundWhereinDssBIt is defended for sun sensor is opposite
The installation matrix of star body coordinate system,S in formulaAX、SAY、SAZRespectively indicate sun sensor survey
Measure coordinate system xssAxis and satellite body coordinate system xb、yb、zbThe included angle cosine value of three axis, SBX、SBY、SBZRespectively indicate sun sensitivity
Device measures coordinate system yssAxis and satellite body coordinate system xb、yb、zbThe included angle cosine value of three axis, SCX、SCY、SCZRespectively indicate the sun
Sensor measures coordinate system zssAxis and satellite body coordinate system xb、yb、zbThe included angle cosine value of three axis；So=[SoX SoY SoZ]T
The three axis component forms for being solar vector under orbital coordinate system, α are that solar direction vector in sun sensor measures coordinate system
yssozssPlane projection and ozssThe angle of axis.
Input picture in the step (1) should meet: cmos imaging device institute is at earth image Horizon profile real arcs
Length accounts for the 50% of total arc length.The method of smoothed image is the median filtering of 3 × 3 convolution masks in the step (2).The step
Suddenly in (3) edge detection method are as follows: all edges are detected using Sobel operator first, are then detected using Hough transform
Efficient frontier out.It is least square method to the method that efficient frontier is fitted in the step (4).
The advantages of the present invention over the prior art are that:
(1) the method for the present invention is combined by visible light earth sensor and sun sensor determines the three axis appearances of satellite over the ground
State angle, compared to the attitude angle over the ground that traditional infrared earth sensor and sun sensor determine satellite, this method eliminates
Time resolution, sweep mechanism angular velocity measurement precision, orbital-elements measurement precision measure satellite to infrared earth sensor
The influence of attitude angle, so that the precision of satellite attitude measurement is improved, and this method realizes and be easy that algorithm is simple, it is flat to satellite
The control of platform is of less demanding, is the determination satellite three-axis attitude method over the ground of complete set；
(2) extraction for focusing on efficient frontier of the method for the present invention, the essence of the quality of edge extracting directly determining method
Degree.The extraction edge method that the present invention takes can not only accurately extract efficient frontier, meet the requirement of calculation amount, and algorithm
Robustness is fine；
(3) the method for the present invention can be integrated into controller, lightweight be collectively formed with other parts, precision is met the requirements
Visible light sensor, for determine the attitude of satellite another solution is provided.
Detailed description of the invention
Fig. 1 is the flow diagram of the method for the present invention；
Fig. 2 is CMOS camera imaging schematic diagram of the present invention；
Fig. 3 is attitude of satellite coordinate system schematic diagram of the present invention；
Fig. 4 is sun sensor coordinate system schematic diagram of the present invention.
Specific embodiment
With the fast development of microminaturization cmos imaging device, the real use of only earth visible light sensor is not mentioned
Very big possibility has been supplied, and has complied with the development trend of satellite.The present invention is based on visible light earth sensor and the sun is sensitive
How device combines determining satellite three-axis attitude over the ground, the method for proposing complete set.
As shown in Fig. 1 process, the method for the present invention is main including the following steps: and first, cmos imaging device ball over the ground
Input picture after imaging；Second, image preprocessing；Third extracts image efficient frontier；4th, using least square fitting figure
As efficient frontier center；5th, under in satellite, attitude angle is low-angle supposed premise over the ground, according to efficient frontier center approximation meter
Satellite is calculated to roll and pitch attitude angle；6th, the roll angle and pitch angle posture obtained using sun sensor and approximate calculation
Information calculates the yaw angle of satellite over the ground.Above-mentioned six steps are discussed in detail individually below.
First, it is seen that cmos imaging device in light earth sensor input picture after ball imaging over the ground.
In order to guarantee the satellite pitch angle finally determined and roll angular accuracy, it is desirable that cmos imaging device institute is at earth image
Horizon profile has sufficiently long radian.Generally the effective arc length of imaging plane is required to account for the 50% of total arc length.
Image captured by cmos imaging device meets input picture after arc length requires, and carries out the processing of next step.
Second, image preprocessing.
Pretreated purpose is carried out to image to be to remove noise, smoothed image.The method of image preprocessing has very much, such as
Median filtering, mean filter etc., the present invention in use median filtering.
Median filtering uses 3 × 3 convolution masks, and specific practice is as follows: by 8 pixels in center pixel and its 3 × 3 neighborhoods
Gray value by sorting from small to large, take the 5th element gray value be center pixel gray value.Median filtering can not only remove
Noise spot reduces calculation amount for later edge extracting, moreover it is possible to which the edge for protecting image obtains relatively satisfactory image preprocessing effect
Fruit.
Third detects efficient frontier.
This step is committed step and core of the method for the present invention by solving the problems, such as above-mentioned proposition.Detection has
The quality of effect edge method directly determines the size of calculation amount needed for the method for the present invention and the precision of determined attitude angle.
At present there are many kinds of the algorithms at detection edge, common are: differential operator, Laplce's Gauss operator, Canny are calculated
Son, fitting process etc..The present invention detects efficient frontier using Sobel differential operator and Hough transform, and Sobel operator is detection institute
Some edges, including efficient frontier and invalid edge, and Hough transform is then to reject invalid edges detected by Sobel operator
Edge retains its efficient frontier.Both methods is introduced separately below.
Sobel operator calculates the shade of gray amplitude at a certain pixel first, then by this value and set threshold values into
Row compares, if shade of gray amplitude is greater than threshold values at the point, retains the point, otherwise rejects the point.The middle thought of this method
Think it is that gray-value variation is more violent in image edge pixels vertex neighborhood, adopting said method determines gray-value variation ratio in neighborhood
Biggish pixel, to detect image border.The size of set threshold values is constant in Sobel operator, and threshold values size needs root
It is determined according to real image, the size of general threshold values is chosen for the maximum positive integer no more than whole image gray value average value.
The method that Sobel operator calculates shade of gray amplitude at pixel is described below.
Gradient magnitude calculates as follows at any point:
Wherein sx, syX, y directional derivative at this point is respectively indicated, and for digital picture, it can be replaced with first-order difference
First differential.Usually calculating shade of gray amplitude at a certain pixel has 3 × 3 templates, 5 × 5 templates etc., and the present invention uses 3 × 3
Template, therefore sx, syIt is writeable are as follows:
sx=(a2+2a3+a4)-(a0+a7+a6) (2)
sy=(a0+2a1+a2)-(a6+2a5+a4) (3)
Wherein a0…a7Represent the gray value of 3 × 3 neighborhood interior elements at the pixel, a0…a7Each pixel gray value is 3
Position is as follows in × 3 neighborhoods:
| a0 | a1 | a2 |
| a7 | [i,j] | a3 |
| a6 | a5 | a4 |
Note: [i, j] represents required shade of gray amplitude pixel coordinate.
As previously mentioned, the edge that Sobel operator detects includes efficient frontier and invalid edge, wherein efficient frontier is ground
The contour edge of ball, invalid edge include earth terminator edge, image background bright celestial body edge, moon edge etc..In order to
The precision for improving determined attitude angle, needs further to detect the edge that Sobel operator detects, rejects invalid edge, protects
Stay efficient frontier.The Hough transform of detection efficient frontier method is described below.
Hough transform has two kinds of known image radius of circle Hough transform and unknown images radius of circle Hough transform.For
This problem would naturally fall within the latter, a kind of because the radius of formed image circle can not be obtained accurately in advance after earth imaging
The Hough transform method of quick unknown images radius of circle is as follows:
Firstly, for round detection, using gradient information, Ballard is detected with following formula to be justified
x0=x ± r cos (h) (4)
y0=y ± r sin (h) (5)
In formula, (x0,y0) it is central coordinate of circle, h is gradient direction angle at the edge of the circle point (x, y), and r is image radius of circle.
Based on above-mentioned formula, detecting step is as follows:
1. first according to priori knowledge or the variation range r ∈ [r of given data estimation parameter radius rmin,rmax] and
Step delta r size.
For this problem, radius r range determination approximate circle on imaging surface is as follows after earth imaging:
Firstly, camera imaging schematic diagram is as shown in Fig. 2, in figure, R is the distinguishable atmosphere layer height arrived of camera, changes model
Enclose R ∈ [6371,6371+120]；H is satellite altitude, it is known that value；For half angle, sizeRE
For earth radius, then it thus can determine half angleRange；F is camera focus, it is known that value；R is imaging surface after earth imaging
Upper approximate circle radius, byIt can determine the range r ∈ [r of rmin,rmax]；Step delta r is by sensor controller
Computing capability and required precision etc. determine that, if the computing capability of controller is strong, that step-length can be selected a little bit smaller, no
Only meet attitude accuracy requirement, while meeting the requirement of posture real-time.
From radius r0=rminStart, recycles radius each time and increase Δ r, until r=rmax, for any circulation, radius ri
=rmin+ Δ r × i, (i=0,1 ... m), and wherein m is the total degree of circulation, expression formulaEach edge is sat
Punctuate (xj, yj) (j=1 ... n) and its gradient direction angle hjIt successively substitutes into above-mentioned (4) and (5) Ballard formula, calculates
Central coordinate of circle (x0ij,y0ij), wherein j indicates that effectively resulting central coordinate of circle four then will be calculated with invalid marginal point total number
House five enters to be rounded, and above-mentioned each marginal point coordinate is to the right x-axis, is upwards y-axis using the image lower left corner as coordinate origin
It is indicated under coordinate system.It is negative to prevent the integer central coordinate of circle being calculated, by original coordinate system respectively along x-axis, y-axis losing side
To 1024 pixels of translation, new positive integer central coordinate of circle value is obtained after translation, using new coordinate value as the seat of two-dimensional array
Scale value, the Data-Statistics of two-dimensional array add up the number of two-dimensional coordinate value appearance, all marginal points are traversed, finally by cumulative time
Number peak value determines this time circulation r=riWhen Circle Parameters ci={ x0i,y0i,ri, it recycles determine one group of parameter c each timei=
{x0i,y0i,riAnd accumulative frequency Ki.Pass through this cycle accumulor number KiWith last cycle accumulor number Ki-1Comparison, note
Record relatively high accumulative frequency KmaxAnd corresponding Circle Parameters cmax={ X0,Y0, R }, work as ri=rmaxWhen, the circle finally determined is joined
Number cmax={ X0,Y0, R } and it is efficient frontier Pixel-level central coordinate of circle { X0,Y0And radius R.
(4) and (5) Ballard formula takes positive sign, and the sine and cosine of shade of gray deflection h calculates as follows:
Wherein sx,syIt can be calculated by (2) and (3).
2. unknown radius Hough transform obtains Pixel-level central coordinate of circle { X0,Y0And radius R, if certain marginal point is to have
Marginal point, i.e. boundary point on circumference are imitated, then the point and Pixel-level central coordinate of circle { X0,Y0Between distance meet:
The Clutter edge point for not meeting relational expression (8) will be deleted, to obtain efficient frontier.
4th, using least square fitting efficient frontier center；
After obtaining efficient frontier, least square fitting efficient frontier center is next applied.
To one group of given data point (xk,yk) (k=1 ... l), l be data point number, at least 3, if with this group of data
The quadratic sum of the distance of point is the equation of the smallest circle are as follows:
(x-x0)2+(y-y0)2=r2 (9)
(x, y) indicates the coordinate value on above-mentioned circumference, (x in formula0,y0) indicate that the central coordinate of circle of above-mentioned circle, r indicate above-mentioned
Round radius.
A is set again1=∑ xk, b1=∑ yk, a2=∑ (xk)2, b2=∑ (yk)2, a3=∑ (xk)3, b3=∑ (yk)3, c11=
∑xk*yk, c12=∑ xk*yk 2, c21=∑ (xk)2*yk。
If f (x0, y0, R) and=∑ ((xk-x0)2+(yk-y0)2-r2)2, enable respectively
I.e.
-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
By a1, b1, a2, b2, a3, b3, c11, c12, c21Three equations above are substituted into arrange:
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
Solve central coordinate of circle and radius are as follows:
(10), (11), (12) calculating acquired results are coordinate values as a result, without actual physical meaning, are needed in the 5th step
Have physical significance as a result, set unit pixel length as d, to above-mentioned (10), (11), (12) calculated result all multiplied by unit picture
Plain length d, acquired results are the central coordinate of circle and radius for having length physical significance:
Cx=x0×d (13)
Cy=y0×d (14)
R=r × d (15)
5th, calculate satellite roll angle over the ground and pitch angle.
Before calculating attitude of satellite angle, attitude of satellite coordinate system O-XYZ (such as attached drawing 3) and sensor coordinate are initially set up
Device O-XTYTZT.Attitude of satellite coordinate system O-XYZ is defined as follows: origin O is centroid of satellite, and OX axis is directed toward satellite flight direction, OZ
Axis is directed toward the earth's core, and OY axis is determined by the right-hand rule.Sensor coordinates device OT-XTYTZTEach reference axis with attitude of satellite coordinate
It is that reference axis is parallel, unique difference origin OTPositioned at the mass center of sensor.
According to the definition of above-mentioned coordinate system, attitude angle has as described below: satellite rotates to be roll angle around X-axis, rotates around Y-axis
For pitch angle, yaw angle is rotated to be about the z axis, with around X-axis, Y-axis, Z axis rotates in the forward direction the attitude angle being positive.
On this basis, the roll angle φ and pitching angle theta of satellite are determined: in attached drawing 3,For half angle, can be obtained by attached drawing 2In formula, R is acquired by (15), and f is it is known that then half angleHave:
If coordinate of the earth's core in imaging surface is (Mx,My), by image geometry center, that is, resulting efficient frontier of least square
Centre coordinate is (Cx,Cy), wherein coordinate (M of the earth's core in imaging surfacex,My) have:
Mx=f*tan θ (17)
My=f*tan φ (18)
Derive efficient frontier centre coordinate (Cx,Cy) calculation expression, for describe the problem for the sake of simplicity, assume satellite roll
Dynamic angle φ=0, pitching angle theta ≠ 0, the then coordinate on imaging surface in X-direction apart from efficient frontier center farthest point and closest approach
RespectivelyThen efficient frontier center X-coordinate are as follows: Cx=Xmax-Xmin,
That is:
It can similarly obtain:
In the case where roll angle φ and pitching angle theta are low-angle, can be obtained by (17) and (19), (18) and (20):
Therefore by (17) and (21), (18) and (22) and half angleCalculation formula (16) can obtain satellite roll angle φ and bow
Elevation angle theta has:
6th, calculate satellite yaw-position angle over the ground.
First as shown in figure 4, establishing the measurement coordinate system o-x of digital sun sensorssysszss, wherein ozssAxis edge
The normal direction of code-disc plane, oxssAxis determines oy according to the right hand is orthogonal along the slit direction of sunlight incidence surelyssAxis, the then sun
Direction vector measurement coordinate system in coordinate be
Sss=[Sssx Sssy Sssz]T (25)
The measured value of measurement model based on sun sensor, sun sensor meets:
α is measurement coordinate system y of the solar direction vector in sun sensorssozssPlane projection and ozssThe angle of axis.
Right ascension and declination of the measurement moment sun on ecliptic can be calculated according to sun ephemeris data, it according to this can be with
Determine that coordinate of the solar vector S in the inertial coodinate system of the earth's core equator is SI, it may be assumed that
SI=[cos δs cosαs cosδs sinαs sinδs]T (27)
α in formulas,δsThe respectively right ascension and declination in solar vector direction, the earth's core equator inertial coodinate system Oe-xIyIzIOrigin
In the earth's core, xIAxis is directed toward the ascending node (being directed toward first point of Aries γ) that equatorial plane intersects nodel line with ecliptic plan, zIOneself of axis and the earth
Spin axis is parallel, yIAxis and xIAnd yIAxis forms right hand orthogonal coordinate system.
Then, orbital coordinate system can be calculated relative to the earth's core equator inertial coordinate according to the orbital tracking at measurement moment
The transition matrix of system, calculation formula are as follows:
I is orbit inclination angle in formula, and Ω is right ascension of ascending node, and u is orbital latitude argument, orbital coordinate system O-xoyozoIt is former
For point in the mass center of satellite, the orbit plane of satellite is coordinate plane, is z by the reference axis that mass center is directed toward the earth's coreoAxis, xoAxis is in-orbit
In road plane with zoVertically, the directional velocity of satellite, y are directed towardoAxis and xo,zoForm right hand orthogonal coordinate system.
Coordinate S of the solar direction vector S in orbital coordinate system can then be exported0Are as follows:
S0=CoiSI (29)
CoiSpin matrix of the expression geocentric inertial coordinate system to orbital coordinate system.
Solar vector is converted by orbital coordinate system the expression formula of the measurement coordinate system of sun sensor are as follows:
Sss=DssBCboS0 (30)
In formula: S0It is acquired by (29), CboFor by the rotation attitude matrix of orbital coordinate system to satellite body coordinate system, if
Z is pressed firstBRotation is then rotated by middle coordinate system y-axis, is finally rotated with the x-axis of new coordinate system, then rotating coordinate system CboTable
It is as follows up to formula；
The roll angle φ and pitching angle theta of formula Satellite are acquired by earth visible light sensor, and only yaw angle ψ is unknown
Amount, and DssBFor the installation matrix of sun sensor relative satellite body coordinate system, it is known quantity, there can be following expression:
S in formulaAX、SAY、SAZRespectively indicate sun sensor xssAxis and celestial body coordinate system xb、yb、zbThe included angle cosine of three axis
Value, SBX、SBY、SBZRespectively indicate sun sensor yssAxis and celestial body coordinate system xb、yb、zbThe included angle cosine value of three axis, SCX、SCY、
SCZRespectively indicate sun sensor zssAxis and celestial body coordinate system xb、yb、zbThe included angle cosine value of three axis.
Solar vector coordinate S under orbital coordinate system0It indicates to become three axis component forms are as follows:
So=[SoX SoY SoZ]T (33)
Remember again
It brings formula (31), (32) into (30), can obtain:
It brings formula (35), (36) into (26), can obtain:
A1cosψ+A2sinψ+A3=0 (37)
Ask this equation can be in the hope of the yaw angle of satellite over the ground are as follows:
A in formula1, A2, A3It can be acquired by following formula:
The earth in space is imaged in the cmos imaging device of the visible light earth sensor of satellite carrying first, then inputs
Secondly image detects gained image border and rejects Clutter edge, so that the final efficient frontier of earth image is obtained,
Furthermore the center for going out efficient frontier using least square fitting, then, according to image center location coordinate and the attitude of satellite
Angle is that low-angle assumes that lower aprons calculate satellite rolling and two axis attitude angle of pitching, finally, the measurement using sun sensor is former
Reason and the satellite being calculated roll the yaw angle that satellite is determined with pitch attitude angle information over the ground.So far, satellite is over the ground
Triaxial attitude angle it has been determined that finally achievable sustainable observation of the satellite to the earth.
The content that description in the present invention is not described in detail belongs to the well-known technique of those skilled in the art.
Claims (2)
1. the satellite attitude determination method based on visible light earth sensor and sun sensor, it is characterised in that including walking as follows
It is rapid:
(1) using the cmos imaging device in visible light earth sensor, ball is imaged over the ground, obtains input picture；Input figure
As should meet: cmos imaging device account for the 50% of total arc length at the effective arc length of earth image Horizon profile；
(2) input picture is pre-processed, removes noise, smoothed image；The method of smoothed image is 3 × 3 convolution
The median filtering of template, specific practice are as follows: by the gray value of 8 pixels in center pixel and its 3 × 3 neighborhoods by from small to large
Sequence, taking the 5th element gray value is the gray value of center pixel；
(3) edge detection is carried out to the image that step (2) obtain, obtains efficient frontier；The method of edge detection are as follows: use first
Sobel operator detects all edges, then detects efficient frontier using Hough transform, Hough transform detection is using unknown
Image radius of circle Hough transform, specific as follows:
A, first according to priori knowledge or the variation range r ∈ [r of given data estimation parameter radius rmin,rmax] and step-length
Δ r size；
H is satellite altitude,For half angle, REFor earth radius, f is camera
Focal length, r are radius approximate circle on imaging surface after the earth is imaged, and thereby determine that the range r ∈ [r of rmin,rmax]；
Step delta r is determined by the computing capability and required precision of sensor controller；
From radius r0=rminStart, recycles radius each time and increase Δ r, until r=rmax, for any circulation, radius ri=
rmin+ Δ r × i, (i=0,1 ... m), and wherein m is the total degree of circulation, expression formulaBy each edge coordinate
Point (xj,yj) (j=1 ... n) and its gradient direction angle hjSuccessively substitute into x0=x ± r cos (h) and y0=y ± r sin (h),
(x0,y0) it is central coordinate of circle, h is gradient direction angle at the edge of the circle point (x, y), and r is image radius of circle, calculates central coordinate of circle
(x0ij,y0ij), wherein j indicates effectively and invalid marginal point total number, then to calculate resulting central coordinate of circle and round up and take
Whole, it is to the right x-axis that above-mentioned each marginal point coordinate, which is using the image lower left corner as coordinate origin, upwards under the coordinate system of y-axis
It indicates；It is negative to prevent the integer central coordinate of circle being calculated, by original coordinate system respectively along x-axis, negative direction of the y-axis translation
1024 pixels obtain new positive integer central coordinate of circle value, using new coordinate value as the coordinate value of two-dimensional array, two after translation
The Data-Statistics of dimension group add up the two-dimensional coordinate value appearance number, all marginal points are traversed, finally by accumulative frequency peak value
Determine this time circulation r=riWhen Circle Parameters ci={ x0i,y0i,ri, it recycles determine one group of parameter c each timei={ x0i,y0i,
riAnd accumulative frequency Ki；Pass through this cycle accumulor number KiWith last cycle accumulor number Ki-1Comparison, record relatively high
Accumulative frequency KmaxAnd corresponding Circle Parameters cmax={ X0,Y0, R }, work as ri=rmaxWhen, final determining Circle Parameters cmax=
{X0,Y0, R } and it is efficient frontier Pixel-level central coordinate of circle { X0,Y0And radius R；
B, unknown radius Hough transform obtains Pixel-level central coordinate of circle { X0,Y0And radius R, if certain marginal point is effective edge
Edge point, i.e. boundary point on circumference, the then point and Pixel-level central coordinate of circle { X0,Y0Between distance meetThe Clutter edge point for not meeting the relational expression will be deleted, to be had
Imitate edge；
(4) efficient frontier is fitted, obtains efficient frontier center and radius；
(5) the roll angle φ obtained satellite over the ground and pitching angle theta are calculated, wherein
(Cx,Cy) it is efficient frontier centre coordinate,R is the radius of efficient frontier, and f is cmos imaging device
Focal length；
(6) it calculates and obtains the yaw angle of satellite over the groundWhereinDssBIt is defended for sun sensor is opposite
The installation matrix of star body coordinate system,S in formulaAX、SAY、SAZRespectively indicate sun sensor survey
Measure coordinate system xssAxis and satellite body coordinate system xb、yb、zbThe included angle cosine value of three axis, SBX、SBY、SBZRespectively indicate sun sensitivity
Device measures coordinate system yssAxis and satellite body coordinate system xb、yb、zbThe included angle cosine value of three axis, SCX、SCY、SCZRespectively indicate the sun
Sensor measures coordinate system zssAxis and satellite body coordinate system xb、yb、zbThe included angle cosine value of three axis；So=[SoX SoY SoZ]T
The three axis component forms for being solar vector under orbital coordinate system, α are that solar direction vector in sun sensor measures coordinate system
yssozssPlane projection and ozssThe angle of axis.
2. the satellite attitude determination method according to claim 1 based on visible light earth sensor and sun sensor,
It is characterized by: being least square method to the method that efficient frontier is fitted in the step (4).
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610262687.5A CN105928527B (en) | 2016-04-25 | 2016-04-25 | Satellite attitude determination method based on visible light earth sensor and sun sensor |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610262687.5A CN105928527B (en) | 2016-04-25 | 2016-04-25 | Satellite attitude determination method based on visible light earth sensor and sun sensor |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105928527A CN105928527A (en) | 2016-09-07 |
| CN105928527B true CN105928527B (en) | 2019-04-30 |
Family
ID=56837225
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610262687.5A Active CN105928527B (en) | 2016-04-25 | 2016-04-25 | Satellite attitude determination method based on visible light earth sensor and sun sensor |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105928527B (en) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106525054B (en) * | 2016-10-27 | 2019-04-09 | 上海航天控制技术研究所 | A kind of above pushed away using star is swept single star of remote sensing images information and independently surveys orbit determination method |
| CN108639385B (en) * | 2018-05-15 | 2021-04-13 | 浙江大学 | Method for realizing rapid and stable control of non-reference simplest attitude control system |
| CN109655080B (en) * | 2018-12-13 | 2020-05-29 | 上海航天控制技术研究所 | On-orbit calibration method for digital sun sensor |
| CN110514187B (en) * | 2019-08-30 | 2023-01-06 | 北京航空航天大学 | Astronomical north finding technical method and device for small field of view camera |
| CN110901956B (en) * | 2019-12-10 | 2021-05-07 | 中国人民解放军国防科技大学 | Satellite stable sun-to-day orientation method with earth orientation deviation as constraint |
| CN111323020B (en) * | 2020-02-25 | 2021-11-16 | 上海航天控制技术研究所 | Independent orbit determination method based on multi-vector observation of Mars edges and centers |
| CN111366149B (en) * | 2020-03-25 | 2022-02-15 | 中国科学院微小卫星创新研究院 | Method for extracting sun vector direction based on lunar phase |
| CN111605739B (en) * | 2020-05-15 | 2021-09-03 | 北京控制工程研究所 | Method for calculating attitude determination angle of rolling pitching yawing sun sensor |
| CN112632756B (en) * | 2020-12-07 | 2022-12-13 | 上海卫星工程研究所 | Sun sensor-based satellite earth shadow autonomous forecasting method and system |
| CN112977889B (en) * | 2021-03-19 | 2022-08-30 | 航天科工空间工程发展有限公司 | Satellite attitude capturing method based on sun sensor and earth sensor |
| CN114252151A (en) * | 2021-11-29 | 2022-03-29 | 中国科学院国家空间科学中心 | On-orbit working mode automatic adjustment method and device for satellite-borne optical instrument |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0589387A1 (en) * | 1992-09-21 | 1994-03-30 | Honeywell Inc. | Method and system for determining 3-axis spacecraft attitude |
| CN101236085A (en) * | 2008-01-31 | 2008-08-06 | 北京控制工程研究所 | Ultraviolet sensor,star sensor based independent rail estimation method for moon tracking |
| CN101556155A (en) * | 2009-05-20 | 2009-10-14 | 上海微小卫星工程中心 | Small satellite attitude determination system and method thereof |
| CN102591349A (en) * | 2012-03-12 | 2012-07-18 | 北京控制工程研究所 | No-gyroscope sun capture control method of high orbit satellite large initial angular rate condition |
| CN102749924A (en) * | 2012-06-18 | 2012-10-24 | 北京控制工程研究所 | Method for identifying reconfigurable weak link of satellite control system |
| CN103439974A (en) * | 2013-08-21 | 2013-12-11 | 南京航空航天大学 | Satellite posture determining system hardware reconfigurability measurement method |
-
2016
- 2016-04-25 CN CN201610262687.5A patent/CN105928527B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0589387A1 (en) * | 1992-09-21 | 1994-03-30 | Honeywell Inc. | Method and system for determining 3-axis spacecraft attitude |
| CN101236085A (en) * | 2008-01-31 | 2008-08-06 | 北京控制工程研究所 | Ultraviolet sensor,star sensor based independent rail estimation method for moon tracking |
| CN101556155A (en) * | 2009-05-20 | 2009-10-14 | 上海微小卫星工程中心 | Small satellite attitude determination system and method thereof |
| CN102591349A (en) * | 2012-03-12 | 2012-07-18 | 北京控制工程研究所 | No-gyroscope sun capture control method of high orbit satellite large initial angular rate condition |
| CN102749924A (en) * | 2012-06-18 | 2012-10-24 | 北京控制工程研究所 | Method for identifying reconfigurable weak link of satellite control system |
| CN103439974A (en) * | 2013-08-21 | 2013-12-11 | 南京航空航天大学 | Satellite posture determining system hardware reconfigurability measurement method |
Non-Patent Citations (3)
| Title |
|---|
| 三轴稳定航天器的姿态确定;北京航空航天大学;《https://wenku.baidu.com/view/9e0467662e3f5727a5e96259.html?rec_flag=default&mark_pay_doc=2&mark_rec_page=1&mark_rec_position=3&mark_rec=view_r_1&clear_uda_param=1&pn=NaN》;20131022;第1-28页 |
| 应用于微小卫星的地球敏感器设计与优化;郭振东;《中国优秀硕士学位论文全文数据库·信息科技辑》;20130615;第I135-171页 |
| 面向微小卫星的可见光地球敏感器设计;郭振东等;《传感技术学报》;20121031;第25卷(第6期);第1400-1405页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105928527A (en) | 2016-09-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105928527B (en) | Satellite attitude determination method based on visible light earth sensor and sun sensor | |
| EP3132231B1 (en) | A method and system for estimating information related to a vehicle pitch and/or roll angle | |
| CN107292927A (en) | A kind of symmetric motion platform's position and pose measuring method based on binocular vision | |
| CN106525001B (en) | It is directed toward calculation method in geostationary orbit remote sensing satellite camera optical axis space | |
| CN107677274A (en) | Unmanned plane independent landing navigation information real-time resolving method based on binocular vision | |
| CN108492333A (en) | Spacecraft attitude method of estimation based on satellite-rocket docking ring image information | |
| CN106197265A (en) | A kind of space free flight simulator precision visual localization method | |
| CN101451844B (en) | Integration star sensor based on biology compound eye structure and use thereof | |
| CN113624231B (en) | Inertial vision integrated navigation positioning method based on heterogeneous image matching and aircraft | |
| CN104864852B (en) | High resolution satellite attitude fluttering detection method based on intensive control points | |
| CN110083177A (en) | A kind of quadrotor and control method of view-based access control model landing | |
| CN107449416B (en) | Fixed star hangover asterism extracting method based on vector accumulation | |
| CN104751451B (en) | Point off density cloud extracting method based on unmanned plane low latitude high resolution image | |
| Jiang et al. | High-accuracy synchronous extraction algorithm of star and celestial body features for optical navigation sensor | |
| CN105160341B (en) | A kind of triangle detection method based on line segment detection and heart property | |
| CN109459043A (en) | A kind of spacecraft Relative Navigation based on production reconstructed image | |
| CN106289156B (en) | The method of photography point solar elevation is obtained when a kind of satellite is imaged with any attitude | |
| CN103487032A (en) | Low earth orbit space camera free-pointing image motion vector calculation method | |
| CN105928526B (en) | A method of the attitude of satellite is determined based on visible light earth sensor | |
| CN107389072B (en) | Small celestial body surface navigation feature region detection method | |
| CN117760404A (en) | All-weather autonomous navigation method, system, computer equipment and medium based on single star orientation | |
| CN113295171A (en) | Monocular vision-based attitude estimation method for rotating rigid body spacecraft | |
| CN106353032A (en) | Method for rapidly detecting celestial body centroid under light deficiency condition | |
| Kikuya et al. | Attitude determination algorithm using Earth sensor images and image recognition | |
| CN114565653B (en) | Heterologous remote sensing image matching method with rotation change and scale difference |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |