CN105300376B - A kind of the earth's core direction vector extracting method based on ultraviolet radioactive model - Google Patents
A kind of the earth's core direction vector extracting method based on ultraviolet radioactive model Download PDFInfo
- Publication number
- CN105300376B CN105300376B CN201510741026.6A CN201510741026A CN105300376B CN 105300376 B CN105300376 B CN 105300376B CN 201510741026 A CN201510741026 A CN 201510741026A CN 105300376 B CN105300376 B CN 105300376B
- Authority
- CN
- China
- Prior art keywords
- msub
- mrow
- msup
- mfrac
- earth
- 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
- 239000013598 vector Substances 0.000 title claims abstract description 30
- 238000000034 method Methods 0.000 title claims abstract description 18
- 230000002285 radioactive effect Effects 0.000 title abstract 2
- 238000000605 extraction Methods 0.000 claims abstract description 6
- 238000004364 calculation method Methods 0.000 claims description 22
- 230000005855 radiation Effects 0.000 claims description 15
- 238000012216 screening Methods 0.000 claims description 6
- 238000013461 design Methods 0.000 claims description 3
- 230000036544 posture Effects 0.000 claims description 3
- 238000005096 rolling process Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000002211 ultraviolet spectrum Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 2
- 230000000630 rising effect Effects 0.000 claims 1
- 238000001514 detection method Methods 0.000 abstract 1
- 238000003384 imaging method Methods 0.000 description 4
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
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
-
- 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/20—Instruments for performing navigational calculations
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Automation & Control Theory (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Nuclear Medicine (AREA)
- Image Analysis (AREA)
Abstract
The invention discloses a kind of the earth's core direction vector extracting methods based on ultraviolet radioactive model, include the following steps：Gather image, by using Sobel differential detections operator calculate pixel in the X direction with Grad in Y-direction, to obtain the gradient map of original image；By the local binary patterns of 5 × 5 neighborhoods of statistical pixel point in gradient map, by choosing 4 kinds of edge point feature candidate patterns in all 32 kinds of patterns of LBP features, the extraction of marginal point is realized；Construct the error function based on ellipsoidal model；By bringing marginal point into, the geocentric coordinates and radius of image plane are obtained using least-squares algorithm iterative solution, so resolve obtain the earth's core direction vector and the earth's core away from.The present invention calculates the earth's core disc centre during low rail and high rail independent navigation in by least square fitting algorithm, can obtain the geocentric position Accuracy extimate of sub-pixel.
Description
Technical Field
The invention relates to a satellite relative navigation technology, in particular to a rapid high-precision geocentric vector direction extraction method based on an ultraviolet radiation model.
Background
The autonomous navigation capability of the spacecraft is vital, the anti-interference performance and the high reliability of the navigation are vital requirements, the ultraviolet radiation characteristic of the earth has extremely stable characteristics at the edge of the earth, the spacecraft is not influenced by landform and atmosphere, obvious edge characteristics can be observed, and the autonomous navigation method based on the ultraviolet radiation model has the characteristics of high reliability, high precision and the like. Through ultraviolet imaging of the earth, the geocentric coordinates and the radius of the earth in an image plane are extracted, and a foundation is provided for subsequent autonomous navigation.
The estimation accuracy of the earth plane geocentric coordinates and the radius is a key factor influencing the autonomous navigation accuracy, and due to the factors such as different sensor track heights, different view fields, different earth ultraviolet radiation distribution at different time and the like, how to effectively process the image information to obtain the geocentric coordinates and the radius of a sub-pixel level and solve the unit geocentric vector direction is a key problem in the autonomous navigation of the ultraviolet sensor.
The application of the ultraviolet sensor in China is mainly in moon exploration engineering, imaging is carried out on Chang' e series by using the ultraviolet sensor of the moon, edge information obtained by annular 8 discontinuous sub-fields of the imaging sensor with a large field of view is effectively integrated to obtain a moon center vector and navigation information with the track height of 200km, and a related algorithm is researched on the basis of the sensors to compensate the navigation precision.
The sensor research aiming at the earth orbit is concentrated on a navigation algorithm, the navigation error depends on the measurement error of the geocentric direction, and the need of accurately estimating the geocentric direction from the earth ultraviolet characteristic image is the key point to be researched. The existing centroid determining method comprises a centroid method and a fitting method, the centroid method is adopted more, and the algorithms do not consider the computing capability and the internal storage capability of the spaceborne computer and are not suitable for on-orbit real-time computing.
Disclosure of Invention
The invention provides an edge extraction method based on gradient statistics based on an ultraviolet navigation sensor for complete once imaging of the earth, which solves the geocentric disk center through a least square fitting algorithm in the process of autonomous navigation of medium and low rails and high rails and can obtain accuracy estimation of the geocentric position of a subpixel level.
The purpose of the invention is realized by the following technical scheme: a geocentric vector direction extraction method based on an ultraviolet radiation model comprises the following steps:
step 1: collecting an earth image on a track through an ultraviolet band camera to obtain an ultraviolet spectrum image of the earth, respectively calculating pixel gradient values from the original collected image and the right image along the X direction and the Y direction, performing gradient map calculation on the image, including Sobel operator gradient calculation, and selecting edge candidate points according to the sequence of the gradient values;
step 2: carrying out LBP statistical calculation on edge-by-edge candidate points of the gradient image obtained by processing, analyzing the edge candidate point gradient statistical quantity, and extracting the obtained earth disc edge points which accord with the earth ultraviolet radiation characteristic; the method comprises the following specific steps:
s21, selecting 5 x 5 neighborhood size to calculate G obtained in step 1E(x, y) gradient statistic features of all pixels in the edge candidate point geometry, including eight unidirectional features of 10000000, 01000000, 00100000, 00010000, 00001000, 00000100, 00000010, 00000001, and 4 continuous features of 10001000, 01000100, 00100010, 00010001;
s22 screening G by the twelve characteristic quantitiesE(x, y) obtaining earth disk edge points on the image.
And step 3: constructing an error function based on an earth disk ellipsoid model, substituting the coordinates of the edge points obtained in the step (2) in a detector coordinate system into the function, performing least square iterative fitting, and solving to obtain sub-pixel-level geocentric vectors;
wherein the sub-pixel level geocentric vector is obtained by the following steps:
when the orbit is a stationary orbit of the earth:
A. when the geocentric vector is estimated on the earth stationary orbit, because the proportion between the orbit height and the change of the earth radius can be ignored, an error function of a sphere model without considering the earth oblateness is established, and the three-axis postures of the ultraviolet band camera are respectively considered as follows: roll angle phiYaw angle phi, pitch angle theta, error functionComprises the following steps:
in the formula, i represents the number of edge points, and the value range is [1, n ]]，xiAnd yiIs the coordinate of the edge point in the sensor coordinate system, xcAnd ycIs the coordinate of the earth's center under the sensor coordinate system, ziIs the actual focal length value, ReRepresenting the size of the earth radius in a sensor coordinate system, wherein a superscript 2 represents an arithmetic square, and a symbol sigma represents summation;
wherein, the coefficients A, B, C, D, E, F are respectively:
A＝-k2sin2θ+cos2θ，
wherein k is the radius of the view angle, and the calculation formula is as follows:
wherein R is the earth radius, R is the geocentric distance, the symbol √ represents the arithmetic root, and the symbol tan represents the tangent trigonometric operation;
when the orbit altitude belongs to the low-medium earth orbit:
B. on the medium and low orbit, an error function of an ellipsoid model considering the earth oblateness is established, and the error function is establishedComprises the following steps:
in the range of [1, n]，xiAnd yiIs the coordinate of the edge point in the sensor coordinate system, xcAnd ycThe earth center is under the sensor coordinate system
Wherein i represents the number of edge points, the coordinate of the value, ziIs the actual focal length value, ReRepresenting the size of the earth radius in a sensor coordinate system, wherein a superscript 2 represents an arithmetic square, and a symbol sigma represents summation;
wherein, the coefficients A, B, C, D, E, F are respectively:
wherein e is1And e2Respectively longitude and latitude flat ratios under the geocentric inertial coordinate system, the superscript 2 represents the arithmetic square, lijThe calculation formula of (2) is as follows:
in the formula, phi and theta are respectively a rolling angle, a yaw angle and a pitch angle of a three-axis attitude angle between a sensor coordinate system and a track coordinate system, omega, I and U are respectively a rising-crossing declination, a track inclination angle and a true-near-point angle between the track coordinate system and a geocentric inertial coordinate system, cos and sin respectively represent trigonometric cosine and sine function operation, and symbol x represents an arithmetic product;
introducing the edges obtained in step 3 into an error functionOrIterative solution is carried out by adopting least square fitting, and when the threshold value of minimum convergence of an error function is 2.2204 multiplied by 10-16When not higher than the threshold, the optimal solution x is obtainedc、ycAnd ReConstructing a geocentric vector;
the formula for calculating the geocentric unit vector R is as follows:
in the formula, xc、ycFor the optimal solution of the geocentric coordinate obtained by calculation, f is the focal length of the sensor, the symbol x represents the arithmetic product, and the symbol V represents the arithmetic evolution.
The formula for calculating the distance between centers of earth is as follows:
in the formula, RearthIs the radius of the earth, ReFor the calculated radius of the earth on the image plane, IFOV is the instantaneous field of view for the sensor design, the symbol sin represents the sine trigonometric operation, and the symbol x represents the arithmetic product.
Wherein, the gradient calculation formula of the Sobel operator in the step 1 is as follows:
wherein G (x, y) represents a gradient value of a pixel having (x, y) as an image coordinate, and GX(X, y) represents the gradient component in the X direction of the point, GY(x, Y) represents the Y-direction gradient component of the point, x, Y are the pixel coordinates in the image coordinate system, and the symbol √ represents the arithmetic power.
The specific steps of selecting the edge candidate points according to the gradient value sorting in the step 1 are as follows:
1) arranging the gradient characteristic values of pixel points on the gradient body according to the size, and taking 0.75-0.85 times of the maximum value as a threshold value;
2) calculating candidate edge points by the following formula:
GE(x,y)≥(0.75～0.85)×Gmax(x,y)
in the formula, Gmax(x, y) is the maximum of the gradient map, GEAnd (x, y) is the candidate edge point set obtained by screening in the step 1, and the symbol x represents a product.
Wherein, the gradient statistic LBP in step 2 is calculated by the following formula:
in the formula,gpis the neighborhood pixel gray value, gcIs the central pixel gray value, (x)c,yc) As the center pixel coordinate, P is the template sampling number, R is the template pixel radius, gl is the difference between the two pixel gray values, the sign Σ represents the summation, and R takes the value of 5 pixels.
Compared with the prior art, the invention has the following advantages:
the geocentric disk center is solved by a least square fitting algorithm in the process of middle-low orbit and high orbit autonomous navigation, and the accuracy estimation of the geocentric position at the sub-pixel level can be obtained.
Drawings
Fig. 1 is a schematic flow chart of a method for rapidly extracting a geocentric vector direction with high precision based on an ultraviolet radiation model according to an embodiment of the present invention.
Detailed Description
The present invention will be described in detail with reference to specific examples. The following examples will assist those skilled in the art in further understanding the invention, but are not intended to limit the invention in any way. It should be noted that variations and modifications can be made by persons skilled in the art without departing from the spirit of the invention. All falling within the scope of the present invention.
As shown in fig. 1, an embodiment of the present invention provides a method for extracting a centroid vector direction quickly and accurately based on an ultraviolet radiation model, including the following steps:
step 1:
the method comprises the following steps that an ultraviolet band camera collects an earth image on a track to obtain an ultraviolet spectrum image of the earth, pixel gradient values are respectively calculated from an original collected image and a right image along the X direction and the Y direction, gradient map calculation is carried out on the image, the gradient calculation comprises Sobel operator gradient calculation, and the calculation formula is as follows:
wherein G (x, y) represents a gradient value of a pixel having (x, y) as an image coordinate, GX(X, y) represents the gradient component in the X direction of the point, GY(x, Y) represents the Y-direction gradient component of the point, x, Y are the pixel coordinates in the image coordinate system, and the symbol √ represents the arithmetic power.
Arranging the gradient characteristic values of pixel points on the gradient body according to the size, taking 0.75-0.85 times of the maximum value as a threshold value, and the candidate edge point calculation method comprises the following steps:
GE(x,y)≥(0.75～0.85)×Gmax(x,y)
wherein G ismax(x, y) is the maximum of the gradient map, GEAnd (x, y) is the candidate edge point set obtained by screening in the step 1, and the symbol x represents a product.
Illustrated in table 1:
step 2: calculating gradient statistics of candidate feature points extracted from a gradient map of an original image, wherein the LBP has the following calculation formula:
wherein,gpis the neighborhood pixel gray value, gcIs the central pixel gray value, (x)c,yc) As the center pixel coordinate, P is the template sampling number, R is the template pixel radius, gl is the difference between the two pixel gray values, the sign Σ represents the summation, and R takes the value of 5 pixels.
The specific steps of the step 2 are as follows:
selecting 5 x 5 neighborhood size to calculate G obtained in step 1E(x, y) gradient statistics for all pixels in the edge candidate point geometry, the gradient statistics comprising: eight unidirectional characteristic quantities of 10000000, 01000000, 00100000, 00010000, 00001000, 00000100, 00000010 and 00000001, and 4 continuous characteristic quantities of 10001000, 01000100, 00100010 and 00010001;
screening G by the twelve characteristic quantitiesE(x, y) obtaining earth disk edge points on the image.
And step 3: and (3) constructing an error function based on the earth disk ellipsoid model, substituting the coordinates of the edge points obtained in the step (2) in a detector coordinate system into the function, performing least square iterative fitting, and solving to obtain sub-pixel-level geocentric vectors.
The sub-pixel level geocentric vector is calculated by the following steps:
when the orbit is a stationary orbit of the earth:
A. when estimating geocentric vectors on stationary orbits of the earth, the orbit altitude and the earthThe proportion between the radius changes can be ignored, an error function of a spherical ball model without considering the oblateness of the earth is established, and the three-axis postures of the ultraviolet band camera are respectively considered as follows: roll angle phi, yaw angle phi, pitch angle theta, error functionComprises the following steps:
wherein i represents the number of edge points, and the value range is [1, n ]]，xiAnd yiIs the coordinate of the edge point in the sensor coordinate system, xcAnd ycIs the coordinate of the earth's center under the sensor coordinate system, ziIs the actual focal length value, ReThe radius of the earth is shown in the sensor coordinate system, the superscript 2 is the arithmetic square, and the symbol Σ is the sum.
Wherein, the coefficients A, B, C, D, E, F are respectively:
A＝-k2sin2θ+cos2θ，
wherein k is the radius of the view angle, and the calculation formula is as follows:
where R is the earth radius, R is the geocentric distance, the symbol √ represents the arithmetic root, and the symbol tan represents the tangent-trigonometric operation.
When the orbit altitude belongs to the low-medium earth orbit:
B. on the medium and low orbit, an error function of an ellipsoid model considering the earth oblateness is established, and the error function is establishedComprises the following steps:
wherein i represents the number of edge points, and the value range is [1, n ]]，xiAnd yiIs the coordinate of the edge point in the sensor coordinate system, xcAnd ycIs the coordinate of the earth's center under the sensor coordinate system, ziIs the actual focal length value, ReThe radius of the earth is shown in the sensor coordinate system, the superscript 2 is the arithmetic square, and the symbol Σ is the sum.
Wherein, the coefficients A, B, C, D, E, F are respectively:
wherein e is1And e2Respectively longitude and latitude flat ratios under the geocentric inertial coordinate system, the superscript 2 represents the arithmetic square, lijThe calculation formula of (2) is as follows:
phi, phi and theta are respectively a rolling angle, a yaw angle and a pitch angle of a three-axis attitude angle between a sensor coordinate system and a track coordinate system, omega, I and U are respectively a rising-crossing declination, a track inclination angle and a true-near-point angle between the track coordinate system and a geocentric inertial coordinate system, cos and sin respectively represent trigonometric cosine and sine function operation, and symbol x represents an arithmetic product.
Introducing the edges obtained in step 3 into an error functionOrIterative solution is carried out by adopting least square fitting, and when the threshold value of minimum convergence of an error function is 2.2204 multiplied by 10-16When not higher than the threshold, the optimal solution x is obtainedc、ycAnd ReAnd constructing the geocentric vector.
The formula for calculating the geocentric unit vector R is as follows:
wherein x isc、ycFor the optimal solution of the geocentric coordinate obtained by calculation, f is the focal length of the sensor, the symbol x represents the arithmetic product, and the symbol V represents the arithmetic evolution.
The formula for calculating the distance between centers of earth is as follows:
wherein R isearthIs the radius of the earth, ReFor the calculated radius of the earth on the image plane, IFOV is the instantaneous field of view for the sensor design, the symbol sin represents the sine trigonometric operation, and the symbol x represents the arithmetic product.
The foregoing description of specific embodiments of the present invention has been presented. It is to be understood that the present invention is not limited to the specific embodiments described above, and that various changes and modifications may be made by one skilled in the art within the scope of the appended claims without departing from the spirit of the invention.
Claims (5)
1. A geocentric vector direction extraction method based on an ultraviolet radiation model is characterized by comprising the following steps:
step 1: collecting an earth image on a track through an ultraviolet band camera to obtain an ultraviolet spectrum image of the earth, respectively calculating pixel gradient values from the original collected image and the right image along the X direction and the Y direction, performing gradient map calculation on the image, including Sobel operator gradient calculation, and selecting edge candidate points according to the sequence of the gradient values;
step 2: carrying out LBP statistical calculation on edge-by-edge candidate points of the gradient image obtained by processing, analyzing the edge candidate point gradient statistical quantity, and extracting the obtained earth disc edge points which accord with the earth ultraviolet radiation characteristic;
and step 3: constructing an error function based on an earth disk ellipsoid model, substituting the coordinates of the edge points obtained in the step (2) in a detector coordinate system into the function, performing least square iterative fitting, and solving to obtain sub-pixel-level geocentric vectors;
the sub-pixel level geocentric vector is obtained by the following steps:
when the orbit is a stationary orbit of the earth:
A. when the geocentric vector is estimated on the earth stationary orbit, because the proportion between the orbit height and the change of the earth radius can be ignored, an error function of a sphere model without considering the earth oblateness is established, and the three-axis postures of the ultraviolet band camera are respectively considered as follows: roll angleYaw angle phi, pitch angle theta, error functionComprises the following steps:
<mrow> <msubsup> <mi>E</mi> <mi>s</mi> <mi>h</mi> </msubsup> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>A</mi> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <mi>B</mi> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msubsup> <mi>Cz</mi> <mi>i</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mn>2</mn> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mn>2</mn> <mi>E</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>+</mo> <mn>2</mn> <mi>F</mi> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>R</mi> <mi>e</mi> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mn>2</mn> </msup> </mrow>
in the formula, i represents the number of edge points, and the value range is [1, n ]]，xiAnd yiIs the coordinate of the edge point in the sensor coordinate system, xcAnd ycIs the coordinate of the earth's center under the sensor coordinate system, ziIs the actual focal length value, ReRepresenting the size of the earth radius in a sensor coordinate system, wherein a superscript 2 represents an arithmetic square, and a symbol sigma represents summation;
wherein, the coefficients A, B, C, D, E, F are respectively:
A＝-k2sin2θ+cos2θ，
wherein k is the radius of the view angle, and the calculation formula is as follows:
<mrow> <mi>k</mi> <mo>=</mo> <mi>t</mi> <mi>a</mi> <mi>n</mi> <mrow> <mo>(</mo> <mi>R</mi> <mo>/</mo> <msqrt> <mrow> <msup> <mi>r</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>R</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
wherein R is the earth radius, R is the geocentric distance, the symbol √ represents the arithmetic root, and the symbol tan represents the tangent trigonometric operation;
when the orbit altitude belongs to the low-medium earth orbit:
B. on the medium and low orbit, an error function of an ellipsoid model considering the earth oblateness is established, and the error function is established
<mrow> <msubsup> <mi>E</mi> <mi>s</mi> <mrow> <mi>l</mi> <mo>.</mo> <mi>m</mi> </mrow> </msubsup> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>A</mi> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <mi>B</mi> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msubsup> <mi>Cz</mi> <mi>i</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mn>2</mn> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mn>2</mn> <mi>E</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>+</mo> <mn>2</mn> <mi>F</mi> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>R</mi> <mi>e</mi> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mn>2</mn> </msup> <mo>;</mo> </mrow>
In the formula, i represents the number of edge points, and the value range is [1, n ]]，ziIs the actual focal length value, ReRepresenting the size of the earth radius in the sensor coordinate system, superscript 2 representing the arithmetic square, symbol Σ representing the sum, xiAnd yiIs the coordinate of the edge point in the sensor coordinate system, xcAnd ycIs the coordinate of the geocenter under the sensor coordinate system;
wherein, the coefficients A, B, C, D, E, F are respectively:
<mrow> <mi>A</mi> <mo>=</mo> <msubsup> <mi>l</mi> <mn>11</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mfrac> <msubsup> <mi>l</mi> <mn>21</mn> <mn>2</mn> </msubsup> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>+</mo> <mfrac> <msubsup> <mi>l</mi> <mn>31</mn> <mn>2</mn> </msubsup> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow>
<mrow> <mi>B</mi> <mo>=</mo> <msubsup> <mi>l</mi> <mn>12</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mfrac> <msubsup> <mi>l</mi> <mn>22</mn> <mn>2</mn> </msubsup> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>+</mo> <mfrac> <msubsup> <mi>l</mi> <mn>32</mn> <mn>2</mn> </msubsup> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow>
<mrow> <mi>C</mi> <mo>=</mo> <msubsup> <mi>l</mi> <mn>13</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mfrac> <msubsup> <mi>l</mi> <mn>23</mn> <mn>2</mn> </msubsup> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>+</mo> <mfrac> <msubsup> <mi>l</mi> <mn>33</mn> <mn>2</mn> </msubsup> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow>
<mrow> <mi>D</mi> <mo>=</mo> <msub> <mi>l</mi> <mn>11</mn> </msub> <msub> <mi>l</mi> <mn>12</mn> </msub> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <msub> <mi>l</mi> <mn>21</mn> </msub> <msub> <mi>l</mi> <mn>22</mn> </msub> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <msub> <mi>l</mi> <mn>31</mn> </msub> <msub> <mi>l</mi> <mn>32</mn> </msub> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow>
<mrow> <mi>E</mi> <mo>=</mo> <msub> <mi>l</mi> <mn>11</mn> </msub> <msub> <mi>l</mi> <mn>13</mn> </msub> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <msub> <mi>l</mi> <mn>21</mn> </msub> <msub> <mi>l</mi> <mn>23</mn> </msub> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <msub> <mi>l</mi> <mn>31</mn> </msub> <msub> <mi>l</mi> <mn>33</mn> </msub> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow>
<mrow> <mi>F</mi> <mo>=</mo> <msub> <mi>l</mi> <mn>12</mn> </msub> <msub> <mi>l</mi> <mn>13</mn> </msub> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <msub> <mi>l</mi> <mn>22</mn> </msub> <msub> <mi>l</mi> <mn>23</mn> </msub> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <msub> <mi>l</mi> <mn>32</mn> </msub> <msub> <mi>l</mi> <mn>33</mn> </msub> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>e</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>;</mo> </mrow>
wherein e is1And e2Respectively longitude and latitude flat ratios under the geocentric inertial coordinate system, the superscript 2 represents the arithmetic square, lijThe calculation formula of (2) is as follows:
in the formula,phi and theta are respectively a rolling angle, a yaw angle and a pitch angle of a three-axis attitude angle between a sensor coordinate system and a track coordinate system, omega, I and U are respectively a rising intersection declination, a track inclination angle and a true anomaly angle between the track coordinate system and a geocentric inertia coordinate system, cos and sin respectively represent trigonometric cosine and sine function operation, and symbol x represents an arithmetic product;
introducing the edges obtained in step 3 into an error functionOrIterative solution is carried out by adopting least square fitting, and when the threshold value of minimum convergence of an error function is 2.2204 multiplied by 10-16When not higher than the threshold, the optimal solution x is obtainedc、ycAnd ReConstructing a geocentric vector;
the formula for calculating the geocentric unit vector R is as follows:
<mrow> <mi>R</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mfrac> <mrow> <mo>-</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> </mrow> <msqrt> <mrow> <msup> <msub> <mi>x</mi> <mi>c</mi> </msub> <mn>2</mn> </msup> <mo>+</mo> <msup> <msub> <mi>y</mi> <mi>c</mi> </msub> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>f</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mtd> </mtr> <mtr> <mtd> <mfrac> <mrow> <mo>-</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> </mrow> <msqrt> <mrow> <msup> <msub> <mi>x</mi> <mi>c</mi> </msub> <mn>2</mn> </msup> <mo>+</mo> <msup> <msub> <mi>y</mi> <mi>c</mi> </msub> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>f</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mtd> </mtr> <mtr> <mtd> <mfrac> <mrow> <mo>-</mo> <mi>f</mi> </mrow> <msqrt> <mrow> <msup> <msub> <mi>x</mi> <mi>c</mi> </msub> <mn>2</mn> </msup> <mo>+</mo> <msup> <msub> <mi>y</mi> <mi>c</mi> </msub> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>f</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow>
in the formula, xc、ycF is the focal length of the sensor, the symbol x represents the arithmetic product, and the symbol x represents the optimal solution of the geocentric coordinateRepresenting an arithmetic evolution;
the formula for calculating the distance between centers of earth is as follows:
<mrow> <mi>r</mi> <mo>=</mo> <mfrac> <msub> <mi>R</mi> <mrow> <mi>e</mi> <mi>a</mi> <mi>r</mi> <mi>t</mi> <mi>h</mi> </mrow> </msub> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>R</mi> <mi>e</mi> </msub> <mo>&times;</mo> <mi>I</mi> <mi>F</mi> <mi>O</mi> <mi>V</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
in the formula, RearthIs the radius of the earth, ReFor the calculated radius of the earth on the image plane, IFOV is the instantaneous field of view for the sensor design, the symbol sin represents the sine trigonometric operation, and the symbol x represents the arithmetic product.
2. The method for extracting the geocentric vector direction based on the ultraviolet radiation model according to claim 1, wherein the gradient calculation formula of the Sobel operator in the step 1 is as follows:
<mrow> <mi>G</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>X</mi> </msub> <mo>(</mo> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>Y</mi> </msub> <mo>(</mo> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mrow>
wherein G (x, y) represents a gradient value of a pixel having (x, y) as an image coordinate, and GX(X, y) represents the gradient component in the X direction of the point, GY(x, Y) represents the Y-direction gradient component of the point, x, Y are the pixel coordinates in the image coordinate system, and the symbol √ represents the arithmetic power.
3. The method for extracting geocentric vector direction based on ultraviolet radiation model as claimed in claim 1, wherein the specific step of selecting edge candidate points according to the gradient value sorting in step 1 is:
1) arranging the gradient characteristic values of pixel points on the gradient body according to the size, and taking 0.75-0.85 times of the maximum value as a threshold value;
2) calculating candidate edge points by the following formula:
GE(x,y)≥(0.75～0.85)×Gmax(x,y)
in the formula, Gmax(x, y) is the maximum of the gradient map, GEAnd (x, y) is the candidate edge point set obtained by screening in the step 1, and the symbol x represents a product.
4. The method for extracting geocentric vector direction based on ultraviolet radiation model of claim 1, wherein the gradient statistic LBP in step 2 is calculated by the formula:
<mrow> <msub> <mi>LBP</mi> <mrow> <mi>P</mi> <mo>,</mo> <mi>R</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>P</mi> </munderover> <mi>s</mi> <mrow> <mo>(</mo> <msub> <mi>g</mi> <mi>p</mi> </msub> <mo>-</mo> <msub> <mi>g</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow>
in the formula,gpis a neighborhood pixel grayValue of gcIs the central pixel gray value, (x)c,yc) As the center pixel coordinate, P is the template sampling number, R is the template pixel radius, gl is the difference between the two pixel gray values, the sign Σ represents the summation, and R takes the value of 5 pixels.
5. The method for extracting geocentric vector direction based on the ultraviolet radiation model according to claim 3, wherein the specific steps of the step 2 are as follows:
s21, selecting 5 x 5 neighborhood size to calculate G obtained in step 1E(x, y) gradient statistic features of all pixels in the edge candidate point geometry, including eight unidirectional features of 10000000, 01000000, 00100000, 00010000, 00001000, 00000100, 00000010, 00000001, and 4 continuous features of 10001000, 01000100, 00100010, 00010001;
s22 screening G by the twelve characteristic quantitiesE(x, y) obtaining earth disk edge points on the image.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510741026.6A CN105300376B (en) | 2015-11-04 | 2015-11-04 | A kind of the earth's core direction vector extracting method based on ultraviolet radioactive model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510741026.6A CN105300376B (en) | 2015-11-04 | 2015-11-04 | A kind of the earth's core direction vector extracting method based on ultraviolet radioactive model |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105300376A CN105300376A (en) | 2016-02-03 |
| CN105300376B true CN105300376B (en) | 2018-05-29 |
Family
ID=55197913
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510741026.6A Active CN105300376B (en) | 2015-11-04 | 2015-11-04 | A kind of the earth's core direction vector extracting method based on ultraviolet radioactive model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105300376B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111695564B (en) * | 2020-06-12 | 2023-11-14 | 上海航天控制技术研究所 | Target identification and navigation method |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0589387B1 (en) * | 1992-09-21 | 1997-04-02 | Honeywell Inc. | Method and system for determining 3-axis spacecraft attitude |
| CN101236092A (en) * | 2008-01-31 | 2008-08-06 | 北京控制工程研究所 | Ultraviolet navigation sensor |
| CN101713655A (en) * | 2009-11-16 | 2010-05-26 | 北京控制工程研究所 | Dual-mode control method of UV navigation sensor |
-
2015
- 2015-11-04 CN CN201510741026.6A patent/CN105300376B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0589387B1 (en) * | 1992-09-21 | 1997-04-02 | Honeywell Inc. | Method and system for determining 3-axis spacecraft attitude |
| CN101236092A (en) * | 2008-01-31 | 2008-08-06 | 北京控制工程研究所 | Ultraviolet navigation sensor |
| CN101713655A (en) * | 2009-11-16 | 2010-05-26 | 北京控制工程研究所 | Dual-mode control method of UV navigation sensor |
Non-Patent Citations (4)
| Title |
|---|
| Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns;Timo Ojala etc.;《IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE》;20020731;第24卷(第7期);第971-987页 * |
| Texture discrimination with multidimensional distributions of signed gray-level differences;Timo Ojala etc.;《Pattern Recognition》;20010331;第34卷(第3期);第727-739页 * |
| 亚像素精度的行星中心定位算法;陈阔等;《光学精密工程》;20130731;第21卷(第7期);第1881-1890页 * |
| 基于圆球与椭球的地球敏感器地心矢量算法分析;刘军等;《航天控制》;20080430;第26卷(第2期);第86-91页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105300376A (en) | 2016-02-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN101839721B (en) | Visual navigation method in autonomous rendezvous and docking | |
| CN102607526B (en) | Target posture measuring method based on binocular vision under double mediums | |
| CN103236064B (en) | A kind of some cloud autoegistration method based on normal vector | |
| CN107451593B (en) | High-precision GPS positioning method based on image feature points | |
| CN103697855B (en) | A kind of hull horizontal attitude measuring method detected based on sea horizon | |
| CN106056605B (en) | A kind of in-orbit high precision image localization method based on images match | |
| CN106780729A (en) | A kind of unmanned plane sequential images batch processing three-dimensional rebuilding method | |
| Salamuniccar et al. | Method for crater detection from martian digital topography data using gradient value/orientation, morphometry, vote analysis, slip tuning, and calibration | |
| CN106556412A (en) | The RGB D visual odometry methods of surface constraints are considered under a kind of indoor environment | |
| CN103632366A (en) | Parameter identification method for elliptical target | |
| CN102116626B (en) | Prediction and correction method of node of star point track image | |
| CN106529587A (en) | Visual course identification method based on target point identification | |
| CN102865859A (en) | Aviation sequence image position estimating method based on SURF (Speeded Up Robust Features) | |
| CN104749598B (en) | A kind of method in generation GNSS occultation path | |
| CN102853835A (en) | Scale invariant feature transform-based unmanned aerial vehicle scene matching positioning method | |
| CN103644918A (en) | Method for performing positioning processing on lunar exploration data by satellite | |
| WO2022104260A1 (en) | Data normalization of aerial images | |
| Zhang et al. | High-accuracy location algorithm of planetary centers for spacecraft autonomous optical navigation | |
| CN107097256A (en) | Model-free method for tracking target of the view-based access control model nonholonomic mobile robot under polar coordinates | |
| CN110927765A (en) | Laser radar and satellite navigation fused target online positioning method | |
| CN105300376B (en) | A kind of the earth's core direction vector extracting method based on ultraviolet radioactive model | |
| CN113436313B (en) | Three-dimensional reconstruction error active correction method based on unmanned aerial vehicle | |
| CN104567879A (en) | Method for extracting geocentric direction of combined view field navigation sensor | |
| CN108921896B (en) | Downward vision compass integrating dotted line characteristics | |
| CN109946682A (en) | GF3 data baseline estimation method based on ICESat/GLAS |
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 |