Secant-Improved State Estimation Method for Moving Target Tracking Under Video Satellite

Abstract

A video satellite has continuous imaging capabilities, which grants it great potential for tracking and monitoring moving targets. The Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF) are commonly used in the above process. However, the accuracy of EKF estimation is low, and the computational complexity of UKF estimation is high. To address the contradiction between estimation accuracy and real-time performance in mobile-target state estimation, in this paper, we propose a new Kalman Filter with a secant-approximating nonlinear function. Firstly, the truncation error mechanism in the EKF is analysed here to illustrate the limitation of the EKF in approximating the nonlinear function. Then, the paper recommended a secant method to approximate the nonlinear function, which improved fitting accuracy without excessively increasing computational complexity. In order to improve the robustness of the proposed method, an adaptive selection strategy for correction elements is designed based on the advantageous range of secant approximation. The simulation results show that, in conventional ship motion scenarios, the computational accuracy is comparable to that of the EKF. In constant-power acceleration scenarios, the target positioning accuracy was 28.6% better than that of the EKF, and the computational speed was an order of magnitude greater than that of the UKF.

1. Introduction

Target motion state estimation plays a crucial role in various fields, such as ship navigation and collision warning. Specially, estimating the motion state of surface targets is indispensable in ship collision avoidance, accident rescue, and key vessel monitoring, making this one of the most important application scenarios for target motion state estimation today [1]. In terms of observation methods, using optical remote sensing satellites for target motion state estimation is one of the most commonly used approaches [2,3], offering advantages such as high stealthiness [4], a wide field of view, and minimal influence from the Earth’s curvature and terrain features. However, the high relative velocity between the satellite and target, coupled with the long imaging intervals of conventional optical satellites [5], makes it challenging to achieve continuous tracking and observation of a target [6]. Video satellites have shorter imaging intervals and continuous tracking imaging capabilities, effectively overcoming the limitations of conventional optical satellites in estimating target motion states [7]. Therefore, this paper primarily focuses on target state estimation under video satellite observation. Additionally, due to the highly dynamic nature of targets, target movement state estimation requires real-time results. That is to say, while there are certain requirements for estimation accuracy, there are also certain requirements for computational speed in target motion state estimation. Therefore, target movement state estimation methods must be designed in consideration of low computational complexity.

In recent years, methods for estimating the states of moving targets based on optical observations have included angle-only positioning [8,9,10], image matching [11], and dynamic estimation [12,13,14,15]. Angle-only positioning typically utilises satellite-observed angle information, employing mathematical transformations and least squares methods to estimate the target position, before calculating other states based on this information. Image matching methods typically use optical image recognition techniques such as convolutional neural networks, in conjunction with geographic information systems, to achieve target state estimation. In practical scenarios, positioning based solely on angular information does not fully utilise all the available information, resulting in lower positioning accuracy. Image matching methods require surrounding scene feature information, but sea surface targets have weak background information, making them difficult to apply effectively. Kinematic estimation methods typically combine target observation information with kinematic models, using filtering algorithms to estimate target states. This method fully utilises known information and can achieve high state estimation accuracy for sea surface targets. Currently, most filtering algorithms are based on the Kalman Filter (KF). The KF is a state estimation method that uses error distribution information to weight and fuse observation information with the target state model, achieving optimal estimation accuracy under linear conditions.

In actual scenarios, targets often move in a non-uniformly accelerating linear pattern. In such cases, both the target’s kinematic and optical observation equations become nonlinear equations, rendering the KF inapplicable. To address filtering problems in nonlinear systems, mathematical methods are typically employed for approximation, encompassing both linear and stochastic model approximation methods. Linear approximation methods include the batch least squares methods and the Extended Kalman Filter (EKF), and stochastic model approximation methods include the Unscented Kalman Filter (UKF), Particle Filter (PF), Cubature Kalman Filter (CKF) [16,17], and Gaussian Sum Filter (GSF) [18,19]. Nonlinear equation processing methods based on stochastic model approximation have the advantage of high equation fitting accuracy; however, they suffer from high computational complexity, which limits their application in practical engineering scenarios such as in on-board computations [13,20,21]. In practical engineering applications, linear approximation methods are more widely used due to their computational simplicity. The batch least squares methods require all obtained data to be used in calculations, which is not conducive to real-time online state estimation. In contrast, the EKF enables real-time system state estimation and is currently the most widely used filtering algorithm in practical applications. Eang C et al. [22] combined EKF with Deep Neural Network algorithms to achieve ultra-wideband localization for distance loss optimisation. Wang S et al. [23] proposed an improved interactive multi model EKF in their research, which improved the accuracy of target tracking. Zhao W et al. [24] designed a collaborative localization method for multiple autonomous underwater robots by combining Federated Learning and EKF.

The EKF linearises the system model by truncating the first-order terms of the Taylor series expansion of the nonlinear system model, before using the KF to perform state estimation [25]. In certain practical application scenarios, the EKF cannot achieve optimal estimation accuracy, primarily due to truncation errors and the presence of non-ideal noise in certain scenarios. To address the issue of non-ideal noise, adaptive correction methods such as the Adaptive EKF [26] and the Augmented EKF can be used to handle different types of noise. To address the issue of truncation error, the authors of [27] employed a Gauss–Newton method to achieve maximum likelihood estimation of the system state, thereby proposing an iterative EKF. And the authors of [28] employed the Lagrange multiplier optimisation technique to propose a Constrained Kalman Filter algorithm with quadratic state equation constraints based on the EKF. While these algorithms effectively reduce truncation errors, they require iterative calculations, thereby losing the EKF’s advantages of low computational complexity and limiting their potential in practical engineering applications.

Considering the above issues, this paper addresses truncated errors in estimation accuracy when using the EKF for motion target positioning and proposes a secant-line-improved Kalman filtering estimation method, which involves re-fitting the nonlinear equations using the secant line method for estimation result correction. This reduces the truncated error of EKF while keeping the computational complexity approximately unchanged, thereby improving estimation accuracy. The contributions of this paper can be summarised as follows:

(1). This paper analysed the causes of truncation errors in the EKF. Based on this, a new filter method was proposed by re-fitting the nonlinear function using secant lines to correct truncation errors in the first-order Taylor expansion of the EKF. This method improves estimation accuracy without losing the advantages of low computational complexity and fast operation speed of EKF.

(2). This paper also analysed the advantages of secant approximation over first-order Taylor series approximation. Based on this, adaptive criteria for determining correction elements were established. These criteria can be used to identify coefficient matrix parameters with superior estimation accuracy, thereby ensuring the robustness of the proposed method.

The rest of this study is structured as follows. Section 2 presents the target motion and observation models in video satellite gaze mode. Section 3 proposes a Kalman filtering method based on secant correction based on the analysis of EKF truncation error. Section 4 presents experimental simulations of the methods presented in this paper. Section 5 conclude the study by summarising the key research findings.

Notation: In this study, scalar values are represented by italicised letters, column vectors by bold lowercase letters, and matrices by bold italicised uppercase letters. The observations of a are denoted as $\tilde{a}$, and the estimates of a are represented as $\widehat{a}$ and $\widehat{a}$. cov(a) is the covariance matrix of a, and ${\left.f\left(x\right)\right|}_{x=x_1}$ denotes the value of the function f(x) when x = x₁.

2. Motion Target State Estimation Modelling

Nowadays, filtering methods usually combine dynamic models and observation data to achieve state estimation of moving targets. Therefore, this paper also comprehensively utilises ship motion dynamics information and video satellite observation data to estimate the position and velocity of moving ships. Therefore, when estimating the state of moving targets, it is first necessary to establish target kinematic and observation models.

Coordinate system: For ease of calculation, this paper uses the local Cartesian coordinate system (East, North, Up) with the initial target position as the origin. That is, the initial position of the target is taken as the origin, the geodetic tangent plane at that point is taken as the xOy plane, the x-axis points eastward in the plane, the y-axis points northward in the plane, and the z-axis is perpendicular to the tangent plane and points towards the sky.

2.1. Target Kinematic Modelling

Given the relatively short interval between video satellite imaging and the limited speed of sea surface targets, their position during observation is minimally affected by the Earth’s curvature. Therefore, this motion can be approximated as the target’s planar motion. The specific tracking process is shown in Figure 1. Let the target motion information be denoted as ${\mathit{\theta }}_{\mathit{h},\mathit{j}}={\left[x_{m,j},y_{m,j},z_{m,j},v_j,{\psi }_j\right]}^{\mathrm{T}}$, where ${\mathit{r}}_{\mathit{m},\mathit{j}}={\left[x_{m,j},y_{m,j},z_{m,j}\right]}^{\mathrm{T}}$ represents the target position at time j, v_j represents the target velocity at time j, and ψ_j represents the target heading angle at time j. The target kinematic equation can be established as follows:

(1)${\mathit{\theta }}_{\mathit{h},\mathit{j}}=f_h\left({\mathit{\theta }}_{\mathit{h},\mathit{j}-1}\right)$

In actual situations, the motion state of the target is subject to noise due to various factors such as waves, wind force, and engine operating conditions. Therefore, the formula needs to be modified to

(2)${\mathit{\theta }}_{\mathit{h},\mathit{j}}=f_h\left({\mathit{\theta }}_{\mathit{h},\mathit{j}-1}\right)+{\mathit{\epsilon }}_{\mathit{h},\mathit{j}}$

where ${\mathit{\epsilon }}_{\mathit{h},\mathit{j}}={\left[{\mathit{\epsilon }}_{\mathit{x},\mathit{j}},{\mathit{\epsilon }}_{\mathit{y},\mathit{j}},{\mathit{\epsilon }}_{\mathit{z},\mathit{j}},{\mathit{\epsilon }}_{\mathit{v},\mathit{j}},{\mathit{\epsilon }}_{\mathit{\psi },\mathit{j}}\right]}^{\mathrm{T}}$ represents the noise at time j in each state parameter, with a mean of 0 and covariance matrix of Q_h,j.

2.2. Video Satellite Observation Modelling

When video satellites track and locate targets, they often use an angle-only positioning and tracking model, which expresses the positional relationship between the target and the satellite through the azimuth and elevation angles, as shown in Figure 2. Given that the coordinates of the observation satellite at time j are ${\mathit{s}}_{\mathit{j}}={\left[s_{x,j},s_{y,j},s_{z,j}\right]}^{\mathrm{T}}$, the azimuth angle of the observation satellite observing the target is α_j, and the elevation angle is β_j, denoted as ${\mathit{\phi }}_{\mathit{j}}={\left[{\alpha }_j,{\beta }_j\right]}^{\mathrm{T}}$. The observation equation can be obtained as follows:

(3)$\left\{\begin{array}{l}{\alpha }_j=\arctan {\displaystyle \frac{y_{m,j}-s_{y,j}}{x_{m,j}-s_{x,j}}}\hfill \\ {\beta }_j=\arctan {\displaystyle \frac{z_{m,j}-s_{z,j}}{\sqrt{{\left(x_{m,j}-s_{x,j}\right)}^2+{\left(y_{m,j}-s_{y,j}\right)}^2}}}\hfill \end{array}\right.$

Considering the actual situation, where various factors such as satellite positioning accuracy and stability have an impact [8], there is an error between the measured value of the line-of-sight (LOS) angle and the true value. The observation equation set that uses the measured LOS angle value to replace the true value is as follows:

(4)$\left\{\begin{array}{l}{\tilde{\alpha }}_j=\arctan {\displaystyle \frac{y_m-s_{y,j}}{x_m-s_{x,j}}}+{\epsilon }_{\alpha ,j}\hfill \\ {\tilde{\beta }}_j=\arctan {\displaystyle \frac{z_m-s_{z,j}}{\sqrt{{\left(x_m-s_{x,j}\right)}^2+{\left(y_m-s_{y,j}\right)}^2}}}+{\epsilon }_{\beta ,j}\hfill \end{array}\right.$

where ${\tilde{\mathit{\phi }}}_j={\left[{\tilde{\alpha }}_j,{\tilde{\beta }}_j\right]}^{\mathrm{T}}$ is the LOS angle measurement at time j, and ${\mathit{\epsilon }}_{\mathit{\phi },\mathit{j}}={\left[{\epsilon }_{\alpha ,j},{\epsilon }_{\beta ,j}\right]}^{\mathrm{T}}$ is the angular error of the satellite at time j, with a mean of 0 and covariance matrix of R_s,j. For convenience of expression, Equation (4) can be written in matrix form as follows:

(5)${\tilde{\mathit{\phi }}}_{\mathit{j}}=h_s\left({\mathit{\theta }}_{\mathit{j}}\right)+{\mathit{\epsilon }}_{\mathit{\phi },\mathit{j}}$

3. Secant Kalman Filter Method

The KF can effectively fuse target state and observation information in linear situations to obtain high-precision target state estimation results, but it cannot effectively solve nonlinear problems. The EKF linearises nonlinear problems by using linear fitting based on the KF, thereby enabling filtering estimation of nonlinear problems. Due to its simplicity and low computational load, the EKF has become the most commonly used filtering method in current engineering applications. However, the linear approximation method used by the EKF is to approximate nonlinear equations using tangent line, which will introduce truncation errors and result in low estimation accuracy. This paper proposes replacing the linear approximation method of the EKF with an appropriate secant approximation to effectively reduce truncation errors and improve estimation accuracy.

3.1. Derivation and Error Analysis of the Extended Kalman Filter

The EKF mainly uses first-order Taylor expansion to perform linear fitting on the nonlinear state and observation equations. The state and observation equations are

(6)$\left\{\begin{array}{l}{\mathit{\theta }}_{\mathit{j}}=f\left({\mathit{\theta }}_{\mathit{j}-1},{\mathit{w}}_{\mathit{j}}\right)\hfill \\ {\mathit{\phi }}_{\mathit{j}}=h\left({\mathit{\theta }}_{\mathit{j}},{\mathit{w}}_{\mathit{j}}\right)+{\mathit{v}}_{\mathit{j}}\hfill \end{array}\right.$

where θ_j is the target motion state, w_j is the process noise, with a covariance matrix of Q_j, ${\tilde{\mathit{\phi }}}_j$ is the measurement, and v_j is the measurement error, with a covariance matrix of R_j, all at time j. By using the state equation, the predicted value of the state at time j can be given as:

(7)${\overline{\mathit{\theta }}}_{\mathit{j}}=f\left({\mathit{\theta }}_{\mathit{j}-1}\right)$

based on the state at time j − 1.

In order to facilitate fusion coefficient calculations, it is necessary to convert the formula system state and observation equations into linear equations. Taylor expansion is a commonly used polynomial function fitting method, and the first-order Taylor expansion results in linear equations. Therefore, using Taylor expansion, the state (Equation (6)) and observation equations are expanded to the first-order at ${\overline{\mathit{\theta }}}_{\mathit{j}}$ as follows:

(8)$\{\begin{array}{c}{\mathit{\theta }}_{\mathit{j}}=f({\mathit{\theta }}_{\mathit{j}-1},{\mathit{w}}_{\mathit{j}})+{\mathit{F}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{j}})\\ {\mathit{\phi }}_{\mathit{j}}=h({\overline{\mathit{\theta }}}_{\mathit{j}},{\mathit{w}}_{\mathit{j}})+{\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{j}})\end{array}$

where

(9)$\left\{\begin{array}{l}{\mathit{F}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}={\left.{\displaystyle \frac{\partial f}{\partial \mathit{\theta }}}\right|}_{\mathit{\theta }={\overline{\mathit{\theta }}}_{\mathit{j}}}\hfill \\ {\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}={\left.{\displaystyle \frac{\partial h}{\partial \mathit{\theta }}}\right|}_{\mathit{\theta }={\overline{\mathit{\theta }}}_{\mathit{j}}}\hfill \end{array}\right.$

constituting the Jacobian matrix of the state equation f and observation equation h with respect to the state variable θ. At this point, the filtering can be converted to a KF based on the expanded state equation and the expanded observation equation, and the covariance matrix of the prior probability error can be calculated as follows:

(10)${\overline{\mathit{P}}}_{\mathit{j}}={\mathit{F}}_{\mathit{E}\mathit{K}\mathit{F}}{\mathit{P}}_{\mathit{j}-1}{\mathit{F}}_{\mathit{E}\mathit{K}\mathit{F}}^{\mathrm{T}}-{\mathit{Q}}_{\mathit{j}}$

The gain matrix can be calculated as follows:

(11)${\mathit{K}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}={\overline{\mathit{P}}}_{\mathit{j}}{\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}^{\mathrm{T}}{\left({\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F}}{\overline{\mathit{P}}}_{\mathit{j}}{\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F},\mathit{j}}^{\mathrm{T}}+{\mathit{R}}_{\mathit{j}}\right)}^{-1}$

The data fusion equation is

(12)${\widehat{\mathit{\theta }}}_{\mathit{j}}={\overline{\mathit{\theta }}}_{\mathit{j}}+{\mathit{K}}_{\mathit{E}\mathit{K}\mathit{F}}\left[{\mathit{\phi }}_{\mathit{j}}-h\left({\overline{\mathit{\theta }}}_{\mathit{j}}\right)\right]$

The covariance matrix of the posterior state estimation error is

(13)${\mathit{P}}_{\mathit{j}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{E}\mathit{K}\mathit{F}}{\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F}}\right){\overline{\mathit{P}}}_{\mathit{j}}$

When linearising the equation, the EKF uses a first-order Taylor expansion, discarding higher-order terms, which introduces a certain amount of truncation error, as shown in Figure 3 (The red dot in Figure 3 is the origin, and Figures 6–8 are the same). To better illustrate the variation in truncation error, Figure 4 presents the error from truncating the first-order Taylor expansion at x = 1 for several typical functions. In the figure, the distance between the true value point and the expansion point is used to characterise the EKF prediction error. As shown in Figure 4, the magnitude of the truncation error is positively correlated with the nonlinearity of the equation and the prediction error. When the distance between the true value and expansion points is 0.9, the truncation error of the cube of the function reaches 5. The presence of a truncation error reduces the filtering accuracy, and this error may accumulate continuously during the filtering process. Additionally, when the equation exhibits strong nonlinearity or significant prediction errors, the truncation error of the EKF may cause the fusion result error to increase further compared to the prediction result, leading to estimation divergence.

3.2. Correction of Extended Kalman Filter Truncation Error

Selecting an appropriate linearisation method to approximate a nonlinear system can effectively reduce truncation errors. Since filtering calculations are performed sequentially for each observation, it is only necessary to separately consider the fitting effect of the approximate equation for each calculation. Taking a univariate function as an example, for any two points θ₁ and θ₂ on the nonlinear function h(θ), there exists a secant line g(θ), such that

(14)$\left\{\begin{array}{l}h\left({\mathit{\theta }}_1\right)=g\left({\mathit{\theta }}_1\right)\hfill \\ h\left({\mathit{\theta }}_2\right)=g\left({\mathit{\theta }}_2\right)\hfill \end{array}\right.$

The analytical expression of the secant function g(θ) is

(15)$g\left(\mathit{\theta }\right)={\displaystyle \frac{h\left({\mathit{\theta }}_2\right)-h\left({\mathit{\theta }}_1\right)}{{\mathit{\theta }}_2-{\mathit{\theta }}_1}}\left(\mathit{\theta }-{\mathit{\theta }}_1\right)+h\left({\mathit{\theta }}_1\right)$

When considering only the calculation of θ₂ through θ₁, there is no truncation error when using the secant g(θ) to approximate h(θ).

We can extend the above content to a multivariate function φ = h(θ), with q being dependent and p independent variables, and perform secant calculations using the ‘separate dependent variables, stepwise independent variables’ method. For the dependent variables, h(θ) is split into equations for each individual dependent variable with respect to all independent variables for further secant calculations, i.e.,

(16)$\left\{\begin{array}{l}{\mathit{\phi }}_{\left(1\right)}=h_1\left(\mathit{\theta }\right)\hfill \\ {\mathit{\phi }}_{\left(2\right)}=h_2\left(\mathit{\theta }\right)\hfill \\ ⋮\hfill \\ {\mathit{\phi }}_{\left(q\right)}=h_q\left(\mathit{\theta }\right)\hfill \end{array}\right.$

Considering the independent variables, we can take h₁(θ) as an example. For the two points θ₁ and θ₂, when calculating the secant line for the i-th component, the first i-1 components of the independent variable are fixed to the first i-1 components of θ₂, and the remaining p-i components are fixed to the remaining p-i components of θ₁. At this point, the function becomes a univariate function separated by the i-th independent variable, and the secant line for the i-th component can be calculated based on Equation (15) as follows:

(17)$\{\begin{array}{l}{g}_1({\phi }_{\left(i\right)})={\displaystyle \frac{h_1\left({\mathit{\theta }}_i^{\prime }\right)-h\left({\mathit{\theta }}_{i-1}^{\prime }\right)}{{\theta }_{2\left(i\right)}-{\theta }_{1\left(i\right)}}}\\ {\mathit{\theta }}_{\mathit{i}}^{\prime }={[{\theta }_{2\left(1\right)},{\theta }_{2\left(2\right)},\cdots ,{\theta }_{2\left(i\right)},{\theta }_{1(i+1)},\cdots ,{\theta }_{1\left(p\right)}]}^{\mathrm{T}}\end{array}$

The calculation process is shown in Figure 5, where the red line indicates the secant calculation of the first component, and the green line represents that of the second component.

Based on the above analysis, the secant operator SecLine[h(θ),θ₁,θ₂] can be defined as follows:

(18)$\{\begin{array}{l} \mathrm{S}\mathrm{e}\mathrm{c}\mathrm{L}\mathrm{i}\mathrm{n}\mathrm{e} [h(\mathit{\theta }),{\mathit{\theta }}_1,{\mathit{\theta }}_2]=\left[\begin{array}{cccc}{\displaystyle \frac{{\phi }_{0\left(1\right)}^{\prime }-{\phi }_{1\left(1\right)}^{\prime }}{{\theta }_{2\left(1\right)}-{\theta }_{1\left(1\right)}}}& {\displaystyle \frac{{\phi }_{1\left(1\right)}^{\prime }-{\phi }_{2\left(1\right)}^{\prime }}{{\theta }_{2\left(2\right)}-{\theta }_{1\left(2\right)}}}& \cdots & {\displaystyle \frac{{\phi }_{p-1\left(1\right)}^{\prime }-{\phi }_{p\left(1\right)}^{\prime }}{{\theta }_{2\left(p\right)}-{\theta }_{1\left(p\right)}}}\\ {\displaystyle \frac{{\phi }_{0\left(2\right)}^{\prime }-{\phi }_{1\left(2\right)}^{\prime }}{{\theta }_{2\left(1\right)}-{\theta }_{1\left(1\right)}}}& {\displaystyle \frac{{\phi }_{1\left(2\right)}^{\prime }-{\phi }_{2\left(2\right)}^{\prime }}{{\theta }_{2\left(2\right)}-{\theta }_{1\left(2\right)}}}& \cdots & {\displaystyle \frac{{\phi }_{p-1\left(2\right)}^{\prime }-{\phi }_{p\left(2\right)}^{\prime }}{{\theta }_{2\left(p\right)}-{\theta }_{1\left(p\right)}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{{\phi }_{0\left(q\right)}^{\prime }-{\phi }_{1\left(q\right)}^{\prime }}{{\theta }_{2\left(1\right)}-{\theta }_{1\left(1\right)}}}& {\displaystyle \frac{{\phi }_{1\left(q\right)}^{\prime }-{\phi }_{2\left(q\right)}^{\prime }}{{\theta }_{2\left(2\right)}-{\theta }_{1\left(2\right)}}}& \cdots & {\displaystyle \frac{{\phi }_{p-1\left(q\right)}^{\prime }-{\phi }_{p\left(q\right)}^{\prime }}{{\theta }_{2\left(p\right)}-{\theta }_{1\left(p\right)}}}\end{array}\right]\\ {\mathit{\phi }}_i^{\prime }=h\left({\mathit{\theta }}_i^{\prime }\right)\\ {\mathit{\phi }}_0^{\prime }=h\left({\mathit{\theta }}_1\right)\end{array}$

For multivariate nonlinear functions, there exists a multivariate linear function g(θ), such that

(19)$\left\{\begin{array}{l}h\left({\mathit{\theta }}_1\right)=g\left({\mathit{\theta }}_1\right)\hfill \\ h\left({\mathit{\theta }}_2\right)=g\left({\mathit{\theta }}_2\right)\hfill \end{array}\right.$

The analytical expression for g(θ) can be calculated as follows:

(20)$g\left(\mathit{\theta }\right)=\mathrm{SecLine}\left[g\left(\mathit{\theta }\right),{\mathit{\theta }}_1,{\mathit{\theta }}_2\right]\left(\mathit{\theta }-{\mathit{\theta }}_1\right)+g\left({\mathit{\theta }}_1\right)$

Applying the above method to expand the target state and observation equations, we obtain

(21)$\left\{\begin{array}{l}{\mathit{\theta }}_{\mathit{j}}=f\left({\mathit{\theta }}_{\mathit{j}-1}\right)+{\mathit{F}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}\left(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{j}}\right)\hfill \\ {\mathit{\phi }}_{\mathit{j}}=h\left({\overline{\mathit{\theta }}}_{\mathit{j}}\right)+{\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}\left(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{j}}\right)\hfill \end{array}\right.$

where F_GKF,j and H_GKF,j are the secant approximation coefficient matrices, which can be expanded as follows:

(22)$\left\{\begin{array}{l}{\mathit{F}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}=\mathrm{SecLine}\left[f\left(\mathit{\theta }\right),{\overline{\mathit{\theta }}}_{\mathit{j}},{\mathit{\theta }}_{\mathit{j}}\right]\hfill \\ {\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}=\mathrm{SecLine}\left[h\left(\mathit{\theta }\right),{\overline{\mathit{\theta }}}_{\mathit{j}},{\mathit{\theta }}_{\mathit{j}}\right]\hfill \end{array}\right.$

At this point, applying the KF algorithm based on Equation (21) at step j yields a filtered result with no truncation error, as shown in Figure 6.

However, the calculation of Equation (22) includes the true value of θ_j, which cannot be obtained. Therefore, it is necessary to use the filter value ${\widehat{\mathit{\theta }}}_j$ obtained by EKF to replace θ_j, that is

(23)$\left\{\begin{array}{l}{\mathit{F}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}=\mathrm{SecLine}\left[f\left(\mathit{\theta }\right),{\overline{\mathit{\theta }}}_{\mathit{j}},{\widehat{\mathit{\theta }}}_{\mathit{j}}\right]\hfill \\ {\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}=\mathrm{SecLine}\left[h\left(\mathit{\theta }\right),{\overline{\mathit{\theta }}}_{\mathit{j}},{\widehat{\mathit{\theta }}}_{\mathit{j}}\right]\hfill \end{array}\right.$

At this point, a certain amount of truncation error will occur, as shown in Figure 7.

On this basis, the new prior probability error can be calculated as follows:

(24)${\overline{\mathit{P}}}_{\mathit{G},\mathit{j}}={\mathit{F}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}{\mathit{P}}_{\mathit{G},\mathit{j}-1}{\mathit{F}}_{\mathit{G}\mathit{K}\mathit{F}}^{\mathrm{T}}-{\mathit{Q}}_{\mathit{j}}$

The gain matrix is

(25)${\mathit{K}}_{\mathit{G}\mathit{K}\mathit{F}}={\overline{\mathit{P}}}_{\mathit{j}}{\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}^{\mathrm{T}}{\left({\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F}}{\overline{\mathit{P}}}_{\mathit{j}}{\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F},\mathit{j}}^{\mathrm{T}}+{\mathit{R}}_{\mathit{j}}\right)}^{-1}$

The data fusion equation is

(26)${\widehat{\mathit{\theta }}}_{\mathit{j}}={\overline{\mathit{\theta }}}_{\mathit{j}}+{\mathit{K}}_{\mathit{G}\mathit{K}\mathit{F}}\left[{\tilde{\mathit{\phi }}}_{\mathit{j}}-h\left({\overline{\mathit{\theta }}}_{\mathit{j}}\right)\right]$

The covariance matrix of the posterior state estimation error is

(27)${\mathit{P}}_{\mathit{G},\mathit{j}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{G}\mathit{K}\mathit{F}}{\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F}}\right){\overline{\mathit{P}}}_{\mathit{G},\mathit{j}}$

3.3. Determining the Adaptive Judgement Conditions for Correction Elements

Using estimated target positions instead of actual values results in a certain amount of truncation error in secant approximation calculations. As shown in Figure 8, this error is generally smaller than that of a first-order Taylor expansion in most cases, but there are also instances where the truncation error exceeds that of the Taylor expansion. In order to ensure stable estimation accuracy improvement of the proposed method, the method must establish adaptive decision conditions to dynamically select correction parameters. Therefore, it is necessary to theoretically analyse the applicable range within which the secant approximation is superior to the first-order Taylor series expansion.

As shown in Figure 8, for scenarios where the secant approximation is better than the first-order Taylor expansion, the function needs to satisfy the following:

(28)$\left|{\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F}}^{-1}\right[{\mathit{\phi }}_{\mathit{i}}-h\left({\overline{\mathit{\theta }}}_{\mathit{i}}\right)+{\mathit{H}}_{\mathit{G}\mathit{K}\mathit{F}}{\overline{\mathit{\theta }}}_{\mathit{i}}]-{\mathit{\theta }}_{\mathit{i}}|<\left|{\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F}}^{-1}\right[{\mathit{\phi }}_{\mathit{i}}-h\left({\overline{\mathit{\theta }}}_{\mathit{i}}\right)+{\mathit{H}}_{\mathit{E}\mathit{K}\mathit{F}}{\overline{\mathit{\theta }}}_{\mathit{i}}]-{\mathit{\theta }}_{\mathit{i}}|$

Considering that the calculation of secant lines for multivariate functions can be decomposed into multiple univariate Equation (28), we can analyse the conditions that satisfy the formula in the case of univariate functions. Furthermore, by generalising to the case of multivariate functions and combining the calculation process of Equation (18), we can treat each element in the secant line correction matrix as a univariate secant line approximation process. Therefore, each element in the secant correction matrix can be analysed independently. When an element satisfies the aforementioned judgement conditions, it adopts the corresponding element from the secant correction matrix; when it does not satisfy the conditions, it adopts the corresponding element from the first-order Taylor expansion matrix. Based on this, considering that directly determining the judgement conditions is relatively difficult, the condition conversion method is adopted to determine whether the secant correction has advantages based on two judgement metrics.

(1). Judgement Metric 1

Condition conversion: Under a univariate function, proving that the secant approximation has an advantage in certain situations is equivalent to proving that this situation may be equivalent to when the calculated secant slope is closer to the true value than the slope given by the first-order Taylor expansion. That is, within the interval, the function satisfies

(29)$\left|\dot{h}\left({\overline{\theta }}_i\right)-{\displaystyle \frac{h\left({\theta }_i\right)-h\left({\overline{\theta }}_i\right)}{{\theta }_i-{\overline{\theta }}_i}}\right|>\left|{\displaystyle \frac{h\left({\widehat{\theta }}_i\right)-h\left({\overline{\theta }}_i\right)}{{\widehat{\theta }}_i-{\overline{\theta }}_i}}-{\displaystyle \frac{h\left({\theta }_i\right)-h\left({\overline{\theta }}_i\right)}{{\theta }_i-{\overline{\theta }}_i}}\right|$

In the case of Equation (29), the secant approximation has an advantage. Considering that, in EKF filtering, the fusion estimate accuracy is higher than the prediction value, the above problem is therefore transformed into a mathematical problem. For θ₁, θ₂, θ₃, satisfying $\left|{\theta }_1-{\theta }_3\right|>\left|{\theta }_2-{\theta }_3\right|$, the function h(θ) in the interval $\left[{\theta }_3-\left|{\theta }_3-{\theta }_1\right|,{\theta }_3+\left|{\theta }_3-{\theta }_1\right|\right]$ satisfies the following:

(30)$\left|\dot{h}\left({\theta }_1\right)-{\displaystyle \frac{h\left({\theta }_3\right)-h\left({\theta }_1\right)}{{\theta }_3-{\theta }_1}}\right|>\left|{\displaystyle \frac{h\left({\theta }_2\right)-h\left({\theta }_1\right)}{{\theta }_2-{\theta }_1}}-{\displaystyle \frac{h\left({\theta }_3\right)-h\left({\theta }_1\right)}{{\theta }_3-{\theta }_1}}\right|$

In the case of Equation (30), the secant correction results can effectively reduce truncation errors.

Condition determination: According to the concavity theorem [29], when h(θ) is a convex function in the interval, the above condition holds [30]. Therefore, when the function is convex, the tangent line correction result has a smaller truncation error.

(2). Judgement Metric 2

Condition conversion: When the function is not convex, the advantage of the secant approximation can be determined according to whether the fusion result has a tendency to approach the true value. Let

(31)${\theta }_{\delta }=\overline{\theta }+\delta$

where δ is a sufficiently small positive real number. It can be considered that

(32)$\dot{h}\left(\overline{\theta }\right)={\displaystyle \frac{h\left({\theta }_{\delta }\right)-h\left(\overline{\theta }\right)}{{\theta }_{\delta }-\overline{\theta }}}$

At this point, the tangent expansion can be regarded as a special secant expansion, with $\overline{\theta }$, $\widehat{\theta }$, $\widehat{\theta }$ regarded as the initial value and the solutions of the previous two steps in a special iterative process.

Condition determination: According to the Lipschitz condition [31], when the Lipschitz constant L_h of the function h(θ) is less than 1, it satisfies the following:

(33)$L_h\left|{\overline{\theta }}_j-{\widehat{\theta }}_j\right|\ge \left|{\widehat{\theta }}_j-{\widehat{\theta }}_j\right|$

At this point, the secant correction results can effectively reduce truncation errors.

Furthermore, even if the above conditions are not met, it is still possible that secant correction may display better fitting accuracy; however, no relevant proof can be provided at this time.

At the same time, common nonlinear state estimation scenarios were analysed based on mathematical analysis. As shown in Table 1, in most cases, tangent correction can satisfy the judgement conditions and effectively reduce truncation errors. In the table, x represents the state variable, θ represents the measurement variable, and the other letters represent parameters. The judgement analysis of the state equation and observation equation in specific situations is shown in Appendix B.

3.4. Secant Kalman Filter Process

This paper refers to the above method as Secant Kalman Filtering (SecKF). Combined with the scenario given in Section 2, the specific process is shown in Figure 9 (The solid box in Figure 9 represents the steps, and the dashed box represents the stages), and the main steps include the following:

(1). Prediction: Calculate the prediction result as ${\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}=f_h\left({\mathit{\theta }}_{\mathit{h},\mathit{j}-1}\right)$.

(2). Expansion: Perform a first-order Taylor expansion of the state equation and observation equation at the location to obtain the linearised state equation ${\mathit{\theta }}_{\mathit{j}}=f\left({\mathit{\theta }}_{\mathit{j}-1}\right)+{\mathit{w}}_{\mathit{j}}+{\mathit{F}}_{\mathit{h},\mathit{j}}\left(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{j}}\right)$ and observation equation ${\mathit{\phi }}_{\mathit{j}}=\mathit{h}\left({\overline{\mathit{\theta }}}_{\mathit{j}}\right)+{\mathit{w}}_{\mathit{j}}+{\mathit{H}}_{\mathit{s},\mathit{j}}\left(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{j}}\right)$, where the expansion of F_h,j is determined based on the actual scenario, the expansion of H_s,j is shown in Appendix A, and the prior probability error is calculated as ${\overline{\mathit{P}}}_{\mathit{h},\mathit{j}}={\mathit{F}}_{\mathit{h}}{\mathit{P}}_{\mathit{s},\mathit{j}-1}{\mathit{F}}_{\mathit{h}}^{\mathrm{T}}-{\mathit{Q}}_{\mathit{h},\mathit{j}}$.

(3). Initial fusion: Calculate the gain matrix as ${\mathit{K}}_{\mathit{h},\mathit{j}}={\overline{\mathit{P}}}_{\mathit{h},\mathit{j}}{\mathit{H}}_{\mathit{s},\mathit{j}}^{\mathrm{T}}{\left({\mathit{H}}_{\mathit{s},\mathit{j}}{\overline{\mathit{P}}}_{\mathit{h},\mathit{j}}{\mathit{H}}_{\mathit{s},\mathit{j}}^{\mathrm{T}}+{\mathit{R}}_{\mathit{s},\mathit{j}}\right)}^{-1}$.

(4). Initial update: Calculate the initial fusion result as ${\widehat{\mathit{\theta }}}_{\mathit{h},\mathit{j}}={\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}+{\mathit{K}}_{\mathit{h},\mathit{j}}\left[{\mathit{\phi }}_{\mathit{h},\mathit{j}}-\mathit{h}\left({\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}\right)\right]$ and the initial fusion posterior probability error as ${\widehat{\mathit{P}}}_{\mathit{h},\mathit{j}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{h}}{\mathit{H}}_{\mathit{s}}\right){\overline{\mathit{P}}}_{\mathit{h},\mathit{j}}$.

(5). Secant calculation: Calculate the approximate secant expansion equations of the state and observation equations, and obtain the linearised state ${\mathit{\theta }}_{\mathit{h},\mathit{j}}=f\left({\mathit{\theta }}_{\mathit{h},\mathit{j}-1}\right)+{\mathit{F}}_{\mathit{s}\mathit{g},\mathit{j}}\left(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}\right)$ and observation equations ${\mathit{\phi }}_{\mathit{h},\mathit{j}}=\mathit{h}\left({\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}\right)+{\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}\left(\mathit{\theta }-{\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}\right)$ of the secant. The expansion of F_hg,j is determined based on the actual scenario, and the expansion of H_sg,j is shown in Appendix A.

(6). Element determination: Using adaptive determination conditions, analyse whether the elements in F_h,j and H_s,j have smaller truncation errors than the corresponding elements in F_hg,j and H_sg,j. If not, replace the corresponding elements with those in F_h,j and H_s,j. In the dynamic equations, the nonlinear components include the equations for x_m and y_m with respect to v and ψ, as well as the equation for v with respect to v. In the observation equations, the nonlinear components include the equations for ɑ and β with respect to x_m and y_m. Among these, the equations for x_m and y_m with respect to ψ in the dynamic equations require replacement with the corresponding elements from F_h,j, while the remaining elements use those from F_hg,j and H_sg,j.

(7). Secondary fusion: Calculate the gain matrix as ${\mathit{K}}_{\mathit{h}\mathit{g},\mathit{j}}={\overline{\mathit{P}}}_{\mathit{j}}{\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}^{\mathrm{T}}{\left({\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}{\overline{\mathit{P}}}_{\mathit{j}}{\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}^{\mathrm{T}}+{\mathit{R}}_{\mathit{s},\mathit{j}}+{\mathit{P}}_{\mathit{H}\mathit{s},\mathit{j}}\right)}^{-1}$.

(8). Secondary update: Calculate the secondary fusion result as ${\widehat{\mathit{\theta }}}_{\mathit{h},\mathit{j}}={\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}+{\mathit{K}}_{\mathit{h}\mathit{g},\mathit{j}}\left[{\tilde{\mathit{\phi }}}_{\mathit{j}}-h\left({\overline{\mathit{\theta }}}_{\mathit{h},\mathit{j}}\right)\right]$, and the secondary fusion posterior probability error as ${\mathit{P}}_{\mathit{h},\mathit{j}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{h}\mathit{g},\mathit{j}}{\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}\right){\overline{\mathit{P}}}_{\mathit{h},\mathit{j}}$.

Figure 9

Basic SecKF process.

[Figure omitted. See PDF]

3.5. Computational Analysis

Furthermore, theoretically analyse the computational complexity of the proposed method. Due to the influence of equations and parameters on actual computational complexity, only the computational complexity is estimated. For ease of comparison, a multiplication operation is used as the benchmark for calculating computational complexity. The dimension of the target state is set to a, the dimension of the observed quantity is set to b, the average computational complexity of each equation in the state equation is c, the average computational complexity of each equation in the observed equation is d, and the computational complexity of inverting the n-dimensional matrix is n²e. Meanwhile, considering that the computational complexity of the function after differentiation is relatively small compared to before differentiation, it is assumed that the average computational complexity of the equation remains unchanged after differentiation. Based on the above settings, the computational complexity of each step in SecKF can be calculated as shown in Table 2. Using the same method to analyse EKF, the computational complexity of UKF is shown in Table 3. It can be concluded that the computational complexity of SecKF varies by countless orders of magnitude compared to EKF and is much lower than that of UKF.

4. Simulation Experiment and Analysis

In order to effectively verify the advantages of the method proposed in this paper, simulation experiments and physical simulation verification were conducted separately.

4.1. Simulation Experiment

The simulation uses a single satellite to continuously observe targets on the sea surface in two different scenarios and estimate their motion states. The number of satellite orbits at the initial observation time is shown in Table 4. The observation interval is 1 s, and the noise distribution for each observation is the same, with a standard deviation of 10⁻⁴ rad. The hardware configuration of the simulation runtime environment is an AMD R7-7840HS CPU and 16 GB of RAM. (Manufacturer: Advanced Micro Devices, Inc. Location: Santa Clara, CA, USA).

4.1.1. Target Motion Condition Settings

In order to better verify the universality of the proposed method, two water surface target motion scenarios (normal and constant power acceleration motions) were set up for simulation.

Scenario 1: Constant acceleration motion. In constant acceleration motion, changes in the motion state of target are mainly reflected in changes in speed and heading angle, so the kinematic equation of target can be approximated as follows:

(34)$\{\begin{array}{l}{x}_{m,j}=x_{m,j-1}+{\displaystyle \frac{a_m\mathrm{s}\mathrm{i}\mathrm{n}({\psi }_{j-1}+k_{\psi }{t}_d)-k_{\psi }(v_{j-1}+a_m{t}_d)\mathrm{c}\mathrm{o}\mathrm{s}({\psi }_{j-1}+k_{\psi }{t}_d)}{k_{\psi }^2}}+{\epsilon }_{x,j-1}\\ y_{m,j}=y_{m,j-1}+{\displaystyle \frac{a_m\mathrm{c}\mathrm{o}\mathrm{s}({\psi }_{j-1}+k_{\psi }{t}_d)+k_{\psi }(v_{j-1}+a_m{t}_d)\mathrm{s}\mathrm{i}\mathrm{n}({\psi }_{j-1}+k_{\psi }{t}_d)}{k_{\psi }^2}}+{\epsilon }_{y,j-1}\\ z_{m,j}=z_{m,j-1}+{\epsilon }_{z,j-1}\\ v_j=v_{j-1}+a_m{t}_d+{\epsilon }_{v,j-1}\\ {\psi }_j={\psi }_{j-1}+k_{\psi }{t}_d+{\epsilon }_{\psi ,j-1}\end{array}$

where t_d is the imaging time interval, a_m is the acceleration parameter, and k_ψ is the heading change rate parameter.

Scenario 2: Constant power acceleration movement. In scenarios where the target needs to accelerate to the maximum extent, the target is generally in a state of constant power acceleration movement. At this time, the target’s movement is mainly affected by factors such as engine acceleration and resistance. Therefore, the target’s kinematic equation is established as

(35)$\{\begin{array}{l}{x}_{m,j}=x_{m,j-1}+\sin {\psi }_{j-1}[v_{j-1}{t}_d+{\displaystyle \frac{(P-\eta v_{j-1}^3)t_d^2}{2m{v}_{j-1}}}]+{\epsilon }_{x,j-1}\\ y_{m,j}=y_{m,j-1}+\cos {\psi }_{j-1}[v_{j-1}{t}_d+{\displaystyle \frac{(P-\eta v_{j-1}^3)t_d^2}{2m{v}_{j-1}}}]+{\epsilon }_{y,j-1}\\ z_{m,j}=z_{m,j-1}+{\epsilon }_{z,j-1}\\ v_j=v_{j-1}+{\displaystyle \frac{(P-\eta v_{j-1}^3)t_d}{m{v}_{j-1}}}+{\epsilon }_{v,j-1}\\ {\psi }_j={\psi }_{j-1}+{\epsilon }_{\psi ,j-1}\end{array}$

where P is power, η is a coefficient related to resistance, and m is mass.

For Scenarios 1 and 2, the dynamic parameters of the ship target are shown in Table 5. The initial position coordinates of each scene are [0, 0, 0] m, the velocity is 10 m/s, and the heading angle is 0.25π rad. During the observation process, the distribution of the target’s motion state noise remains unchanged, with a standard deviation of 1 m for position noise, 0.1 m/s for velocity noise, and 10⁻⁵ rad for heading angle noise ($\mathit{Q}=\mathrm{diag}\left(1 {\mathrm{m}}^2,1 {\mathrm{m}}^2,1 {\mathrm{m}}^2,0.01 {\mathrm{m}}^2/{\mathrm{s}}^2,{10}^{-10} {\mathrm{rad}}^2\right)$). The initial target estimation distribution is ${\mathit{P}}_0=\mathrm{diag}\left(41 {\mathrm{m}}^2,41 {\mathrm{m}}^2,1 {\mathrm{m}}^2,1 {\mathrm{m}}^2/{\mathrm{s}}^2,{10}^{-4} {\mathrm{rad}}^2\right)$. The standard deviation of satellite angle observation error is 10⁻⁴ rad ($\mathit{R}=\mathrm{diag}\left({10}^{-8} {\mathrm{rad}}^2,{10}^{-8} {\mathrm{rad}}^2\right)$).

4.1.2. Analysis of Simulation Results

To validate the estimation performance of the proposed method, we compared it with the linearisation solution method [8], EKF, and UKF. Due to the significant randomness of single noise, 100,000 Monte Carlo simulations were performed to better reflect the method’s overall estimation performance. The simulations used mean absolute error (MAE) to characterise the positioning, speed, and heading angle determination accuracy, as calculated using the following formula:

(36)$\left\{\begin{array}{l}MA{E}_r={\displaystyle \frac{{\displaystyle \sum _{i=1}^{m_c}\left|{\widehat{\mathit{r}}}_{\mathit{m}\mathit{c},\mathit{i}}-{\mathit{r}}_{\mathit{m}}\right|}}{m_c}}\hfill \\ MA{E}_v={\displaystyle \frac{{\displaystyle \sum _{i=1}^{m_c}\left|{\widehat{v}}_{mc,i}-v_m\right|}}{m_c}}\hfill \\ MA{E}_{\psi }={\displaystyle \frac{{\displaystyle \sum _{i=1}^{m_c}\left|{\widehat{\psi }}_{mc,i}-{\psi }_m\right|}}{m_c}}\hfill \end{array}\right.$

where m_c is the number of Monte Carlo simulations, and ${\widehat{\mathit{r}}}_{\mathit{m}\mathit{c},\mathit{i}}$ and ${\widehat{v}}_{mc,i}$, ${\widehat{\psi }}_{mc,i}$, are the position, velocity, and heading angle estimates obtained from the i-th simulation.

In the simulation results, the green line with a circle indicates the estimation accuracy of the linearisation solution method, the blue line with a triangle indicates that of the EKF, the blue line with an asterisk indicates that of the UKF, and the red line with a five-pointed star indicates that of the SecKF proposed in this paper.

The simulation results for Scenario 1 are shown in Figure 10, Figure 11 and Figure 12 (the linearisation solution method cannot calculate velocity and heading angle, so its results are not plotted in the comparison plots). Figure 10 shows the error curve for estimating the target position, Figure 11 shows the error curve for estimating target velocity, and Figure 12 shows the error curve for estimating target heading angle estimation accuracies with observation time. Based on the figures and simulation data, it can be observed that as observation time increases (i.e., the number of observations increases), the positioning, velocity estimation, and heading angle determination accuracy of the EKF, UKF, and SecKF gradually improve, while the positioning error of the linearisation solution method remains almost unchanged. After 20 observations, the positioning accuracy of SecKF improves by approximately 57.1%, the velocity estimation accuracy by approximately 10.6%, and the heading angle determination accuracy by approximately 2.89%. Additionally, in this scenario, the positioning, velocity estimation, and heading angle determination accuracies of SecKF are similar to those of the EKF but inferior to those of the UKF.

In this simulation scenario, SecKF maintains estimation accuracy similar to EKF due to the fact that the target state equation does not undergo secant correction based on the Judgement Metrics, and the observation equation has a low degree of nonlinearity.

The simulation results for Scenario 2 are shown in Figure 13, Figure 14 and Figure 15. Figure 13 shows error curve for estimating target position, Figure 14 shows error curve for estimating target velocity, and Figure 15 shows terror curve for estimating target heading angle estimation accuracies with observation time. Based on Figure 13 and the simulation data, it can be observed that UKF dramatically improving after only a few observations, then showing a slower improvement after that, EKF and SecKF show a steady improvement throughout., while that of the linearisation solution method remain largely unchanged. After 20 observations, the positioning accuracy of the SecKF improves by approximately 28.6%. Additionally, the SecKF positioning accuracy is superior to the EKF but inferior to the UKF. More specifically, at the 20th observation, the SecKF positioning accuracy is approximately 12.5% higher than that of the EKF. Based on Figure 14 and simulation data, the error in target velocity estimation for both the SecKF and EKF stabilises around 0.1 m/s after more than seven observations, exhibiting similar estimation accuracy, though slightly lower than that of the UKF. As shown in Figure 15 and the simulation data, the error in target heading angle estimation for both SecKF and EKF gradually improves from 8 × 10⁻³ to 7.05 × 10⁻³ after 20 observations, achieving similar estimation accuracy, though slightly lower than the UKF.

In this simulation scenario, SecKF effectively improves the estimation accuracy compared to EKF because the target state equation can be corrected by secant based on the Judgement Metrics.

To further validate the applicability of the method proposed in this paper, we conducted simulations to evaluate positioning accuracy by randomly generating 50 satellite trajectories at orbital heights of 500–800 km within the target’s visible range, based on Scenario 2. All other conditions remained unchanged. Due to the significant randomness of single-shot noise, to better demonstrate the effectiveness of the method, each satellite trajectory was subjected to 100,000 Monte Carlo simulations. The simulation results are shown in Figure 16, Figure 17 and Figure 18. Based on the simulation results, it can be seen that under different satellite and target position relationships, the SecKF proposed in this paper can still significantly improve target positioning accuracy compared to the EKF, indicating its applicability to a wide range of scenarios.

Furthermore, we compared and analysed the computational speeds of various state estimation methods, with the results shown in Table 6 (using the Scenario 2 simulation). From this table, it can be seen that the linearisation solution method has the fastest computational speed. Among the filtering methods, the EKF has the fastest speed. The SecKF presented in this paper is slightly slower, but it has the same order of magnitude, and is about one order of magnitude faster than the UKF. At the same time, in order to further demonstrate the differences in computation time between different methods, embedded systems were used to test the running time of the methods. Simulate using the STM32F407 Kirin development board from Prechin to run the positioning algorithm. (Manufacturer: Shenzhen Puzhong Technology Co., Ltd. Location: Shenzhen, Guangdong, China). The results are shown in Table 6. According to the results, it can be concluded that in embedded systems, the method proposed in this paper improves the computational speed by approximately 0.03 s compared to UKF each step.

The main reason why SecKF and EKF have similar computation times and are significantly faster than UKF is that SecKF does not have an excessive increase in computational complexity compared to EKF, while UKF has a larger computational complexity.

By comparing Figure 10 and Figure 13 comprehensively, and combining with the analysis in the Appendix B, it can be shown that SecKF can improve the filtering estimation accuracy relative to EKF in the face of certain nonlinear models. The reason for this is that it uses secant approximation instead of first-order Taylor expansion, which improves the fitting accuracy of nonlinear models and reduces truncation errors. Meanwhile, according to the comparison of calculation time, it can be concluded that SecKF does not increase the computational complexity significantly compared to EKF. SecKF effectively achieved the first contribution mentioned in the introduction section. In addition, for some scenarios where the secant approximation does not have better fitting accuracy, SecKF’s adaptive judgement can effectively determine the coefficients, ensuring that the estimation accuracy of any model is not lower than EKF. SecKF effectively achieved the second contribution mentioned in the introduction section.

4.2. Physical Simulation Verification

In real situations, it is difficult to obtain satellite tracking observation data for ships, so an equivalent experiment is designed to verify the above method. The experiment uses a camera to simulate satellite observations and obtains camera attitude information through an attitude sensor. A robotic arm is used to drive the target point to simulate the motion of the target and obtain the true position of the target through measurements by the robotic arm. By observing the target point driven by the robotic arm through a camera, simulate the observation process of a moving target by a video satellite. The equivalent experimental scenario is shown in Figure 19.

The equivalent experimental camera and target position are shown in Table 5, and the target kinematic equation is

(37)$\{\begin{array}{l}{x}_{m,j}=x_{m,j-1}+t_d\sin {\psi }_{j-1}{v}_{j-1}+{\displaystyle \frac{k_1{t}_d^2\mathrm{s}\mathrm{i}\mathrm{n}{\psi }_{j-1}}{v_{j-1}}}-k_2{t}_d^2\sin {\psi }_{j-1}{v}_{j-1}^2+{\epsilon }_{x,j-1}\\ y_{m,j}=y_{m,j-1}+{\epsilon }_{z,j-1}\\ z_{m,j}=z_{m,j-1}+t_d\cos {\psi }_{j-1}{v}_{j-1}+{\displaystyle \frac{k_1{t}_d^2}{v_{j-1}}}-k_2{t}_d^2{v}_{j-1}^2+{\epsilon }_{y,j-1}\\ v_j=v_{j-1}+{\displaystyle \frac{2k_1{t}_d^2}{v_{j-1}}}-2k_2{t}_d^2{v}_{j-1}^2+{\epsilon }_{v,j-1}\\ {\psi }_j={\psi }_{j-1}+{\epsilon }_{\psi ,j-1}\end{array}$

The initial conditions and error distribution are shown in Table 7. The target image in the first second of the first experiment is shown in Figure 20. The position coordinates of the observation camera are [−0.2, 2, 1] m, with an observation error of 10⁻⁴ rad. The initial position coordinates are [0, 0, 0] m, the velocity is 0.1 m/s, and the heading angle is 0.25π rad. The noise distribution of the target motion state remains unchanged during the observation process, with a standard deviation of 0.01 m for position noise, 0.001 m/s for velocity noise, and 10⁻⁷ rad for heading angle noise. Considering the randomness of errors, 10 repeated experiments were conducted, and the average errors of each method are shown in Table 8.

Based on actual data analysis, it can be concluded that the positioning accuracy of SecKF in this article is higher than EKF and linearization methods, slightly lower than UKF, and consistent with simulation results, indicating that SecKF has achieved the desired effect.

5. Conclusions

This article proposes a Secant Kalman filtering method (SecKF) based on secant calculation to address the issue of truncation error in using EKF for state estimation of moving targets. This method improves estimation accuracy without significantly increasing computational complexity compared to EKF. Additionally, based on an analysis of the advantages of secant correction, this paper provides a method for selecting correction elements. The simulation results show that under certain circumstances, the method proposed in this paper can improve the positioning accuracy by 28.6% compared to EKF, and at the same time, it can increase the computation speed by 10 times compared to UKF. The method proposed in this article can be applied to tasks such as tracking and positioning of moving targets, autonomous navigation of spacecrafts, etc., that require satellite onboard computing. Meanwhile, for some other state estimation tasks that require high precision, good timeliness, and limited computational power, the method proposed in this paper can also be widely applied. In future research, we plan to further analyse the range of correction elements to expand the model range where SecKF has higher accuracy compared to EKF.

Author Contributions

Conceptualization, X.B. and H.S.; methodology, X.B.; software, X.B. and H.S.; validation, X.B. and J.Z.; formal analysis, X.B. and H.H.; investigation, X.B. and C.F.; resources, C.F. and H.S.; data curation, H.S.; writing—original draft preparation, X.B.; writing—review and editing, H.S. and C.F.; visualisation, X.B. and J.Z.; supervision, H.S.; project administration, C.F.; funding acquisition, C.F. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

KF	Kalman Filter
EKF	Extended Kalman Filter
UKF	Unscented Kalman Filter
PF	Particle Filter
CKF	Cubature Kalman Filter
GSF	Gaussian Sum Filter
LOS	line-of-sight
SecKF	Secant Kalman Filtering

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Sea surface target tracking observation process.

Figure 2 Schematic diagram of the relative position between the target and satellite.

Figure 3 Schematic diagram of EKF approximation process.

Figure 4 Truncation error of first-order expansion for some typical functions.

Figure 5 Schematic diagram of multidimensional function secant calculation process.

Figure 6 Schematic diagram of theoretical secant line approximation process.

Figure 7 Schematic diagram of truncation error of actual secant line approximation process.

Figure 8 Compare with secant approximation and first-order Taylor series expansion approximation truncation errors.

Figure 10 Scenario 1 target MAE_r chart showing changes with observation time.

Figure 11 Scenario 1 target MAE_v chart showing changes with observation time.

Figure 12 Scenario 1 target MAE_ψ chart showing changes with observation time.

Figure 13 Scenario 2 target MAE_r chart showing changes with observation time.

Figure 14 Scenario 2 target MAE_v chart showing changes with observation time.

Figure 15 Scenario 2 target MAE_ψ chart showing changes with observation time.

Figure 16 Graph showing changes in MAE_r accuracy for random satellite targets with observation time.

Figure 17 Graph showing changes in MAE_v accuracy for random satellite targets with observation time.

Figure 18 Graph showing changes in MAE_ψ accuracy for random satellite targets with observation time.

Figure 19 Equivalent experimental scenario.

Figure 20 Initial camera image.

Table 1

Common nonlinear models.

Model Type	Basic Formula	Does Secant Correction Reduce Truncation Error
Kinematic model under fluid resistance	$\left\{\begin{array}{l}x={\displaystyle \frac{\ln \left(C_1-kt\right)}{k}}+C_2\hfill \\ \dot{x}={\displaystyle \frac{1}{C_1-kt}}\hfill \end{array}\right.$	Positive
Kinematic model of re-entry velocity variation	$\left\{\begin{array}{l}x=k\ln \left[\cosh \left(\sqrt{{\displaystyle \frac{k}{g}}}t\right)\right]\hfill \\ \dot{x}=\sqrt{{\displaystyle \frac{g}{k}}}\tanh \left(\sqrt{{\displaystyle \frac{k}{g}}}t\right)\hfill \end{array}\right.$	Positive
Kinematic model of constant thrust within the atmosphere	$\left\{\begin{array}{l}x=k\ln \left[\cosh \left(\sqrt{{\displaystyle \frac{k}{a}}}t\right)\right]\hfill \\ \dot{x}=\sqrt{{\displaystyle \frac{a}{k}}}\tanh \left(\sqrt{{\displaystyle \frac{k}{a}}}t\right)\hfill \end{array}\right.$	Positive
Kinematic model of simple harmonic motion	$\left\{\begin{array}{l}x=A\cos \left(\omega t+\psi \right)\hfill \\ \dot{x}=-A\sin \left(\omega t+\psi \right)\hfill \end{array}\right.$	Negative
Angle measurement observation model	$\theta =\arctan {\displaystyle \frac{y}{x}}$	Positive
Pseudo-range observation model	$\left\{\begin{array}{l}{\theta }_1={\displaystyle \frac{\sqrt{{\left(r_{x1}-x_1\right)}^2+{\left(r_{y1}-x_2\right)}^2+{\left(r_{z1}-x_3\right)}^2}}{c}}\hfill \\ {\theta }_2={\displaystyle \frac{\sqrt{{\left(r_{x2}-x_1\right)}^2+{\left(r_{y2}-x_2\right)}^2+{\left(r_{z2}-x_3\right)}^2}}{c}}\hfill \\ {\theta }_3={\displaystyle \frac{\sqrt{{\left(r_{x3}-x_1\right)}^2+{\left(r_{y3}-x_2\right)}^2+{\left(r_{z3}-x_3\right)}^2}}{c}}\hfill \end{array}\right.$	Positive
Time difference in arrival observation model	$\left\{\begin{array}{l}{\theta }_1={\tau }_1-{\tau }_2\hfill \\ {\theta }_2={\tau }_2-{\tau }_3\hfill \\ {\theta }_3={\tau }_3-{\tau }_1\hfill \\ {\tau }_1={\displaystyle \frac{\sqrt{{\left(r_{x1}-x_1\right)}^2+{\left(r_{y1}-x_2\right)}^2+{\left(r_{z1}-x_3\right)}^2}}{c}}\hfill \\ {\tau }_2={\displaystyle \frac{\sqrt{{\left(r_{x2}-x_1\right)}^2+{\left(r_{y2}-x_2\right)}^2+{\left(r_{z2}-x_3\right)}^2}}{c}}\hfill \\ {\tau }_3={\displaystyle \frac{\sqrt{{\left(r_{x3}-x_1\right)}^2+{\left(r_{y3}-x_2\right)}^2+{\left(r_{z3}-x_3\right)}^2}}{c}}\hfill \end{array}\right.$	Positive

Table 2

Calculation volume of each step in SecKF.

Step	Computational
Prediction	ac
Expansion	a2c + b2d + 2a3
Initial fusion	2a2b + 2ab2 + b2e
Initial update	ab + a 2 b + a 3
Secant calculation	2ac + 2bd
Element determination	2ac + 2bd
Secondary fusion	2a2b + 2ab2 + b2e
Secondary update	ab + a 2 b + a 3

Table 3

Comparison of computational complexity of filtering methods.

Method	Computational
EKF	ab + ac + 3a3 + 3a2b + a2c + 2ab2 + b2d + b2e
UKF	24a + ab + a3 + 3a2b + 2a2c + 2ab2 + 2a4 + 2ab2d + 2a3c + b2e
SecKF	2ab + 5ac + 4bd + 4a3 + 6a2b + 4ab2 + a2c + b2d + 2b2e

Table 4

Number of satellite initial moment orbits.

Parameter	Value
Semi-major axis	7071.004 km
Inclination	97.57 deg
Eccentricity	0.000001
RAAN	300 deg
Perigee	203 deg
Mean anomaly	157 deg

Table 5

Parameters of ship target dynamics.

Parameter	Value
am	0.5 m/s2
kψ	1 × 10−5 rad/s2
P	2.4 × 105 w
η	10 Ns2/m2
m	1 × 104 kg

Table 6

Comparison of calculation speed performance.

Method	Computational Time (PC) (s)	Computational Time (Embedded System) (s)
Linearisation solution method	2.578 × 10−6	1.203 × 10−4
UKF	7.675 × 10−4	0.0381
EKF	2.405 × 10−5	1.562 × 10−3
SecKF	7.367 × 10−5	3.021 × 10−3

Table 7

Parameters of ship target dynamics.

Parameter	Value
k 1	0.01 m2/s3
k 2	0.0005 m−1
td	1 s

Table 8

Comparison of calculation speed performance.

Method	Positioning Accuracy (m)
Linearisation solution method	0.036
UKF	0.017
EKF	0.023
SecKF	0.020

Appendix A

The expansion of H_s,j is(A1)$${\mathit{H}}_{\mathit{s},\mathit{j}}=\left[\begin{array}{ccccc}{h}_{11}& h_{12}& 0& 0& 0\\ h_{21}& h_{22}& h_{23}& 0& 0\end{array}\right]$$ where(A2)$$\{\begin{array}{l}{h}_{11}={\displaystyle \frac{s_{y,j}-{\overline{y}}_{m,j}}{q_1^2}}\\ h_{12}=-{\displaystyle \frac{s_{x,j}-{\overline{x}}_{m,j}}{q_1^2}}\\ h_{21}=-{\displaystyle \frac{(s_{x,j}-{\overline{x}}_{m,j})(s_{z,j}-{\overline{z}}_{m,j})}{q_1{q}_2}}\\ h_{22}=-{\displaystyle \frac{(s_{y,j}-{\overline{y}}_{m,j})(s_{z,j}-{\overline{z}}_{m,j})}{q_1{q}_2}}\\ h_{23}={\displaystyle \frac{q_1}{{(s_{z,j}-{\overline{z}}_{m,j})}^2+q_1}}\\ q_1=\sqrt{{(s_{x,j}-{\overline{x}}_{m,j})}^2+{(s_{y,j}-{\overline{y}}_{m,j})}^2}\\ q_2={(s_{z,j}-{\overline{z}}_{m,j})}^2+{(s_{x,j}-{\overline{x}}_{m,j})}^2+{(s_{y,j}-{\overline{y}}_{m,j})}^2\end{array}$$

The expansion of H_sg,j is(A3)$${\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}=\mathrm{SecLine}\left[h_s\left(\mathit{\theta }\right),{\overline{\mathit{\theta }}}_{\mathit{j}},{\widehat{\mathit{\theta }}}_{\mathit{j}}\right]$$

Additionally, it can be further expanded to(A4)$${\mathit{H}}_{\mathit{s}\mathit{g},\mathit{j}}=\left[\begin{array}{ccccc}{g}_{h, 11}& g_{h, 12}& 0& 0& 0\\ g_{h, 21}& g_{h, 22}& g_{h, 23}& 0& 0\end{array}\right]$$ where(A5)$$\{\begin{array}{l}{g}_{h, 11}={\displaystyle \frac{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{y}}_{m,j}-s_{y,j}}{{\hat{x}}_{m,j}-s_{x,j}}-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{y}}_{m,j}-s_{y,j}}{{\overline{x}}_{m,j}-s_{x,j}}}{{\hat{x}}_{m,j}-{\overline{x}}_{m,j}}}\\ g_{h, 12}={\displaystyle \frac{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\hat{y}}_{m,j}-s_{y,j}}{{\hat{x}}_{m,j}-s_{x,j}}-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{y}}_{m,j}-s_{y,j}}{{\hat{x}}_{m,j}-s_{x,j}}}{{\hat{y}}_{m,j}-{\overline{y}}_{m,j}}}\\ g_{h, 21}={\displaystyle \frac{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{z}}_{m,j}-s_{z,j}}{\sqrt{{({\hat{x}}_{m,j}-s_{x,j})}^2+{({\overline{y}}_{m,j}-s_{y,j})}^2}}-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{z}}_{m,j}-s_{z,j}}{\sqrt{{({\overline{x}}_{m,j}-s_{x,j})}^2+{({\overline{y}}_{m,j}-s_{y,j})}^2}}}{{\hat{x}}_{m,j}-{\overline{x}}_{m,j}}}\\ g_{h, 22}={\displaystyle \frac{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{z}}_{m,j}-s_{z,j}}{\sqrt{{({\hat{x}}_{m,j}-s_{x,j})}^2+{({\hat{y}}_{m,j}-s_{y,j})}^2}}-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{z}}_{m,j}-s_{z,j}}{\sqrt{{({\hat{x}}_{m,j}-s_{x,j})}^2+{({\overline{y}}_{m,j}-s_{y,j})}^2}}}{{\hat{y}}_{m,j}-{\overline{y}}_{m,j}}}\\ g_{h, 23}={\displaystyle \frac{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\hat{z}}_{m,j}-s_{z,j}}{\sqrt{{({\hat{x}}_{m,j}-s_{x,j})}^2+{({\hat{y}}_{m,j}-s_{y,j})}^2}}-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\overline{z}}_{m,j}-s_{z,j}}{\sqrt{{({\hat{x}}_{m,j}-s_{x,j})}^2+{({\hat{y}}_{m,j}-s_{y,j})}^2}}}{{\hat{z}}_{m,j}-{\overline{z}}_{m,j}}}\end{array}$$

Appendix B

Analyse the state equation in Scenario 2. Among them, x_m, y_m, z_m, and ψ are all linear functions themselves, so there is no need to consider the fitting situation. The second derivative of x_m and y_m with respect to v is(A6)$${\displaystyle \frac{{\partial }^2{x}_m}{\partial v^2}}={\displaystyle \frac{{\partial }^2{y}_m}{\partial v^2}}={\displaystyle \frac{P}{m{v}^3}}-{\displaystyle \frac{\eta t_d^2}{m}}$$

By combining Equation (A6) with simulation parameters, it can be concluded that the second derivative of x_m and y_m with respect to v in the interval is non-negative, meaning that the functions of x_m and y_m with respect to v in the interval are convex. According to Judgement Metric 1, it is determined that secant fitting correction can be performed on the part of x_m and y_m with respect to v to achieve higher fitting accuracy. The second derivative of v with respect to v is(A7)$${\displaystyle \frac{{\partial }^2{x}_m}{\partial v^2}}={\displaystyle \frac{{\partial }^2{y}_m}{\partial v^2}}={\displaystyle \frac{P}{m{v}^3}}-{\displaystyle \frac{\eta t_d^2}{m}}$$

By analogy with Equations (A6) and (A7), it can be concluded that the function of v with respect to v within the interval is convex. Based on Judgement Metric 1, it can be determined that secant fitting correction can be performed on the part of v with respect to v to achieve higher fitting accuracy. The second derivative of x_m and y_m with respect to ψ is(A8)$$\left\{\begin{array}{l}{\displaystyle \frac{{\partial }^2{x}_m}{\partial {\psi }^2}}=-v_{j-1}{t}_d\sin {\psi }_{j-1}\hfill \\ {\displaystyle \frac{{\partial }^2{y}_m}{\partial {\psi }^2}}=-v_{j-1}{t}_d\cos {\psi }_{j-1}\hfill \end{array}\right.$$

By combining Equation (A8) with simulation parameters, it can be concluded that the second derivative of x_m and y_m with respect to ψ within the interval is negative, meaning that the functions of x_m and y_m with respect to ψ within the interval are concave functions and cannot be determined using Judgement Metric 1. Furthermore, considering Judgement Metric 2, calculate that the Lipschitz constants of x_m and y_m with respect to the ψ function are greater than or equal to 1. Therefore, criterion 2 is not met, and the corresponding parts cannot be corrected by secant fitting. Therefore, the modified state transition matrix can be obtained as(A9)$${\mathit{F}}_{\mathit{j}}=\left[\begin{array}{ccccc}1& 0& 0& g_{f, 14}& f_{15}\\ 0& 1& 0& g_{f, 24}& f_{25}\\ 0& 0& 1& 0& 0\\ 0& 0& 0& g_{f, 44}& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$ where(A10)$$\{\begin{array}{l}{g}_{f, 14}={\displaystyle \frac{f_{xv}({\overline{x}}_{m,j-1},{\hat{v}}_{j-1},{\overline{\psi }}_{j-1})-f_{xv}({\overline{x}}_{m,j-1},{\overline{v}}_{j-1},{\overline{\psi }}_{j-1})}{{\hat{v}}_{j-1}-{\overline{v}}_{j-1}}}\\ g_{f, 24}={\displaystyle \frac{f_{yv}({\hat{x}}_{m,j-1},{\hat{v}}_{j-1},{\overline{\psi }}_{j-1})-f_{yv}({\hat{x}}_{m,j-1},{\overline{v}}_{j-1},{\overline{\psi }}_{j-1})}{{\hat{v}}_{j-1}-{\overline{v}}_{j-1}}}\\ g_{f, 44}={\displaystyle \frac{f_{vv}\left({\hat{v}}_{j-1}\right)-f_{vv}\left({\overline{v}}_{j-1}\right)}{{\hat{v}}_{j-1}-{\overline{v}}_{j-1}}}\\ f_{15}=\cos {\hat{\psi }}_{j-1}[{\hat{v}}_{j-1}{t}_d+{\displaystyle \frac{(P-\eta {\hat{v}}_{j-1}^3)t_d^2}{2m{\hat{v}}_{j-1}}}]\\ f_{25}=-\sin {\hat{\psi }}_{j-1}[{\hat{v}}_{j-1}{t}_d+{\displaystyle \frac{(P-\eta {\hat{v}}_{j-1}^3)t_d^2}{2m{\hat{v}}_{j-1}}}]\\ f_{xv}({\theta }_1,{\theta }_2,{\theta }_3)={\theta }_1+\sin {\theta }_3[{\theta }_2{t}_d+{\displaystyle \frac{(P-\eta {\theta }_2^3)t_d^2}{2m{\theta }_2}}]\\ f_{yv}({\theta }_1,{\theta }_2,{\theta }_3)={\theta }_1+\cos {\theta }_3[{\theta }_2{t}_d+{\displaystyle \frac{(P-\eta {\theta }_2^3)t_d^2}{2m{\theta }_2}}]\\ f_{vv}\left({\theta }_1\right)={\theta }_1+{\displaystyle \frac{(P-\eta {\theta }_2^3)t_d}{m{v}_{j-1}}}\end{array}$$

Similarly, for the observation equation in Scenario 2, according to the results given in Table 1, the equations for α_j and β_j are convex functions with respect to x_m, y_m, and z_m. Based on Judgement Metric 1, secant fitting correction can be performed.

References

1. Zhang, H.; Liu, W.; Zhang, L.; Meng, Y.; Han, W.; Song, T.; Yang, R. An allocation strategy integrated power, bandwidth, and subchannel in a RCC network. Def. Technol.; 2025; [DOI: https://dx.doi.org/10.1016/j.dt.2025.08.010]

2. Cao, H.; Gao, W.; Zhang, X.; Liu, X.; Fan, B.; Li, S. Overview of ZY-3 satellite research and application. Proceedings of the 63rd IAC (International Astronautical Congress); Naples, Italy, 1–5 October 2012; pp. 1-5.

3. Wang, Z.; Luo, J.-A.; Zhang, X.-P. A novel location-penalized maximum likelihood estimator for bearing-only target localization. IEEE Trans. Signal Process.; 2012; 60, pp. 6166-6181. [DOI: https://dx.doi.org/10.1109/TSP.2012.2218809]

4. Liu, J. Research on Single-Base Angle Measurement Passive Positioning Technology. Master’s Thesis; University of Electronic Science and Technology of China: Chengdu, China, 2020.

5. Fan, L.J.; Wang, Y.; Yang, W.T.; Yu, L.J.; Zhang, G.B. GFDM-1 Satellite System Design and Technical Characteristics. Spacecr. Eng.; 2021; 30, pp. 10-19. [DOI: https://dx.doi.org/10.3969/j.issn.1673-8748.2021.03.002]

6. Zhang, H.; Liu, W.; Zhang, Q.; Liu, B. Joint customer assignment, power allocation, and subchannel allocation in a UAV-based joint radar and communication network. IEEE Internet Things J.; 2024; 11, pp. 29643-29660. [DOI: https://dx.doi.org/10.1109/JIOT.2024.3397029]

7. Zhang, X. Study on Moving Objects Intelligent Sensing and Tracking Control for Video Satellite. Ph.D. Thesis; National University of Defense Technology: Changsha, China, 2017.

8. Bai, X.; Fan, C.; Song, H.; Zhang, Y. Space Target Positioning Under Gaze Mode of Collaborative Distributed Video Satellite. Proceedings of the 2024 International Symposium on Intelligent Robotics and Systems (ISoIRS); Changsha, China, 14–16 June 2024; pp. 227-231.

9. Bai, X.; Song, H.; Fan, C.; Hao, L.; Yang, Y. A System Error Self-Correction Target-Positioning Method in Video Satellite Observation. Remote Sens.; 2025; 17, 2935. [DOI: https://dx.doi.org/10.3390/rs17172935]

10. Bai, X.; Song, H.; Fan, C. A Space Target Positioning Method Under Distributed Video Satellite with LOS Angle and Orbit Position Error. Adv. Space Res.; 2025; [DOI: https://dx.doi.org/10.1016/j.asr.2025.10.034]

11. Lei, T.; Guan, B.; Liang, M.; Liu, Z.; Liu, J.; Shang, Y.; Yu, Q. Motion measurements of explosive shock waves based on an event camera. Opt. Express; 2024; 32, pp. 15390-15409. [DOI: https://dx.doi.org/10.1364/OE.506662]

12. Gong, B.C. Research on Angles-only Relative Orbit Determination Algorithms for Spacecraft Autonomous Rendezvous. Ph.D. Thesis; Northwestern Polytechnical University: Xi’an, China, 2016.

13. Wang, D.Y.; Hou, B.W.; Wang, J.Q.; Ge, D.M.; Li, M.D.; Xu, C.; Zhou, H.Y. State estimation method for spacecraft autonomous navigation: Review. Acta Aeronaut. Astronaut. Sin.; 2021; 42, pp. 72-89. [DOI: https://dx.doi.org/10.7527/S1000-6893.2020.24310]

14. Garg, S.K. Initial Relative-Orbit Determination Using Second-Order Dynamics and Line-of-Sight Measurements. Master’s Thesis; Auburn University: Auburn, AL, USA, 2015.

15. Gong, B.; Wang, S.; Li, S.; Li, X. Review of space relative navigation based on angles-only measurements. Astrodynamics; 2023; 7, pp. 131-152. [DOI: https://dx.doi.org/10.1007/s42064-022-0152-2]

16. Ge, D.; Wang, D.; Zou, Y.; Shi, J. Motion and inertial parameter estimation of non-cooperative target on orbit using stereo vision. Adv. Space Res.; 2020; 66, pp. 1475-1484. [DOI: https://dx.doi.org/10.1016/j.asr.2020.05.029]

17. Leong, P.H.; Arulampalam, S.; Lamahewa, T.A.; Abhayapala, T.D. A Gaussian-sum based cubature Kalman filter for bearings-only tracking. IEEE Trans. Aerosp. Electron. Syst.; 2013; 49, pp. 1161-1176. [DOI: https://dx.doi.org/10.1109/TAES.2013.6494405]

18. Yin, J.; Zhang, J.; Zhuang, Z. Gaussian sum PHD filtering algorithm for nonlinear non-Gaussian models. Chin. J. Aeronaut.; 2008; 21, pp. 341-351. [DOI: https://dx.doi.org/10.1016/s1000-9361(08)60045-x]

19. Kotecha, J.H.; Djuric, P.M. Gaussian sum particle filtering. IEEE Trans. Signal Process.; 2003; 51, pp. 2602-2612. [DOI: https://dx.doi.org/10.1109/TSP.2003.816754]

20. Liu, Y. The Robust and High Accurate Fusion Filter for Trajectory Target Realtime Tracking. Ph.D. Thesis; National University of Defense Technology: Changsha, China, 2011.

21. Lanz, O. Approximate bayesian multibody tracking. IEEE Trans. Pattern Anal. Mach. Intell.; 2006; 28, pp. 1436-1449. [DOI: https://dx.doi.org/10.1109/TPAMI.2006.177] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16929730]

22. Eang, C.; Lee, S. An integration of deep neural network-based extended Kalman filter (DNN-EKF) method in ultra-wideband (UWB) localization for distance loss optimization. Sensors; 2024; 24, 7643. [DOI: https://dx.doi.org/10.3390/s24237643]

23. Wang, S.; Men, C.; Li, R.; Yeo, T.-S. A maneuvering extended target tracking IMM algorithm based on second-order EKF. IEEE Trans. Instrum. Meas.; 2024; 73, 8505011. [DOI: https://dx.doi.org/10.1109/TIM.2024.3418076]

24. Zhao, W.; Zhao, S.; Zhang, G.; Liu, G.; Meng, W. FL-EKF-based cooperative localization method for multi-AUVs. IEEE Internet Things J.; 2024; 11, pp. 30742-30753. [DOI: https://dx.doi.org/10.1109/JIOT.2024.3414501]

25. Simon, D. Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches; John Wiley & Sons: Hoboken, NJ, USA, 2006.

26. Wang, D.; Zhang, L.; Zhu, W.; Shi, J.; Huang, M.; Zou, Y. Autonomous relative navigation for noncooperative target using stereo vision measurement. Sci. Sin. Phys. Mech. Astron.; 2019; 49, 024509. [DOI: https://dx.doi.org/10.1360/SSPMA2018-00294]

27. Bell, B.M.; Cathey, F.W. The iterated Kalman filter update as a Gauss-Newton method. IEEE Trans. Autom. Control; 1993; 38, pp. 294-297. [DOI: https://dx.doi.org/10.1109/9.250476]

28. Wang, D.; Li, M.; Huang, X.; Li, J. Kalman filtering for a quadratic form state equality constraint. J. Guid. Control. Dyn.; 2014; 37, pp. 951-958. [DOI: https://dx.doi.org/10.2514/1.62890]

29. School of Mathematical Sciences, East China Normal University. Mathematical Analysis; Higher Education Press: Beijing, China, 2019.

30. Zhang, M. The Secant Algorithm-A Hardware Algorithm of theComputer to Evaluate Approximation of the Functionsat High Speed. J. Natl. Univ. Def. Technol.; 1983; pp. 77-81.

31. Cobzaş, Ş.; Miculescu, R.; Nicolae, A. Lipschitz Functions; Springer: Berlin/Heidelberg, Germany, 2019; Volume 2241.