On-orbit calibration and performance validation of the Yunyao polarimetric radio occultation system

Polarimetric radio occultation (PRO) extends the capability of standard radio occultation (RO) by providing not only the conventional thermodynamic profiles but also information on clouds and precipitation. The differential phase (ΔΦ) is the cumulative phase shift between horizontal and vertical polarizations observed from PRO caused by aspherical hydrometeors along the propagation path, typically measured in millimeters. In early 2025, Yunyao Aerospace Technology Co., Ltd. successfully launched the first Chinese low Earth orbit satellite equipped with a PRO payload, generating over 500 measurements per day. Based on this mission, we established an end-to-end PRO data processing chain incorporating on-orbit calibration and tailored for operational applications. We analysed approximately 53 000 events collected between March and June 2025, in conjunction with the Integrated Multi-satellite Retrievals for Global Precipitation Measurement (GPM) precipitation product (IMERG). The results show that ΔΦ remains close to zero under non-precipitating conditions but exhibits distinct peaks at 3–5 km altitude when traversing precipitation layers, with amplitudes strongly correlated with path-averaged rainfall rates. Thresholds of 1, 2, and 5 mm h⁻¹ are proposed as indicators of precipitation sensitivity, detection confidence, and heavy-rain events, respectively, and a ΔΦ-to-rainfall intensity mapping table is derived to quantify this relationship. Yunyao PRO data preserve the thermodynamic retrieval quality of conventional RO while enabling effective precipitation detection, thereby providing important data support for the theoretical, technical and data research on the transition of meteorological observations from “temperature, humidity and pressure” observations to new types of observations such as precipitation.

1 Introduction

Precipitation events exert long-term and widespread impacts on human society, particularly in the context of intense rainfall. Sudden episodes of high-intensity precipitation and persistent extreme wet spells pose substantial threats to human health and critical infrastructure, while also triggering secondary hazards such as flash floods, landslides, and urban waterlogging. The mesoscale processes that govern the formation and evolution of precipitation are fundamentally controlled by atmospheric convection: upward motions transport near-surface water vapor into the upper troposphere, where decreasing temperatures induce condensation, freezing, and coalescence. These microphysical transformations ultimately manifest as various hydrometeors (e.g., liquid raindrops, supercooled droplets, ice crystals, graupel, and snow), undergoing phase changes and sedimentation that result in precipitation at the surface. Despite notable advances in convective theory and numerical modeling over the past decades, significant uncertainties remain regarding convective initiation, organizational patterns, microphysical processes, and land–atmosphere interactions. These uncertainties continue to constrain the ability of operational forecasting systems and climate models to accurately characterize the spatiotemporal distribution, intensity, and persistence of precipitation.

Figure 1 Illustration of the PRO concept.

Standard radio occultation (RO) retrieves scalar refractivity profiles for temperature, pressure, and water vapor, with the advantages of all-weather capability, high vertical resolution, and global coverage (Spilker et al., 1996; Jin et al., 2014; Teunissen and Montenbruck, 2017). Polarimetric radio occultation (PRO) enhances the standard RO by receiving the GNSS (Global Navigation Satellite System) signals in two orthogonal linear (horizontal and vertical (H/V)) polarizations. The induced cross polarization (denoted by ΔΦ) represents a cumulative phase delay between the horizontal and vertical measurements, arising from the presence of aspherical hydrometeors such as oblate-shaped raindrops and complex ice crystal shapes along each ray path (Turk et al., 2024), as shown in Fig. 1. Operational techniques capable of systematically resolving the vertical structure of convection remain lacking, limiting our understanding of the three-dimensional evolution of water vapor and hydrometeors during precipitation (Cardellach et al., 2010). Passive microwave radiometers (e.g., the Advanced Technology Microwave Sounder, ATMS) are sensitive to upper-level ice clouds and provide clues to cloud–ice processes, while imaging radiometers (e.g., the Global Precipitation Measurement Microwave Imager, GMI) retrieve column-integrated quantities such as total precipitable water and cloud liquid water. Yet, passive observations suffer from weak vertical resolution, strong dependence on prior information, and non-unique emission–scattering signals. Active precipitation radars (e.g., the Dual-frequency Precipitation Radar, DPR) infer microphysical profiles from backscattering, but cannot separate water vapor from hydrometeors and are strongly affected by absorption and attenuation (An et al., 2019). Furthermore, the limited spatiotemporal coverage of low Earth orbit sensors restricts continuous monitoring of rapidly evolving convection. The PRO extends GNSS RO for precipitation-sensitive sensing by measuring the ΔΦ (Cardellach et al., 2018). This polarization-resolved information enhances sensitivity to cloud–precipitation microphysics and convective organization.

With respect to instrumentation and methodology, Cardellach al. (2018) and Padullés et al. (2018) have systematically documented PRO antenna designs, receiver chains, calibration strategies, and retrieval algorithms, and clarified how signal-to-noise ratio, geometric configuration, and frequency plan govern ΔΦ detection sensitivity (Padullés et al., 2018). Theoretical analyses and early in-orbit experiments show that ΔΦ during precipitation markedly exceeds the expected noise floor, and that the differential-measurement design helps suppress certain systematic errors and ionospheric common-mode residuals. The Spanish PAZ satellite, launched on 22 February 2018, began radio-occultation and heavy-precipitation experiments on 10 May 2018, acquiring on the order of ∼ 300 GPS (Global Positioning System) profiles per day and confirming PRO's pronounced sensitivity to precipitation (Cardellach et al., 2019). At the mechanistic and application levels, Turk et al. (2022) revealed correlations between the vertical structure of ΔΦ profiles in convective systems and lower-tropospheric moisture profiles – implying diagnostic value for moisture–hydrometeor coupling in deep convection (Turk et al., 2019) – while Padullés et al. (2023) attributed strong ΔΦ peaks to electromagnetic waves traversing layers of oriented frozen hydrometeors (e.g., ice crystals and aggregates). Paz et al. (2024a) used Next Generation Weather Radars (NEXRAD) data to compare the specific differential phase with ΔΦ from PRO observations, and found good agreement on the PRO ΔΦ and collocated with ground-based NEXRAD radars. Retrieving Level-2 variables along the rays remains challenging due to geometric effects – long path integration and vertical superposition – which induce non-uniqueness: different configurations of particle location, amount, and type along the ray path can produce the same observable (integrated polarimetric phase shift). Lookup-table retrievals were proposed and validated with synthetic experiments (Cardellach et al., 2018). It is worth noting that PRO does not degrade traditional thermodynamics measurements (Paz et al., 2024b).

To fully exploit PRO's spatiotemporal sampling and microphysical sensitivity, Turk et al. (2024) proposed deploying small, tightly clustered PRO constellations to retrieve quasi-simultaneous water-vapor fields for diagnosing inter-model differences in the precipitation–moisture relationship, and they further examined synergistic retrieval frameworks that combine PRO with microwave imagers/sounders to quantify PRO's marginal contribution to cloud-ice characterization and forecast gain. Because PRO entails only limited hardware and processing-chain modifications to conventional RO payloads, it represents a cost-effective, technically continuous, and scalable evolution of GNSS radio occultation (Turk et al., 2024). Recent demonstrations include the ESA–Spire–Spanish research team integration of a novel triplet polarimetric antenna on the LEMUR-2 platform, validated on orbit in 2023 (Talpe et al., 2025), and PlanetiQ's launch of its first PRO satellite on 16 August 2024 (Kursinski et al., 2023; Padullés et al., 2025). The Chinese government is similarly promoting the adoption of GNSS polarimetric radio occultation payloads on the Fengyun series, thereby driving payload development and data-processing capability advances at the China Commercial Meteorological Satellite Corporation.

At present, Chinese commercial aerospace enterprises are actively laying out in this field. Yunyao Aerospace Technology Co., Ltd. is committed to integrating polarized occultation technology into its high-timeliness meteorological constellation. The point is to make up for the deficiency of traditional occultation observations in detecting the microphysical processes of water vapor and precipitation. By capturing ΔΦ caused by aspherical hydrometeors within the cloud and rain area, PRO can significantly enhance the constraints on water vapor condensation and phase transition paths during severe convective weather processes, thereby improving the accuracy of numerical models in short-term and imminent precipitation prediction.

Yunyao Aerospace Technology Co., Ltd., established in March 2019, is developing the “Yunyao Meteorological Constellation” comprising 90 high–temporal-resolution meteorological satellites (Xu et al., 2025). The constellation's core payloads include GNSS radio occultation and GNSS-Reflectometry instruments. The initiative is designed to deliver an end-to-end capability – spanning spaceborne observation, data processing, and service provision – to furnish industry users with high-cadence, reliable, and scalable meteorological space-data products. The schematic in Fig. 2 illustrates the planned network topology, highlighting on-orbit regional redundancy and enhanced temporal sampling that together support monitoring and nowcasting of rapidly evolving meso- and small-scale weather systems.

Figure 2 The 90-satellite Yunyao Meteorological Constellation.

Figure 3 The ground track of Yunyao polarimetric radio occultation satellite on 21 March 2025.

On 21 March 2025, Yunyao launched its first radio-occultation satellite carrying a GNSS Polarimetric Radio Occultation Instrument (GPROI) payload into a sun-synchronous, polar orbit (inclination 97.71°, nominal altitude ≈ 540 km). The satellite's sub-satellite ground track on the launch date is shown in Fig. 3. From initial insertion, the spacecraft has continuously acquired and downlinked occultation measurements, and the Yunyao Data Processing Center promptly established the “Yunyao PRO” dataset. Following geometric and polarimetric calibration and quality control, the center initiated a GNSS-PRO heavy-precipitation experiment and a suite of sensitivity analyses. To date, nearly three months of observations have been used to evaluate PRO sensitivity to precipitation, including: statistical characterization of ΔΦ profiles across distinct weather regimes; co-consistency tests between ΔΦ and independent ground/spaceborne references (reanalyses, ground-based radar reflectivity, and passive-microwave retrievals); and assessment of ΔΦ as an indicator of near-tropospheric moisture stratification and layers of oriented frozen hydrometeors. With the feasibility and cost-effectiveness assessments, Yunyao is considering expansion of the PRO constellation and optimization of orbital configurations to enhance the spatiotemporal coverage and robustness of near-real-time hydro-meteorological profiling.

This paper systematically presents the technical background of the Yunyao meteorological constellation and its PRO payload, the polarimetric radio occultation data-processing workflow and its key elements, and the sensitivity tests and experimental results assessing PRO's response to precipitation. Section 2 describes the Yunyao occultation and polarimetric-occultation data sources, the end-to-end processing chain and retrieval methods, and the reference products and co-registration strategies used for cross-validation. Section 3 details the heavy-precipitation experiment design and case-selection criteria, presents detection performance results and integrated analyses, and discusses sources of uncertainty, methodological limitations, and potential pathways for subsequent assimilation studies.

2 Data and Method

2.1 Data collection

The satellite mission was conducted as a technology demonstration payload. Building upon our existing, standardized occultation platform, we implemented a moderate upgrade to validate our capability in tracking and processing polarized radio occultation signals. Accordingly, the forward-facing antenna, signal channels, and processing methods retained the conventional occultation configuration. In contrast, the backward-facing antenna, channels, and data processing procedures were specifically designed to accommodate the characteristics of polarimetric occultation. Between 21 March and 21 June 2025, we acquired approximately 53 000 PRO events from Yunyao PRO payload. The global distribution of these events collocated with Integrated Multi-satellitE Retrievals for GPM (IMERG) is shown in Fig. 4, with overall spatial coverage that is relatively balanced due to the limb-grazing geometry of GNSS radio occultation. However, parts of Eastern Europe and the Middle East exhibit pronounced data sparsity and gaps. These localized deficiencies arise not only from complex radio-frequency environments (including elevated radio frequency interference, RFI) but also from the inherently lower signal-to-noise ratio of the polarimetric RO configuration relative to conventional RO, which increases PRO's vulnerability to such adverse radio conditions. Figure 5a shows the two-dimensional coverage of latitude and local time. Events exhibit a typical dual peak at dawn and dusk (approximately 04:00–06:00 and 17:00–19:00 UTC), with near symmetry across the northern and southern hemispheres. This suggests that orbital/geometric factors, rather than regional quality issues, are the primary controlling factor. Figure 5b shows that the daily counts for the three major constellations are generally stable, but with a natural fluctuation of approximately 10 %. The BDS (BeiDou Navigation Satellite System) averages approximately 330 daily counts, GPS approximately 230, and GLONASS (Global Navigation Satellite System) approximately 110. Data reception was suspended in early April and mid-June due to a combination of constellation maintenance and receive link status. Because channel resources were not initially allocated perfectly evenly across the constellations, the actual distribution is consistent with expectations, as shown in Fig. 5c.

Figure 4 Geographic distribution of Yunyao polarimetric radio occultation events from 21 March to 21 June 2025. Events are collocated with the IMERG precipitation product. Black dots denote non-precipitating events. Colored markers (blue to dark red) indicate events where IMERG precipitation > 0 mm h⁻¹; warmer colors and larger sizes indicate heavier rain.

Figure 5 (a) 1° × 30 min latitude – local time bin of PRO profiles. (b) Time series of PRO profile counts by GNSS constellation. (c) The proportion of PRO profiles of GPS, Galileo and GLONASS.

2.2 Standard RO data-processing workflow

Figure 6 summarizes the standard RO processing workflow at the Yunyao Data Processing Center. The workflow takes the “synthesized signal” – the complex baseband signal obtained by coherently combining the horizontal and vertical polarizations – as the scalar-channel input, sequentially deriving the excess phase, bending angle, and refractivity, and then inverting these to produce high-quality profiles of temperature, humidity, and pressure.

Figure 6 Standard GNSS radio occultation processing workflow.

2.2.1 Precise Orbit Determination

Precise orbit determination (POD) for low Earth orbit (LEO) satellites is a key driver of the accuracy of radio occultation retrievals – including excess phase, bending angle, and refractivity. Accordingly, the Yunyao Data Processing Center, in collaboration with Wuhan University, devised a POD scheme tailored to Yunyao LEO meteorological satellites (Shi et al., 2008), accounting for orbital altitude, refined force models, and spacecraft characteristics such as area-to-mass ratio, volume, and solar-array design (see Table 1). Over the experimental period (21 March–21 June 2025), the orbit of POD was evaluated via the overlapping-arc method in the radial (R), transverse (T), and normal (N) components, as summarized in Table 2. Because the Yunyao PRO receiver hosts both GPS and GLONASS tracking on the same receiver board, the POD accuracies for the two constellations are identical (Yue et al., 2025).

Table 1 Precision orbit determination strategy for Yunyao LEO satellite.

Table 2 Root Mean Square (RMS) difference for the overlap comparison in radial, transverse, normal and 3D directions for GPS GLONASS and BDS.

2.2.2 Atmospheric Excess Phase Processing

The observed GNSS carrier phase measurement on frequency can be expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{L}_f=\mathit{\rho }+c\left(\mathit{\delta }{t}_{\mathrm{r}}-\mathit{\delta }{t}^{\mathrm{s}}\right)+T+I_f+{\mathit{\lambda }}_f{N}_f+{\mathit{\epsilon }}_f\text{(2)}& {\displaystyle }\mathit{\rho }=\sqrt{(x^{\mathrm{s}}-x_{\mathrm{r}}{)}^{\mathrm{2}}+(y^{\mathrm{s}}-y_{\mathrm{r}}{)}^{\mathrm{2}}+(z^{\mathrm{s}}-z_{\mathrm{r}}{)}^{\mathrm{2}}}\end{array}$$

where ρ denotes the geometric range between the satellite and the receiver antennas; ($x^{\mathrm{s}},y^{\mathrm{s}},z^{\mathrm{s}}$) and ($x_{\mathrm{r}},y_{\mathrm{r}},z_{\mathrm{r}}$) are the Earth-centered inertial (ECI) coordinate of the GNSS satellite and LEO receiver in the True of Date (TOD) frame, respectively; c denotes the speed of light in vacuum; δt_r and δt^s are the receiver clock offset and satellite clock offset, respectively; T denotes the neutral-atmosphere tropospheric delay; I_f denotes the ionospheric delay at frequency f; λ_f and N_f are the carrier wavelength and the integer carrier-phase ambiguity at frequency f, respectively; ε_f contains other errors (multipath, phase center variations, etc.).

The atmospheric excess phase is then defined as

$$\begin{array}{}\text{(3)}& {\mathrm{\Phi }}_{\mathrm{atm}}=L_f-\mathit{\rho }-c\left(\mathit{\delta }{t}_{\mathrm{r}}-\mathit{\delta }{t}^{\mathrm{s}}\right)-{\mathit{\lambda }}_f{N}_f-{\mathit{\epsilon }}_f=T+I_f\end{array}$$

The single-difference phase at frequency f is used by taking the phase measured to the occulting GNSS satellite and subtracting the ionosphere-free phase in Eq. (4) measured from a high-elevation, non-occulting reference satellite, both tracked by the same receiver. This subtraction cancels the common receiver clock error and leaves only the differential contributions from geometric range and satellite clocks, together with the neutral-atmosphere and ionospheric delays, the integer phase ambiguity, and residual errors.

$$\begin{array}{}\text{(4)}& L_{\mathrm{IF}}^{\mathrm{ref}}={\displaystyle \frac{f_{\mathrm{1}}^{\mathrm{2}}}{f_{\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{2}}^{\mathrm{2}}}}{L}_{\mathrm{1}}^{\mathrm{ref}}-{\displaystyle \frac{f_{\mathrm{2}}^{\mathrm{2}}}{f_{\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{2}}^{\mathrm{2}}}}{L}_{\mathrm{2}}^{\mathrm{ref}}\end{array}$$

Thus, $\mathrm{\Delta }{L}_f=L_f-L_f^{\mathrm{ref}}$preserves the atmospheric/ionospheric information along the occulting path while removing the receiver clock term, providing the starting point for excess-phase processing. We obtain the final excess phase for each frequency as:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{\Phi }}_{\mathrm{atm},f}=& L_f-L_{\mathrm{IF}}^{\mathrm{ref}}\\ & -\left[\mathrm{\Delta }\mathit{\rho }+c\mathrm{\Delta }\mathit{\delta }{t}^{\mathrm{s}}+{\mathit{\lambda }}_f{N}_f-{\mathit{\lambda }}_{\mathrm{IF}}{N}_{\mathrm{IF}}^{\mathrm{ref}}+\mathrm{\Delta }{\mathit{\epsilon }}_f\right],\end{array}\end{array}$$

where L_f is the occulting-satellite carrier phase; $L_{\mathrm{IF}}^{\mathrm{ref}}$ is ionosphere-free combination from the reference satellite; Δρ and Δδt^s are geometric-range and satellite-clock differences.

In the lower troposphere, GNSS receivers often switch between closed-loop and open-loop tracking to maintain lock under severe multipath. Specifically, closed-loop tracking employs a feedback-controlled architecture where the receiver continuously adjusts its carrier and code replicas based on discriminator outputs from received signals. This phase-locked loop mechanism enables real-time compensation for Doppler shifts and dynamic trajectory variations, providing optimal tracking precision under stable signal conditions. In contrast, open-loop tracking operates without phase feedback, instead utilizing predicted ephemeris, atmospheric models, and pre-calibrated trajectory data to directly generate local signal replicas at anticipated signal phases. This forward-prediction approach eliminates loop latency to rapid signal dynamics, making it particularly robust in the lower troposphere where refraction and intense multipath frequently cause deep signal fades and phase discontinuities.

Processing coherently integrates In-phase and Quad-phase ($I/Q$) over defined intervals, extracts phase via $\arctan (Q/I)$ and amplitude via $\sqrt{I^{\mathrm{2}}+Q^{\mathrm{2}}}$, unwraps the phase to form a continuous time series, and then merges it with closed-loop data using dedicated transition logic (Chuang et al., 2013; Sokolovskiy et al., 2006). By performing independent or combined phase reconstruction on the open-loop observations $I/Q$ of vertical and horizontal polarization signals, we can obtain the vertical-/horizontal-polarization excess phases or the combined ones. The former is used to obtain polarization phase differences or polarization observation occultation profiles based on the traditional RO retrieval algorithm, and the latter is processed similarly to standard RO.

2.2.3 Atmospheric Profile Inversion

Our atmospheric-profile inversion adopts a dual-method strategy tailored to altitude. Above 20 km, bending angles are retrieved with geometric optics (GO) from Doppler-derived excess phase; below 20 km, full-spectrum inversion (FSI) transforms the complex RO signal into impact-parameter space (Gorbunov and Lauritsen, 2004), resolves multipath, and outputs Local Spectral Width (LSW) sequences used for quality control.

Ionospheric correction is altitude dependent. For regions above 20 km, the first-order ionosphere is removed with the ionosphere-free combination.

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_{\mathrm{corrected}}\left(a\right)={\mathit{\alpha }}_{\mathrm{L}\mathrm{1}}\left(a\right)+{\displaystyle \frac{f_{\mathrm{1}}^{\mathrm{2}}}{f_{\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{2}}^{\mathrm{2}}}}\left({\mathit{\alpha }}_{\mathrm{L}\mathrm{1}}\left(a\right)-{\mathit{\alpha }}_{\mathrm{L}\mathrm{2}}\left(a\right)\right),\end{array}$$

where α_L1(a) and α_L2(a) are bending angles derived from L1 and L2 frequencies. For regions below 20 km, we employ a fixed correction application strategy.

$$\begin{array}{}\text{(7)}& {\mathit{\alpha }}_{\mathrm{corrected}}\left(a\right)={\mathit{\alpha }}_{\mathrm{L}\mathrm{1}}\left(a\right)+\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{corrected}}\left(a=a_{\mathrm{20}\mathrm{km}}\right),\end{array}$$

where Δα_corrected(a=a_20 km) is the ionospheric correction at the 20 km impact parameter.

To obtain a physically consistent bending-angle profile, we merge GO- and FSI-retrieved bending angles from the geometric optics method (upper/middle atmosphere) based on Doppler shift analysis and the Full Spectrum Inversion (FSI) method (lower atmosphere) with a smooth window around 20 km. The weight formula is:

$$\begin{array}{}\text{(8)}& {\mathit{\alpha }}_{\mathrm{combined}}\left(a\right)=W\left(a\right){\mathit{\alpha }}_{\mathrm{GO}}\left(a\right)+(\mathrm{1}-W(a\left)\right){\mathit{\alpha }}_{\mathrm{FSI}}\left(a\right),\end{array}$$

where a is the impact parameter; W(a) is a smooth window increasing from 0 to 1 near ∼ 20 km; α_GO(a), α_FSI(a) are the GO- and FSI-retrieved bending angles; α_combined(a) is the fused bending.

LSW from FSI is used for threshold-based quality control to flag multipath. The merged profile is statistically optimized against a background using covariance weighting.

$$\begin{array}{}\text{(9)}& {\mathit{\alpha }}_{\mathrm{opt}}={\mathit{\alpha }}_{\mathrm{obs}}+C_{\mathrm{obs}}{\left(C_{\mathrm{obs}}+C_{\mathrm{bg}}\right)}^{-\mathrm{1}}({\mathit{\alpha }}_{\mathrm{bg}}-{\mathit{\alpha }}_{\mathrm{obs}}),\end{array}$$

where α_obs is the observation; α_bg is the background bending; C_obs, C_bg are the corresponding error covariances; α_opt is the optimized bending used for inversion.

The optimized bending angle is converted into refractive index by top-down integration of the inverse Abel transformation corrected for Earth ellipticity. The optimized bending angles are inverted to refractivity via the Abel transform with top-down integration and Earth-ellipticity correction:

$$\begin{array}{}\text{(10)}& N\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{x}{\overset{\propto }{\int }}{\displaystyle \frac{\mathit{\alpha }\left(a\right)}{\sqrt{a^{\mathrm{2}}-r^{\mathrm{2}}}}}\mathrm{d}a,x=nr,\end{array}$$

where N(x) is refractivity multiplied by radius; n is refractive index; r is geocentric radius. Dry temperature and pressure profiles are derived from refractivity using the Smith-Weintraub relation in combination with the hydrostatic balance and the ideal gas law, integrated downward from a prescribed top-of-atmosphere boundary condition. To obtain full profiles of temperature, pressure, and humidity, we apply a one-dimensional variational (1D-Var) retrieval, in which GNSS-RO refractivity serves as the observation constraint and the Yunyao Meteorological Numerical Forecast Model provides the background state. The minimum-variance (optimal-estimation) solution yields vertically resolved fields of temperature, pressure, and humidity.

2.3 PRO data-processing workflow

YunYao in-house GNSS-RO receiver onboard the satellite continuously acquires RO data at 100 Hz. In the YunYao configuration, only the occulting link is tracked, with closed-loop and open-loop channels running in parallel to allow flexible handover control in post-processing. Each RO event is tracked independently by two polarization ports (H and V), which operate independently yet are time-synchronized and output in parallel. The processing strategy and product convention for PRO follow PAZ (Padullés et al., 2024), and the overall workflow is illustrated in Fig. 7.

Figure 7 GNSS Polarimetric radio occultation processing workflow.

The processing chain comprises four stages: (i) $I/Q$ synchronization and geometric co-registration; (ii) generation of standard RO L1B products; (iii) construction and calibration of the polarimetric differential phase; (iv) cycle-slip removal, smoothing, and detrending to produce L2 products.

The H- and V-channel $I/Q$ streams are time- and frequency-aligned and coherently combined into a complex baseband as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }{I}_H+j{Q}_H\cdot \mathrm{conj}\left(I_V+j{Q}_V\right)=\left(I_{H-V}+j{Q}_{H-V}\right),\text{(12)}& {\displaystyle }\mathrm{\Delta }{\mathit{\phi }}_{H-V}=\mathrm{\angle }\left(I_{H-V}+j{Q}_{H-V}\right),\text{(13)}& {\displaystyle }{S}_{\mathrm{combined}}=\left(I_H+j{Q}_H\right)+\left(I_V+j{Q}_V\right)\cdot e^{-j\mathrm{\Delta }\mathit{\phi }},\end{array}$$

where I_H and I_V represents horizontal and vertical polarization in-phase component, respectively; Q_H and Q_V represents horizontal and vertical polarization orthogonal component, respectively; Δφ_H−V represents the relative phase difference between horizontal and vertical polarized signals; S_combined represents the combined complex signal. Signal-to-noise ratio (SNR) is computed and harmonized, and epochs and sampling rates are unified to ensure strict consistency for differencing and parallel processing. Standard RO processing is then applied independently to the H, V, and “synthesized” branches to produce the L1B atmospheric excess phases HatmPhs, VatmPhs, and atmPhs, as well as the polarimetric differential phase defined as the ΔΦ between the atmospheric excess phases from the two polarization ports:

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }\mathit{\varphi }={\mathit{\varphi }}_H-{\mathit{\varphi }}_V,\end{array}$$

This includes single differencing against a reference link to remove receiver clock and hardware delays, the dual-frequency ionosphere-free combination to suppress first-order ionospheric effects, and detection/repair of carrier cycle slips with lock-loss detection and phase reconstruction.

Figure 8 (a) L1 signal-to-noise ratio time series for the H- and V-polarization channels; (b) time series of the raw differential phase ΔΦ (blue) and the cycle-slip–corrected ΔΦ (red); (c) time series of the corrected ΔΦ (blue) and its 1 s moving-average smoothed series (1 s Smoothed). Gray dashed lines indicate the closed-loop–open-loop transition in each panel.

The PRO observable is defined as the ΔΦ between the atmospheric excess phases from the two polarization ports: In the absence of polarizing media along the path, ΔΦ should remain constant in time. Under ideal RHCP illumination the absolute H–V phase offset is $\mathit{\pi }/\mathrm{2}$; however, because H and V are tracked independently, only a random constant remains in the measurements. We estimate this constant over a stable stratospheric segment and remove it to align the zero level of ΔΦ. Even after conventional per-port cycle-slip repair, forming ΔΦ can expose residual jumps. We therefore adopt the polarimetric phase calibration scheme from Padullés (Padullés et al., 2020; Ao et al., 2009): half-cycle slips are corrected during CL using Eq. (15) and full-cycle slips during OL using Eq. (16):

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{\Delta }\mathrm{\Phi }\left(t\right)=\arctan \left(\tan \left(\mathrm{\Delta }\mathrm{\Phi }\right(t\left)\right)\right)\text{(16)}& {\displaystyle }\mathrm{\Delta }\mathrm{\Phi }\left(t\right)=\arctan \mathrm{2}\left(\sin \left(\mathrm{\Delta }\mathrm{\Phi }\right(t\left)\right)\cos \left(\mathrm{\Delta }\mathrm{\Phi }\right(t\left)\right)\right).\end{array}$$

This procedure removes remaining half- and full-cycle slips. Figure 8 shows an example for a GPS event (PRN-31) on 29 April 2025: Fig. 8a presents calibrated L1 SNR for H and V, decreasing as the tangent point descends and showing stronger jitter at the CL–OL transition (dashed lines); Fig. 8b compares the raw and de-slipped ΔΦ, where the corrected series becomes continuous and nearly constant aloft; Fig. 8c plots the corrected ΔΦ converted to path length together with a 1 s moving average.

Figure 9 Example of linear trend correction applied to the differential phase.

Afterwards, ΔΦ is smoothed with a 1 s moving window to suppress high-frequency noise and linearly detrended along the full profile to produce the polPhs file (Padullés et al., 2024). In parallel, the excess phase from the synthesized signal undergoes standard RO retrieval to derive dry and wet profiles, which are interpolated to a 0.1 km grid to form the resPrf file (Padullés et al., 2024) for collocation of ΔΦ with thermodynamic fields. For calibration, IMERG-identified rain-free events are used to derive the in-orbit antenna pattern and to remove any residual ionospheric imprint on ΔΦ (Padullés et al., 2020). Because the ΔΦ originates from hydrometeors, we set ΔΦ to zero at 30 km under water-free conditions to remove a profile-wide offset; all subsequent ΔΦ values are referenced to this level. The calibrated phase is again smoothed with a 1 s filter and linearly detrended; the resulting calibrated and smoothed ΔΦ is the primary PRO observable. After these steps, remaining ΔΦ variability can be attributed to differences in H- and V-component propagation induced by non-spherical, preferentially oriented hydrometeors along the path. Figure 9 illustrates the smoothing and detrending: in Fig. 9a the light-blue curve is ΔΦ after de-slipping, the blue curve is the 1 s smoothed ΔΦ, and the red dashed line is the linear trend; in Fig. 9b the detrended ΔΦ after smoothing is shown.

2.4 Global Precipitation Measurement (GPM)

In the calibration and validation of PRO, the correlation between polarimetric occultation data and precipitation products provides an independent assessment of rainfall influence. Rain-free events serve as references to calibrate systematic biases, while rainy events test ΔΦ sensitivity. IMERG, which fuses multi-satellite passive microwave and infrared observations, offers reliable, global precipitation estimates at ∼ 0.1° × 0.1° spatial and 30 min temporal resolution, with both near-real-time and final products (Katona et al., 2025; Watters et al., 2024). Its coverage and resolution make it suitable for independent co-location and sample selection.

Unlike prior approaches that fixed circular regions (Katona et al., 2025), we project the actual occultation ray paths – tilted elongated bands due to GNSS–LEO relative motion – onto the surface, considering only the troposphere below 6 km where precipitation primarily affects propagation. Polygonal footprints thus reflect both along-ray elongation and lateral extent, reducing geometric misclassification. Each occultation event is temporally matched to the nearest IMERG 30 min interval, and the surface-projected ray points are weighted to obtain an event-level mean precipitation. Rain-free events are used for calibration of the differential phase baseline, while rainy events assess PRO sensitivity and discriminative capability. Although only the lower-troposphere projection is used, the long-term dataset provides sufficient statistics for robust validation.

2.5 ERA5

ERA5 is the fifth-generation global atmospheric reanalysis dataset developed and maintained by the European Centre for Medium-Range Weather Forecasts (ECMWF), with a horizontal resolution of approximately 30 km and 137 vertical levels spanning the surface to 80 km. Covering the period from 1979 to the present, ERA5 undergoes rigorous quality control, with near-real-time products released within three months and reanalysis updates available with a 6 d lag. ERA5 data can be accessed and downloaded via the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/datasets, last access: 10 September 2025) (Hersbach et al., 2020). For this study, the dataset corresponding to the temporal span of the polarimetric occultation observations (21 March–22 June 2025) was retrieved and stored on the servers of Yunyao in netCDF format.

3 The correlation experiment of PRO profile and precipitation

3.1 Quality evaluation of PRO profile

H-polarized, V-polarized, and their combined circularly polarized observation data can all be used to retrieve meteorological elements, and the retrieved profiles are nearly identical. Compared to H- and V-polarized observations, the synthesized data exhibit a higher signal-to-noise ratio and a greater success rate in retrieval. Therefore, the accuracy of the profiles derived from the synthesized observations is representative of the performance of polarimetric occultation. Using ERA5 as a reference, a comparative evaluation was performed on the profiles of bending angle, refractivity, temperature and pressure retrieved from the polarimetric occultation.

The formula for calculating refractivity is given:

$$\begin{array}{}\text{(17)}& N=k_{\mathrm{1}}\times {\displaystyle \frac{P-e}{T}}+k_{\mathrm{2}}\times {\displaystyle \frac{e}{T}}+k_{\mathrm{3}}\times {\displaystyle \frac{e}{T^{\mathrm{2}}}},\end{array}$$

where N is the refractivity (in N-units), P is pressure (in hPa), T is the temperature (in K), and e is the water vapor pressure (in hPa). The coefficients are k₁=77.604 K hPa⁻¹, k₂=64.79 K hPa⁻¹, and k₃=377 600.0 K² hPa⁻¹.

The profile data corresponding to PRO profile is calculated by applying logarithmic cubic spline interpolation in the vertical dimension and bilinear interpolation in the horizontal plane to the ERA5 data. The relative deviations of the retrieved profiles (including bending angle, refractivity, temperature, and pressure) and the ERA5 values are then computed using Eq. (10), as shown below.

$$\begin{array}{}\text{(18)}& R={\displaystyle \frac{O-B}{B}}\times \mathrm{100}\mathit{\%},\end{array}$$

where R represents the relative deviation, O denotes the observed value, and B represents the reference value.

Figure 10 Comparison between PRO profile and ERA5 data.

Figure 10 shows the accuracy of PRO profile between 21 March and 21 June 2025. Statistical analyses were performed on the deviations and standard deviations of bending angle, refractivity, temperature, and pressure profiles at various altitudes for the BDS, GPS, and GLONASS systems (solid lines and dashed lines represent the deviations and standard deviations of each system respectively). Overall, all three systems exhibit minimal bias and stable errors from the mid-upper troposphere to the lower stratosphere (approximately 2–30 km). In contrast, the uncertainty increases significantly in the boundary layer below 1 km and the upper stratosphere above 30 km, which is generally consistent with the characteristic error profile of conventional radio occultation retrievals.

The accuracy of bending angle and refractivity profiles is shown in Fig. 10a and b, respectively. Both exhibit similar error variation characteristics. Taking the impact height as the y-axis in Fig. 10a, the bending angle accuracy shows its most pronounced negative bias at 2–4 km, where the mean reaches about −2.9 % and the standard deviation peaks near 9.6 %, consistent with multipath and super-refraction. From 4 to 30 km the bias is slightly positive and generally decreases with height; above 30 km both the mean and the spread increase modestly and the curves begin to diverge. In Fig. 10b, which uses altitude on the y-axis, refractivity exhibits a modest positive near-surface bias and large variability due to sharp boundary-layer moisture gradients and residual multipath; between 2 and 10 km, absolute mean differences are typically below 0.5 % and decrease with height; from 10 to 30 km, means remain within 0 %–0.2 % with a spread of about 1 %, followed by slight increases and divergence above 30 km.

Figure 10c shows the temperature profile accuracy. A slight negative bias appears between 2–3 km, transitions to a positive bias at 10–18 km, and stabilizes near zero from 20–30 km. This pattern reflects the propagation of multipath and humidity-induced errors from the lower troposphere into the retrieval. Above 30 km, the negative bias is dominated by higher-order ionospheric residuals. Figure 10d presents the pressure profile accuracy. A small negative bias persists throughout the column, reaching about −0.07 % near the tropopause. The bias increases slightly between 30–40 km due to error amplification from sequential integration, but remains small (2–30 km).

In summary, the thermodynamic elements retrieved from the synthesized PRO data demonstrate robust performance between 2–30 km, with both bias and dispersion remaining small.

3.2 Rainfall detection capability of PRO profile

This section presents the differential phase observed by the PRO payload in relation to precipitation. The detection capabilities are similar across different GNSS constellations.

Figure 11 ΔΦ not induced by precipitation. (a) The IMERG precipitation product collocated with the Yunyao PRO profile on the southern coastline of Africa. (b) The brightness temperature data collocated with the Yunyao PRO profile. (c) The corresponding ΔΦ values for each PRO tangent point. Each line represents the ΔΦ along a single occultation ray path, with color denoting the altitude of its tangent point. The gray dashed line represents the freezing level.

Figure 12 ΔΦ induced by precipitation. (a) The IMERG precipitation product collocated with the Yunyao PRO profile over the Pacific. (b) The brightness temperature data collocated with the Yunyao PRO profile. (c) The corresponding ΔΦ values for each PRO tangent point. The presence of precipitation induces an increase in ΔΦ values, which is particularly substantial at altitudes below 7 km. The gray dashed line represents the freezing level.

To evaluate the precipitation detection capability of the Yunyao PRO profile, we matched the ΔΦ along the ray path with the IMERG precipitation rate product from GPM. The geometric relationship was constrained using the tangent point trajectory in the longitude-latitude plane. Figures 11 and 12 present two representative cases, illustrating scenarios “ΔΦ without precipitation” and “ΔΦ with precipitation”, respectively, to verify the correlation between precipitation and ΔΦ. Meanwhile, we used 4 km pixel-resolution infrared brightness temperature data to match the PRO, so as to be able to observe the difference between ice and liquid water on the ray path.

In the non-precipitation case, Fig. 11a shows that the ray path remains entirely within the IMERG zero-precipitation region. Figure 11b shows that the infrared brightness temperature along the ray path. The corresponding vertical profile of ΔΦ in Fig. 11c exhibits minor fluctuations around zero across all altitude layers, with mean ΔΦ values of −0.05 and −0.09 mm for the 0–6 and 0–12 km ranges, respectively. These magnitudes are consistent with the slight negative bias observed in the profile and can be attributed to noise and residual correction errors rather than hydrometeor signals. Under the same processing chain and geometric conditions, the sub-millimeter mean ΔΦ can be regarded as the noise background of the instrument and retrieval system. There is no significant systematic deviation, indicating that the calibration is accurate.

In contrast, as shown in Fig. 12a, the ray path traverses multiple intense precipitation regions starting from the near-surface layer. Figure 12b shows that the infrared brightness temperature along the ray path. Figure 12c reveals a significant increase in ΔΦ within the troposphere, with rapid attenuation observed near the ground. A peak value of 15–16 mm occurs around 3–4 km, decreasing to approximately 6–8 mm at 6 km, further declining to 1–2 mm between 8–10 km, and approaching zero at 14–18 km. The mean ΔΦ values over the 0–6 and 0–12 km layers are 12.25 and 6.96 mm, respectively, with the matched IMERG precipitation rates being 2.15 and 1.13 mm h⁻¹. This vertical structure, which is characterized by strong signals in the lower atmosphere and rapid attenuation with height, is consistent with the distribution of raindrop-dominated liquid water and its geometric anisotropy: large, oblate, and oriented raindrops in the near-surface layer cause significant ΔΦ accumulation at L-band, while the liquid water content decreases rapidly with altitude and the contribution from ice-phase particles remains weak, leading to ΔΦ approaching zero. The correlation between precipitation and ΔΦ is also confirmed: tangent points at lower altitudes are closer to the main precipitation core, corresponding to higher ΔΦ peaks; as the tangent points move away from the precipitation core with increasing height, ΔΦ decreases accordingly.

It should be noted that ΔΦ is an integrated quantity along the oblique path, influenced by factors such as the horizontal drift of tangent points and the inhomogeneity of the precipitation volume. The spatiotemporal representativeness error of IMERG may also affect the matched precipitation rate. The consistency in geometry, vertical structure, and external precipitation fields between these two cases clearly demonstrates that the Yunyao PRO profile exhibits a distinct and reproducible response characteristic to liquid precipitation.

3.3 Precipitation sensitivity experiment of PRO profile

Based on the registration statistics between ΔΦ of Yunyao PRO profile and precipitation of GPM IMERG, only the segment below 6 km altitude was used. The ΔΦ profiles smoothed with a 1 s window were interpolated onto a 100 m equidistant vertical grid from 0 to 40 km. The profiles were then grouped according to precipitation intensity: no precipitation (0 mm h⁻¹), light precipitation (0–0.1, 0.1–0.2, 0.2–0.3, 0.3–0.5, 0.5–0.7 mm h⁻¹ ), and moderate-to-heavy precipitation (0.7–1, 1–2, 2–3, 3–5, 5–7, ≥ 7 mm h⁻¹). For each group, the mean and standard deviation of ΔΦ were computed at each altitude layer. Figure 13a–c show the sample size, mean ΔΦ, and the standard deviation profile under non-precipitating conditions for each precipitation interval. This method follows previous studies (Teunissen and Montenbruck, 2017; Cardellach et al., 2018; Gorbunov et al., 2011), enabling direct comparison with published results.

Figure 13 Statistical analysis of results over 0–40 km. (a) Sample size for each precipitation group; (b) mean ΔΦ per precipitation group, with the black curve indicating the bias under non-precipitating conditions; (c) standard deviation of the ΔΦ for the non-precipitation (0 mm h⁻¹) group.

Figure 13a shows the distribution of sample sizes across altitude and precipitation intensity groups. The non-precipitation group contains approximately 13 000 samples per layer between 6–15 km, while several light precipitation groups (0–0.5 mm h⁻¹) maintain sample sizes on the order of several thousand. In contrast, heavier precipitation groups (2–7 mm h⁻¹) exhibit an order-of-magnitude decrease in sample count. This sample size structure reflects both the spatiotemporal sparsity of global precipitation and ensures the robustness of mean and variance estimates for each group.

Figure 13b presents the mean ΔΦ profiles for different precipitation intensity intervals. The non-precipitation group shows ΔΦ values close to 0 mm across all altitude layers, as expected, and is consistent with the results reported by Cardellach et al. (2014). As precipitation intensifies, the magnitude of ΔΦ increases and exhibits significant altitude dependence. At a representative altitude of 3 km, the mean ΔΦ values for each group – from light to heavier precipitation – are 0.07, 0.48, 0.89, 1.23, 1.54, 2.07, 2.81, 4.27, 6.69, 9.54, and 9.03 mm, respectively. These results indicate an approximately linear response of ΔΦ to path-averaged precipitation rate, with higher sensitivity in the lower troposphere. It should be noted that the ΔΦ values of the heavier precipitation groups peak around 3–5 km and then decrease rapidly with height, falling below 1–2 mm above 10–12 km. This vertical structure which characterized by strong signals at lower levels and rapid attenuation aloft reflects both the optical anisotropy dominated by liquid raindrops and geometric factors: in the mid-troposphere, a few kilometers in altitude, certain cloud layers such as altocumulus and altostratus are commonly present. Radio occultation rays that intersect these layers may traverse an extended optical path through moist/cloudy regions, potentially influencing attenuation or refractivity retrievals depending on local water/ice content and particle phase distribution. Such behavior of the ΔΦ reveals the coupling between strong liquid-phase contributions in the lower atmosphere and weak ice-phase contributions at higher altitudes, and is consistent with the geometric mechanism of “ray path length maximization.”

Figure 13c shows the altitude-dependent standard deviation of ΔΦ under the 0 mm h⁻¹ condition (non-precipitation regions), which can be regarded as the “noise background” of the instrument and retrieval system. This provides a statistical basis for batch quality control and threshold setting. The standard deviation decreases with increasing altitude: it is approximately 1.5 mm near the surface£¬the result generally consistent with Padullés et al. (2020), though the near-surface trend differs, likely due to factors such as the Yunyao PRO receiver performance or the open-loop tracking algorithm. The value decreases to around 1.40 mm at 2 km, 0.70 mm at 5 km, and remains below 0.30 mm above 10 km, stabilizing between 0.12–0.16 mm in the 20–30 km range. This magnitude is consistent with the results reported by Padullés et al. (2020), further validating the stability of calibration and correction.

It should be noted that in this study, the precipitation-based grouping using the IMERG product was performed based on the path-averaged precipitation rate below 6 km altitude. Therefore, the correspondence between ΔΦ and precipitation reflects a coupling of path-integrated and area-averaged relationships. Inhomogeneity in the precipitation structure and beam-filling effects may introduce a certain degree of dispersion, particularly in the lower atmosphere. Nevertheless, the consistency in geometry, vertical structure, and magnitude of the grouped profiles along with their good agreement with previous literature collectively demonstrates that the ΔΦ of Yunyao PRO profile exhibits a distinct, reproducible, and dynamically responsive detection capability for liquid precipitation.

To further quantify the detection performance of ΔΦ in 0–6 km, we computed

$$\begin{array}{}\text{(19)}& P(\mathrm{\Delta }\mathrm{\Phi }>\mathit{\tau }|R),\end{array}$$

where R is the path-averaged precipitation rate from IMERG, and τ denotes the phase difference threshold (0.5, 1, 2, and 5 mm). As shown in Table 3, the threshold exceedance rate increases monotonically with R. Under non-precipitating conditions, only 7.85 % and 2.25 % of the cases exceed 1 and 2 mm, respectively, while exceedance of 5 mm is nearly negligible (0.09 %). However, approximately 17 % of cases still show ΔΦ exceeding 0.5 mm, indicating that additional factors influencing ΔΦ correction need to be considered. Under light precipitation, the exceedance ratio rises. When R ≥ 1 mm h⁻¹, the rates exceeding 1 and 2 mm reach 75.40 % and 52.44 %, respectively. For R = 2–5 mm h⁻¹, these values increase further to 89.37 % and 76.29 %. When R > 5 mm h⁻¹, the corresponding rates are 97.91 % and 95.20 %, with 77.61 % of cases exceeding 5 mm. Therefore, a threshold of 1 mm can be used as a highly sensitive indicator for “precipitation presence”, 2 mm as a conservative and highly specific threshold for “confident detection”, and 5 mm as an indicator for heavy precipitation/deep convection.

Table 3 Statistics on some ΔΦ relative thresholds for each elevation layer.

Precipitation is assessed using inverse conditional probability:

$$\begin{array}{}\text{(20)}& P(R>\mathit{\rho }|\mathrm{\Delta }\mathrm{\Phi }).\end{array}$$

As shown in Table 4, when ΔΦ < 0.1 mm, the probability of R > 0.1 mm h⁻¹ is only 1.85 %, and R > 1 mm h⁻¹ is nearly zero. This indicates that small phase differences are highly indicative of non-precipitation scenarios. For ΔΦ in the range of 0.1–1 mm, the probability of R falling within 0.1–1 mm h⁻¹ is 61.71 % (with 16.34 % for the sub-interval), demonstrating moderate predictive capability. When ΔΦ = 1–2 mm, the probability for R = 1–2 mm h⁻¹ increases to 72.17 % (41.36 %). For ΔΦ = 2–5 mm, the probabilities of R > 1 mm h⁻¹ and R > 2 mm h⁻¹ reach 87.90 % and 74.53 %, respectively. When ΔΦ > 5 mm, the probabilities for R > 1 mm h⁻¹, R > 2 mm h⁻¹, and R > 5 mm h⁻¹ are 97.18 %, 94.34 %, and 75.63 %, respectively, indicating a high likelihood of significant heavy precipitation.

Table 4 Statistics on some rain rates relative thresholds for each elevation layer.

Figure 14 The IMERG precipitation intensity-phase difference statistics of the 0–6 and 0–10 km datasets.

Figure 14 shows boxplots of precipitation intensity versus ΔΦ, with the path-averaged precipitation rate R along the occultation ray path as the independent variable. The mean and variability of the layer-averaged ΔΦ over 0–6 and 0–10 km are presented as functions of R. Both exhibit pronounced monotonic and convex growth curves. Taking the 0–6 km precipitation intensity–phase difference relationship as an example, ΔΦ fluctuates slightly around zero under near-zero precipitation conditions. When R = 0.1–0.3 mm h⁻¹, ΔΦ increases to approximately 0.3–0.8 mm. Beyond R ≥ 1 mm h⁻¹, the curve rises sharply, with R ≈ 3, 5, 7, and 10 mm h⁻¹ corresponding to ΔΦ values of approximately 3, 5, 7, and 9 mm, respectively. The expansion of error bars with increasing R reflects stronger spatial inhomogeneity in heavy precipitation and the variability introduced by geometric path weighting.

The difference between the 0–6 and 0–10 km statistical results reflects the influence of altitude on the correlation. The ΔΦ values in the 0–10 km are systematically lower than those in the 0–6 km panel at the same R, indicating that the inclusion of higher altitudes (> 6 km) dilutes the lower-layer polarimetric signal dominated by liquid precipitation, which is a direct consequence of path-integration. This result is consistent with the physical picture established in previous literature, where liquid water dominates and lower-layer contributions prevail: the shape/orientation anisotropy of raindrops in the near-surface layer, combined with longer cloud-path lengths, amplifies the ΔΦ accumulation. In contrast, the decrease in liquid water content and the weaker contribution from ice-phase particles at higher altitudes result in a continued yet diminishing marginal increase in the curve.

Table 5 The relationship of Polarization phase difference and rain rate.

Meanwhile, we also performed binning of ΔΦ based on the 0–6 km layer average (e.g., 0–0.1, 0.5–0.7, …, ≥ 18 mm). Within each ΔΦ bin, only samples whose R values fell within the range of mean ±1σ were retained for quality control, thereby obtaining a conditional mapping of $\stackrel{\mathrm{‾}}{R}\left(\mathrm{\Delta }\mathrm{\Phi }\right)\pm \mathit{\sigma }\left(R\right)$. This mapping provides the most probable precipitation rate and its uncertainty when a specific ΔΦ is observed, as shown in Table 5.

4 Conclusions

This study establishes the first end-to-end framework for PRO observations from the Yunyao mission. Quality assessment confirms that the synthesized H–V signals achieve low noise and bias in bending angle and refractivity within 2–30 km, with stable temperature and pressure retrievals; boundary-layer and upper-stratospheric error patterns are consistent with previous findings, demonstrating that the conventional RO capability meets operational standards and provides a reliable baseline for interpreting ΔΦ. For precipitation sensitivity, we introduce a collocation method based on realistic ray geometry, in which only the ray segments below 6 km are matched with IMERG using temporal and spatial weighting, effectively reducing the geometric mismatches of fixed circular windows. Case studies further reveal a distinct ΔΦ vertical structure: near-zero under rain-free conditions, but peaking at several to over ten millimeters at 3–5 km within precipitation, followed by rapid attenuation to 1–2 mm above 10–12 km. These results, together with statistical analyses, support the use of ΔΦ as a robust precipitation-sensitive observable and motivate further development of forward operators, assimilation schemes, and constellation-based applications.

Beyond precipitation detection, PRO observations have broader potential in atmospheric science. The sensitivity of ΔΦ to aspherical hydrometeors enables its use in discriminating precipitation types and in identifying mixed-phase and ice-dominated cloud regions. Combined with conventional RO profiles, PRO can constrain not only thermodynamic structures but also microphysical processes aloft, providing a pathway to improve cloud parameterizations in weather and climate models. Furthermore, the high spatiotemporal sampling of PRO constellations supports the analysis of moist processes in data assimilation systems, potentially enhancing the accuracy of short-term precipitation forecasts and the representation of latent heating in tropical cyclones and mesoscale convective systems. As a cost-effective extension of existing RO infrastructure, PRO is poised to bridge gaps between thermodynamic sounding and precipitation observation, advancing the integrated profiling of the moist atmosphere.