Three-Dimensional Extended Target Tracking and Shape Learning Based on Double Fourier Series and Expectation Maximization

TL;DR

This paper introduces a new method for tracking 3D extended targets using shape learning and optimization techniques.

Contribution

The novel approach uses ECM framework and DFS to jointly estimate kinematics, extent, and orientation without predefined shape evolution.

Findings

The proposed method models 3D shape using radial function via DFS expansion.

ECM algorithm effectively estimates kinematics and shape parameters in experiments.

Abstract

This paper investigates the problem of tracking targets with unknown but fixed 3D star-convex shapes using point cloud measurements. While existing methods typically model shape parameters as random variables evolving according to predefined prior models, this evolution process is often unknown in practice. We propose a particular approach within the Expectation Conditional Maximization (ECM) framework that circumvents this limitation by treating shape-defining quantities as parameters estimated directly via optimization. The objective is the joint estimation of target kinematics, extent, and orientation in 3D space. Specifically, the 3D shape is modeled using a radial function estimated via double Fourier series (DFS) expansion, and orientation is represented using the compact, singularity-free axis-angle method. The ECM algorithm facilitates this joint estimation: an Unscented Kalman Smoother infers kinematics in the E-step, while the M-step estimates DFS shape parameters and rotation angles by minimizing regularized cost functions, promoting robustness and smoothness. The effectiveness of the proposed algorithm is substantiated through two experimental evaluations.

Figures

Full text

1. Introduction

The accurate perception and tracking of targets in three-dimensional (3D) environments are a fundamental requirement for many modern intelligent systems, including autonomous vehicles, robotic platforms, and advanced surveillance systems [1,2,3,4]. While traditional point target tracking algorithms assume that one target only gives rise to one measurement, they are insufficient when dealing with objects possessing significant spatial extent relative to sensor resolution, known as extended target [5,6]. The effective tracking of these extended targets demands algorithms that can estimate not only their kinematic state (position, velocity) but also their physical characteristics, including shape and orientation, within a 3D space.

A substantial body of work considers the extended target tracking (ETT) problem, especially in 2D space [5,6]. Corresponding methods process 2D measurements based on diverse extent models. One common approach employs random matrix to model the target’s shape as an ellipse [7,8]. These extent models often assume that the extent state follows a specific probability distribution, such as a Gaussian distribution. However, this underlying distribution is frequently unknow in practice. Addressing this limitation, paper [9] presents an alternative method where extent characteristics, including semi-axis lengths, are treated as parameters for direct estimation. For complex star-convex targets, more sophisticated shape modeling methods, such as B-splines [10] and Gaussian processes (GPs) [11,12], have been developed. Notably, the random hyper-surface model (RHM) is a typical approach for representing arbitrary star-convex shapes within the frequency domain [13,14]. This model defines the target’s spatial extent using a Fourier series expansion, with its coefficients jointly estimated alongside the target’s kinematics. RHM offers a parameterization particularly well-suited for star-convex extents, effectively balancing representational accuracy with computational complexity. Numerous adaptations of such Fourier-based approaches have subsequently been explored across diverse application settings [15,16,17,18,19,20,21].

When 3D surface measurements from devices like multilayer LiDAR sensor or depth cameras are available, it becomes possible to model the target’s 3D shape. A few approaches to 3D extended target tracking involved measurement models using basic geometric shapes [22] or restricted shape estimation [23,24,25]. The methods are suitable for when the detailed shape information is not critical for the application. For instance, where the target’s shape is unknown beforehand, random matrix models can represent the target using forms like ellipsoid, or rectangular cuboids. However, the representation power of RM models remains limited to these relatively simple geometric shapes. To address more complex geometries, Gaussian Processes have been developed for 3D space, allowing the modeling of arbitrary shaped star-convex targets [26,27]. While GPs effectively capture complex 3D shapes, their representation (often via a radius function) typically results in a high-dimensional state vector. This high dimensionality poses challenges for filtering algorithms and increases the computational burden. In contrast, the double Fourier series (DFS) shape model proposed in [28] utilizes a shape vector of significantly lower dimensionality compared to GP models, while also enabling the estimation of both shape and pose. It’s worth noting that a common characteristic across many of these 3D shape measurement models is the modeling of the shape-characterizing state as a random variable, typically assumed to follow a Gaussian distribution.

Differing from traditional 3D methods that rely on prior shape evolution models, our approach models the parameters defining the extent state directly, treating them as quantities with clear physical meaning. Within the Expectation Conditional Maximization (ECM) framework, we propose an extended target tracking method for the joint estimation of kinematics, extent, and orientation in 3D space. The main contributions are listed as follows.

(1)The unknown but fixed 3D shape is represented via a radial function in spherical coordinates, with a Double Fourier Series (DFS) expansion employed for its modeling. This approach converts the continuous radial function into a parametric representation characterized by a finite set of coefficients.(2)The axis-angle representation constitutes a compact and singularity-free approach for parameterizing 3D orientation. It describes a rotation by specifying a unit vector, which defines the axis of rotation, and a scalar angle, which denotes the magnitude of rotation about this axis. Furthermore, this method exhibits both geometric intuitiveness and notable computational efficiency.(3)Joint estimation of the target’s kinematic state and shape parameters is achieved via the Expectation Conditional Maximization (ECM) framework. The E-step infers the kinematic state using an unscented Kalman smoother filter, while the M-step estimates shape and rotation parameters by minimizing separate cost functions with added regularization for robustness and smoothness.

The remainder of this paper is organized as follows. Section 2 formulates the problem and presents the system, including the motion model, shape representation, and measurement model. Section 3 derives an ECM-based algorithm for 3D ETT problem. Section 4 presents simulation results to validate the proposed approach across two distinct scenarios. Finally, Section 5 summarizes the main findings and concludes the article.

Notation: Throughout this paper, the superscripts “−1” and “ $[eqn]$ ” represent the inverse and transpose of a matrix, respectively; $[eqn]$ represents the Gaussian distribution with mean $[eqn]$ and covariance $[eqn]$ ; $[eqn]$ and $[eqn]$ is mathematical expectation and conditional expectation. The notation $[eqn]$ refers to the Kronecker product. $[eqn]$ denote the identity matrix with proper dimensions; $[eqn]$ and $[eqn]$ denote the determinant and trace of matrices, respectively; $[eqn]$ denotes the Euclidean norm.

2. Problem Formulation

This section first introduces a model for the spatial extent of 3D target, which employs a radial function parameterized via a Double Fourier Series (DFS). The target’s orientation is then described using the axis–angle representation, a method characterized by its compactness and singularity-free nature. Following this, two kinematic state evolution models specifically designed for 3D extended target tracking are presented. Finally, the paper details the measurement model that integrates this DFS-based shape representation with the selected axis–angle orientation.

2.1. DFS-Based Shape Model for 3D Target

Advanced sensors, such as depth cameras and LiDAR, can generate detailed three-dimensional (3D) point cloud data, extracting precise shape estimates from such data is crucial for object classification and predicting future behavior. To effectively learn shapes from these point clouds, it is essential to formulate a suitable object extent description method that possesses high representational power for diverse 3D shapes while also being sufficiently compact to support efficient online tracking algorithms.

This section details the modeling of the target shape in spherical coordinates via a radial function $[eqn]$ . Such functions relate angular coordinate pairs $[eqn]$ to a radial distance in each orientation and are adept at representing star-convex geometries. In our methodology, this radial function is modeled using a Double Fourier Series (DFS) expansion, the specific form of which is given by:

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$

Here, $[eqn]$ is the base radial parameter, $[eqn]$ are the Fourier coefficients, and $[eqn]$ are the number of angle $[eqn]$ and $[eqn]$ . $[eqn]$ is azimuth and $[eqn]$ is elevation. The shape parameter $[eqn]$ collects $[eqn]$ and all $[eqn]$ coefficients, resulting in a dimension of $[eqn]$ .

2.2. Orientation Representation (Axis-Angle)

To represent the target’s three-dimensional (3D) orientation, this work utilizes the axis-angle parameterization, often expressed compactly as an orientation vector $[eqn]$ . In this formulation, $[eqn]$ is a unit vector denoting the axis of orientation, and $[eqn]$ is the scalar angle defining the magnitude of rotation orientation according to the right-hand rule.

This choice of representation is underpinned by several compelling advantages pertinent to tracking applications. Firstly, the axis–angle representation (particularly in its three-parameter rotation vector form $[eqn]$ ) provides a geometrically intuitive and unique mapping for any 3D rotation (up to 2π periodicity in angle), without requiring the satisfaction of complex inter-parameter constraints for its validity—unlike, for instance, the nine-parameter rotation matrix which must adhere to strict orthogonality and determinant conditions. Secondly, the practical application and conversion of this representation are significantly facilitated by Rodrigues’ rotation formula. This formula not only provides a straightforward and computationally efficient method for determining the effect of the rotation on a vector (i.e., rotating a vector) but also, more fundamentally, offers an explicit algorithm to compute the exponential map from the Lie algebra so (3) (the space of 3 × 3 skew-symmetric matrices, representing infinitesimal rotations) to the Special Orthogonal group SO (3) (the Lie group of finite 3D rotations) [29]. The corresponding rotation matrix $[eqn]$ is computed using Rodrigues’ rotation formula:

$$[eqn]$$

where

$$[eqn]$$

$[eqn]$ is the skew-symmetric cross-product matrix of $[eqn]$ . The axis-angle coordinates $[eqn]$ are a natural and compact rotation representation in terms of its geometric building blocks.

2.3. Kinematic State Model

A.line kinematic model

The kinematic state $[eqn]$ at time $[eqn]$ includes the center position $[eqn]$ and velocity components $[eqn]$ . The dynamic model of extended target is modeled by

$$[eqn]$$

where $[eqn]$ is the state transition matrix, the process noise $[eqn]$ .

For constant velocity model (CV) [27]:

$$[eqn]$$

where $[eqn]$ is sampling time, $[eqn]$ is the process noise variance for the centroid, $[eqn]$ is the Kronecker product, and $[eqn]$ denotes an identity matrix.

B.nonlinear kinematic model

While the target moves with a constant acceleration (CA) between time step, the acceleration state vector $[eqn]$ is augmented the kinematic state. Combining these into the 9-dimensional state vector $[eqn]$ , the state evolution model is given as follows [30],

$$[eqn]$$

where $[eqn]$ is the state transition matrix, the process noise $[eqn]$ .

$$[eqn]$$

The scalar $[eqn]$ represents the intensity of the underlying noise process and is a crucial tuning parameter for the CA model.

2.4. Measurement Model

Sensor data at time $[eqn]$ consists of a set of $[eqn]$ noisy 3D measurements $[eqn]$ originating from the target’s surface. Each measurement $[eqn]$ is related to the state and parameters via the non-line measurement function plus noise. A single measurement can be expressed as

$$[eqn]$$

where $[eqn]$ is the center of the target, $[eqn]$ is the radial function representation which describes the target’s shape in local coordinate. However, the measurement $[eqn]$ is obtained in global coordinates. Consequently, a global-to-local transformation is essential. This transformation yields the following local representation,

$$[eqn]$$

$[eqn]$ represents the radial function form in global coordinate, which modeled using DFS parameterized by the coefficient $[eqn]$ , and $[eqn]$ is zero-mean Gaussian measurement noise with covariance $[eqn]$ . $[eqn]$ is the angle pair of the measurement source on the target surface that originates $[eqn]$ .

$[eqn]$ is the rotation matrix that converts the global coordinates to local coordinates for the corresponding orientation parameter $[eqn]$ . This parameters $[eqn]$ utilize the axis-angle representation, expressed as $[eqn]$ , where $[eqn]$ is a unit vector ( $[eqn]$ ) indicating the axis of rotation orientation, and $[eqn]$ is the angle of rotation about the axis $[eqn]$ (typically in radial, following the right-hand rule). This parameter $[eqn]$ compactly stores both axis and angle information. Specifically, its magnitude $[eqn]$ give the rotation angle, while its direction $[eqn]$ give the axis for rotation orientation (for $[eqn]$ ). This axis-angle representation uses only 3 parameters (the components of $[eqn]$ ) and advantageously avoids the gimbal lock singularities associated with Euler angles.

The parameter $[eqn]$ is the set of coefficients for the DFS that defines the entire 3D shape of the target. The core assumption of the model is that the target is a rigid body. Therefore, this single, unified shape parameter $[eqn]$ is used to calculate the surface radius for any point $[eqn]$ on the object’s surface. When process the $[eqn]$ th measurement, and we use its specific projected angles $[eqn]$ as inputs to the universal shape function defined by $[eqn]$ to get the predicted radius at that exact spot.

Let $[eqn]$ and $[eqn]$ be expressed as,

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ .

As the measurement $[eqn]$ cannot be isolated as an explicit function of the state vector and noise, Equation (8) defines an implicit measurement model. It can be rewritten as:

$$[eqn]$$

The zero vector appearing in Equation (12) can be interpreted as representing a pseudo-measurement. This pseudo-measurement is effectively a nonlinear function of the state and the original measurement, assumed to be corrupted by additive Gaussian noise.

3. ECM Tracking Algorithm Based on DFS

This section details the DFS-ECM framework for joint motion state, shape, and orientation estimation. Its core, the Expectation Conditional Maximization (ECM) algorithm, uses a batch method involving iterative Kalman smoothing over a sliding data window. This incorporates information from previous and subsequent scans. For a batch window from time 1 to $[eqn]$ , the ECM framework defines these variables:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ denotes the latent variables within the sliding window; $[eqn]$ represents the set of measurements within the sliding window, including $[eqn]$ ; $[eqn]$ denotes the complete sequence of unknown parameters, where each $[eqn]$ of the sliding window at time step $[eqn]$ . ECM is an iterative optimization method whose objective is to maximize the log-likelihood function of the incomplete data (latent variable) $[eqn]$ with respect to the parameter $[eqn]$ .

The main terms used in this algorithm can be found in Table 1.

We define the complete-data log-likelihood function $[eqn]$ and its corresponding conditional expectation function $[eqn]$ of the $[eqn]$ th iteration as follows:

$$[eqn]$$ $$[eqn]$$

As the $[eqn]$ function involves multiple unknown parameters, a direct application of the traditional maximization step is not feasible. Instead, we utilize multiple conditional maximization steps to iteratively estimate the parameters.

3.1. E-Step

The primary purpose of the E-step is to estimate the state $[eqn]$ and its corresponding error covariance matrix $[eqn]$ , a process that involves an Unscented Kalman forward filter (UKF) [31,32] and backward Rauch–Tung–Striebel Smoother smoother [33,34]. The forward propagation is performed as,

Let $[eqn]$ be the dimension of the kinematic state $[eqn]$ , we generate $[eqn]$ sigma points $[eqn]$ based on the previous state estimate $[eqn]$ and error covariance $[eqn]$ (state initialized as $[eqn]$ , $[eqn]$ ):

$$[eqn]$$

The weights for mean $[eqn]$ and covariance $[eqn]$ are:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is a scaling parameter, $[eqn]$ controls the distance of the Sigma points from mean, $[eqn]$ is a secondary parameter $[eqn]$ is used to incorporate covariance information.

Prediction step of UKF:

(1) Propagate Sigma Points: Pass each state sigma point through the process model $[eqn]$ ,

$$[eqn]$$

(2) Compute the predicted state $[eqn]$ and predict error covariance $[eqn]$ :

$$[eqn]$$ $$[eqn]$$

Update step of UKF: the update proceeds sequentially for each measurement. Initialize the updated state and error covariance for this time step with predicted values,

$$[eqn]$$ $$[eqn]$$

the measurements $[eqn]$ are processed using the following steps,

(1) Predict the $[eqn]$ th measurement $[eqn]$ for each sigma point $[eqn]$ using the parameter $[eqn]$ ,

$$[eqn]$$

(2) Compute the $[eqn]$ th weighted average measurement $[eqn]$ of the predicted sigma point,

$$[eqn]$$

(3) Estimate the $[eqn]$ th innovation covariance $[eqn]$ and the cross-covariance $[eqn]$ ,

$$[eqn]$$ $$[eqn]$$

(4) Calculate the $[eqn]$ th Kalman Gain $[eqn]$ ,

$$[eqn]$$

(5) Update the $[eqn]$ th state $[eqn]$ and its error covariance $[eqn]$ ,

$$[eqn]$$ $$[eqn]$$

When processing all $[eqn]$ measurements, the UKF state $[eqn]$ estimation and its error covariance $[eqn]$ are:

$$[eqn]$$ $$[eqn]$$

After sequentially processing the data from time $[eqn]$ to $[eqn]$ with an Unscented Kalman Filter (UKF), we obtain the series of filter states $[eqn]$ and error covariances $[eqn]$ . Subsequently, the Rauch–Tung–Striebel smoothing algorithm is applied to compute the smoothed state $[eqn]$ and its error covariance $[eqn]$ as follows,

The backward pass begins at the final time step $[eqn]$ . (1) Initialize the smoothed state $[eqn]$ and error covariance $[eqn]$ for this time step with UKF estimation,

$$[eqn]$$ $$[eqn]$$

The algorithm then iterates backward in time, from $[eqn]$ down to $[eqn]$ . For each time step $[eqn]$ , the following computations with the measurements $[eqn]$ are performed:

(2) Compute the $[eqn]$ th Smoothed Gain $[eqn]$ ,

$$[eqn]$$

where $[eqn]$ is the $[eqn]$ th cross covariance and $[eqn]$ is the predict error covariance of forward UKF, both obtained from the forward UKF pass.

(3) Compute the $[eqn]$ th Smoothed state $[eqn]$ ,

$$[eqn]$$

(4) Compute the $[eqn]$ th Smoothed error Covariance $[eqn]$ ,

$$[eqn]$$

When processing all $[eqn]$ measurements, the smoothed state $[eqn]$ estimation and its error covariance $[eqn]$ are,

$$[eqn]$$ $$[eqn]$$

Considering the Markov property of state propagation and the independence of measurement sources, the log-likelihood function $[eqn]$ is written as,

$$[eqn]$$

So,

$$[eqn]$$

Under the assumption of Gaussian measurement noise, the measurement likelihood becomes,

$$[eqn]$$

Insert Equation (41) into (40), $[eqn]$ is,

$$[eqn]$$

where $[eqn]$ denotes the term same as the former one.

3.2. M-Step

After obtaining $[eqn]$ in the E-step of the algorithm, the task of the M-step is to find $[eqn]$ , thereby updating the target’s shape and axis-angle parameters from the previous iteration. Since $[eqn]$ is a non-linear function of the parameters, its analytical solution is generally difficult to obtain. Here, parameters are estimated by minimizing regularized cost function which equivalent to maxing $[eqn]$ estimation.

(1)Shape parameter optimization

The shape parameter $[eqn]$ are estimated by minimizing the cost function $[eqn]$ over a sliding window of $[eqn]$ time steps ( $[eqn]$ to $[eqn]$ ).

$$[eqn]$$

where, $[eqn]$ , $[eqn]$ are the estimates of state and orientation at time $[eqn]$ within the window. $[eqn]$ denotes the squared Mahalanobis distance. The first term is the average squared Mahalanobis distance over all valid measurements in the window. $[eqn]$ is the L2 regularization parameter applied to the coefficients (excluding $[eqn]$ ) to prevent overfitting and improving conditioning. $[eqn]$ is the smoothness penalty term.

Remark 1. Using the original cost function $[eqn]$  can lead to artificial bulges or distortions near the poles ( $[eqn]$ ). This occurs because the spherical coordinate system used by the DFS is degenerate at the poles, and measurement near these regions can exert disproportionate influence during optimization, leading to overfitting localized at the poles.

To mitigate these artifacts, we introduce a latitude-dependent weight factor $[eqn]$ that down-weights measurements near the poles. A common choice is

$$[eqn]$$

where $[eqn]$ (e.g., $[eqn]$ or $[eqn]$ ) is a tuning parameter controlling the strength of the suppression.

The modified cost function $[eqn]$ is,

$$[eqn]$$

Remark 2. By multiplying the squared Mahalanobis distance of each measurement by $[eqn]$ , we reduces the contribution of measurements where $[eqn]$  is close to $[eqn]$ . This is mathematically equivalent to assuming a higher effective measurement noise variance ( $[eqn]$ ) near the poles. It prevents the optimizer from overly fitting to potentially unreliable projections in these geometrically sensitive regions. Allowing the regularization terms to enforce a smoother, more plausible shape. The tuning parameter $[eqn]$  control how aggressively measurements near the poles are down-weighted.

Improved Smoothness Penalty: A physically meaningful smoothness penalty for the DFS on a sphere penalizes higher spatial frequencies more heavily. The indices $[eqn]$ and $[eqn]$ in the DFS correspond to the frequencies in the azimuth and elevation directions, respectively. Higher values of $[eqn]$ and $[eqn]$ represent finer details or oscillations on the surface. A suitable penalty terms $[eqn]$ is given by:

$$[eqn]$$

where $[eqn]$ is the smoothness regularization strength parameter. The term ( $[eqn]$ ) acts as a weighting factor that increases with the frequency indices $[eqn]$ and $[eqn]$ . A common is $[eqn]$ , making the penalty proportional to ( $[eqn]$ ). The sum of squares involves all four coefficient types for the given ( $[eqn]$ , $[eqn]$ ) pair.

(2)Orientation parameter optimization

Minimize $[eqn]$ at step $[eqn]$ , regularizing towards the previous estimate $[eqn]$ :

$$[eqn]$$

where $[eqn]$ are the current estimates of state and shape. $[eqn]$ is the estimated orientation parameter from the previous time step. $[eqn]$ is the regularization parameter controlling the penalty for deviating from the previous orientation estimation, promoting temporal smoothness.

(3)Optimization Constraints

Constraints can be applied during optimization:

• A.Shape parameter constraints:

$[eqn]$ (small positive value)

B.Rotation parameter constraints:

angle limit: $[eqn]$ .

Iteration of the EM algorithm ceases upon convergence or reaching a limit. Convergence is typically determined when the change in the likelihood function value between successive iterations is sufficiently small. Alternatively, termination occurs if the iteration count reaches a predefined upper bound.

The overall algorithm flow is summarized in Algorithm 1. Algorithm 1. DFS-ECM Initial parameters: Measurement set batch $[eqn]$ , State batch $[eqn]$ , orientation and extent parameters $[eqn]$ , Maximum Iterations $[eqn]$ begin   Setup $[eqn]$ 1       while not converged and $[eqn]$ do 2              expectation: 3                  for $[eqn]$ do 4                       for $[eqn]$ do 5                           Calculate $[eqn]$ and $[eqn]$ according to (28–29) 6                        end for 7                    end for 8                    for $[eqn]$ do 9                        Calculate smoothed $[eqn]$ and $[eqn]$  according to (37–38) 10                    end for

11               maximization: 12                   Calculate $[eqn]$ , $[eqn]$ according to (45,47) 13                $[eqn]$ ; 14         end while end begin

4. Simulation Results

Our proposed method, referred to as DFS-ECM in this section, are evaluated in two distinct scenarios. The simulations are run on an Intel Core™ i7-7700 processor at 3.6 GHz with 32 GB of RAM using MATLAB 2021a. Root Mean Square Error (RMSE) and the Intersection-Over-Union (IOU) are used to assess overall performance.

4.1. Evaluation Index

Root-Mean-Square Error (RMSE) is defined below:

$$[eqn]$$

where $[eqn]$ represents the RMSE of the parameter $[eqn]$ , $[eqn]$ represents the true and $[eqn]$ represents the estimated value.

The shape estimation performance is evaluated by Intersection-Over-Union (IOU). IOU is defined as follows:

$$[eqn]$$

where volumes $[eqn]$ is the true target shape, and $[eqn]$ represents the estimate. $[eqn]$ , $[eqn]$ corresponds to a perfect match, while $[eqn]$ indicates no intersection between the true and the estimated extent.

4.2. Scenario I

In this scenario, we test the DFS-ECM algorithm with ellipsoid and cubic target. Every target moves in a linear trajectory with the constant velocity (CV) model, and with a sensor located at the origin. The trajectory of target motion is presented in Figure 1. The initial kinematic state is set to [0 m, 0 m, 0 m, 10 m/s, 10 m/s, 0 m/s]. At each instant, 100 points measurements are originated from random sources which are sample from a uniform distribution defined over the target surface. The sampling time of all algorithms is set to $[eqn]$ . The orientation is represented by the axis angle parameter $[eqn]$ .

The process noise standard deviations are set to $[eqn]$ , For the sigma points, a scaling parameter $[eqn]$ , $[eqn]$ , $[eqn]$ . For the extent parameter, $[eqn]$ and $[eqn]$ , the covariance matrix is set to $[eqn]$ . The degree of double Fourier series is chosen as $[eqn]$ .

Figure 2 illustrates the typical results and tracking trajectories in one Monte Carlo run. Figure 2a,b depict the tracking scene for the ellipsoid, cube, respectively. For enhanced clarity regarding the reconstruction outcomes, Figure 3 illustrates the reconstruction results for 3D shapes obtained without incorporating the weight. A comparison between Figure 2 and Figure 3 reveals that optimizing shape reconstruction with the weight yields superior accuracy.

We evaluated centroid estimation for linearly moving ellipsoidal and cubic targets using 100 Monte Carlo trials (randomizing noise and sample points). Figure 4 shows the centroid RMSE for sliding window lengths $[eqn]$ , $[eqn]$ , and $[eqn]$ . For the ellipsoid (Figure 4a), $[eqn]$ yields the minimum RMSE; performance worsens for $[eqn]$ and $[eqn]$ , suggesting potential overfitting. For the cube (Figure 4b), $[eqn]$ is optimal. Generally, errors decrease over the simulation (indicating stabilization), and the relative performance of different window sizes remains consistent.

Figure 5 illustrates the RMSE values of the velocity estimates with the sliding window length $[eqn]$ , $[eqn]$ , and $[eqn]$ . This indicates that while there’s a large initial uncertainty in velocity, the estimation becomes highly accurate over time, and the choice of window size (within this range) has minimal impact on the long-term velocity estimation accuracy.

Figure 6 presents the Intersection over Union (IOU) values for the shape estimates, obtained using sliding window lengths of $[eqn]$ , $[eqn]$ , and $[eqn]$ . For both the ellipsoidal and cubic targets, the analysis indicates that the highest IOU score (representing the best shape estimate accuracy) is achieved with the window length $[eqn]$ .

As detailed in Table 2, the errors associated with the estimated orientataion parameter components remain consistently below 0.1 for both ellipsoidal and cubic target types, indicating a high degree of estimation accuracy.

For comparison, we use a standard random matrix-based extended target tracker [35], henceforce denoted as RM. Table 3 summarizes the comparative performance of the RM and DFS-ECM approaches. The results consistently show that the DFS-ECM method outperforms the RM model in terms of centroid estimation accuracy, velocity estimation accuracy, and shape estimation quality (as measured by IOU) across all tested scenarios.

As presented in Table 4, increasing the number of measurements positively impacts overall tracking performance. Specifically, a higher measurement count leads to a lower Centroid RMSE, indicating improved positioning accuracy. Concurrently, the shape estimation quality, measured by IOU, steadily improves for both the ellipsoid and the cube, confirming that more data yields more accurate reconstructions.

4.3. Scenario II

The performance of the DFS-ECM algorithm is further evaluated using a nonlinear scenario. In this experiment, a cubic target moves along a curved trajectory in three-dimensional space, as depicted in Figure 7. The target’s motion is governed by the constant acceleration (CA) model, which is defined in Equation (6). This simulation utilizes the same DFS model hyperparameters and EM parameters as those employed in Scenario Ⅰ. Within the CA model, the scalar parameter $[eqn]$ is set to 0.001. At each instant, 100 measurement points are sampled uniformly from the target’s surface. The orientation is represented by the axis angle parameter $[eqn]$ .

To better illustrate the reconstruction quality, Figure 8 presents magnified views of typical reconstruction examples for the cubic target. These examples demonstrate satisfactory reconstruction performance.

Figure 9 presents the centroid RMSE over time for the cubic target, comparing window lengths $[eqn]$ , $[eqn]$ , and $[eqn]$ . Optimal performance (lowest RMSE) is achieved with the shortest window, $[eqn]$ . Conversely, the larger window lengths ( $[eqn]$ and $[eqn]$ ) yield higher errors, potentially indicating underfitting in this nonlinear scenario. However, all three RMSE curves exhibit convergence towards stable values over time.

The velocity RMSE for the cubic target is presented in Figure 10. Notably, while the larger window sizes ( $[eqn]$ and $[eqn]$ ) show similar performance with lower peak errors and better long-term stability (lower RMSE), the smallest window ( $[eqn]$ ) suffers from a higher peak error (occurring later) and maintains a higher error level beyond approximately 18 s.

Figure 11 plots the volume Intersection over Union (IOU), representing shape estimation accuracy, over time for three different window sizes. Towards the end of the observation period (after 40 s), the performance stabilizes for all window sizes. They converge to similar and relatively high IOU values, approximately in the range of 0.86–0.88. While these results are strong, the accuracy is slightly below that achieved for the cube in the linear scenario. We attribute this minor performance difference to the inherent challenges of the nonlinear model. The nonlinearity likely introduces small errors into the centroid position estimate, which subsequently impacts the overall accuracy of the shape reconstruction.

The performance of orientation parameter estimation is detailed in Figure 12. The error dynamics for individual rotation components within a single Monte Carlo run are illustrated in Figure 12a–c. The overall accuracy, depicted by the average Root Mean Square Error (RMSE) across all Monte Carlo runs, is shown in Figure 12d. The low error values demonstrate the algorithm’s capability for accurate angle estimation.

Figure 13 illustrates the computational cost, specifically the average processing time per step required by the algorithm, as a function of different window lengths ( $[eqn]$ , $[eqn]$ , and $[eqn]$ ). A clear positive correlation is evident: longer window lengths necessitate incNreased computation time per step.

5. Conclusions

This paper proposed a method for joint 3D extended object tracking with unknown but fixed shape in the framework of Expectation Conditional Maximization (ECM). By introducing a shape representation based on a Double Fourier Series (DFS) radial function and utilizing the axis-angle representation for orientation, this method enables the direct estimation of target shape and orientation parameters, overcoming the reliance on prior shape evolution models. Leveraging the ECM algorithm, kinematic state inference is performed using an Unscented Kalman Filter (UKF) in the E-step, while shape and orientation parameters are estimated via minimization of regularized cost functions in the M-step. This ultimately achieves the joint, robust, and smooth estimation of kinematics, shape, and orientation. The methods provide a full expression for the target’s shape, this detailed information can be used in 3D applications such as navigation, identification and classification. Future work will focus on adaptive tuning and joint estimation schemes to further advance the capability of 3D target tracking.