## 1 Introduction  
Global Navigation Satellite Systems (GNSS) offer positioning services through a common method known as the two-step procedure (2 SP): satellite signals are independently processed to derive synchronization parameters, which are subsequently used for receiver multilateration to estimate the receiver's position, velocity, and time (PVT) state. While effective in open-sky environments, 2 SP-based receivers encounter accuracy and stability challenges in complex scenarios, including jamming, multipath, and channel fading. These limitations stem from the inability of the 2 SP to ensure that positioning parameters correspond to the same PVT state, thereby hindering its overall performance. Recently, the growing demand for high-precision positioning and navigation services in challenging environments has highlighted the limitations of traditional positioning methods. With advancements in data transmission and computational capabilities, there has been increasing interest in Direct Position Estimation (DPE), also referred to as the one-step procedure, which involves obtaining the Maximum Likelihood estimate (MLE) of the PVT by directly solving the problem in the navigation domain, where signals from all visible satellites are processed jointly [1; 2; 3]. By eliminating intermediate estimation steps, DPE can partially mitigate some of the shortcomings of 2 SP, such as degradation in position accuracy due to near-far resistance and severe channel fading conditions [4]. Despite its benefits, DPE presents a non-convex optimization problem, which makes it unsuitable for gradient-based algorithms and highly dependent on the optimization algorithm employed [5]. One such method, the Space Alternating Generalised Expectation-Maximization (SAGE) algorithm, decomposes a multi-dimensional search into a series of uni-dimensional searches [6]. However, as the clock dimension is not separable, an incorrect time estimate in the initial iteration can cause SAGE to converge to an incorrect result [7]. Another method uses the Accelerated Random Search (ARS) algorithm [8]. As a variant of the random search algorithm, its mechanism of adaptively adjusting the search radius enables dynamic coordination of the balance between global and local searches [9]. However, the performance of ARS highly depends on the setting of parameters, such as the initial search radius, contraction factor, and iterations. Currently, grid-based search algorithms remain the primary solution for DPE. Although assisted GPS (A-GPS) can provide a priori position-time information to effectively reduce the search space, the computational burden of performing signal correlation across a large number of candidate points remains a significant challenge. To overcome these limitations, various approaches have sought to improve search efficiency. One of the most effective methods is the multi-resolution grid (MRG) method, which reduces the number of computational points by adjusting the grid size at each stage [10; 11; 12]. This method also uses an averaged correlogram at the search initiation to accelerate the process. However, generating the correlogram requires calculating the correlation value for each candidate point, which is computationally intensive. Another approach involves using a pseudorange correlogram-based algorithm to significantly reduce the correlation process computational load [13], although it does not account for the impact of the satellite elevation on correlation values. While these methods can accelerate the correlation process, they do not alleviate one of the most significant challenges of DPE: the need to evaluate a large number of candidates. Swarm intelligence optimizers have been widely used to solve non-convex optimization problems in engineering fields, including positioning systems. An expanded and contracted Pigeon-Inspired Optimization (PIO) algorithm was developed for Collective Detection, a combined GNSS signal acquisition method [14].