A Rapid Geolocation Algorithm for Three-axis Stabilization Geostationary Remote Sensing Satellite

1. Introduction

Geostationary orbit has become a contested strategic resource, with major powers vying for control. Geostationary orbit remote sensing satellites possess the unique advantage of being able to "see far and stay high," making them optimal for continuous and dynamic observation of targets. This has led to an increasing number of geostationary orbit satellite applications being used worldwide. Some notable examples include the US GOES-R [1], Japan's "Himawari 8/9 [2]", China's "GF-4 [3]", "Fengyun-4A [4]", and "Fengyun-4B [5]". Furthermore, the third-generation imaging satellite MTG-I1 of the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) was successfully launched on December 13, 2022.

Georeferencing of geostationary orbit remote sensing images is the cornerstone of subsequent applications of remote sensing data. Georeferencing algorithms for satellite images can be classified into non-parametric and parametric methods [6]. Non-parametric methods can be applied to all types of geometric distortions but require multiple ground control points to establish a mapping relationship between remote sensing data and geographic coordinates. However, they are unsuitable for georeferencing automated operational geostationary satellite imagery, such as the Fengyun series, due to the time- and resource-intensive process of selecting high-precision ground control points. Moreover, it remains challenging to identify usable ground control points in cloudy conditions, which fails to meet the requirements of high timeliness and stability for meteorological satellite operations [7,8]. In contrast, parametric methods establish models based on the observation geometry, spatial position, and pointing of satellite sensors to calculate the georeferencing information of remote sensing data [9]. Ground control points are used to eliminate georeferencing errors caused by uncertainties in model param [10,11,12], and terrain information is utilized to remove georeferencing errors caused by surface topography. These errors depend on satellite altitude, terrain, and the distance between observation and nadir points [13]. This paper focuses on parametric methods for georeferencing, particularly those tailored for geostationary orbit meteorological satellite imagery.

Geostationary orbit remote sensing satellites have two primary modes of observation: full disk and regional. Full-disk observation is primarily utilized to support the tracking of large-scale weather systems, numerical weather prediction applications, and climate dataset construction. Regional observation is mainly conducted for weather systems on a scale of 1,000–2,000 km, specifically to monitor small to medium-scale convective systems and typhoons, providing weather analysis and warning services [14].

As shown in Table 1, the current full disk observation time efficiency is 5 min for the US GOES-R series, Japan's Himawari-8/9 satellites, 10 min for MTG-I1, and 15 min for Fengyun-4A and Fengyun-4B satellites [15,16]. In terms of medium-scale regional scan time, Himawari-8/9 requires 2.5 min for typhoon-scale regional observation, Fengyun-4A requires 3 min for typhoon-scale regional observation, and Fengyun-4B's rapid imaging instrument requires 1 min for a 2000×2000 km regional observation. In terms of spatial resolution, GEOS-R, Himawari-8/9, and Fengyun-4A are similar, while Fengyun-4B can achieve the highest resolution of 0.25 km. Remote sensing applications such as weather forecasting and environmental monitoring increasingly demand real-time processing of remote sensing data, even at a minute level. With the improvement of the spatial resolution of remote sensing satellites, the volume of acquired data increases significantly, requiring higher efficiency in data processing and georeferencing calculations. For example, the spatial resolution of Fengyun-4B's rapid imaging instrument is 250 m, which is twice the resolution of existing imagers. Efficient georeferencing algorithms can meet the requirements of regional observation but cannot satisfy the needs of full-disk observation with rapid imaging instruments. Therefore, to achieve efficient processing of geostationary orbit remote sensing satellite data, it is crucial to research more efficient georeferencing algorithms that have practical significance and economic value in addition to hardware improvements.

This paper centers on the param-based geolocation process commonly utilized for three-axis stabilized geostationary orbit remote sensing satellites. Based on our analysis of the standard data release format of geostationary orbit satellite products, we propose a rapid geolocation algorithm and confirm its feasibility through experiments.

2. Algorithm Principle

2.1. Conventional Geolocation Process Based on Intersection of Viewing Vectors

Three-axis stabilized geostationary orbit remote sensing satellites often adopt a dual scanner drive mode, where scanning mirrors are driven in an east-west direction for cross-track scanning and in a north-south direction for along-track stepping, enabling large-area coverage. Figure 1 presents the scanning mechanism and ideal optical path of the satellite payload.

The geolocation algorithm utilizes the scanning mirror angles from the sensor's Level 0 data to computationally determine the corresponding ground latitude and longitude. This conventional geolocation algorithm commences by modeling the instrument's observation geometry to acquire the initial viewing vectors for each pixel. Through a series of coordinate transformations, the image-plane coordinates are projected onto the Earth-fixed coordinate system. The algorithm then computes the intersection coordinates of the sensor's viewing vector with the ground, which are subsequently transformed into latitude and longitude. Finally, these coordinates are mapped onto the coordination group for meteorological satellites (CGMS) grid (a nominal grid with a specific projection defined by the coordination group for meteorological satellites), one of the international standards for meteorological remote sensing. The specific coordinate system transformations involved in this process are graphically depicted in Figure 2.

The initial viewing vector, denoted as $\stackrel{⃑}{L_0}$, varies with the changes in the angles of the two scanning mirrors. Assuming the outgoing ray vector in the payload coordinate system is ${\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{in}$, the resulting viewing vector after two reflections from the scanning mirrors, denoted as $\stackrel{⃑}{L_0}$, was expressed as:

$${\stackrel{⃑}{L_0}=\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{in}=F_1·F_2·{\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_1$$

(1)

Here, 𝐹1 and 𝐹2 are the reflection matrices for the north-south and east-west mirrors, respectively. These matrices are determined by the angles of rotation 𝛼 for the north-south mirror and 𝛽 for the east-west mirror. The specific derivation process can be found in references [17,18,19].

$${\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{in}=F_1·F_2·{\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_1=F_1·F_2·\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)=\left(\begin{array}{c}-\mathrm{c}\mathrm{o}\mathrm{s}\left(2\alpha \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(2\beta \right)\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left(2\alpha \right)\\ \mathrm{c}\mathrm{o}\mathrm{s}\left(2\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(2\beta \right)\end{array}\right)$$

(2)

After the viewing vector exits the satellite body, it is transformed from the satellite body coordinate system to the Earth-fixed coordinate system based on the satellite attitude information, orbit information, and param such as nutation, precession, Greenwich sidereal time, and polar motion. For a given time $t$, the geolocation process for the viewing vector $\stackrel{⃑}{L_0}$ in the payload coordinate system can be divided into the following eight steps:

1) calculating the viewing vector $\stackrel{⃑}{L_0}$ in the payload coordinate system based on the mirror angles and the optical path geometry ;

2) transforming the viewing vector $\stackrel{⃑}{L_0}$ from the payload coordinate system to the body coordinate system by using the installation matrix, obtaining $\stackrel{⃑}{L_1}$;

3) based on the orbit attitude, transforming the viewing vector from the body coordinate system to the orbit coordinate system, obtaining $\stackrel{⃑}{L_2}$ .The satellite body coordinate system was rotated in the order of $\mathsf{\phi }，\mathsf{\theta }$ and $\mathsf{\omega }$ to obtain the orbit coordinate system.

$$\stackrel{⃑}{L_2}=T_{orb}\stackrel{⃑}{L_1}=T_z\left(-\phi \right)T_x\left(-\theta \right)T_y\left(-\omega \right)\stackrel{⃑}{L_1}$$

(3)

4) based on the satellite's position and velocity, transforming the viewing vector from the orbit coordinate system to the Earth-centered inertial coordinate system.

The transformation from the orbit coordinate system to the Earth-centered inertial coordinate system is determined by the satellite's instantaneous inertial frame of reference, which considers the position and velocity. This transformation can be obtained by constructing the orbit coordinate system in the inertial frame of reference. The Earth-centered inertial coordinate system was $\left(X_{eci},Y_{eci},Z_{eci}\right)$with unit vectors as ${\overrightarrow{u}}_{eci}=\left({\overrightarrow{u}}_x,{\overrightarrow{u}}_y,{\overrightarrow{u}}_z\right)$. The orbit coordinate system was$\left(X_{orb},Y_{orb},Z_{orb}\right)$with unit vector $\stackrel{⃑}{{\mathrm{L}}_2}=\left({\overrightarrow{\mathrm{w}}}_{\mathrm{x}},{\overrightarrow{\mathrm{w}}}_{\mathrm{y}},{\overrightarrow{\mathrm{w}}}_{\mathrm{z}}\right)$For any given satellite position vector$\overrightarrow{P}\left(x,y,z\right)$ in the Earth-centered inertial coordinate system, the coordinates of the orbit coordinate system's axes were obtained using the following formula:

$${\overrightarrow{u}}_{eci}=T_{ec/orb}\stackrel{⃑}{L_2}=\left[\begin{array}{ccc}{\overrightarrow{w}}_x& {\overrightarrow{w}}_y& {\overrightarrow{w}}_z\end{array}\right]·\stackrel{⃑}{L_2}$$

(4)

According to the definition of the orbit coordinate system, we derived:

$${\overrightarrow{w}}_z=-\frac{\overrightarrow{P}}{‖\overrightarrow{P}‖}$$

(5)

The unit vector ${\overrightarrow{w}}_y$, which is orthogonal to both ${\overrightarrow{w}}_z$ and the satellite's position vector $\stackrel{⃑}{r}$, were determined by the following equation:

$${\overrightarrow{w}}_y=-\frac{{\overrightarrow{v}}_p\times {\overrightarrow{w}}_z}{‖{\overrightarrow{v}}_p\times {\overrightarrow{w}}_z‖}$$

(6)

where ${\overrightarrow{v}}_p$ is the velocity vector in the Earth-centered inertial coordinate system. Therefore, the unit vector ${\overrightarrow{w}}_x$ was determined by the following equation:

$${\overrightarrow{w}}_x={\overrightarrow{w}}_y\times {\overrightarrow{w}}_z$$

(7)

By performing the calculations mentioned above, the transformation matrix from the orbit coordinate system to the Earth-centered inertial coordinate system was obtained as follows:

$$T_{ec/orb}=\left[\begin{array}{ccc}{\overrightarrow{w}}_x& {\overrightarrow{w}}_y& {\overrightarrow{w}}_z\end{array}\right]$$

(8)

5) transforming the viewing vector from the Earth-centered inertial coordinate system to the Earth-fixed coordinate system based on time information.

This step mainly relies on information such as precession, nutation, Greenwich sidereal time, and polar motion to transform the viewing vector and satellite position vector from the Earth-centered inertial coordinate system to the Earth-fixed coordinate system. As a result, the viewing vector in the Earth-fixed coordinate system was obtained as follows:

$${\stackrel{⃑}{u}}_{WGS}=T_{WGS}{\overrightarrow{u}}_{eci}=T_4\left(t\right)T_3\left(t\right)T_2\left(t\right)T_1\left(t\right){\overrightarrow{u}}_{eci}$$

(9)

Satellite position vector in the Earth-fixed coordinate system was calculated as:

$$r_{WGS}=T_{WGS}{r}_i=T_4\left(t\right)T_3\left(t\right)T_2\left(t\right)T_1\left(t\right)\stackrel{⃑}{r}$$

(10)

Here, the subscript WGS represents the Earth-fixed coordinate system. $T_1\left(t\right)$ denotes the precession correction matrix, $T_2\left(t\right)$ represents the nutation correction matrix, $T_3\left(t\right)$ signifies the Greenwich sidereal time rotation matrix, $T_4\left(t\right)$ denotes the polar motion correction matrix, and $T_{WGS}$ represents the direction cosine matrix from the Earth-centered inertial coordinate system to the Earth-fixed coordinate system.

6) computing the position vector of the intersection point between the viewing vector and the WGS84 ellipsoid in the Earth-fixed coordinate system.

In the Earth-fixed coordinate system, the coordinates of the intersection point between the viewing vector and the surface of the Earth represented by the WGS-84 geodetic reference ellipsoid were calculated. The equation for the WGS-84 geodetic reference ellipsoid was represented as:

$$\frac{{\mathrm{x}}^2+{\mathrm{y}}^2}{{\mathrm{a}}^2}+\frac{{\mathrm{z}}^2}{{\mathrm{b}}^2}=1$$

(11)

where $a$ represents the major axis and $b$ represents the minor axis. The intersection point $P_{WGS}$ between the viewing vector and the Earth's surface was calculated using formula (12) as follows:

$$P_{WGS}=r_{WGS}+s·V_{WGS}$$

(12)

$s$ represents the distance between the satellite and the intersection point:

$$s=\frac{-\left(V^{\prime }{r}^{\prime }\right)-\sqrt{{\left(V^{\prime }{r}^{\prime }\right)}^2-{|V\prime |}^2(\left|r^{\prime }\right|-1)}}{{|V\prime |}^2},V^{\prime }=\left[\begin{array}{c}{V}_{WGSx}/a\\ V_{WGSy}/a\\ V_{WGSz}/b\end{array}\right],r^{\prime }=\left[\begin{array}{c}{r}_{WGSx}/a\\ r_{WGSy}/a\\ r_{WGSz}/b\end{array}\right]$$

7) converting the coordinates of the intersection point in the Earth-fixed coordinate system to the coordinates in the geodetic coordinate system.

Based on the coordinates of the intersection point in the Earth-fixed coordinate system, the geographic latitude and longitude of the intersection point were calculated using formula (13) as follows:

$$lat={tan}^{-1}\frac{P_{WGSz}}{\sqrt{P_{WGSx}^2+P_{WGSy}^2}},lon={2tan}^{-1}\frac{P_{WGSy}}{P_{WGSx}+\sqrt{P_{WGSx}^2+P_{WGSy}^2}}$$

(13)

8) Based on the latitude and longitude of the intersection point, its row and column number (I,J) were calculated on the CGMS grid.

By employing the latitude and longitude coordinates, in conjunction with the north-south scan angle, east-west scan angle, satellite's geolocation, and the Earth model employed, it is possible to computationally determine the row and column number of the point on the nominal grid. For specific details regarding the calculation methodology, please refer to reference [16].

The computation time and proportion of each step for the conventional geolocation algorithm were computed using simulated data. The results are presented in Table 2.

2.2. Rapid Geolocation Algorithm for CGMS Nominal Grid

The nominal data format is the most commonly employed format for applying and disseminating geostationary remote sensing data. For instance, China's FY-2 and FY-4 L1A level data, along with the United States' GOES-R satellite L1B level data, are all published in nominal data formats [14]. The nominal grid comprises fixed observation angles with equal intervals, ensuring that identical data points in all products correspond to the same location on Earth. This provides a basis for further applications of remote sensing imagery. Currently, most geostationary meteorological satellites employ image geolocation methods based on the nominal grid. This allows each pixel in the remote sensing image to be associated with a specific latitude and longitude on Earth. Therefore, there is no need to store latitude and longitude information within the image data itself, significantly reducing the requirements for data transmission and storage.

The geolocation algorithm has a strict correspondence between the initial ray vectors of the remote sensing sensor's pixels, latitude, longitude, and the CGMS nominal grid. The row-column numbers (I, J) of the nominal grid have a unique correspondence with the geographic latitude and longitude (B, L). Therefore, during data geolocation, it is possible to first calculate the row-column numbers (I, J) of the ray vectors on the nominal grid and then obtain the corresponding geographic latitude and longitude (B, L) through a quick query. The ultimate goal of geolocation is to compute the position of the ray on the CGMS nominal grid. By analyzing the transformation process of the ray vector in different coordinate systems, this paper proposes a rapid algorithm that eliminates the time-consuming process of computing latitude and longitude. Instead, it directly establishes the correspondence between the ray vectors of the remote sensing sensor's pixels in the satellite body coordinate system and the CGMS nominal grid's row-column numbers (I, J), thus achieving the final geolocation.

After converting the viewing vector from the Earth-centered inertial coordinate system to the Earth-fixed coordinate system, errors caused by attitude, orbit, precession, nutation, Greenwich sidereal time, and polar motion are eliminated. It can be considered an idealized viewing vector in the satellite body coordinate system.

At this point, the viewing vector in the Earth-fixed coordinate system can be directly transformed into the equivalent ideal state in the satellite body coordinate system. The representation of this viewing vector was given as:

$$\stackrel{⃑}{L_1^{\prime }}={\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{ideal}=T_5·{\stackrel{⃑}{u}}_{WGS}$$

(14)

where $T_5$ represents the transformation matrix from the Earth-fixed coordinate system to the satellite body coordinate system. The specific method for numerical conversion between the geostationary sensor imaging grid and the CGMS grid can be found in reference [16,17].

Based on the viewing vector in the satellite body coordinate system, the relationship between the corresponding CGMS mirror angles $\left({\alpha }^{\prime },{\beta }^{\prime }\right)$ was established as follows:

$${\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)}_{ideal}=\left(\begin{array}{c}-\mathrm{c}\mathrm{o}\mathrm{s}\left(2{\alpha }^{\prime }\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(2{\beta }^{\prime }\right)\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left(2{\alpha }^{\prime }\right)\\ \mathrm{c}\mathrm{o}\mathrm{s}\left(2{\alpha }^{\prime }\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(2{\beta }^{\prime }\right)\end{array}\right)$$

(15)

The $\left({\alpha }^{\prime }，{\beta }^{\prime }\right)$ values were calculated as follows:

$$\left\{\begin{array}{c}{\alpha }^{\prime }=-\frac{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(X_{ideal}\right)}{2}\\ {\beta }^{\prime }=\frac{arc\left(\raisebox{1ex}{$-Y_{ideal}$}\!\left/ \!\raisebox{-1ex}{$Z_{ideal}$}\right.\right)}{2}\end{array}\right.$$

(16)

Subsequently, the corresponding $\left(I，J\right)$ for $\left({\alpha }^{\prime }，{\beta }^{\prime }\right)$ were calculated as follows:

$$\left\{\begin{array}{c}I=Int\left(\frac{m}{2}-\frac{\alpha }{∆\theta }+1\right)\\ J=Int\left(\frac{\beta }{∆\theta }+\frac{n}{2}+1\right)\end{array}\right.$$

(17)

where $m\times n$ represents the size of the nominal grid, and $∆\theta$ is the angular size of a pixel in this band.

Through the aforementioned refinements, the computationally burdensome process of finding the intersection point between the viewing vector and the Earth's ellipsoid was obviated. Instead, a direct mapping relationship was established between the original initial vector and the nominal grid. In the geolocation process, after transforming the ray vector into the Earth-fixed coordinate system, the positioning result was obtained directly through the conversion relationship between the ray vectors in the sensor imaging grid and the CGMS grid, without the need for iterative calculations to solve for the intersection point between the ray and the Earth's ellipsoidal surface. The flowchart of the improved rapid geolocation algorithm can be found in Figure 3.

Figure 3 illustrates that in the rapid algorithm flow, the calculations prior to the Earth-fixed coordinate system are consistent with the conventional algorithm flow. In the Earth-fixed coordinate system, the viewing vector no longer intersects with the geodetic ellipsoid. Instead, the viewing vector in the Earth-fixed coordinate system is transformed into the ideal state in the satellite body coordinate system. This is equivalent to the viewing vector of the actual satellite's downlink data being transformed into the idealized viewing vector (i.e., the corrected viewing vector after considering various error factors) through positioning and compensation calculations. Then, the mirror angles $\left({\alpha }^{\prime }，{\beta }^{\prime }\right)$ corresponding to the idealized viewing vector are calculated. By using the relationship between the angles and the initial position, along with the angle step size, the corresponding row-column numbers $\left(I，J\right)$ can be computed.

3. Experiment and Analysis

To validate the feasibility of the proposed algorithm, programming experiments were conducted using C++. The efficiency and accuracy of the geolocation calculations were compared. The computer environment used was an 11th Gen Intel(R) Core(TM) i7-1165G7@2.80GHz processor with 16 GB of memory, running on a 64-bit Windows 10 operating system.

3.1. Efficiency Comparison

To compare the computational efficiency of each step, we took a sample of geospatial observation data with a spatial resolution of 250 m. We simulated 10 million uniformly distributed nominal grid points within the full disk range, corresponding to 10 million sets of mirror angles. Employing both algorithms, we performed geolocation calculations and recorded the time consumption for each step, as shown in Table 2.

The percentage increase in computational efficiency ($∆r$) can be calculated as follows:

$$∆r=\frac{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$24143$}\right.-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$44260$}\right.}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$24143$}\right.}\times 100\mathrm{\%}=83.3\mathrm{\%}$$

(17)

We observed that the rapid geolocation algorithm eliminated the time-consuming Earth-fixed to Geodetic coordinate transformation process compared to the conventional geolocation algorithm. This results in an improvement in computational efficiency of 83.3%.

3.2. Accuracy Comparison of Algorithms

To verify the accuracy of the two algorithms, representative regions were selected for statistical analysis of the geolocation results. These regions include the satellite subpoint area, the region towards the western edge of the Earth, and the region towards the northern edge of the Earth. The geolocation results in these three regions are presented in magnified display maps in Figure 4, Figure 5 and Figure 6. The line segment in the middle represents the east-west and north-south directions of the satellite coverage. The dots represent the IJ positions of the rapid algorithm on the CGMS nominal grid, while the circles represent the IJ positions of the conventional algorithm on the CGMS nominal grid. The difference between the positions of the dots and circles represents the difference in calculation results for the viewing vector corresponding to the same set of mirror angles. In the satellite subpoint area, the dots and circles coincide very well, but there is a slight difference as we move towards the edges of the Earth. To analyze the magnitude of this difference, the differences between the positions of the dots and circles in the three regions were statistically analyzed, as shown in Table 3. It can be observed that the mean difference in the geolocation results between the two algorithms is in the order of 10^-5 pixels in the satellite subpoint area and in the order of 10^-2 pixels in the regions towards the western and northern edges of the Earth. Therefore, it can be concluded that the computational accuracy of the two geolocation algorithms is comparable.

Table 3 lists the statistical values of the geolocation differences in the three test regions.

In summary, through simulation experiments, the proposed rapid geolocation algorithm for CGMS has been verified to achieve an 83.3% improvement in computational efficiency compared to that observed in the conventional geolocation algorithm without any loss in computational accuracy.

4. Conclusion

With the development of geospatial observation technology in geostationary orbit remote sensing, enhancing the timeliness of remote sensing data processing is crucial for improving practical and economic value. In this paper, by analyzing the conventional geolocation methods for geostationary orbit remote sensing satellite data and the data format of CGMS, we propose a fast geolocation algorithm for three-axis stabilized geostationary scanning remote sensing satellite images, specifically designed for CGMS nominal grids. This algorithm eliminates the computationally intensive process of calculating the intersection points between the viewing vectors and the Earth ellipsoid in the Earth-fixed coordinate system, thus improving the efficiency of geolocation calculations. Through experimental analysis with simulated data, the proposed rapid geolocation algorithm achieves improved timeliness in the geolocation processing of geostationary orbit remote sensing data without sacrificing accuracy when compared to conventional geolocation methods. This research provides a reference for efficient geolocation processing of geospatial observation data in geostationary orbit remote sensing.