Fast Determination of Satellite-to-Moon Visibility Using an Adaptive Interpolation Method Based on Vertex Protection

Abstract

Fast determination of satellite visibility with respect to a target area is important for satellite navigation and positioning. In this paper, we propose an adaptive interpolation algorithm based on vertex protection to solve the satellite visibility period problem more accurately and quickly, where “vertex” refers to the local extremum point. The algorithm can avoid the error in the visibility period calculation caused by skimming the vertices when fitting the multi-hump visibility function under certain fitting accuracy requirements with the traditional adaptive interpolation method. The algorithm does not need to construct a cubic polynomial in each subinterval to determine whether the satellite is visible or not; it only constructs a cubic polynomial to solve the problem if the visibility function of that subinterval is judged to have a solution from the existence theorem of zero points, which can improve the computational efficiency. For the lunar navigation problem, a solution to satellite–Moon visibility calculations based on a vertex-protected adaptive interpolation is given, and the experimental results show that the computation time of the algorithm can be reduced by approximately 98% compared with the brute force method and by approximately 30% compared with the traditional adaptive interpolation algorithm.

1. Introduction

Artificial satellites provide important support for global observation services, such as global navigation and remote sensing. As data must be transmitted within a certain visibility range [1,2], the critical foundation for providing services is that the region of interest must be visible from the viewpoint of the satellite. Satellite orbit analysis and constellation coverage analysis are directly related to satellite-to-site visibility [3,4]. The accuracy and speed of satellite coverage are determined by calculating the visible period of the satellite. Task scheduling optimizations, which can reduce service costs by avoiding excessive rescans, require satellite-to-ground visibility predictions [5,6]. At least four satellites need to be used for positioning using navigation satellites, so satellite-to-site visibility is an important indicator of the ability to use satellites for navigation and positioning. Satellite visibility is very important in all aspects of satellite analysis.

The most primitive way to calculate the satellite-to-site visibility is to track the satellite’s trajectory and then to determine the visibility at each moment [5], which is the traditional brute force method. This method is often used as a comparison standard for improving subsequent algorithms. Although the brute force method is accurate, it is very time-consuming. Lawton et al. [7] developed an algorithm for calculating the visible period of a low eccentricity satellite orbit using iteration and a Fourier transform, which greatly improved the calculation speed compared with the brute force method. Alfano et al. [8] used the parabolic blending technique to construct the waveform of the visibility function at equal intervals and used the root of the local cubic polynomial to represent the rise and set times of the satellite. The algorithm is applicable to all satellite orbit types. Sun et al. [9] improved Alfano’s method and proposed an adaptive Hermite interpolation method in order to solve the satellite visibility problem. This method determines the interpolation step by checking the consistency of the second derivative and comparing the extreme value variance. Based on this work, Han et al. [5] derived an adaptive Hermite interpolation technology in a strict sense, which uses the fourth derivative to control the approximation error and achieve a better balance between the accuracy and efficiency. Han et al. [10] used a radial basis function to fit the satellite visibility function and accelerated the calculation of satellite visibility according to the interval contraction strategy. Wan et al. [11] designed a metamodeling framework based on adaptive interpolation and used different metamodel technologies, such as the radial basis function, Kriging, and support vector regression, to solve the satellite visibility problem.

The return to the Moon program has been proposed [12], and it has been realized that precise Moon position information is of great importance for an in-depth exploration of the Moon. In addition, the Magnetospheric Multiscale (MMS) mission proposed by NASA has also confirmed that satellite signals can acquire the position information of spacecraft located beyond the Earth’s orbit [13]. In recent years, the application of Global Navigation Satellite System (GNSS) technology to lunar navigation has attracted the attention of scientists. Kar-Ming et al. [14] simulated lunar navigation with GPS, GLONASS, Galileo, etc., using weak satellite signals for positioning with an accuracy of 200–300 m. The European Space Agency’s proposed Moonlight initiative will carry an advanced satellite signal receiver and perform the first satellite navigation and positioning mission in lunar orbit [15]. In order to carry out navigation and positioning on the Moon, the visibility of the Moon from satellites is an issue that must be analyzed and studied. Both the ESA and NASA have performed detailed analyses of the expected visibility of GNSS signals at Moon altitude [16]. A high-precision and efficient method for determining satellite-to-Moon visibility can provide a theoretical basis for achieving navigation and positioning on the Moon.

Since there are almost no existing algorithms for satellite-to-Moon visibility, this paper draws on the idea of satellite-to-site visibility calculations and identifies a problem when applying the adaptive interpolation algorithm to the satellite-to-Moon visibility calculation (explained in detail in the next section). To avoid this problem, an adaptive interpolation algorithm based on vertex protection is proposed, which is applicable not only to satellite-to-site visibility, but also to satellite-to-Moon visibility. At present, this paper might be the first to provide a fast solution for determining satellite-to-Moon visibility data. In Section 2, the satellite-to-Moon visibility model, which includes the elevation angle model and the Earth occultation model, highlights the drawbacks that arise from directly using the traditional adaptive interpolation algorithm to solve the satellite-to-Moon visibility problem. In Section 3, an adaptive interpolation algorithm based on vertex protection is introduced, and a specific scheme for calculating the satellite-to-Moon visibility using an adaptive interpolation algorithm based on vertex protection is described. In Section 4, the algorithm given in Section 3 is used to experimentally analyze the BDS satellite data. Section 5 concludes this paper.

2. Mathematical Models

2.1. Elevation Angle Function

Since the radius of the Moon is small, approximately 1738 km on average, and the distance from the satellite to the Moon can reach more than 350,000 km, the Moon is abstracted as a point for analysis in the following model.

Due to the ionospheric effect, tropospheric effect, etc., the satellite must first be considered under the condition of meeting a certain elevation angle with the station in order to achieve a good positioning effect. The elevation angle criterion is thus given as follows:

$$V(t)=\frac{\Delta r·r_s}{\Delta r·r_s}>\sin {\theta }_0$$

(1)

where $\Delta r=r_m-r_s$ and $r_s$ represent the coordinate vectors of the satellite and $r_m$ represents the coordinate vector of the Moon. The calculation of the satellite visibility at a certain elevation angle can also be approximated as the problem of solving $V(t)-\chi =0$, where $\chi$ is the cosine of elevation angle threshold ${\theta }_0$. The solution to this equation can be used for the time points of the satellite’s rise and set with respect to the Moon.

2.2. Earth Occultation

Since satellites are designed to serve the Earth, they rotate around the Earth and are oriented towards the Earth, so when applying satellites to the Moon, the Earth occultation problem needs to be considered. As shown in Figure 1, when the vertical distance between the satellite and the Moon is less than the radius $R$ of the Earth, the satellite’s line of sight to the Moon is blocked by the Earth, and the satellite is not visible to the Moon. The above problem can be described mathematically as follows:

$$L(t)=\frac{||r_s\times \Delta r||}{||\Delta r||R}<1$$

(2)

where $r_s$ is the coordinate vector of the satellite in an Earth-centered, Earth-fixed coordinate system, and $\Delta r=r_m-r_s$. When $L<R$, the line of sight from the satellite to the Moon is unobstructed at this time, and it is obscured by the Earth in all other cases.

2.3. Adaptive Interpolation Algorithm

As seen in Figure 2, where C01 represents a BDS GEO satellite, both the rise and set functions of the satellite-to-Moon visibility and the period function of the Earth occultation between the Earth and the Moon are multi-hump functions, which are irregular and do not have specific analytical expressions. It is known from previous work [5,10] that the solution to a similar function can be approximated using the idea of segment interpolation.

For the satellite-to-Moon visibility, this paper utilizes the idea of segmented Hermite interpolation to construct a cubic polynomial between the subintervals $t_0$ to $t_1$ and calculates the maximum step from the second-order derivatives of the corresponding visibility functions, as well as the visibility function values in satisfying a certain error limit. According to [17], the function of segmented cubic Hermite interpolation can be written as:

$$\begin{array}{l}S(t)=f(t_0)+f^{\prime }(t_0)(t-t_0)+a{(t-t_0)}^2+b{(t-t_0)}^3\\ a=(3\frac{f(t_1)-f(t_0)}{h^2}-\frac{f^{\prime }(t_1)+2f^{\prime }(t_0)}{h})\\ b=(\frac{f^{\prime }(t_1)+f^{\prime }(t_0)}{h^2}-2\frac{f(t_1)-f(t_0)}{h^3})\end{array}$$

(3)

where $t_0$ denotes the starting point for constructing the cubic polynomial, $t_1$ denotes the end point, and $h=t_1-t_0$. When the fourth-order derivative exists, the approximation error can be described by:

$$R(t)=V(t)-S(t)=\frac{V^{(4)}(\eta )}{4!}{((t-t_0)(t-t_1))}^2$$

(4)

where $t,\eta \in \left[t_0,t_1\right]$, $R(t)$, denotes the interpolation error residual term, $V(t)$ denotes the true function, $V^{(4)}(\eta )$ is the fourth-order derivative of $V(t)$ and can be obtained using higher-order Hermite interpolation, and $S(t)$ denotes the constructed cubic polynomial. The quadratic polynomial of the function can be approximated by the following equation:

$$\begin{array}{l}{V}^{(4)}(t)\approx 120{\kappa }_{5}t+24{\kappa }_4\\ {\kappa }_4=\frac{4}{h^4}\alpha -\frac{4}{h^4}\beta -\frac{24}{h^5}\chi \\ {\kappa }_5=\frac{24}{h^5}\left[V(t_0)-V(t_1)\right]+\frac{4}{h^4}\left[\dot{V}(t_0)-\dot{V}(t_0+\frac{h}{2})+\dot{V}(t_1)\right]\\ \alpha =\left[V(t_0)+4V(t_0+\frac{h}{2})+V(t_1)\right]\\ \beta =\left[\dot{V}(t_0)(2t_0+3t_1)+10\dot{V}(t_0+\frac{h}{2})(t_0+t_1)+3\dot{V}(t_1)(3t_0+2t_1)\right]\\ \chi =\left[V(t_0)(2t_0+3t_1)-V(t_1)(3t_0+2t_1)\right]\end{array}$$

(5)

Using $V^{(4)}(t)$ to place limits on the errors, the following is obtained:

$$|R(t){|}_{\max }\le |\frac{V^{(4)}{(n)}_{\max }}{4!}|{\left[(t-t_0){(t-t_1)}^2\right]}_{\max }\le |5a_{5}n+a_4{|}_{\max }{(\frac{h^2}{4})}^2$$

(6)

Then, the step is:

$$\widehat{h}={(\frac{16\epsilon }{|5{\kappa }_5\eta +{\kappa }_4{|}_{\max }})}^{1/4}$$

(7)

The solution for the step is an iterative process, and the termination condition of the iterative solution is:

$$\frac{|{\widehat{h}}_k-{\widehat{h}}_{k-1}|}{{\widehat{h}}_{k-1}}\le u$$

(8)

From the above equation, it is clear that the first-order derivative of the original function is required to perform the segmented Hermite interpolation, where the first-order derivative of the satellite elevation angle is:

$$\dot{V}(t)=\frac{1}{\Vert \Delta r\Vert }(\Delta \dot{r}·r_s+\Delta r·{\dot{r}}_s)-\frac{1}{\Vert \Delta r\Vert }(\Delta r·\Delta \dot{r})\Delta r·r_s$$

(9)

The first-order derivative of the Earth occultation model is:

$$\dot{L}(t)=\frac{(r_s\times \Delta r)·({\dot{r}}_s\times \Delta r+r_s\times \Delta \dot{r})}{\Vert r_s\times \Delta r\Vert ·\Vert \Delta \Vert }-\frac{\Delta ·\dot{\Delta }·\Vert r_s\times \Delta r\Vert }{\Vert \Delta {\Vert }^3}$$

(10)

where $\Delta \dot{r}=v_m-v_s$.

3. Algorithm Design

3.1. Vertex Protection Algorithm

The satellite-to-Moon visibility consists of two aspects: the availability of the satellite signal is limited by the elevation angle (as deduced in Section 2.1), and the Earth may obscure the line of sight from the satellite to the Moon (as deduced in Section 2.2). The algorithm provided in Section 2.3 can calculate the maximum interpolation step while guaranteeing the fitting accuracy. As shown in Figure 3, under the low fitting accuracy requirement, when calculating whether the Earth obscures the line of sight between the satellite and the Moon, an intersection exists between the original visibility function curve and the threshold, but the fitted curve has no intersection with the threshold, which causes errors when calculating the Earth occultation period in one cycle (28 days). Such a problem is related to the threshold selection and fitting accuracy requirement.

To avoid this problem, in this paper, the basic idea is to try to make the local extremal points of the multi-hump function the endpoints of the subinterval, where the local extremal points are the “vertices”. This method not only avoids the visibility calculation error, but also has the advantage that it can directly determine whether there is a solution to the original function within that subinterval without constructing a cubic polynomial in each subinterval. The interpolation method must know the original function values of the two endpoints. Then, according to the existence theorem of zero points, only when the two endpoints’ values corresponding to the function are different must the subinterval exist within the zero solution. Then, a cubic polynomial is constructed and solved. The vertex protection algorithm is designed as follows:

(1) Calculate the adaptive step $h_0$ according to the algorithm in Section 2.3, determine the termination time $t_e$ according to step $h_0$ and start time $t_s$, and introduce the search factor $\sigma$.

(2) Calculate the values $f_1$, $f_2$, $f_3$, and $f_4$ of the function $F(t)$ at $t_1=t_s$, $t_2=t_s+\sigma$, $t_3=t_e-\sigma$, and $t_4=t_e$ based on the original function (satellite rise and set function or occultation function, respectively).

(3) Compare the magnitudes of $f_1$, $f_2$, and $f_4$ or $f_1$, $f_3$, and $f_4$. When $f_2$ is smaller than both $f_1$ and $f_4$, go to step (4), and when $f_3$ is larger than both $f_1$ and $f_4$, go to step (5). When $f_2$ is greater than both $f_1$ and $f_4$, go to step (6), and when $f_3$ is smaller than both $f_4$ and $f_4$, go to step (7).

(4) Let $t_4=t_2$, and calculate the corresponding $f_4$. Calculate the corresponding $f_2$ when $t_2=t_2+\sigma$. Compare $f_4$ and $f_2$. If $f_2<f_4$, repeat step (4); otherwise, go to step (8).

(5) Let $t_4=t_3$, and calculate the corresponding $f_4$. Calculate the corresponding $f_3$ when $t_3=t_3-\sigma$. Compare $f_4$ and $f_3$. If $f_3>f_4$, repeat step (5); otherwise, go to step (8).

(6) Let $t_4=t_2$, and calculate the corresponding $f_4$. Calculate the corresponding $f_2$ when $t_2=t_2+\sigma$. Compare $f_2$ and $f_4$. If $f_2>f_4$, repeat step (4); otherwise, go to step (8).

(7) Let $t_4=t_3$, and calculate the corresponding $f_4$. Calculate the corresponding $f_3$ when $t_3=t_3-\sigma$. Compare $f_4$ and $f_3$. If $f_3<f_4$, repeat step (5); otherwise, go to step (8).

(8) Calculate when $t_s=t_1$ and $t_e=t_4$ correspond to $f_s=f_1$ and $f_e=f_4$ and obtain $h=t_e-t_s$.

Here, $\sigma$ denotes the step for searching for the local maxima of the multi-hump function, which determines the search accuracy and efficiency. The smaller $\sigma$ is, the more accurate and slower the local maxima of the multi-hump function will be as the endpoint of the interpolated subinterval; the larger $\sigma$ is, the coarser and faster the local maxima of the multi-hump function will be as the endpoint of the interpolated subinterval.

3.2. Adaptive Interpolation Based on Vertex Protection

The basic idea of adaptive interpolation based on the vertex protection method is as follows: according to the interpolation error function and the fourth-order derivative constraints, calculate the interpolation step $h_0$, and then use the obtained step $h_0$ and the starting time $t_s$ in the vertex protection algorithm to obtain $h$. Next, try to make the vertex of the multi-hump function the interpolation endpoint, and then, according to step $h$ and the starting time point $t_0$, determine the termination time point $t_1$ of the subinterval. The values $f_0$ and $f_1$ corresponding to the original function $F(t)$ at times $t_0$ and $t_1$, respectively, are calculated to determine whether the original function in this region has a solution. This shortens the process of constructing a cubic polynomial function for each subinterval in order to determine whether there is a solution and then solving for it, thereby effectively reducing the amount of computation. For the satellite-to-Moon visibility problem, the step calculation is based on the mathematical model in Section 2.3, and the vertex protection algorithm is the method proposed in Section 3.1. The rise and set periods of the satellite relative to the Moon are quickly solved using the adaptive interpolation method based on vertex protection, and then the Earth’s occultation for the line of sight between the satellite and the Moon is calculated within such periods, which can shorten the Earth occultation calculated time. The specific computational flow chart of the adaptive interpolation based on vertex protection is shown in Figure 4.

4. Experiment and Analysis

Since the period of the Moon’s revolution around the Earth is approximately 27.32 days, the period of a geosynchronous Earth orbit (GEO) satellite is approximately 24 h, the period of an inclined geosynchronous orbit (IGSO) satellite is approximately 24 h, and the period of a medium Earth orbit (MEO) satellite is approximately 12 h. Satellite data with a period of 28 days for the Moon’s revolution are used for the analysis. The data used here are from the multisystem precision ephemeris file provided by the Analysis Centre of Wuhan University for 31 December–30 January 2021, which contains the coordinates of each BDS, GLONASS, and Galileo satellite at five-minute intervals in the geocentric geostationary coordinate system. The lunar coordinates were calculated using a simplified model [18].

Figure 5 and Figure 6 show a comparison between the adaptive interpolation algorithm based on vertex protection and the traditional adaptive interpolation algorithm [5] in fitting the elevation angle function and the Earth occlusion function from the C01 satellite to the Moon when $\epsilon =0.1$ and $\sigma =0.2$. From the two figures, it can be seen that the adaptive interpolation algorithm based on vertex protection can make the local extreme points of the multi-hump functions, such as the elevation angle function and the Earth occlusion function, the endpoints of the interpolation sub-interval (to the greatest extent possible), and the fit is good.

The period of visibility from the satellite to the Moon is calculated as the true value using the brute force method, including the solution for the rise and set moments of the satellite with respect to the Moon and the period of the Earth’s line-of-sight obscuration between the satellite and the Moon.

Table 1. Satellite visibility times.

Satellite	Methods	$\mathit{\epsilon }$	Rise and Set Times	Occlusion Times
C01 GEO Satellite	Brute force method	0.1	28	6
		0.01	28	6
		0.001	28	6
	Self-adaptive interpolation technique without vertex protection	0.1	28	4
		0.01	28	5
		0.001	28	6
	Self-adaptive interpolation technique with vertex protection	0.1	28	6
		0.01	28	6
		0.001	28	6
C06 IGSO Satellite	Brute force method	0.1	28	3
		0.01	28	3
		0.001	28	3
	Self-adaptive interpolation technique without vertex protection	0.1	27	1
		0.01	28	3
		0.001	28	3
	Self-adaptive interpolation technique with vertex protection	0.1	28	3
		0.01	28	3
		0.001	28	3
C11 MEO Satellite	Brute force method	0.1	53	13
		0.01	53	13
		0.001	53	13
	Self-adaptive interpolation technique without vertex protection	0.1	53	12
		0.01	53	13
		0.001	53	13
	Self-adaptive interpolation technique with vertex protection	0.1	53	13
		0.01	53	13
		0.001	53	13

shows a comparison between the adaptive interpolation based on the vertex protection algorithm, the conventional adaptive interpolation algorithm, and the brute force method in solving the number of periods of visibility between the satellite and the Moon; that is, the times of the satellite’s rise and set with respect to the Moon and the times of the Earth’s occultation. For example, for the C01 satellite at fitting errors $\epsilon =0.1$ and $\epsilon =0.01$, the C06 satellite at fitting error $\epsilon =0.1$, and the C11 satellite at fitting error $\epsilon =0.1$, the solved Earth occlusion times do not match the real times, which are errors that cannot occur in practical applications. Adaptive interpolation based on vertex protection can accurately solve the Earth occlusion count for satellites C01, C06, and C11 when the search factor in the protection vertex algorithm is $\sigma =0.2$ for the fitting errors $\epsilon =0.1$, $\epsilon =0.01$, and $\epsilon =0.001$.

The efficiency of the adaptive interpolation method based on vertex protection proposed in this paper is compared to that of the traditional adaptive interpolation method when applied to the satellite-to-Moon visibility calculation in

Table 2. Comparison of the C01 satellite visibility solving efficiencies.

Case	$\mathit{\epsilon }=0.1$ After/Before		$\mathit{\epsilon }=0.01$ After/Before		$\mathit{\epsilon }=0.001$ After/Before
Step calculation times	248	295	404	539	1104	1611
Resolution times	71	295	69	539	70	1611
Run time	98.34%	98.8%	98.03%	98.2%	96.5%	94.6%

,

Table 3. Comparison of the C06 satellite visibility solving efficiencies.

Case	$\mathit{\epsilon }=0.1$ After/Before		$\mathit{\epsilon }=0.01$ After/Before		$\mathit{\epsilon }=0.001$ After/Before
Step calculation times	218	239	472	745	2990	6083
Resolution times	79	239	63	745	63	6083
Run time	98.6%	99.2%	98.1%	97.6%	93.5%	84.6%

and

Table 4. Comparison of the C11 satellite visibility solving efficiencies.

Case	$\mathit{\epsilon }=0.1$ After/Before		$\mathit{\epsilon }=0.01$ After/Before		$\mathit{\epsilon }=0.001$ After/Before
Step calculation times	467	579	1318	2196	10,539	19,282
Resolution times	135	579	135	2196	135	19,282
Run time	97.5%	98.2%	95.0%	93.6%	76.2%	47.6%

. For satellites C01, C06, and C11, under the limits of interpolation accuracies of 0.1, 0.01, and 0.001, respectively, the number of step calculations, the number of times the algorithm actually constructs a cubic polynomial for the solution, and the time efficiency improvement (compared with the brute force method) are compared.

As seen from these tables, the number of step calculations is always lower than that of the traditional adaptive interpolation algorithm, and the difference in the number of step calculations is more obvious as the fitting accuracy decreases. The algorithm proposed in this paper constructs a cubic polynomial for the number of solutions, and a cubic polynomial is constructed for when there is a solution in that subinterval; however, the traditional adaptive interpolation algorithm requires the construction of cubic polynomials in each subinterval to determine whether the interval has a solution and then calculate it. These tables give the percentage improvements in the time efficiency for adaptive interpolation based on the vertex protection algorithm and traditional adaptive interpolation compared with the brute force method. Under the requirement of a low fitting accuracy $\epsilon =0.1$, the computational efficiency of the algorithm proposed in this paper and the traditional algorithm are essentially comparable, and both improve by approximately 98% compared with the brute force method. The higher the requirement is, the more obvious the improvement in efficiency. When $\epsilon =0.001$ and $\sigma =0.2$, it can be seen from the following table that, for the C11 MEO satellite, the algorithm proposed in this paper improves the computational efficiency by approximately 30% compared with the traditional adaptive computation. This affects the accuracy of the vertex protection algorithm when seeking the vertex as the endpoint of the fitted subinterval, where the larger $\sigma$ is, the faster the algorithm is computed under the condition that the overall cycle is computed error-free.

To represent the accuracy of the adaptive interpolation algorithm based on vertex protection, the percentage normalization error defined by [19] is listed, where the percentage normalization is specifically defined as follows:

$$PNE=\frac{|Predicted satellite rise/set time-Acutal satellite rise/set time|}{Actual in-view period}\times 100$$

(11)

Figure 7 and Figure 8 show the $PNE$ of the satellite rise and set times and the $PNE$ of the Earth occultation start and end times, respectively; the red squares indicate the $PNE$ of the start time and the green stars indicate the $PNE$ of the end time. The subgraph C01-0.1 in Figure 7 shows the $PNE$ of the C01 satellite at the rise and set times relative to the Moon when $\epsilon =0.1$. Other subgraphs can also be interpreted in this way. The subgraph C01-0.1 in Figure 8 indicates the $PNE$ values at the start and end times of the Earth occultation for the line of sight from satellite C01 to the Moon when $\epsilon =0.1$. Other subgraphs are also expressed in this way.

5. Conclusions

Due to the existence of the traditional adaptive interpolation de-fitting error, the solution to the satellite visibility problem results in a cycle calculation error. In this paper, we introduced an adaptive interpolation algorithm based on vertex protection that makes the vertices of the multi-hump function the endpoints of the adaptive interpolation subinterval to the greatest extent possible in order to avoid the period calculation error in the satellite visibility problem, making the adaptive interpolation algorithm applicable to a wider range of satellite visibility problems. The first solution for quickly determining the satellite-to-Moon visibility problem was given, and experiments were conducted using data from three BDS orbiting satellites to analyze its efficiency and accuracy. The experiments showed that the method reduces the computation time by approximately 98% compared to the brute force method, and the computational accuracy and efficiency of the algorithm proposed in this paper are better than those of the traditional adaptive interpolation method, with a reduction in computation time of up to 30%.