The Earth–Moon libration point orbit (LPO) offers unique advantages in terms of location and dynamical characteristics, providing new opportunities for designing communication and navigation constellations. However, LPOs cannot be described by Keplerian elements and require significant computational resources for initial value searches, limiting their optimization potential. This paper proposes a two-step optimization algorithm based on a two-layer initial value library. The first layer represents orbit families, while the second layer contains orbits with varying amplitudes within each family. The first step of the algorithm quickly filters orbit families, and the second step selects the orbits that form the constellation, significantly reducing both the optimization scope and the number of parameters. We introduce a novel grid-based division of key cislunar regions and expand the constellation service area using a three-phase construction strategy. To evaluate the optimized constellation’s navigation capabilities, we test its performance with typical orbits, including Earth–Moon transfer orbit (EMTO), elliptical lunar orbits (ELO), and geosynchronous orbit (GEO). Experimental results show that a single near rectilinear halo orbit (NRHO) satellite provides single coverage of the Earth–Moon transfer critical region at sampling times, with an orbit determination (OD) accuracy of 475.7 m for the EMTO. Seven satellites, positioned on the L1, L2, L4, L5, and NRHO orbits, achieve quadruple coverage for both the Earth–Moon transfer and near-Moon regions, with OD accuracies of 25.7 and 17.3 m for the EMTO and ELO, respectively. In the third phase, adding two L3 LPO satellites forms a nine-satellite constellation, extending coverage across the entire cislunar space, achieving OD accuracies of 15.1, 13.4, and 2.2 m for the EMTO, ELO, and GEO, respectively. This study provides valuable insights for the future design and deployment of cislunar communication and navigation constellations.

Formation flying Low Earth Orbiters (LEOs) are important for implementing new and advanced concepts in Earth observation missions. Precise Baseline Determination (PBD) is a prerequisite for LEOs to complete specified mission targets. PBD is usually performed based on space-borne GNSS data, the relative corrections of phase center and code residual variations play crucial roles in achieving the best relative orbit accuracy. Herein, the influences of antenna Relative Phase Centre Variations (RPCVs) and Single-Difference (SD) Melbourne-Wübbena (MW) Combination Residuals Variations (SD MWVs) on PBD are studied. The methods were tested using flight data from Gravity Recovery And Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO). Results showed that the maximum values for RPCVs and SD MWVs were 14 mm and 0.32 cycles, respectively. Then, the RPCVs correction significantly enhanced the baseline accuracy; the K-Band Ranging (KBR) measurement consistency improved by 30.1% and 37.5% for GRACE and GRACE-FO, respectively. The application of SD MWVs further improved the accuracy and reliability of PBD results. For GRACE, the ambiguities fixing success rate increased from 85.1% to 97.9% and a baseline consistency of 0.57 mm was achieved for the KBR measurements. It was found that the correction of both RPCVs and SD MWVs reduced the carrier phase observation minus computation residuals from double-difference ionosphere-free combination. In addition, in-flight data processing demonstrated that RPCVs and SD MWVs estimations for the current period could be used for the previous and subsequent periods.