Methods and systems for atmospheric refraction compensation in geolocation of optical Earth remote sensing satellite imagery

By constructing a multi-layer atmospheric refraction model and iterative calculation of line-of-sight intersections, the problem of atmospheric refraction in high-resolution optical Earth remote sensing satellite images was solved, and high-precision control-point-free geolocation was achieved.

CN119688649BActive Publication Date: 2025-10-28SHANGHAI SATELLITE ENG INST
View PDF 5 Cites 0 Cited by

Patent Information

Application Number
CN202411626789.1
Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
Filing Date
2024-11-14
Publication Date
2025-10-28
Estimated Expiration
2044-11-14

AI Technical Summary

Technical Problem

Existing technologies have failed to effectively compensate for the effects of atmospheric refraction in geolocation using high-resolution optical Earth remote sensing satellite imagery, resulting in insufficient positioning accuracy.

Method used

A multi-layer atmospheric refraction model is constructed. By obtaining the line-of-sight vector and atmospheric layers at the time of satellite exposure, the average refractive index of each layer is calculated, the position of the line-of-sight intersection is calculated iteratively, and high-precision geolocation is achieved by combining it with a digital elevation model.

Benefits of technology

It improves the control point-free geometric positioning accuracy of high-resolution optical Earth remote sensing satellite imagery, simplifies the calculation process, and enhances positioning accuracy.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN119688649B_ABST
    Figure CN119688649B_ABST
Patent Text Reader

Abstract

This invention provides a method and system for atmospheric refraction compensation in optical Earth remote sensing satellite image geolocation, comprising: acquiring the satellite's position in the Earth-fixed coordinate system at the corresponding exposure time, and the line-of-sight vectors corresponding to each detector; establishing a multi-layer atmospheric refraction model, dividing the atmosphere into layers according to elevation; calculating the average atmospheric refractive index of each layer; performing calculations sequentially starting from layer 0; calculating the intersection point of the line-of-sight vectors based on the incident point, the incident line-of-sight vector, and the elevation. If the elevation of the intersection point is less than the maximum elevation of the region, determining whether the digital elevation corresponding to the intersection point is greater than or equal to the intersection point elevation; if so, calculating the intersection point of the incident line-of-sight vector and the digital elevation surface, and outputting the intersection point coordinates; if not, calculating the refracted line-of-sight vector based on the refractive index of the incident layer and the refractive index of the exit layer, and incrementing the layer number by 1 for repetition. This invention is reasonable, computationally simple, and easy to implement, and can be effectively applied to atmospheric refraction compensation in optical Earth remote sensing images.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of geolocation technology, specifically to a method and system for atmospheric refraction compensation in geolocation of optical remote sensing satellite images. Background Technology

[0002] A crucial aspect of optical Earth remote sensing applications is the geolocation of remote sensing impacts. Traditional methods, based on rigorous geometric imaging models, rely on collinearity equations for mapping calculations. However, with increasing demands for geolocation accuracy, more and more non-ideal factors are being considered and specifically corrected. These include geometric distortion of optical lenses, the accuracy of satellite platform ephemeris measurements, attitude measurement accuracy, and the effects of thermal deformation on the satellite platform. In recent years, with the application of 10-meter, 1-meter, and sub-meter-level optical Earth remote sensing satellites, the demand for geolocation accuracy has further increased.

[0003] Atmospheric refraction alters the rectilinear propagation model of light, causing inaccuracies in the collinearity model. According to existing literature, for low-orbit satellites observing a zenith angle of 55°, atmospheric refraction can cause positioning errors on the order of ten meters. Atmospheric refraction affects optical Earth remote sensing and also impacts ground-based laser strikes, space-based laser altimetry, ground-based telemetry and control equipment, and ground-based astronomical telescope observations. Reference 1 (Xi Hua, et al., Analysis of Factors Affecting Atmospheric Refraction Correction, Value Engineering, 2012, No. 22, pp. 313-314) analyzes the impact of atmospheric dispersion on the deviation of the laser's propagation path from a straight line, leading to inconsistencies between aiming and striking paths. It focuses on the influence of different latitudes and carbon dioxide concentrations on the elevation angle error caused by atmospheric refraction. Reference 2 (Wang Taiping, Research on Waveform Processing and Positioning Accuracy Optimization Technology of Gaofen-7 Satellite-borne Laser Altimeter Data, Wuhan University, Master's Thesis, 2021) describes the deviation of laser footpoint position caused by distance measurement delay due to atmospheric refraction effect for satellite-borne laser altimeters, and presents a correction process for atmospheric interference delay for satellite-borne laser altimeters. Reference 3 (Wu Pengfei, Low Elevation Angle Atmospheric Refraction Correction Method Based on Gridded Atmospheric Parameter Profile Model, Acta Optica Sinica, 2017, Vol. 37, No. 6, 0601004) establishes a gridded atmospheric parameter profile model for low elevation angle tracking below 5° and analyzes the refraction correction results under typical atmospheric conditions. Reference 4 (Liang Xiaobo, Research on Telescope Pointing Correction Based on Astronomical Orientation Technology, Doctoral Dissertation, University of Chinese Academy of Sciences, 2020) describes the impact of atmospheric refraction error (ambient fog) on ​​ground-based astronomical telescopes and uses the atmospheric refraction model in IAU SOFA to correct for ambient fog.

[0004] In the field of high-resolution optical Earth remote sensing, reference 5 (Cheng Yufeng, Research on High-Precision On-Orbit Autonomous Geometric Calibration Method for High-Resolution Optical Remote Sensing Satellites, Doctoral Dissertation, Wuhan University, 2019) analyzes the impact of atmospheric refraction on the on-orbit geometric calibration of optical satellites and presents a geometric calibration method that takes into account atmospheric conditions 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, Tsinghua University, 2009) analyzes the influence of atmospheric refraction. Based on SPOT-5 experimental data, it selects control points to fit the radius that needs compensation, and then uses this as a feature library to perform atmospheric refraction compensation correction on another SPOT-5 remote sensing data. Reference 7 (Man Yiyun, et al., Error Analysis of Planar Positioning Accuracy of Optical Remote Sensing Satellites, *Aerospace Return and Remote Sensing*, Vol. 42, No. 1, 2021, pp. 135-144) addresses the impact of optical path difference and atmospheric refraction offset on the planar positioning of optical remote sensing satellites. It obtains the RPC model coefficients based on virtual control point fitting through the conversion between a physical geometric model and a general geometric RPC model, and conducts positioning accuracy comparison and verification analysis. 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) systematically analyzes the differences in atmospheric refraction error in full-spectrum multi-angle imaging at different regions and times using camera parameters and atmospheric parameters from the Gaofen-5 full-spectrum spectral imager and an optical satellite atmospheric refraction correction model based on ellipsoidal layering. Reference 9 (Zang Wenzhi, 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 influence of factors such as line of sight angle of medium and low orbit satellites, atmospheric layer thickness, target altitude and satellite orbital altitude on atmospheric refraction line of sight errors, and proposes an atmospheric refraction correction method for target positioning scenarios of two satellites.

[0005] A method and system for atmospheric refraction compensation in geolocation of optical satellite remote sensing data (CN102346252A) discloses an atmospheric refraction compensation method and system for geolocation of optical satellite remote sensing data. This method includes an Earth radius compensation step, used to calculate the coordinates of a ground point corresponding to a pixel in a satellite image by assuming the actual Earth radius plus an Earth radius compensation amount. A method for all-weather starlight refraction satellite autonomous positioning (CN104236553A) discloses a starlight refraction correction method for satellite observation sidereal time applied to all-weather high-precision astronomical autonomous navigation. A method for target positioning under atmospheric refraction suitable for airborne electro-optical observation systems (CN108535715A) discloses a target positioning method under atmospheric refraction suitable for airborne electro-optical observation systems. Summary of the Invention

[0006] To address the shortcomings of existing technologies, the purpose of this invention is to provide a method and system for atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite images.

[0007] An atmospheric refraction compensation method for geolocation of optical Earth remote sensing satellite imagery according to the present invention includes:

[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 divide the atmosphere into layers according to elevation;

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

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

[0012] Step S5: Calculate the position of the intersection of the lines of sight based on the incident starting point, the incident line 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 numerical elevation corresponding to the intersection point is greater than or equal to the elevation of the intersection point.

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

[0015] If the digital elevation corresponding to the intersection point is less than the intersection point elevation, calculate the line-of-sight vector after refraction based on the refractive index of the incident layer and the refractive index of the exit layer, increment the layer number by 1, and repeat steps S5 to S6.

[0016] Preferably, in step S1:

[0017] Obtain the satellite's position in the Earth-fixed coordinate system ECEF at the corresponding exposure time. Line of sight for each detector

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

[0019] The line-of-sight vectors corresponding to each detector are accurate directions after optical distortion correction and on-board thermal deformation correction.

[0020] Preferably, in step S2:

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

[0022] The number of atmospheric layers is set according to the angle of deviation of satellite observations from the nadir point. The larger the deviation from the nadir point, the larger the number of layers M. Conversely, when the angle of deviation of satellite observations from the nadir point is small, the number of layers is reduced to reduce the amount of computation.

[0023] Preferably, in step S3:

[0024] Calculate the average atmospheric refractive index n of 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 k-th layer k The calculation method is as follows:

[0026]

[0027] Where k = 1, 2, ..., M, the function f(h) represents the variation of atmospheric refractive index with elevation, h is the elevation, dh is the elevation differential, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively.

[0028] The atmospheric refractive index as a function of elevation, f(h), is approximately calculated using the following formula.

[0029]

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

[0031] Preferably, in step S5:

[0032] Based on the point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point

[0033] Step A.1: Iteration starting 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; location Update once;

[0036]

[0037] in, Represent two unit vectors and Calculate the inner product. Indicates position The outward normal direction;

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

[0039] Preferably, in step S6:

[0040] If the elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? k ;

[0041] If the intersection point is The numerical elevation corresponding to latitude and longitude is greater than or equal to h. k Calculate the incident line-of-sight vector Intersection with digital elevation surface Output the coordinates of the intersection point;

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

[0043] Step C.1: If the intersection point is located The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint

[0044]

[0045] Step C.2: Calculate the intermediate point Digital elevation η corresponding to latitude and longitude middle If η middle -h middle If the calculation threshold is not reached, then 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 It is less than the set threshold.

[0046] Preferably, in step S6:

[0047] If the intersection point is The numerical elevation corresponding to latitude and longitude is less than h. k Based on the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction

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

[0049]

[0050] in, express The outward normal direction, Represent two unit vectors and dot product operation;

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

[0052]

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

[0054]

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

[0056] Preferably, except for step S3, where the empirical function of atmospheric refractive index variation with elevation corresponds to an optical spectral band of 0.4–0.9 μm, the other calculation steps have no spectral band restrictions.

[0057] An optical Earth remote sensing satellite image geolocation atmospheric refraction compensation system according to the present invention includes:

[0058] Module M1: Obtains the satellite's position in the Earth-fixed coordinate system at the corresponding exposure time, and the line-of-sight vector for each detector;

[0059] Module M2: Establishes a multi-layer atmospheric refraction model, dividing the atmosphere into layers according to elevation;

[0060] Module M3: Calculates the average atmospheric refractive index for each layer;

[0061] Module M4: Calculations are performed sequentially starting from layer 0;

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

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

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

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

[0066] Preferably, in module M1:

[0067] Obtain the satellite's position in the Earth-fixed coordinate system ECEF at the corresponding exposure time. Line-of-sight vectors corresponding to each detector is a unit vector, where the subscript j indicates the detector number;

[0068] The line-of-sight vectors for each detector are accurate directions after optical distortion correction and on-board thermal deformation correction.

[0069] In module M2:

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

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

[0072] In module M3:

[0073] Calculate the average atmospheric refractive index n of 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 method for calculating the average atmospheric refractive index nk corresponding to the k-th layer is as follows:

[0075]

[0076] Where k = 1, 2, ..., M, the function f(h) represents the variation of atmospheric refractive index with elevation, h is the elevation, dh is the elevation differential, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively.

[0077] The atmospheric refractive index as a function of elevation, f(h), is approximately calculated using the following formula.

[0078]

[0079] Where h is the elevation, in meters. tropop This represents the average elevation of the tropopause.

[0080] In module M5:

[0081] Based on the point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point

[0082] Step A.1: Iteration starting 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; location Update once;

[0085]

[0086] in, Represent two unit vectors and Calculate the inner product. Indicates position The outward normal direction;

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

[0088] In module M6:

[0089] If the elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? k ;

[0090] If the intersection point is The numerical elevation corresponding to latitude and longitude is greater than or equal to h. k Calculate the incident line-of-sight vector Intersection with digital elevation surface Output the coordinates of the intersection point;

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

[0092] Step C.1: If the intersection point is located The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint

[0093]

[0094] Step C.2: Calculate the intermediate point Digital elevation η corresponding to latitude and longitude middle If η middle -h middle If the calculation threshold is not reached, then 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;

[0095] In module M6:

[0096] If the intersection point is The numerical elevation corresponding to latitude and longitude is less than h. k Based on the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction

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

[0098]

[0099] in, express The outward normal direction, Represent two unit vectors and dot product operation;

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

[0101]

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

[0103]

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

[0105] Except for the empirical function calculation of atmospheric refractive index with elevation in module M3, which corresponds to the optical spectral band of 0.4–0.9 μm, the other calculation steps have no spectral band restrictions.

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

[0107] 1. This invention addresses the need for high-precision geographic positioning in high-resolution optical remote sensing. It models atmospheric refraction, which leads to non-rigid geometric positioning models, and constructs a multi-layer atmospheric refraction model that can be adjusted according to requirements. It also provides the process and method for geographic positioning under this model. Combined with 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 atmospheric refraction compensation of optical remote sensing images of the earth. Attached Figure Description

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

[0110] Figure 1 This is a flowchart of the present invention;

[0111] Figure 2 This is a schematic diagram of light refraction under a layered atmospheric model.

[0112] Figure 3 The result of multi-layer atmospheric refraction modeling for this invention is shown in the figure.

[0113] Figure 4 This refers to the difference between the result of geolocation of a low-orbit satellite using the method of this invention and the result without considering atmospheric refraction. Detailed Implementation

[0114] The present invention will be described in detail below with reference to specific embodiments. The following examples 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 skilled in the art, several changes and improvements can be made without departing from the scope of the present invention. These all fall within the scope of protection of the present invention.

[0115] Example 1:

[0116] This invention comprehensively considers the effects of optical distortion correction, on-board thermal deformation correction, atmospheric refraction, and Earth's elevation, aiming to provide a universally applicable atmospheric refraction method that can improve the geolocation accuracy of high-resolution optical Earth remote sensing satellite imagery. This invention is simple to implement, and atmospheric refraction compensation can improve the geolocation accuracy of high-resolution optical Earth remote sensing satellite imagery.

[0117] According to the present invention, a method for atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite imagery is provided, such as... Figure 1 As shown, it includes:

[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] Obtain the satellite's position in the Earth-fixed coordinate system ECEF at the corresponding exposure time. Line-of-sight vectors corresponding to each detector is a unit vector, where the subscript j indicates the detector number;

[0121] The line-of-sight vectors corresponding to each detector are accurate directions after optical distortion correction and on-board thermal deformation correction.

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

[0123] Specifically, in step S2:

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

[0125] The number of atmospheric layers is set according to the angle of deviation of satellite observations from the nadir point. The larger the deviation from the nadir point, the larger the number of layers M. Conversely, when the angle of deviation of satellite observations from the nadir point is small, the number of layers is reduced to reduce the amount of computation.

[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 of 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 k-th layer k The calculation method is as follows:

[0130]

[0131] Where k = 1, 2, ..., M, the function f(h) represents the variation of atmospheric refractive index with elevation, h is the elevation, dh is the elevation differential, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively.

[0132] The atmospheric refractive index as a function of elevation, f(h), is approximately calculated using the following formula.

[0133]

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

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

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

[0137] Specifically, in step S5:

[0138] Based on the point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point

[0139] Step A.1: Iteration starting 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; location Update once;

[0142]

[0143] in, Represent two unit vectors and Calculate the inner product. Indicates position The outward normal direction;

[0144] Step A.4, repeat steps A.2 to A.3 until h t -h k It is 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 numerical elevation corresponding to the intersection point is greater than or equal to the elevation of the intersection point.

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

[0147] If the digital elevation corresponding to the intersection point is less than the intersection point elevation, calculate the line-of-sight vector after refraction based on the refractive index of the incident layer and the refractive index of the exit layer, increment the layer number by 1, and repeat steps S5 to S6.

[0148] Specifically, in step S6:

[0149] If the elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? k ;

[0150] If the intersection point is The numerical elevation corresponding to latitude and longitude is greater than or equal to h. k Calculate the incident line-of-sight vector Intersection with digital elevation surface Output the coordinates of the intersection point;

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

[0152] Step C.1: If the intersection point is located The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint

[0153]

[0154] Step C.2: Calculate the intermediate point Digital elevation η corresponding to latitude and longitude middle If η middle -h middle If the calculation threshold is not reached, then 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 It is less than the set threshold.

[0155] Specifically, in step S6:

[0156] If the intersection point is The numerical elevation corresponding to latitude and longitude is less than h. k Based on the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction

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

[0158]

[0159] in, express The outward normal direction, Represent two unit vectors and dot product operation;

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

[0161]

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

[0163]

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

[0165] Specifically, except for step S3, where the empirical function of atmospheric refractive index variation with elevation corresponds to an optical spectral band of 0.4–0.9 μm, the other calculation steps have no spectral band restrictions.

[0166] Example 2:

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

[0168] The present invention also provides an atmospheric refraction compensation system for geolocation of optical remote sensing satellite images. The atmospheric refraction compensation system for geolocation of optical remote sensing satellite images can be implemented by executing the process steps of the atmospheric refraction compensation method for geolocation of optical remote sensing satellite images. That is, those skilled in the art can understand the atmospheric refraction compensation method for geolocation of optical remote sensing satellite images as a preferred embodiment of the atmospheric refraction compensation system for geolocation of optical remote sensing satellite images.

[0169] An optical Earth remote sensing satellite image geolocation atmospheric refraction compensation system according to the present invention includes:

[0170] Module M1: Obtains the satellite's position in the Earth-fixed coordinate system at the corresponding exposure time, and the line-of-sight vector for each detector;

[0171] Specifically, in module M1:

[0172] Obtain the satellite's position in the Earth-fixed coordinate system ECEF at the corresponding exposure time. Line-of-sight vectors corresponding to each detector is a unit vector, where the subscript j indicates the detector number;

[0173] The line-of-sight vectors for each detector are accurate directions after optical distortion correction and on-board thermal deformation correction.

[0174] In module M2:

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

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

[0177] Module M2: Establishes a multi-layer atmospheric refraction model, dividing the atmosphere into layers according to elevation;

[0178] Module M3: Calculates the average atmospheric refractive index for each layer;

[0179] In module M3:

[0180] Calculate the average atmospheric refractive index n of 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 k-th layer k The calculation method is as follows:

[0182]

[0183] Where k = 1, 2, ..., M, the function f(h) represents the variation of atmospheric refractive index with elevation, h is the elevation, dh is the elevation differential, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively.

[0184] The atmospheric refractive index as a function of elevation, f(h), is approximately calculated using the following formula.

[0185]

[0186] Where h is the elevation, in meters. tropop This represents the average elevation of the tropopause.

[0187] Module M4: Calculations are performed sequentially starting from layer 0;

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

[0189] In module M5:

[0190] Based on the point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point

[0191] Step A.1: Iteration starting 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; location Update once;

[0194]

[0195] in, Represent two unit vectors and Calculate the inner product. Indicates position The outward normal direction;

[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 point is less than the maximum elevation of the area, determine whether the numerical elevation corresponding to the intersection point is greater than or equal to the elevation of the intersection point.

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

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

[0200] In module M6:

[0201] If the elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? k ;

[0202] If the intersection point is The numerical elevation corresponding to latitude and longitude is greater than or equal to h. k Calculate the incident line-of-sight vector Intersection with digital elevation surface Output the coordinates of the intersection point;

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

[0204] Step C.1: If the intersection point is located The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint

[0205]

[0206] Step C.2: Calculate the intermediate point Digital elevation η corresponding to latitude and longitude middle If η middle -h middle If the calculation threshold is not reached, then 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;

[0207] In module M6:

[0208] If the intersection point is The numerical elevation corresponding to latitude and longitude is less than h. k Based on the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction

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

[0210]

[0211] in, express The outward normal direction, Represent two unit vectors and dot product operation;

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

[0213]

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

[0215]

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

[0217] Except for the empirical function calculation of atmospheric refractive index with elevation in module M3, which corresponds to the optical spectral band of 0.4–0.9 μm, the other calculation steps have no spectral band restrictions.

[0218] Example 3:

[0219] Example 3 is a preferred example of Example 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 a method for atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite images, thereby improving the control point-free geometric positioning accuracy of high-resolution optical Earth remote sensing satellite images.

[0221] When geolocating the satellites based on the impact of optical Earth remote sensing satellites, the first step is to calculate the satellite's position in the Earth-fixed coordinate system (ECEF) at the corresponding exposure time based on a rigorous geometric imaging model. Line-of-sight vectors corresponding to each detector (Unit vector, where the subscript j represents the detector number). When atmospheric refraction is neglected, i.e., the position is determined through the Earth-fixed coordinate system. Line of sight vector And the Digital Elevation Model of the Earth (DEM) can calculate pixel geolocation without control points.

[0222] However, with the improvement of spatial resolution in Earth remote sensing, the accuracy requirements for geolocation also increase. Control point-free geolocation is the foundation of accurate geolocation. Currently, many methods incorporate factors such as optical geometric distortion and satellite thermal deformation during the line-of-sight determination process, thereby improving the accuracy of control point-free geolocation to some extent. Building upon this, this invention further optimizes the rigorous geometric imaging model and models atmospheric refraction.

[0223] Atmospheric refraction affects optical curvature and is related to the atmospheric refractive index and the angle of incidence. The atmospheric refractive index is influenced by multiple factors, including spectral wavelength, temperature, humidity, and air pressure. Atmospheric temperature, humidity, and air pressure are all related to altitude, but not solely by it. To simplify the model, this invention approximates the atmospheric refractive index as a function solely dependent on altitude, ignoring differences in the atmosphere at different latitudes and longitudes, and different seasons. These differences will cause variations in the refractive index, but generally the magnitude is small.

[0224] According to Hoyt and Storey's theory, 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 mean temperature lapse rate r is 0.0065 K / m, then the density factor α in the 0.4–0.9 μm spectral range of the troposphere can be approximated as:

[0227]

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

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

[0230]

[0231] Where α tropop T represents the density factor at the tropopause. tropop h represents the absolute temperature of the tropopause. tropop This represents the average elevation of the tropopause.

[0232] Substituting formulas 9 and 10 into formula 8, we get:

[0233]

[0234] This invention constructs a multi-layer atmospheric refraction model, dividing the atmosphere into M layers according to elevation. The elevations corresponding to the tops of each layer from top to bottom are h0, h1, ..., h0, h0, h1, h2, h3, h4, h5, h6, h7, h8, h9, h1, h1, h2, h9 ...1, h2, h2, h9, h1, h2, h1, h2, h2, h1, h2, h2, h1, h2, h2, h1, M-1 The average atmospheric refractive index of each layer is represented by n. k Let n represent the average atmospheric refractive index n of the k-th layer (k = 1, 2, ..., M), where k represents the layer number, k = 0, 1, 2, ..., M, and n0 is the refractive index in vacuum, n0 = 1. k The calculation method is as follows

[0235]

[0236] Wherein, the function f(h) represents the function of atmospheric refractive index as a function of elevation, h is the elevation, dh is the elevation derivative, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively.

[0237] Under the same conditions, the greater the incident angle, the greater the light refraction. Therefore, the number of atmospheric layers can be set according to the angle of deviation of satellite observation from the nadir. When the deviation from the nadir is greater, the number of layers M is greater. Conversely, when the angle of deviation of satellite observation from the nadir is smaller, the number of layers can be reduced to reduce the amount of computation.

[0238] Outside the atmosphere, the refractive index can be considered to be 1, and light travels in a straight line. In a layered atmosphere, refraction occurs between layers, while light travels in a straight line within a single layer. (According to the appendix...) Figure 2 The schematic diagram shown starts from k=0 and proceeds sequentially from the refractive index n of the incident layer. k , refractive index n of the exit layer k+1 Point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point With the refracted line of sight vector

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

[0240] If the current h t -h k Less than the calculation threshold, currently The output is the coordinates of the intersection point. If the calculation threshold is not reached, the position... Iterative updates:

[0241]

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

[0243] When calculating to lower altitudes, the ground elevation must be considered. If the current layer elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? kIf yes, calculate the intersection point with the surface of the digital elevation model; otherwise, use the already calculated intersection points. and the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction

[0244] First, calculate the incident angle θ. k .

[0245]

[0246] in, express The outward normal direction, Represent two unit vectors and Inner product operation.

[0247] Calculate the angle of incidence using Snell's law.

[0248]

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

[0250]

[0251] Calculate the incident line-of-sight vector Intersection with digital elevation surface It is also done iteratively, if the intersection point position The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint

[0252]

[0253] Calculate the midpoint Digital elevation η corresponding to latitude and longitude middle If η middle -h middle If the calculation threshold is not reached, then if η middle h middle ,Will Updated to If η middle <h middle ,Will Updated to Repeat the iteration until η middle -h middleThe elevation is less than the set threshold. Two-dimensional interpolation is required when calculating the elevation of any point in the digital elevation model (DEM) within the Earth-fixed coordinate system. The DEM only provides elevation data for fixed grid points; the elevation data for any point can be obtained through bilinear interpolation.

[0254] Based on the reversibility of the optical path, the position can be determined. The corresponding light source, after passing through the atmosphere, reaches the corresponding detector, which is the location corresponding to the detector's geolocation result.

[0255] The effectiveness of the method of this invention is illustrated below using data from a low-orbit Earth-based optical remote sensing satellite. The atmospheric refractive index at different elevation locations obtained according to the method of this invention is shown in the attached figure. Figure 3 As shown by the solid line, after performing step 2 to stratify the atmosphere at every 1km elevation, the average atmospheric refractive index of each layer is shown. Figure 3 As shown by discrete points. The satellite's position at a certain moment is... (Unit: m), the line-of-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 nadir point. Combined with DEM data, and disregarding atmospheric refraction, the resulting geographic location is: latitude 30.532304°, longitude 95.181367°, and elevation 4882.2m. The final geographic location obtained using the method of this invention is: latitude 30.532274°, longitude 95.182102°, and elevation 4880.2m. Considering atmospheric refraction, the differences in positioning results at different elevations are shown in the attached figure. Figure 4 As shown, in this embodiment, under a large side-view angle, atmospheric refraction causes an error of nearly 100 meters.

[0256] Those skilled in the art will appreciate that, in addition to implementing the system and its various devices, modules, and units provided by the present invention in purely computer-readable program code, it is entirely possible to implement the same functions of the system and its various devices, modules, and units provided by the present invention 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 a hardware component, and the devices, modules, and units included therein for implementing various functions can also be considered as structures within the hardware component; the devices, modules, and units for implementing various functions can also be considered as both software modules implementing the method and structures within the hardware component.

[0257] The above describes specific embodiments of the present invention. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art may make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. The embodiments of this application and the features in the embodiments may be combined with each other in any manner unless there is a conflict.

Claims

1. A method for atmospheric refraction compensation in optical Earth remote sensing satellite image geolocation, 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 of each detector; Step S2: Establish a multi-layer atmospheric refraction model and divide the atmosphere into layers according to elevation; Step S3: Calculate the average atmospheric refractive index of each layer; Step S4: Perform calculations sequentially starting from layer 0; Step S5: Calculate the position of the intersection of the lines of sight based on the point of incidence, the vector of the line of sight, and the elevation; Step S6: If the elevation of the intersection point is less than the maximum elevation of the area, determine whether the numerical elevation corresponding to the intersection point is greater than or equal to the elevation of the intersection point. If the digital elevation corresponding to the intersection point is greater than or equal to the intersection point elevation, calculate the intersection point of the incident line-of-sight vector and the digital elevation surface, and output the coordinates of the intersection point; If the digital elevation corresponding to the intersection point is less than the intersection point elevation, calculate the line-of-sight vector after refraction based on the refractive index of the incident layer and the refractive index of the exit layer, increment the layer number by 1, and repeat steps S5 to S6.

2. The method for atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that, In step S1: Obtain the satellite's position in the Earth-fixed coordinate system ECEF at the corresponding exposure time. Line-of-sight vectors corresponding to each detector is a unit vector, where the subscript j indicates the detector number; The line-of-sight vectors corresponding to each detector are accurate directions after optical distortion correction and on-board thermal deformation correction.

3. The atmospheric refraction compensation method for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that, In step S2: A multi-layer atmospheric refraction model is established, dividing the atmosphere into M layers according to elevation. The elevations corresponding to the tops of each layer from top to bottom are h0, h1, ..., h. M-1 ; The number of atmospheric layers is set according to the angle of deviation of satellite observations from the nadir point. The larger the deviation from the nadir point, the larger the number of layers M. Conversely, when the angle of deviation of satellite observations from the nadir point is small, the number of layers is reduced to reduce the amount of computation.

4. The atmospheric refraction compensation method for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that, In step S3: Calculate the average atmospheric refractive index n of 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 k-th layer k The calculation method is as follows: Where k = 1, 2, ..., M, the function f(h) represents the variation of atmospheric refractive index with elevation, h is the elevation, dh is the elevation differential, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively. The atmospheric refractive index as a function of elevation, f(h), is approximately calculated using the following formula. Where h is the elevation, in meters. tropop This represents the average elevation of the tropopause.

5. The method for atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that, In step S5: Based on the point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point Step A.1: Iteration starting 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 has not been reached; location... Update once; in, Represent two unit vectors and Calculate the inner product. Indicates position The outward normal direction; Step A.4, repeat steps A.2 to A.3 until |h t -h k | Less than the set threshold.

6. The method for atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that, In step S6: If the elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? k ; If the intersection point is The numerical elevation corresponding to latitude and longitude is greater than or equal to h. k Calculate the incident line-of-sight vector Intersection with digital elevation surface Output the coordinates of the intersection point; Incident line-of-sight vector Intersection with digital elevation surface The calculation includes the following steps: Step C.1: If the intersection point is located The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint Step C.2: Calculate the intermediate point Digital elevation η corresponding to latitude and longitude middle If |η middle -h middle |The calculation threshold has not been reached; 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 atmospheric refraction compensation for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that, In step S6: If the intersection point is The numerical elevation corresponding to latitude and longitude is less than h. k Based on the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction Step B.1: Calculate the incident angle θ k : in, express The outward normal direction, Represent two unit vectors and dot product operation; Step B.2, calculate the angle of incidence. Step B.3, calculate the line-of-sight vector after refraction. Set k = k + 1 and repeat steps S5 to S6.

8. The atmospheric refraction compensation method for geolocation of optical Earth remote sensing satellite imagery according to claim 1, characterized in that: Except for step S3, where the empirical function of atmospheric refractive index variation with elevation corresponds to an optical spectral band of 0.4–0.9 μm, the other calculation steps have no spectral band restrictions.

9. An atmospheric refraction compensation system for geolocation of optical Earth remote sensing satellite imagery, characterized in that, include: Module M1: Obtains the satellite's position in the Earth-fixed coordinate system at the corresponding exposure time, and the line-of-sight vector for each detector; Module M2: Establishes a multi-layer atmospheric refraction model, dividing the atmosphere into layers according to elevation; Module M3: Calculates the average atmospheric refractive index for each layer; Module M4: Calculations are performed sequentially starting from layer 0; Module M5: Calculates the position of the intersection of the lines of sight based on the point of incidence, the vector of the incident line of sight, and the elevation; Module M6: If the elevation of the intersection point is less than the maximum elevation of the area, determine whether the numerical elevation corresponding to the intersection point is greater than or equal to the elevation of the intersection point. If the digital elevation corresponding to the intersection point is greater than or equal to the intersection point elevation, calculate the intersection point of the incident line-of-sight vector and the digital elevation surface, and output the coordinates of the intersection point; If the digital elevation corresponding to the intersection point is less than the intersection point elevation, calculate the line-of-sight vector after refraction based on the refractive index of the incident layer and the refractive index of the exit layer, and increment the layer number by 1 to repeat modules M5 to M6.

10. The optical Earth remote sensing satellite image geolocation atmospheric refraction compensation system according to claim 9, characterized in that: In module M1: Obtain the satellite's position in the Earth-fixed coordinate system ECEF at the corresponding exposure time. Line-of-sight vectors corresponding to each detector is a unit vector, where the subscript j indicates the detector number; The line-of-sight vectors for each detector are accurate directions after optical distortion correction and on-board thermal deformation correction. In module M2: A multi-layer atmospheric refraction model is established, dividing the atmosphere into M layers according to elevation. The elevations corresponding to the tops of each layer from top to bottom are h0, h1, ..., h. M-1 ; The number of atmospheric layers is set according to the angle of deviation of satellite observation from the nadir point. When the deviation from the nadir point is greater, the number of layers M is greater. Conversely, when the angle of deviation of satellite observation from the nadir point is smaller, the number of layers is reduced to reduce the amount of computation. In module M3: Calculate the average atmospheric refractive index n of 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 k-th layer k The calculation method is as follows: Where k = 1, 2, ..., M, the function f(h) represents the variation of atmospheric refractive index with elevation, h is the elevation, dh is the elevation differential, and h k-1 h k These are the elevations corresponding to the top and bottom of the k-th atmospheric layer, respectively. The atmospheric refractive index as a function of elevation, f(h), is approximately calculated using the following formula. Where h is the elevation, in meters. tropop This represents the average elevation of the tropopause. In module M5: Based on the point of incidence Incident line-of-sight vector Elevation h k Calculate the intersection point Step A.1: Iteration starting 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 has not been reached; location... Update once; in, Represent two unit vectors and Calculate the inner product. Indicates position The outward normal direction; Step A.4, repeat steps A.2 to A.3 until |h t -h k | Less than the set threshold; In module M6: If the elevation h k If the elevation is less than the maximum value of the area, determine the location of the intersection point. Are the numerical elevations corresponding to latitude and longitude greater than or equal to h? k ; If the intersection point is The numerical elevation corresponding to latitude and longitude is greater than or equal to h. k Calculate the incident line-of-sight vector Intersection with digital elevation surface Output the coordinates of the intersection point; Incident line-of-sight vector Intersection with digital elevation surface The calculation includes the following steps: Step C.1: If the intersection point is located The digital elevation corresponding to latitude and longitude is η, derived from the incident line-of-sight vector. The starting point Intersection position Calculate the midpoint Step C.2: Calculate the intermediate point Digital elevation η corresponding to latitude and longitude middle If |η middle -h middle If the calculation threshold is not reached, then 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 module M6: If the intersection point is The numerical elevation corresponding to latitude and longitude is less than h. k Based on the refractive index n of the incident layer k , refractive index n of the exit layer k+1 Calculate the line-of-sight vector after refraction Step B.1: Calculate the incident angle θ k : in, express The outward normal direction, Represent two unit vectors and dot product operation; Step B.2, calculate the angle of incidence. Step B.3, calculate the line-of-sight vector after refraction. Set k = k + 1, and repeat modules M5 to M6; Except for the empirical function calculation of atmospheric refractive index with elevation in module M3, which corresponds to the optical spectral band of 0.4–0.9 μm, the other calculation steps have no spectral band restrictions.

Citation Information

Patent Citations

  • Method and system for compensating atmospheric refraction in optical satellite remote sensing data geographic positioning

    CN102346252A

  • Autonomous all-weather stellar refraction satellite location method

    CN104236553A

  • Target positioning method applicable to vehicle-mounted photoelectric watching-aiming system with atmospheric refraction

    CN108535715A

  • Ellipsoidal layered atmospheric refraction-based optical satellite image accurate ground positioning method

    CN110046430A

  • Atmospheric refraction positioning error correction method for optical remote sensing satellite image in Qinghai-Tibet plateau region

    CN113960642A