Atmospheric refraction compensation method and system for optical earth remote sensing satellite image geographic positioning

Abstract

The invention provides an optical earth remote sensing satellite image geographic positioning atmospheric refraction compensation method and system. The method comprises the following steps: acquiring the position of a satellite in an earth fixed connection coordinate system at a corresponding exposure moment and a sight vector corresponding to each detector; establishing a multi-layer atmospheric refraction model, and layering the atmosphere according to the elevation; calculating the average atmospheric refractive index of each layer; sequentially calculating from the 0th layer; and calculating the sight line intersection point position according to the incidence starting point, the incidence sight line vector and the elevation. If the intersection elevation is smaller than the maximum value of the regional elevation, judging whether the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation; if yes, calculating an intersection point of the incident sight line vector and the digital elevation curved surface, and outputting coordinates of the intersection point; and if not, calculating the sight vector after refraction according to the refractive index of the incident layer and the refractive index of the emergent layer, and adding 1 to the layer number for repeating. The method is reasonable, simple in calculation and easy to implement, and can be effectively applied to atmospheric refraction compensation of optical earth remote sensing images.

Description

Technical Field

[0001] The present invention relates to the field of geographic positioning technology, and in particular to an atmospheric refraction compensation method and system for optical earth remote sensing satellite image geographic positioning. Background Art

[0002] An important part of optical earth remote sensing applications is the geolocation of remote sensing effects. On the basis of a rigorous geometric imaging model, traditional methods use collinear equations for mapping calculations. With the increase in the requirements for geolocation accuracy, more and more non-ideal factors are taken into account and targeted corrections are made, such as the geometric distortion of optical lenses, the ephemeris measurement accuracy of satellite platforms, the attitude measurement accuracy, and the influence of thermal deformation of satellite platforms. This year, with the application of 10-meter, meter, and sub-meter optical earth remote sensing satellites, the demand for geolocation accuracy has also been further improved.

[0003] Atmospheric refraction will change the straight-line propagation model of light, causing the collinear model to be non-rigorous. According to existing literature, when observing a low-orbit satellite at a zenith angle of 55°, atmospheric refraction can cause a positioning error of the order of ten meters. Atmospheric refraction affects optical remote sensing of the earth, and will also affect ground-based laser strikes, space-based laser altimetry, ground-based measurement and control equipment, and ground-based astronomical telescope observations. Literature 1 (Xi Hua, et al., Analysis of Influencing Factors of Atmospheric Refraction Correction, Value Engineering, 2012, No. 22, 313-314) analyzes the effect of atmospheric dispersion effects on the deviation of the laser transmission path in the atmosphere from the straight line, resulting in inconsistent aiming and strike paths, and focuses on the influence of different latitudes and carbon dioxide concentrations on atmospheric refraction elevation angle errors. Reference 2 (Wang Taiping, Research on Waveform Processing and Positioning Accuracy Optimization Technology of GF-7 Satellite Laser Altimeter Data, Wuhan University, Master's Degree Thesis, 2021) describes that for satellite-borne laser altimeters, the ranging delay caused by the atmospheric refraction effect leads to the deviation of the laser foot position, and gives the atmospheric dry delay correction process for satellite-borne laser altimeters. Reference 3 (Wu Pengfei, Low-elevation Atmospheric Refraction Correction Method Based on Gridded Atmospheric Parameter Profile Model, Acta Optica Sinica, Vol. 37, No. 6, 2017, 0601004) For the low-elevation angle measurement and control tracking below 5°, a gridded atmospheric parameter profile model is established, and the refraction correction results under typical atmospheric conditions are analyzed. Reference 4 (Liang Xiaobo, Research on Telescope Pointing Correction Based on Astronomical Orientation Technology, Doctoral Dissertation of University of Chinese Academy of Sciences, 2020) describes the influence of atmospheric refraction error (gas error) on ground-based astronomical telescopes, and uses the atmospheric refraction model in IAU SOFA to correct the gas error.

[0004] In the field of high-resolution optical earth remote sensing, Reference 5 (Cheng Yufeng, Research on high-precision on-orbit autonomous geometric calibration methods for high-resolution optical remote sensing satellites, Ph.D. dissertation of Wuhan University, 2019) analyzes the impact of atmospheric refraction on on-orbit geometric calibration of optical satellites, and proposes a geometric calibration method that takes into account the atmosphere using the RPC model. Reference 6 (Zhang Lin, Research on key technologies for high-precision automatic positioning of optical satellite remote sensing data, Master's thesis of Tsinghua University, 2009) analyzes the impact of atmospheric refraction, and based on SPOT-5 experimental data, selects control points to fit the radius that needs to be compensated, and then uses this as a feature library to compensate for atmospheric refraction on another SPOT-5 remote sensing data. Reference 7 (Man Yiyun, et al., Analysis of Planar Positioning Accuracy Error of Optical Remote Sensing Satellite, "Space Return and Remote Sensing", Vol. 42, No. 1, 2021, pp. 135-144) focuses on the influence of optical path difference and atmospheric refraction offset on the planar positioning of optical remote sensing satellites. Through the conversion of physical geometry model and general geometry RPC model, the RPC model coefficients are obtained based on the virtual control point fitting solution, and the positioning accuracy is compared and verified. Reference 8 (Wang Yanli, Research on Atmospheric Refraction Error Analysis of Full-spectrum Multi-angle Optical Imaging, Journal of Wuhan University (Information Science Edition), 2022, https: / / doi.org / 10.13203 / j.whugis20210080) Based on the atmospheric refraction error characteristics of full-spectrum multi-angle imaging, the camera parameters and atmospheric parameters of the GF-5 full-spectrum spectral imager are used, and the atmospheric refraction correction model of optical satellite based on ellipsoid layering is used to systematically analyze the differences in atmospheric refraction errors of full-spectrum multi-angle imaging in different regions and at different times. Reference 9 (Zang Wenchi, Research on the Main Errors in Line of Sight Determination of Medium and Low Orbit Optical Satellites, Master’s Thesis of National University of Defense Technology, 2018) analyzes the specific effects of factors such as line of sight angle of medium and low orbit satellites, atmospheric stratification thickness, target altitude, and satellite orbit altitude on the atmospheric refraction line of sight error, and proposes an atmospheric refraction correction method for binary satellite target positioning scenarios.

[0005] Atmospheric refraction compensation method and system in optical satellite remote sensing data geolocation (CN102346252A) discloses an atmospheric refraction compensation method and system in optical satellite remote sensing data geolocation, the atmospheric refraction compensation method in optical satellite remote sensing data geolocation includes an earth radius compensation step, which is used to calculate the coordinates of the ground point corresponding to a certain pixel point in the satellite image by assuming that the actual radius of the earth plus the earth radius compensation amount is the earth radius. A starlight refraction satellite autonomous positioning method for all-day (CN104236553A) discloses a starlight refraction correction method for satellite observation of stars for all-day high-precision astronomical autonomous navigation of satellites. A target positioning method under atmospheric refraction applicable to airborne optoelectronic sighting system (CN108535715A) discloses a target positioning method under atmospheric refraction applicable to airborne optoelectronic sighting system. Summary of the invention

[0006] In view of the defects in the prior art, the object of the present invention is to provide a method and system for atmospheric refraction compensation for optical earth remote sensing satellite image geolocation.

[0007] According to the present invention, a method for compensating for atmospheric refraction in geo-positioning of optical earth remote sensing satellite images comprises:

[0008] Step S1: Obtain the position of the satellite in the earth-fixed coordinate system at the corresponding exposure time and the line of sight vector corresponding to each detector;

[0009] Step S2: Establish a multi-layer atmospheric refraction model and stratify the atmosphere according to altitude;

[0010] Step S3: Calculate the average atmospheric refractive index of each layer;

[0011] Step S4: Calculate sequentially starting from layer 0;

[0012] Step S5: Calculate the line of sight intersection position based on the incident starting point, the incident line of sight vector, and the elevation.

[0013] Step S6: If the elevation of the intersection point is less than the maximum elevation of the area, determine whether the digital elevation corresponding to the intersection point position is greater than or equal to the elevation of the intersection point;

[0014] If the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation, calculate the intersection of the incident sight vector and the digital elevation surface, and output the intersection coordinates;

[0015] If the digital elevation corresponding to the intersection position is less than the intersection elevation, the line of sight vector after refraction is calculated according to the refractive index of the incident layer and the refractive index of the exit layer, and the layer number is increased by 1 and steps S5 to S6 are repeated.

[0016] Preferably, in step S1:

[0017] Get the satellite position in the Earth-fixed coordinate system ECEF at the corresponding exposure time The line of sight of each detector

[0018] Vector is a unit vector, where the subscript j represents the detector number;

[0019] The line of sight vector corresponding to each detector is accurately pointed after being corrected for optical distortion and on-board thermal deformation.

[0020] Preferably, in step S2:

[0021] A multi-layer atmospheric refraction model is established, and the atmosphere is divided into M layers according to the altitude. The altitudes corresponding to the tops of the layers from top to bottom are h0, h1, ..., h M-1 ;

[0022] The number of atmospheric stratification layers is set according to the angle of the satellite observation deviating from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation deviating from the sub-satellite point is smaller, the number of layers is reduced to reduce the amount of calculation.

[0023] Preferably, in step S3:

[0024] Calculate the average atmospheric refractive index n for each layer k , where k represents the layer number, k = 0, 1, 2, ..., M; n0 is the refractive index under vacuum, n0 = 1;

[0025] The average atmospheric refractive index n corresponding to the kth layer k The calculation method is:

[0026]

[0027] Where k = 1, 2, ..., M, function f(h) represents the function of atmospheric refractive index changing with altitude, h is altitude, dh is altitude differential, h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively;

[0028] The atmospheric refractive index function with altitude f(h) is approximately calculated according to the following formula

[0029]

[0030] Where h is the elevation in meters, h tropop is the average height of the tropopause.

[0031] Preferably, in step S5:

[0032] Based on the starting point of incidence Incident sight vector Height h k , calculate the intersection position

[0033] Step A.1: Iteration start value t = 0, position

[0034] Step A.2: Calculate the position Corresponding geographical latitude Longitude λ t 、Elevation h t ;

[0035] Step A.3: If the current h t -h k The calculation threshold was not reached, the position Update once;

[0036]

[0037] in, Represents two unit vectors and Calculate the inner product, Indicates location The external normal direction at ;

[0038] Step A.4, repeat steps A.2 to A.3 until h t -h k Less than the set threshold.

[0039] Preferably, in step S6:

[0040] If the height h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? k ;

[0041] If the intersection position The digital elevation corresponding to the latitude and longitude is greater than or equal to h k , calculate the incident sight vector Intersection with digital elevation surface Output intersection coordinates;

[0042] Incident sight vector Intersection with digital elevation surface The calculation includes the following steps:

[0043] Step C.1: If the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint

[0044]

[0045] Step C.2: Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle , if η middle -h middle The calculation threshold is not reached. If η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat steps C.1 to C.2 until n middle -h middle Less than the set threshold.

[0046] Preferably, in step S6:

[0047] If the intersection position The digital elevation corresponding to the latitude and longitude is less than h k , according to the refractive index n of the incident layer k , refractive index of the exit layer n k+1 Calculate the view vector after refraction

[0048] Step B.1: Calculate the angle of incidence θ k :

[0049]

[0050] in, express The outward normal direction at Represents two unit vectors and The inner product operation of ;

[0051] Step B.2, calculate the angle of incidence

[0052]

[0053] Step B.3, calculate the sight vector after refraction

[0054]

[0055] Set k=k+1 and repeat steps S5 to S6.

[0056] Preferably, except for the empirical function calculation of the atmospheric refractive index varying with altitude in step S3, which corresponds to an optical spectrum of 0.4 to 0.9 μm, the remaining calculation steps have no spectrum limitation.

[0057] According to the present invention, an atmospheric refraction compensation system for optical earth remote sensing satellite image geolocation is provided, comprising:

[0058] Module M1: Obtain the position of the satellite in the earth-fixed coordinate system at the corresponding exposure time and the line of sight vector corresponding to each detector;

[0059] Module M2: Establish a multi-layer atmospheric refraction model and stratify the atmosphere by altitude;

[0060] Module M3: Calculate the average atmospheric refractive index of each layer;

[0061] Module M4: Calculate sequentially starting from layer 0;

[0062] Module M5: Calculate the line of sight intersection position based on the incident starting point, incident line of sight vector and elevation.

[0063] Module M6: If the elevation of the intersection is less than the maximum elevation of the area, determine whether the digital elevation corresponding to the intersection position is greater than or equal to the elevation of the intersection;

[0064] If the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation, calculate the intersection of the incident sight vector and the digital elevation surface, and output the intersection coordinates;

[0065] If the digital elevation corresponding to the intersection position is less than the intersection elevation, the sight line vector after refraction is calculated based on the refractive index of the incident layer and the refractive index of the exit layer, and the layer number is increased by 1 to repeat modules M5 to M6.

[0066] Preferably, in the module M1:

[0067] Get the satellite position in the Earth-fixed coordinate system ECEF at the corresponding exposure time The line of sight vector corresponding to each detector is a unit vector, where the subscript j represents the detector number;

[0068] The line of sight vectors corresponding to each detector are accurately pointed after optical distortion correction and on-board thermal deformation correction;

[0069] In the module M2:

[0070] A multi-layer atmospheric refraction model is established, and the atmosphere is divided into M layers according to the altitude. The altitudes corresponding to the tops of the layers from top to bottom are h0, h1, ..., h M-1 ;

[0071] The number of atmospheric stratification layers is set according to the angle of the satellite observation from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation from the sub-satellite point is smaller, the number of layers is reduced to reduce the amount of calculation;

[0072] In the module M3:

[0073] Calculate the average atmospheric refractive index n for each layer k , where k represents the layer number, k = 0, 1, 2, ..., M; n0 is the refractive index under vacuum, n0 = 1;

[0074] The calculation method of the average atmospheric refractive index nk corresponding to the kth layer is:

[0075]

[0076] Where k = 1, 2, ..., M, function f(h) represents the function of atmospheric refractive index changing with altitude, h is altitude, dh is altitude differential, h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively;

[0077] The atmospheric refractive index function with altitude f(h) is approximately calculated according to the following formula

[0078]

[0079] Where h is the elevation in meters, h tropop is the average height of the tropopause;

[0080] In the module M5:

[0081] Based on the starting point of incidence Incident sight vector Height h k , calculate the intersection position

[0082] Step A.1: Iteration start value t = 0, position

[0083] Step A.2: Calculate the position Corresponding geographical latitude Longitude λ t 、Elevation h t ;

[0084] Step A.3: If the current h t -h k The calculation threshold was not reached, the position Update once;

[0085]

[0086] in, Represents two unit vectors and Calculate the inner product, Indicates location The external normal direction at ;

[0087] Step A.4, repeat steps A.2 to A.3 until h t -h k Less than the set threshold;

[0088] In the module M6:

[0089] If the height h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? k ;

[0090] If the intersection position The digital elevation corresponding to the latitude and longitude is greater than or equal to h k , calculate the incident sight vector Intersection with digital elevation surface Output intersection coordinates;

[0091] Incident sight vector Intersection with digital elevation surface The calculation includes the following steps:

[0092] Step C.1: If the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint

[0093]

[0094] Step C.2: Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle , if η middle -h middle The calculation threshold is not reached. If η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat steps C.1 to C.2 until n middle -h middle Less than the set threshold;

[0095] In the module M6:

[0096] If the intersection position The digital elevation corresponding to the latitude and longitude is less than h k , according to the refractive index n of the incident layer k , refractive index of the exit layer n k+1 Calculate the view vector after refraction

[0097] Step B.1: Calculate the angle of incidence θ k :

[0098]

[0099] in, express The outward normal direction at Represents two unit vectors and The inner product operation of ;

[0100] Step B.2, calculate the angle of incidence

[0101]

[0102] Step B.3, calculate the sight vector after refraction

[0103]

[0104] Set k=k+1 and repeat modules M5-M6;

[0105] Except for the empirical function calculation of atmospheric refractive index varying with altitude in module M3, which corresponds to the optical spectrum of 0.4-0.9 μm, the remaining calculation steps have no spectrum restrictions.

[0106] Compared with the prior art, the present invention has the following beneficial effects:

[0107] 1. Aiming at the demand for high-precision geographic positioning in high-resolution optical earth remote sensing, the present invention models the atmospheric refraction that causes the non-rigorous geometric positioning model, constructs a multi-layer atmospheric refraction model that can be adjusted according to the needs, and provides the process and method of geolocation under this model. Combined with the digital elevation model, it can provide high-precision geometric positioning results without control points.

[0108] 2. The method of the present invention is reasonable, simple to calculate, and easy to implement, and can be effectively applied to the atmospheric refraction compensation of optical earth remote sensing images. BRIEF DESCRIPTION OF THE DRAWINGS

[0109] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings:

[0110] Figure 1 is a flow chart of the present invention;

[0111] Figure 2 This is a schematic diagram of light refraction under the layered atmosphere model;

[0112] Figure 3 A result diagram of multi-layer atmospheric refraction modeling performed for the present invention;

[0113] Figure 4 This is the difference between the result of geo-positioning a certain low-orbit satellite using the method of the present invention and that without considering atmospheric refraction. DETAILED DESCRIPTION

[0114] The present invention is described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those of ordinary skill in the art, several changes and improvements can also be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0115] Embodiment 1:

[0116] The present invention comprehensively considers the effects of optical distortion correction, on-board thermal deformation correction, atmospheric refraction, earth elevation, etc., aiming to provide a universal atmospheric refraction method that can improve the geolocation accuracy of high-resolution optical earth remote sensing satellite images. The present invention is easy to implement, and can improve the geolocation accuracy of high-resolution optical earth remote sensing satellite images through atmospheric refraction compensation.

[0117] According to the present invention, a method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images is provided. Figure 1 As shown, including:

[0118] Step S1: Obtain the position of the satellite in the earth-fixed coordinate system at the corresponding exposure time and the line of sight vector corresponding to each detector;

[0119] Specifically, in step S1:

[0120] Get the satellite position in the Earth-fixed coordinate system ECEF at the corresponding exposure time The line of sight vector corresponding to each detector is a unit vector, where the subscript j represents the detector number;

[0121] The line of sight vector corresponding to each detector is accurately pointed after being corrected for optical distortion and on-board thermal deformation.

[0122] Step S2: Establish a multi-layer atmospheric refraction model and divide the atmosphere into layers according to altitude;

[0123] Specifically, in step S2:

[0124] A multi-layer atmospheric refraction model is established, and the atmosphere is divided into M layers according to the altitude. The altitudes corresponding to the tops of the layers from top to bottom are h0, h1, ..., h M-1 ;

[0125] The number of atmospheric stratification layers is set according to the angle of the satellite observation deviating from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation deviating from the sub-satellite point is smaller, the number of layers is reduced to reduce the amount of calculation.

[0126] Step S3: Calculate the average atmospheric refractive index of each layer;

[0127] Specifically, in step S3:

[0128] Calculate the average atmospheric refractive index n for each layer k , where k represents the layer number, k = 0, 1, 2, ..., M; n0 is the refractive index under vacuum, n0 = 1;

[0129] The average atmospheric refractive index n corresponding to the kth layer k The calculation method is:

[0130]

[0131] Where k = 1, 2, ..., M, function f(h) represents the function of atmospheric refractive index changing with altitude, h is altitude, dh is altitude differential, h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively;

[0132] The atmospheric refractive index function with altitude f(h) is approximately calculated according to the following formula

[0133]

[0134] Where h is the elevation in meters, h tropop is the average height of the tropopause.

[0135] Step S4: Calculate sequentially starting from layer 0;

[0136] Step S5: Calculate the line of sight intersection position based on the incident starting point, the incident line of sight vector, and the elevation.

[0137] Specifically, in step S5:

[0138] Based on the starting point of incidence Incident sight vector Height h k , calculate the intersection position

[0139] Step A.1: Iteration start value t = 0, position

[0140] Step A.2: Calculate the position Corresponding geographical latitude Longitude λ t 、Elevation h t ;

[0141] Step A.3: If the current h t -h k The calculation threshold was not reached, the position Update once;

[0142]

[0143] in, Represents two unit vectors and Calculate the inner product, Indicates location The external normal direction at ;

[0144] Step A.4, repeat steps A.2 to A.3 until h t -h k Less than the set threshold.

[0145] Step S6: If the elevation of the intersection point is less than the maximum elevation of the area, determine whether the digital elevation corresponding to the intersection point position is greater than or equal to the elevation of the intersection point;

[0146] If the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation, calculate the intersection of the incident sight vector and the digital elevation surface, and output the intersection coordinates;

[0147] If the digital elevation corresponding to the intersection position is less than the intersection elevation, the line of sight vector after refraction is calculated according to the refractive index of the incident layer and the refractive index of the exit layer, and the layer number is increased by 1 and steps S5 to S6 are repeated.

[0148] Specifically, in step S6:

[0149] If the height h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? k ;

[0150] If the intersection position The digital elevation corresponding to the latitude and longitude is greater than or equal to h k , calculate the incident sight vector Intersection with digital elevation surface Output intersection coordinates;

[0151] Incident sight vector Intersection with digital elevation surface The calculation includes the following steps:

[0152] Step C.1: If the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint

[0153]

[0154] Step C.2: Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle , if η middle -h middle The calculation threshold is not reached. If η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat steps C.1 to C.2 until n middle -h middle Less than the set threshold.

[0155] Specifically, in step S6:

[0156] If the intersection position The digital elevation corresponding to the latitude and longitude is less than h k , according to the refractive index n of the incident layer k , refractive index of the exit layer n k+1 Calculate the view vector after refraction

[0157] Step B.1: Calculate the angle of incidence θ k :

[0158]

[0159] in, express The outward normal direction at Represents two unit vectors and The inner product operation of ;

[0160] Step B.2, calculate the angle of incidence

[0161]

[0162] Step B.3, calculate the sight vector after refraction

[0163]

[0164] Set k=k+1 and repeat steps S5 to S6.

[0165] Specifically, except for the empirical function calculation of the atmospheric refractive index varying with altitude in step S3, which corresponds to an optical spectrum of 0.4 to 0.9 μm, the remaining calculation steps have no spectrum limitation.

[0166] Embodiment 2:

[0167] Embodiment 2 is a preferred example of Embodiment 1, and is used to illustrate the present invention in more detail.

[0168] The present invention also provides an atmospheric refraction compensation system for geo-positioning of optical earth remote sensing satellite images. The atmospheric refraction compensation system for geo-positioning of optical earth remote sensing satellite images can be realized by executing the process steps of the atmospheric refraction compensation method for geo-positioning of optical earth remote sensing satellite images, that is, those skilled in the art can understand the atmospheric refraction compensation method for geo-positioning of optical earth remote sensing satellite images as a preferred implementation of the atmospheric refraction compensation system for geo-positioning of optical earth remote sensing satellite images.

[0169] According to the present invention, an atmospheric refraction compensation system for optical earth remote sensing satellite image geolocation is provided, comprising:

[0170] Module M1: Obtain the position of the satellite in the earth-fixed coordinate system at the corresponding exposure time and the line of sight vector corresponding to each detector;

[0171] Specifically, in the module M1:

[0172] Get the satellite position in the Earth-fixed coordinate system ECEF at the corresponding exposure time The line of sight vector corresponding to each detector is a unit vector, where the subscript j represents the detector number;

[0173] The line of sight vectors corresponding to each detector are accurately pointed after optical distortion correction and on-board thermal deformation correction;

[0174] In the module M2:

[0175] A multi-layer atmospheric refraction model is established, and the atmosphere is divided into M layers according to the altitude. The altitudes corresponding to the tops of the layers from top to bottom are h0, h1, ..., h M-1 ;

[0176] The number of atmospheric stratification layers is set according to the angle of the satellite observation from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation from the sub-satellite point is smaller, the number of layers is reduced to reduce the amount of calculation;

[0177] Module M2: Establish a multi-layer atmospheric refraction model and stratify the atmosphere by altitude;

[0178] Module M3: Calculate the average atmospheric refractive index of each layer;

[0179] In the module M3:

[0180] Calculate the average atmospheric refractive index n for each layer k , where k represents the layer number, k = 0, 1, 2, ..., M; n0 is the refractive index under vacuum, n0 = 1;

[0181] The average atmospheric refractive index n corresponding to the kth layer k The calculation method is:

[0182]

[0183] Where k = 1, 2, ..., M, function f(h) represents the function of atmospheric refractive index changing with altitude, h is altitude, dh is altitude differential, h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively;

[0184] The atmospheric refractive index function with altitude f(h) is approximately calculated according to the following formula

[0185]

[0186] Where h is the elevation in meters, h tropop is the average height of the tropopause;

[0187] Module M4: Calculate sequentially starting from layer 0;

[0188] Module M5: Calculate the line of sight intersection position based on the incident starting point, incident line of sight vector and elevation.

[0189] In the module M5:

[0190] Based on the starting point of incidence Incident sight vector Height h k , calculate the intersection position

[0191] Step A.1: Iteration start value t = 0, position

[0192] Step A.2: Calculate the position Corresponding geographical latitude Longitude λ t 、Elevation h t ;

[0193] Step A.3: If the current h t -h k The calculation threshold was not reached, the position Update once;

[0194]

[0195] in, Represents two unit vectors and Calculate the inner product, Indicates location The external normal direction at ;

[0196] Step A.4, repeat steps A.2 to A.3 until h t -h k Less than the set threshold;

[0197] Module M6: If the elevation of the intersection is less than the maximum elevation of the area, determine whether the digital elevation corresponding to the intersection position is greater than or equal to the elevation of the intersection;

[0198] If the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation, calculate the intersection of the incident sight vector and the digital elevation surface, and output the intersection coordinates;

[0199] If the digital elevation corresponding to the intersection position is less than the intersection elevation, the sight line vector after refraction is calculated based on the refractive index of the incident layer and the refractive index of the exit layer, and the layer number is increased by 1 to repeat modules M5 to M6.

[0200] In the module M6:

[0201] If the height h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? k ;

[0202] If the intersection position The digital elevation corresponding to the latitude and longitude is greater than or equal to h k , calculate the incident sight vector Intersection with digital elevation surface Output intersection coordinates;

[0203] Incident sight vector Intersection with digital elevation surface The calculation includes the following steps:

[0204] Step C.1: If the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint

[0205]

[0206] Step C.2: Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle , if η middle -h middle The calculation threshold is not reached. If η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat steps C.1 to C.2 until n middle -h middle Less than the set threshold;

[0207] In the module M6:

[0208] If the intersection position The digital elevation corresponding to the latitude and longitude is less than h k , according to the refractive index n of the incident layer k , refractive index of the exit layer n k+1 Calculate the view vector after refraction

[0209] Step B.1: Calculate the angle of incidence θ k :

[0210]

[0211] in, express The outward normal direction at Represents two unit vectors and The inner product operation of ;

[0212] Step B.2, calculate the angle of incidence

[0213]

[0214] Step B.3, calculate the sight vector after refraction

[0215]

[0216] Set k=k+1 and repeat modules M5-M6;

[0217] Except for the empirical function calculation of atmospheric refractive index varying with altitude in module M3, which corresponds to the optical spectrum of 0.4-0.9 μm, the remaining calculation steps have no spectrum restrictions.

[0218] Embodiment 3:

[0219] Embodiment 3 is a preferred example of Embodiment 1, and is used to illustrate the present invention in more detail.

[0220] The technical problem to be solved by the present invention is to provide an atmospheric refraction compensation method for geographic positioning of optical earth remote sensing satellite images, so as to improve the geometric positioning accuracy of high-resolution optical earth remote sensing satellite images without control points.

[0221] When geo-positioning the impact of optical remote sensing satellites, the position of the satellite in the Earth Fixed Coordinate System (ECEF) at the corresponding exposure time is first calculated based on a rigorous geometric imaging model. The line of sight vector corresponding to each detector (unit vector, where the subscript j represents the detector number). When the influence of atmospheric refraction is not considered, that is, the position of the Earth-fixed coordinate system Line of sight vector And the Earth's Digital Elevation Model (DEM) can calculate the geographic location of pixels without control points.

[0222] However, with the improvement of the spatial resolution of earth remote sensing, the accuracy requirements for geolocation are also increasing. Geolocation without control points is the basis of accurate geolocation. Currently, there are many methods that add factors such as optical geometric distortion and satellite thermal deformation into the process of line of sight determination, thereby improving the accuracy of geolocation without control points to a certain extent. On this basis, the present invention further optimizes the rigorous geometric imaging model and models atmospheric refraction.

[0223] The influence of atmospheric refraction is manifested as deflection of optics, which is related to the atmospheric refractive index and the angle of incidence, and the atmospheric refractive index is related to multiple factors such as spectral wavelength, temperature, humidity, and air pressure. Atmospheric temperature, humidity, and air pressure are all related to altitude, but not only depend on altitude. To simplify the model, the present invention approximates the atmospheric refractive index as a function related only to altitude, ignoring the differences in the atmosphere at different longitudes and latitudes and different seasons of the earth. These differences will cause the refractive index to change, but generally speaking, the magnitude is small.

[0224] According to the theory of Hoyt and Storey, the atmospheric refractive index f can be approximately given by the following empirical formula

[0225] f=1+0.000295α (Formula 8)

[0226] Where α is the density factor. In the troposphere, if the global average temperature lapse rate r is 0.0065K / m, then the density factor α in the spectral range of 0.4 to 0.9 μm in the troposphere can be approximated as

[0227]

[0228] Where h represents the elevation, T sealevel is the average sea level temperature (unit K), M mean is the average molecular weight of tropospheric gas (about 28.825), g0 is the average gravitational acceleration at sea level (9.805 m / s 2 ), R gas is the ideal gas constant (8314.3 J (k mol) -1 K -1 ).

[0229] Above the troposphere, the density factor is considered to decrease exponentially.

[0230]

[0231] where α tropop represents the density factor of the tropopause, T tropop represents the absolute temperature of the tropopause, h tropop is the average height of the tropopause.

[0232] Substituting formula 9 and formula 10 into formula 8, we have

[0233]

[0234] The present invention constructs a multi-layer atmospheric refraction model, and divides the atmosphere into M layers according to the altitude. The altitudes corresponding to the tops of the layers from top to bottom are h0, h1, ..., h M-1 The average atmospheric refractive index of each layer is expressed as n k The average atmospheric refractive index n corresponding to the kth layer (k = 1, 2, ..., M) is represented by (where k represents the layer number, k = 0, 1, 2, ..., M, n0 is the refractive index under vacuum, n0 = 1). k The calculation method is

[0235]

[0236] Among them, the function f(h) represents the change function of atmospheric refractive index with altitude, h is the altitude, dh is the altitude differential, and h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively.

[0237] Under the same conditions, the larger the incident angle, the greater the light deflection. Therefore, the number of atmospheric stratification layers can be set according to the angle of the satellite observation deviating from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation deviating from the sub-satellite point is smaller, the number of layers can be reduced to reduce the amount of calculation.

[0238] Outside the atmosphere, the atmospheric refractive index can be considered to be 1, and light travels in a straight line. In a stratified atmosphere, refraction occurs between layers, and light travels in a straight line within a layer. Figure 2 The schematic diagram shown in the figure starts from k = 0 and is composed of the incident layer refractive index n k , refractive index of the exit layer n k+1 , starting point of incidence Incident sight vector Height h k , calculate the intersection position and the sight vector after refraction

[0239] Intersection location The calculation is performed iteratively, first at t = 0, position Calculate position Corresponding geographical latitude Longitude λ t 、Elevation h t .

[0240] If the current h t -h k Less than the calculated threshold, current That is, the output is the intersection coordinate value. If the calculation threshold is not reached, the position Iterative updates:

[0241]

[0242] in, Represents two unit vectors and Calculate the inner product, Indicates location The outward normal direction at .

[0243] When calculating to a lower height, the surface elevation must be considered. If the current layer elevation is h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? kIf yes, calculate the intersection with the digital elevation model surface; otherwise, use the already calculated intersection and the incident layer refractive index n k , refractive index of the exit layer n k+1 Calculate the view vector after refraction

[0244] First calculate the incident angle θ k .

[0245]

[0246] in, express The outward normal direction at Represents two unit vectors and The inner product operation of .

[0247] Calculate the angle of incidence according to Snell's law

[0248]

[0249] Further calculate the line of sight vector after refraction in the three-dimensional coordinate space

[0250]

[0251] Calculate the incident sight vector Intersection with digital elevation surface Also through iteration, if the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint

[0252]

[0253] Calculate midpoint Digital elevation η corresponding to latitude and longitude middle , if η middle -h middle The calculation threshold is not reached. If η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat the iteration until η middle -h middleLess than the set threshold. When calculating the elevation of a digital elevation model corresponding to any point in the earth's fixed coordinate system, a two-dimensional interpolation operation is required. The digital elevation model only provides elevation data of fixed grid points, and the elevation data of any point can be obtained through bilinear interpolation.

[0254] According to the reversibility of the optical path, the position can be known The corresponding light source reaches the corresponding detector after passing through the atmosphere, that is, the corresponding position of the detector's geolocation result

[0255] The effectiveness of the method of the present invention is illustrated below with reference to data from a low-orbit optical remote sensing satellite. Figure 3 As shown by the solid line, after executing step 2 to stratify the atmosphere at every 1 km altitude, the average atmospheric refractive index of each layer is Figure 3 The satellite position at a certain moment is (Unit: m), the sight vector corresponding to a certain pixel of the detector is The line of sight vector corresponds to a deviation of 58.6879° from the subsatellite point. Combined with the DEM data, if the influence of atmospheric refraction is not considered, the geolocation result is: latitude 30.532304°, longitude 95.181367°, and elevation 4882.2m. The final geolocation result obtained by the method of the present invention is: latitude 30.532274°, longitude 95.182102°, and elevation 4880.2m. After considering atmospheric refraction, the differences in positioning results corresponding to different elevations are shown in the attached figure. Figure 4 As shown, in the scenario with a large side viewing angle in this embodiment, the error caused by atmospheric refraction is close to 100 meters.

[0256] Those skilled in the art know that, in addition to realizing the system and its various devices, modules, and units provided by the present invention in a purely computer-readable program code, it is entirely possible to realize the same functions in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps. Therefore, the system and its various devices, modules, and units provided by the present invention can be considered as a hardware component, and the devices, modules, and units included therein for realizing various functions can also be regarded as structures within the hardware component; the devices, modules, and units for realizing various functions can also be regarded as both software modules for realizing the method and structures within the hardware component.

[0257] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which does not affect the essence of the present invention. In the absence of conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. A method for compensating atmospheric refraction for geo-positioning of optical earth remote sensing satellite images, characterized in that: include: Step S1: Obtain the position of the satellite in the earth-fixed coordinate system at the corresponding exposure time and the line of sight vector corresponding to each detector; Step S2: Establish a multi-layer atmospheric refraction model and stratify the atmosphere according to altitude; Step S3: Calculate the average atmospheric refractive index of each layer; Step S4: Calculate sequentially starting from layer 0; Step S5: Calculate the position of the sight line intersection point according to the incident starting point, the incident sight line vector, and the elevation; Step S6: If the elevation of the intersection point is less than the maximum elevation of the area, determine whether the digital elevation corresponding to the intersection point position is greater than or equal to the elevation of the intersection point; If the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation, calculate the intersection of the incident sight vector and the digital elevation surface, and output the intersection coordinates; If the digital elevation corresponding to the intersection position is less than the intersection elevation, the line of sight vector after refraction is calculated according to the refractive index of the incident layer and the refractive index of the exit layer, and the layer number is increased by 1 and steps S5 to S6 are repeated.

2. The method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images according to claim 1, characterized in that: In step S1: Get the satellite position in the Earth-fixed coordinate system ECEF at the corresponding exposure time The line of sight vector corresponding to each detector is a unit vector, where the subscript j represents the detector number; The line of sight vector corresponding to each detector is accurately pointed after being corrected for optical distortion and on-board thermal deformation.

3. The method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images according to claim 1, characterized in that: In step S2: A multi-layer atmospheric refraction model is established, and the atmosphere is divided into M layers according to the altitude. The altitudes of the tops of the layers from top to bottom are h0, h1, ..., h M-1 ; The number of atmospheric stratification layers is set according to the angle of the satellite observation deviating from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation deviating from the sub-satellite point is smaller, the number of layers is reduced to reduce the amount of calculation.

4. The method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images according to claim 1, characterized in that: In step S3: Calculate the average atmospheric refractive index n for each layer k , where k represents the layer number, k = 0, 1, 2, ..., M; n0 is the refractive index under vacuum, n0 = 1; The average atmospheric refractive index n corresponding to the kth layer k The calculation method is: Where k = 1, 2, ..., M, function f(h) represents the function of atmospheric refractive index changing with altitude, h is altitude, dh is altitude differential, h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively; The atmospheric refractive index function with altitude f(h) is approximately calculated according to the following formula Where h is the elevation in meters, h tropop is the average height of the tropopause.

5. The method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images according to claim 1, characterized in that: In step S5: Based on the starting point of incidence Incident sight vector Height h k , calculate the intersection position Step A.1: Iteration start value t = 0, position Step A.2: Calculate the position Corresponding geographical latitude Longitude λ t 、Elevation h t ; Step A.3: If the current |h t -h k |The calculation threshold was not reached, position Update once; in, Represents two unit vectors and Calculate the inner product, Indicates location The external normal direction at ; Step A.4, repeat steps A.2 to A.3 until |h t -h k |Less than the set threshold.

6. The method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images according to claim 1, characterized in that: In step S6: If the height h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? k ; If the intersection position The digital elevation corresponding to the latitude and longitude is greater than or equal to h k , calculate the incident sight vector Intersection with digital elevation surface Output intersection coordinates; Incident sight vector Intersection with digital elevation surface The calculation includes the following steps: Step C.1: If the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint Step C.2: Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle , if |η middle -h middle |The calculation threshold is not reached. At this time, if η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat steps C.1 to C.2 until |η middle -h middle |Less than the set threshold.

7. The method for compensating atmospheric refraction for geographic positioning of optical earth remote sensing satellite images according to claim 1, characterized in that: In step S6: If the intersection position The digital elevation corresponding to the latitude and longitude is less than h k , according to the refractive index n of the incident layer k , refractive index of the exit layer n k+1 Calculate the view vector after refraction Step B.1: Calculate the angle of incidence θ k : in, express The outward normal direction at Represents two unit vectors and The inner product operation of ; Step B.2, calculate the angle of incidence Step B.3, calculate the sight vector after refraction Set k=k+1 and repeat steps S5 to S6.

8. The method for compensating atmospheric refraction for optical earth remote sensing satellite image geolocation according to claim 1, characterized in that: Except for the empirical function calculation of atmospheric refractive index varying with altitude in step S3, which corresponds to an optical spectrum of 0.4 to 0.9 μm, the remaining calculation steps have no spectrum restrictions.

9. An optical earth remote sensing satellite image geolocation atmospheric refraction compensation system, characterized in that: include: Module M1: obtain the position of the satellite in the earth-fixed coordinate system at the corresponding exposure time and the line of sight vector corresponding to each detector; Module M2: Establish a multi-layer atmospheric refraction model and stratify the atmosphere by altitude; Module M3: Calculate the average atmospheric refractive index of each layer; Module M4: Calculate sequentially starting from layer 0; Module M5: Calculate the line of sight intersection position based on the incident starting point, incident line of sight vector and elevation; Module M6: If the elevation of the intersection is less than the maximum elevation of the area, determine whether the digital elevation corresponding to the intersection position is greater than or equal to the elevation of the intersection; If the digital elevation corresponding to the intersection position is greater than or equal to the intersection elevation, calculate the intersection of the incident sight vector and the digital elevation surface, and output the intersection coordinates; If the digital elevation corresponding to the intersection position is less than the intersection elevation, the sight line vector after refraction is calculated based on the refractive index of the incident layer and the refractive index of the exit layer, and the layer number is increased by 1 to repeat modules M5 to M6.

10. The optical earth remote sensing satellite image geo-positioning atmospheric refraction compensation system according to claim 9, characterized in that: In the module M1: Get the satellite position in the Earth-fixed coordinate system ECEF at the corresponding exposure time The line of sight vector corresponding to each detector is a unit vector, where the subscript j represents the detector number; The line of sight vectors corresponding to each detector are accurately pointed after optical distortion correction and on-board thermal deformation correction; In the module M2: A multi-layer atmospheric refraction model is established, and the atmosphere is divided into M layers according to the altitude. The altitudes of the tops of the layers from top to bottom are h0, h1, ..., h M-1 ; The number of atmospheric stratification layers is set according to the angle of the satellite observation from the sub-satellite point. When the deviation from the sub-satellite point is greater, the number of layers M is greater. Conversely, when the angle of the satellite observation from the sub-satellite point is smaller, the number of layers is reduced to reduce the amount of calculation; In the module M3: Calculate the average atmospheric refractive index n for each layer k , where k represents the layer number, k = 0, 1, 2, ..., M; n0 is the refractive index under vacuum, n0 = 1; The average atmospheric refractive index n corresponding to the kth layer k The calculation method is: Where k = 1, 2, ..., M, function f(h) represents the function of atmospheric refractive index changing with altitude, h is altitude, dh is altitude differential, h k-1 、h k are the corresponding elevations of the top and bottom of the kth atmospheric layer respectively; The atmospheric refractive index function with altitude f(h) is approximately calculated according to the following formula Where h is the elevation in meters, h tropop is the average height of the tropopause; In the module M5: Based on the starting point of incidence Incident sight vector Height h k , calculate the intersection position Step A.1: Iteration start value t = 0, position Step A.2: Calculate the position Corresponding geographical latitude Longitude λ t 、Elevation h t ; Step A.3: If the current |h t -h k |The calculation threshold was not reached, position Update once; in, Represents two unit vectors and Calculate the inner product, Indicates location The external normal direction at ; Step A.4, repeat steps A.2 to A.3 until |h t -h k |Less than the set threshold; In the module M6: If the height h k Less than the maximum elevation of the area, determine the intersection position Is the digital elevation corresponding to the latitude and longitude greater than or equal to h? k ; If the intersection position The digital elevation corresponding to the latitude and longitude is greater than or equal to h k , calculate the incident sight vector Intersection with digital elevation surface Output intersection coordinates; Incident sight vector Intersection with digital elevation surface The calculation includes the following steps: Step C.1: If the intersection position The digital elevation corresponding to the latitude and longitude is η, which is determined by the incident sight vector Starting point Intersection point Calculate midpoint Step C.2: Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle , if |η middle -h middle |The calculation threshold is not reached. At this time, if η middle >h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat steps C.1 to C.2 until |η middle -h middle |Less than the set threshold; In the module M6: If the intersection position The digital elevation corresponding to the latitude and longitude is less than h k , according to the refractive index n of the incident layer k , refractive index of the exit layer n k+1 Calculate the view vector after refraction Step B.1: Calculate the angle of incidence θ k : in, express The outward normal direction at Represents two unit vectors and The inner product operation of ; Step B.2, calculate the angle of incidence Step B.3, calculate the sight vector after refraction Set k=k+1, repeat modules M5 to M6; Except for the empirical function calculation of atmospheric refractive index varying with altitude in module M3, which corresponds to the optical spectrum of 0.4-0.9 μm, the remaining calculation steps have no spectrum restrictions.

Images