Dual-Prior-Constrained Temporal Inversion and FNO Surrogate-Model-Driven Coordinated Optimization of Drilling Engineering Parameters

1. Introduction

With the continuous expansion of oil and gas exploration toward deep, ultra-deep, and structurally complex formations, modern drilling operations are confronted with severe challenges including high-temperature and high-pressure downhole environments, narrow equivalent-circulating-density (ECD) safety windows, and frequent wellbore instability risks [1]. Precise real-time hydraulic monitoring and dynamic parameter optimization are essential to guarantee drilling safety and operational efficiency. Traditional drilling hydraulic analysis and parameter adjustment strategies rely on static pre-drilling design parameters and empirical mechanism models, which fail to capture the dynamic temporal variation of drilling-fluid rheology and annular cuttings transport characteristics during actual drilling processes. Such model errors easily induce downhole complex accidents such as lost circulation, well kick, and borehole collapse [2].

Accurate inversion of downhole hydraulic and cuttings distribution parameters serves as the fundamental premise of real-time while-drilling hydraulic monitoring [3]. Most conventional approaches to managed pressure drilling rely on regulating down-hole pressure to predetermined set points, rather than explicitly estimating the magnitude and location of in-/out-flux events through model-based observers [4]. Conventional drilling hydraulics models are typically calibrated offline and become increasingly inaccurate as downhole conditions evolve during drilling, leading to erroneous downhole predictions and elevated operational risk [5]. Although transient models integrated into real-time digital twins have been proposed to monitor downhole cuttings transport, steady-state models still rely on static parameters that fail to capture the dynamic evolution of cuttings distribution [6]. Furthermore, existing inversion frameworks underutilize multi-source prior information: pre-drilling mechanism simulation baselines and regional offset-well statistical data, which contain abundant formation and fluid characteristic information, are rarely integrated into inversion constraints, resulting in severe fluctuations and abnormal jumps in inversion results under high measurement noise.

The low computational efficiency of traditional hydraulic forward models constitutes another critical bottleneck restricting real-time drilling intelligent decision-making [7]. Conventional computational fluid dynamics (CFD) methods can achieve high-precision simulation of annular two-phase flow, yet a single working condition requires hours of iterative calculation, which is incapable of supporting the massive repeated forward evaluations required by online parameter optimization [7]. In recent years, physics-informed neural networks (PINNs) have been widely applied in fluid-mechanic prediction by embedding physical governing equations into loss functions, enabling physics-driven training that reduces reliance on large labeled datasets [9]. Nevertheless, PINNs are inherently condition-specific, remaining less accurate than traditional solvers for forward problems [10]. For complex eccentric rotating annular two-phase flow in drilling engineering, PINNs suffer from difficult training convergence, prominent local prediction errors, and poor adaptability to time-varying drilling conditions, which greatly restrict their practical application in real-time hydraulic prediction.

As an emerging neural-operator architecture, the Fourier neural operator (FNO) implements parametric modeling of integral kernels in Fourier space and directly learns mapping relationships between function spaces, exhibiting advantages of zero-shot super-resolution capability, high computational efficiency, and the ability to learn an entire family of PDEs [11]. Li et al. originally proposed the FNO framework in their arXiv preprint ‘Fourier Neural Operator for Parametric Partial Differential Equations’ [11]. Kovachki et al. further established the complete theoretical system of neural operators in a full-length journal article [12]. FNO has demonstrated superior performance in solving parametric PDEs, including complex turbulent flow prediction problems. Li et al. (2021) showed that FNO achieves higher accuracy than previous learning-based solvers and is up to three orders of magnitude faster than traditional PDE solvers [11]. Nevertheless, the application of FNO in drilling-fluid hydraulic surrogate modeling and drilling-engineering parameter optimization remains underexplored, and its superiority in dynamic drilling prediction and intelligent optimization requires systematic verification [14].

One of the present authors (Feng Ni) previously developed a closed-loop machine learning framework integrating real-time lithology identification with drilling parameter optimization, in which optimal WOB and RPM setpoints are derived by inverting the classical Teale MSE model upon formation-change detection [15]. In terms of drilling parameter optimization, existing studies have explored the optimization of WOB and RPM based on mechanical specific energy (MSE) minimization to improve rock-breaking efficiency [16]. In practical drilling processes, WOB, RPM, and flow rate present strong hydraulic–mechanical coupling and mutually restrictive relationships: WOB and RPM dominate the rate of penetration (ROP) and the cuttings generation rate, whereas the flow rate determines annular pressure loss, ECD variation amplitude, and cuttings carrying efficiency. Isolated single-parameter optimization inevitably deteriorates other drilling performance indicators and fails to achieve a global optimal operational status. Although multi-objective optimization algorithms have been introduced into drilling parameter optimization [18], most existing strategies adopt offline static optimization modes, which cannot dynamically update with real-time variations of downhole formation, fluid, and hydraulic parameters, leading to poor dynamic adaptability and low engineering practicability [19].

To address the aforementioned defects in conventional inversion and optimization methods, this paper proposes a novel three-layer stacked intelligent decision framework integrating dual-prior-constrained temporal inversion, FNO hydraulic surrogate modeling, and multi-objective coordinated optimization. The main innovative contributions are summarized as follows:

1.

A dual-prior-constrained temporal dynamic inversion method is proposed. By integrating pre-drilling mechanism baseline prior and offset-well statistical prior information, combined with standpipe pressure–rotary torque dual observational coupling constraints and sliding-window temporal smooth regularization, a four-term composite loss function is constructed to realize robust real-time inversion of multiple key hydraulic parameters, significantly improving the anti-noise performance and temporal stability of the inversion results.

2.

An FNO-based eccentric rotating-annular two-phase hydraulic surrogate model is established to replace conventional PINNs. The proposed model realizes millisecond-level full-field hydraulic prediction with high precision and strong generalization, breaking the efficiency and generalization bottlenecks of traditional numerical simulation and intelligent surrogate models in drilling hydraulic prediction.

3.

A hydraulic–mechanical coupled three-parameter coordinated real-time optimization framework is constructed to synchronously optimize WOB, RPM, and flow rate. Combined with online Bayesian multi-objective optimization and sliding-window rolling update strategies, dynamic adaptive parameter optimization under time-varying drilling conditions is realized. Probabilistic safety constraints based on inversion-parameter confidence intervals are established to avoid optimization risks induced by parameter uncertainty, achieving a global optimal balance of drilling safety, efficiency, and energy consumption.

The remainder of this paper is organized as follows. Section 2 first establishes the eccentric rotating-annular two-phase transient hydraulic forward model, and then elaborates the dual-prior-constrained temporal inversion module, the FNO hydraulic surrogate model, and the multi-objective coordinated optimization framework. Section 3 validates the proposed method through numerical simulations, ablation experiments, and comparative analyses, and discusses the underlying mechanisms, limitations, and future directions. Section 4 summarizes the core conclusions.

1.1. Drilling Fluid Hydraulic Parameter Inversion

Drilling-fluid hydraulic parameter inversion is a typical ill-posed inverse problem in drilling engineering, aiming to infer unknown downhole fluid rheological properties and annular-flow parameters through measurable surface or downhole observational data. Traditional hydraulic inversion methods primarily rely on steady-state mechanism models and least-squares fitting algorithms, which invert rheological parameters merely based on single standpipe-pressure data. Single-source observational information cannot satisfy the identification requirements of multi-parameter coupled hydraulic systems, resulting in low inversion accuracy and poor stability.

To improve inversion performance, scholars have introduced advanced filtering algorithms and multi-source data constraints into hydraulic-inversion research. He et al. (2022) developed an inversion-based multi-phase-flow interpretation model to realize real-time dynamic identification of downhole flow parameters during managed-pressure drilling [2]. Kaasa et al. (2012) constructed a simplified hydraulic model for high-timeliness downhole-pressure estimation, which does not support multi-parameter joint inversion [3]. Hauge et al. (2012) proposed a model-based estimation and control scheme for in/out-flux during managed-pressure drilling [4]. In terms of multi-source-data-fusion monitoring, Arévalo et al. (2022) integrated a transient hole-cleaning model with along-string measurements into a real-time digital twin to track cuttings distribution and reduce borehole-cleaning risk [6].

For the real-time determination of drilling-fluid rheology, Vajargah and van Oort (2015) proposed an approach that estimates downhole rheological properties from distributed pressure measurements, improving the physical rationality of rheology estimation [1]. For online hydraulic-model calibration, Altindal et al. (2025) proposed an online parameter calibration approach that dynamically updates drilling hydraulics models from real-time sensor data by integrating physics-based governing equations with data-driven techniques [5]. Habib et al. (2021) reviewed model-based and data-driven methods for early kick detection and estimation during managed-pressure drilling, including unscented-Kalman-filter-based approaches for detecting downhole abnormalities such as lost circulation and gas kick [8].

Nevertheless, existing inversion methods still possess two critical deficiencies. First, most studies exclusively adopt pressure measurements and fail to fully utilize torque signals containing rich rheological and flow-field information, leading to insufficient observational-constraint capability. Second, current prior-constraint strategies are unitary, lacking synchronous integration of mechanism-baseline prior and regional statistical prior information, which results in poor anti-noise performance and temporal stability of inversion results. Different from existing studies, this paper constructs a dual-prior coupled-constraint system based on dual observational data, comprehensively improving the accuracy, stability, and robustness of temporal inversion.

1.2. Fluid Mechanics Surrogate Models and Neural Operators

High-precision CFD numerical simulation is the mainstream approach for drilling-hydraulic-mechanism analysis, whereas its low computational efficiency cannot meet the iterative requirements of online real-time optimization [7]. Surrogate models replace high-cost numerical calculations by constructing input–output mapping relationships of physical systems, serving as an effective solution for efficient hydraulic prediction.

Traditional surrogate models such as response-surface methods and Kriging interpolation are only applicable to low-dimensional parameter-space prediction, with poor performance in high-dimensional full-flow-field fitting. In recent years, PINNs have become a research hotspot in fluid intelligent prediction owing to their physical-constraint characteristics. Raissi et al. (2019) pioneered the application of PINNs in forward and inverse partial-differential-equation (PDE) solving, providing a novel paradigm for fluid-mechanic intelligent computation [9]. Subsequently, Mao et al. (2020) applied PINNs to high-speed aerodynamic flow modeled by the Euler equations, showing that PINNs perform well on inverse problems but are less accurate than traditional numerical solvers for forward problems [10]. However, PINNs belong to discrete point-to-point-mapping models without resolution invariance, resulting in limited cross-condition generalization and frequent prediction failure under time-varying and boundary working conditions.

Neural operators represented by FNO break through the inherent limitations of traditional neural networks. Li et al. (2021) originally proposed the primitive FNO architecture [11]. Kovachki et al. (2023) systematically proposed neural-operator theories and universal-approximation theorems, verifying the superior performance of FNO in infinite-dimensional-function-space mapping [12]. On this basis, Geo-FNO was developed to adapt complex geometric domains, improving physical consistency for complex-domain computation [13], and physics-augmented variants have been proposed to improve prediction accuracy, particularly within the subject region [20]. At present, FNO has been successfully applied in hydraulic-tomography inversion and subsurface-flow-field prediction [14], yet its application in drilling-fluid annular two-phase hydraulic surrogate modeling remains unreported. This paper introduces FNO into drilling-hydraulic prediction and constructs a high-efficiency, high-generalization surrogate model to support real-time drilling optimization.

1.3. Multi-Objective Optimization of Drilling Parameters

Drilling-parameter optimization is the core approach to coordinate drilling efficiency, operational safety, and energy consumption. Mechanical-specific-energy (MSE) theory provides a quantitative evaluation index for rock-breaking efficiency and is widely applied in WOB and RPM optimization. Nystad et al. (2021) proposed an extremum-seeking-control-based real-time MSE-minimization strategy to realize automatic optimization of mechanical drilling parameters [16].

To balance multiple conflicting operational objectives, multi-objective-optimization algorithms have been gradually introduced into drilling-parameter optimization. Song et al. (2022) established a constrained Bayesian multi-objective optimization model that minimizes mechanical specific energy and drilling cost for real-time drilling-parameter optimization [18]. Peng et al. (2023) combined clustering analysis with a deep residual neural network to accurately predict the rate of penetration (ROP) in ultra-deep wells, providing a basis for drilling-parameter optimization [17]. Boukredera et al. (2023) integrated machine-learning models with an optimization algorithm to enhance drilling efficiency and mitigate drill-string vibrations [19].

Existing optimization studies still have prominent limitations. On the one hand, most optimization strategies only optimize single mechanical or hydraulic parameters, ignoring the strong-coupling interactions among WOB, RPM, and flow rate, which cannot achieve global-optimal drilling performance. On the other hand, most optimization methods adopt offline static-solving modes, which fail to dynamically adapt to real-time variations of downhole hydraulic parameters, resulting in poor field adaptability. Targeting the above problems, this paper constructs a hydraulic–mechanical-coupled three-parameter-coordinated online-optimization framework to realize adaptive real-time parameter optimization under time-varying drilling conditions.

2. Methodology

2.1. Basic Assumptions and Governing Equations

This study investigates the drilling-fluid–cuttings two-phase rotating flow in eccentric annuli. The governing equations are established based on the following basic assumptions:

1.

The drilling fluid is a non-Newtonian fluid conforming to the Herschel–Bulkley (H–B) rheological model.

2.

Cuttings are treated as a dispersed phase with particle diameters far smaller than the annular clearance, and the Euler–Euler two-fluid model is adopted for description.

3.

The flow is fully developed axial flow, and the circumferential shear effect induced by drill-string rotation is taken into account.

4.

The annular eccentricity varies gently along the well depth and is locally regarded as a constant eccentric annular space.

5.

The cuttings bed forms on the low side of the annulus, and its height is estimated by an empirical model.

The governing equations for annular two-phase flow consist of continuity equations and momentum equations. The continuity equations for the liquid phase (drilling fluid) and the solid phase (cuttings) are, respectively,

$${\displaystyle \frac{\partial \left({\alpha }_l{\rho }_l\right)}{\partial t}}+\nabla ·\left({\alpha }_l{\rho }_l{\mathbf{u}}_l\right)=0,$$

(1)

$${\displaystyle \frac{\partial \left({\alpha }_s{\rho }_s\right)}{\partial t}}+\nabla ·\left({\alpha }_s{\rho }_s{\mathbf{u}}_s\right)=0,$$

(2)

where ${\alpha }_l$ and ${\alpha }_s$ are the volume fractions of the liquid and solid phases, satisfying ${\alpha }_l+{\alpha }_s=1$; ${\rho }_l$ and ${\rho }_s$ are the densities of the liquid and solid phases; and ${\mathbf{u}}_l$ and ${\mathbf{u}}_s$ are the velocity vectors of the liquid and solid phases.

The axial momentum equations are

$${\displaystyle \frac{\partial \left({\alpha }_l{\rho }_l{u}_{lz}\right)}{\partial t}}+\nabla ·\left({\alpha }_l{\rho }_l{\mathbf{u}}_l{u}_{lz}\right)=-{\alpha }_l{\displaystyle \frac{\partial p}{\partial z}}+\nabla ·\left({\alpha }_l{\mathbf{\tau }}_l\right)+{\alpha }_l{\rho }_l{g}_z+M_{ls,z},$$

(3)

$${\displaystyle \frac{\partial \left({\alpha }_s{\rho }_s{u}_{sz}\right)}{\partial t}}+\nabla ·\left({\alpha }_s{\rho }_s{\mathbf{u}}_s{u}_{sz}\right)=-{\alpha }_s{\displaystyle \frac{\partial p}{\partial z}}+\nabla ·\left({\alpha }_s{\mathbf{\tau }}_s\right)+{\alpha }_s{\rho }_s{g}_z+M_{sl,z},$$

(4)

where p is the pressure; ${\mathbf{\tau }}_l$ and ${\mathbf{\tau }}_s$ are the viscous stress tensors of the liquid and solid phases; $u_{lz}$ and $u_{sz}$ are the axial velocity components; $g_z$ is the axial component of the gravitational acceleration; and $M_{ls,z}$ and $M_{sl,z}$ are the interphase momentum-exchange terms, satisfying $M_{ls,z}=-M_{sl,z}$.

2.2. Rheological Constitutive Relationship

This study adopts the Herschel–Bulkley (H–B) yield power-law model to describe the constitutive relationship between shear stress and shear rate of the drilling fluid, which serves as the core theoretical basis for hydraulic calculations throughout this paper:

$$\tau ={\tau }_y+K{\dot{\gamma }}^n,$$

(5)

where $\tau$ is the shear stress (Pa); ${\tau }_y$ is the fluid yield stress (Pa, corresponding to the engineering yield point $\mathrm{YP}$); K is the consistency index; $\dot{\gamma }$ is the shear rate; and n is the flow-behavior index (dimensionless).

In the H–B rheological system, the apparent viscosity is a dynamically variable quantity dependent on shear rate, rather than the plastic viscosity defined in drilling engineering. To resolve the mismatch between theoretical rheological parameters and engineering measured parameters, and to align with field standard six-speed rotary viscometer testing, this study introduces the industry-standard reference shear rate ${\dot{\gamma }}_0$ and defines the equivalent plastic viscosity:

$$\mathrm{PV}=K{\dot{\gamma }}_0^{n-1}.$$

(6)

In this study, the flow-behavior index n is pre-calibrated by laboratory rheological experiments as a fixed constant and does not participate in the while-drilling inversion calculation. The reference shear rate ${\dot{\gamma }}_0$ is set to the shear rate corresponding to the highest standard viscometer speed of the six-speed rotary viscometer (600 rpm, approximately $1021.7{\mathrm{s}}^{-1}$), so that the equivalent plastic viscosity $\mathrm{PV}$ defined by Equation (6) directly matches the plastic viscosity read from the field six-speed test; in particular, it equals the difference between the dial readings at 600 and 300 rpm for Bingham fluids. The inversion module solves for the two engineering parameters, equivalent plastic viscosity $\mathrm{PV}$ and yield point $\mathrm{YP}$. The consistency index K is then obtained through forward mapping via Equation (6), and combined with the fixed n and the inverted ${\tau }_y$, the three parameters of the H–B model are fully determined. For Bingham fluids ($n=1$), Equation (6) reduces to the classical engineering model, where $\mathrm{PV}=K$ and $\mathrm{YP}={\tau }_y$, fully consistent with drilling industry standard definitions.

Drill-string rotation significantly affects the annular velocity profile through shear-thinning effects and centrifugal forces. The circumferential shear induced by rotation reduces the apparent viscosity, and the equivalent rheological parameters are corrected as

$${\mathrm{PV}}_{\mathrm{eff}}=\mathrm{PV}·f_{\mathrm{rot}}(\mathrm{RPM},e),$$

(7)

$${\mathrm{YP}}_{\mathrm{eff}}=\mathrm{YP}·g_{\mathrm{rot}}(\mathrm{RPM},e),$$

(8)

where $f_{\mathrm{rot}}$ and $g_{\mathrm{rot}}$ are rotation-correction functions related to the rotary speed $\mathrm{RPM}$ and the annular eccentricity e, which are obtained by fitting numerical experiments.

2.3. Correction of Eccentric Annular Velocity Distribution

To obtain an explicit velocity-distribution expression containing spatial coordinates under eccentric geometry, while balancing the nonlinearity of the H–B model and computational efficiency, this study adopts the narrow-slot approximation. Within the local gap, the axial velocity profile is derived using the power-law fluid model (i.e., setting the yield stress to zero); the equivalent consistency index is converted from the inverted $\mathrm{PV}$ through Equation (6). This approximation yields controllable prediction errors for pressure loss and velocity distribution within the engineering-relevant shear-rate range.

The velocity distribution in eccentric annuli differs significantly from that in concentric annuli: the flow velocity is high at the wide gap and low at the narrow gap, causing cuttings to deposit readily on the low side. The local cross-sectional average axial velocity adopts an engineering expression based on the narrow-slot approximation:

$${\overline{u}}_z\left(\theta \right)={\displaystyle \frac{n}{2n+1}}{\left({\displaystyle \frac{\Delta P}{KL}}\right)}^{1/n}{\left(R_2-R_1\right)}^{1+1/n}·h(e,\theta ),$$

(9)

where $R_1$ is the outer radius of the drill string, $R_2$ is the borehole radius, $\Delta P/L$ is the pressure drop per unit length, and $h(e,\theta )$ is the eccentricity-correction function. The eccentricity e is defined as

$$e={\displaystyle \frac{\delta }{R_2-R_1}},$$

(10)

where $\delta$ is the offset distance between the drill-string and wellbore centers. The eccentricity ratio satisfies $e\in [0,1]$; $e=0$ represents a concentric configuration, and $e=1$ represents the drill string contacting the wellbore wall.

2.4. Cuttings Two-Phase Additional Friction and Cuttings Bed Model

The presence of cuttings increases annular flow resistance, and the additional pressure loss is related to the cuttings concentration. A modified two-phase-flow pressure-loss model is adopted:

$${\displaystyle \frac{\Delta P}{L}}={\left({\displaystyle \frac{\Delta P}{L}}\right)}_l·F\left(C_s\right),$$

(11)

where ${(\Delta P/L)}_l$ is the pressure-loss gradient of the pure drilling fluid, $F\left(C_s\right)$ is the cuttings-concentration correction coefficient, and $C_s$ is the average annular cuttings volume concentration.

The cuttings-bed height $h_b$ is estimated through the cuttings-transport equilibrium relationship:

$$h_b=f(C_s,\mathrm{ROP},Q,\mathrm{PV},\mathrm{YP},d_p,e,{\theta }_{\mathrm{dev}}),$$

(12)

where $d_p$ is the average cuttings particle diameter and ${\theta }_{\mathrm{dev}}$ is the well deviation angle. The presence of the cuttings bed further changes the annular effective flow area and velocity distribution, forming a positive-feedback effect.

2.5. Standpipe Pressure and Rotary Torque Calculation Model

Standpipe pressure is the measured circulating pressure at the drill-floor standpipe during drilling-fluid circulation. The circulating standpipe pressure consists of three dynamic components: the frictional pressure loss inside the drill string, the pressure drop across the bit nozzles due to sudden velocity increase, and the frictional pressure loss in the annulus as the fluid returns to the surface. Because the hydrostatic pressure of the drilling-fluid column on the drill-string side is largely balanced by that of the returning column in the annulus, the net static-pressure difference is negligibly small compared with the dynamic friction components and is therefore not included in the circulating standpipe pressure. Accordingly, the expression for standpipe pressure is

$$P_{\mathrm{standpipe}}=\Delta P_{\mathrm{pipe}}+\Delta P_{\mathrm{bit}}+\Delta P_{\mathrm{annulus}},$$

(13)

where $\Delta P_{\mathrm{pipe}}$ is the frictional pressure loss inside the drill string, $\Delta P_{\mathrm{bit}}$ is the pressure drop across the bit nozzles, and $\Delta P_{\mathrm{annulus}}$ is the frictional pressure loss in the annulus. The standpipe-pressure gauge reading primarily reflects the real-time variations of these three dynamic friction components, which are the quantities relevant to the subsequent hydraulic-parameter inversion.

The rotary torque originates from viscous friction between the drill string and fluid, cuttings-bed friction, and bit–formation interaction. The hydraulics-related torque component can be expressed as

$$T_{\mathrm{hyd}}=2\pi {\int }_0^L{\tau }_w(R_1\left(z\right),z)R_1{\left(z\right)}^2\mathrm{d}z,$$

(14)

where L denotes the total length of the drill string; $R_1\left(z\right)$ represents the outer radius of the drill string, which is piecewise-assigned along the wellbore depth to account for the different diameters of the drill pipe and drill collars; and ${\tau }_w$ denotes the circumferential wall shear stress acting on the outer wall of the drill string, which is closely related to the rheological parameters, the velocity distribution, and the eccentricity.

The above forward model constitutes the physical basis for the inversion and optimization in this study. Among the model parameters, the plastic viscosity $\mathrm{PV}$, yield point $\mathrm{YP}$, average annular cuttings concentration $C_s$, and equivalent eccentricity e are the time-varying parameters to be inverted while drilling, whereas the remaining parameters such as wellbore structure, drill-string dimensions, and cuttings density are known quantities.

2.6. Dual-Prior Constrained Temporal Inversion Module

2.6.1. Problem Formulation and Sliding-Window Mechanism

The objective of while-drilling hydraulic parameter inversion is as follows: given the observed standpipe-pressure sequence $P_{\mathrm{obs}}\left(t\right)$, the observed rotary-torque sequence $T_{\mathrm{obs}}\left(t\right)$, and the flow-rate sequence $Q\left(t\right)$, estimate in real time the four time-varying parameters—plastic viscosity $\mathrm{PV}\left(t\right)$, yield point $\mathrm{YP}\left(t\right)$, average annular cuttings concentration $C_s\left(t\right)$, and equivalent eccentricity $e\left(t\right)$.

Define the parameter vector to be inverted as $\mathbf{\theta }\left(t\right)={[\mathrm{PV}\left(t\right),\mathrm{YP}\left(t\right),C_s\left(t\right),e\left(t\right)]}^T$. Because the fluid properties and annular state change relatively slowly during drilling, they can be regarded as quasi-static within a short time window. Therefore, this study adopts a sliding-window mechanism, dividing the continuous time series into fixed windows of length L, within each of which the parameters are treated as constants, and the estimated values are updated window by window.

Let the k-th time window contain N sampling points; the observation vector is

$${\mathbf{y}}_k={[P_{\mathrm{obs},1},\dots ,P_{\mathrm{obs},N},T_{\mathrm{obs},1},\dots ,T_{\mathrm{obs},N}]}^T.$$

(15)

The forward model is denoted as $\mathcal{F}({\mathbf{\theta }}_k,{\mathbf{Q}}_k)$, where ${\mathbf{Q}}_k$ is the flow-rate sequence within the window, and its output is the corresponding predicted pressure and torque values. The inversion problem is transformed into an optimization problem within each window:

$${\mathbf{\theta }}_k^{*}=arg\underset{{\mathbf{\theta }}_k}{min}\mathcal{L}\left({\mathbf{\theta }}_k\right),$$

(16)

where $\mathcal{L}$ is the total loss function.

2.6.2. Dual-Prior Constraint Design

This study introduces two classes of prior information to constrain the inversion process, constituting a dual-prior constraint mechanism.

First class: pre-drilling mechanism baseline prior. During the drilling design phase, hydraulic simulations are usually performed based on formation prediction and drilling-fluid formulation design to obtain a set of baseline parameters ${\mathbf{\theta }}_{\mathrm{base}}$, referred to as the mechanism baseline prior. This prior reflects the physically expected parameter range from the design and serves as the initial baseline and regularization constraint for inversion. The baseline prior adopts an $L_2$ regularization form:

$${\mathcal{L}}_{\mathrm{baseline}}={\parallel {\mathbf{\theta }}_k-{\mathbf{\theta }}_{\mathrm{base}}\parallel }_2^2.$$

(17)

Second class: offset-well statistical prior. Data from multiple offset wells within the same block contain statistical laws of formation and drilling-fluid characteristics in that region. The measured inversion results from offset wells are collected, and a multivariate Gaussian distribution is fitted:

$$\mathbf{\theta }\sim \mathcal{N}({\mathbf{\mu }}_{\mathrm{well}},{\Sigma }_{\mathrm{well}}),$$

(18)

where ${\mathbf{\mu }}_{\mathrm{well}}$ is the mean vector of offset-well parameters and ${\Sigma }_{\mathrm{well}}$ is the covariance matrix. The offset-well prior adopts a regularization term in the form of the Mahalanobis distance, corresponding to the negative log-likelihood of the Gaussian distribution:

$${\mathcal{L}}_{\mathrm{neighbor}}={({\mathbf{\theta }}_k-{\mathbf{\mu }}_{\mathrm{well}})}^T{\Sigma }_{\mathrm{well}}^{-1}({\mathbf{\theta }}_k-{\mathbf{\mu }}_{\mathrm{well}}).$$

(19)

The Mahalanobis distance accounts for correlations among parameters, enabling a more accurate measure of the degree to which the parameters deviate from the block statistical laws.

2.6.3. Four-Term Loss Function

The total loss function consists of four terms, balancing data fitting and prior constraints:

$$\mathcal{L}={\mathcal{L}}_{\mathrm{fit}}+{\lambda }_1{\mathcal{L}}_{\mathrm{baseline}}+{\lambda }_2{\mathcal{L}}_{\mathrm{neighbor}}+{\lambda }_3{\mathcal{L}}_{\mathrm{smooth}}.$$

(20)

Term 1: observation fitting loss. Both standpipe pressure and rotary torque observations are fitted simultaneously, utilizing dual-observation constraints to improve parameter identifiability:

$${\mathcal{L}}_{\mathrm{fit}}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\left[w_p{\left({\displaystyle \frac{P_{\mathrm{pred},i}-P_{\mathrm{obs},i}}{P_{\mathrm{obs},i}}}\right)}^2+w_t{\left({\displaystyle \frac{T_{\mathrm{pred},i}-T_{\mathrm{obs},i}}{T_{\mathrm{obs},i}}}\right)}^2\right],$$

(21)

where $w_p$ and $w_t$ are the weighting coefficients for pressure and torque, respectively, adjusted according to the observation signal-to-noise ratio. The relative-error form is adopted to eliminate dimensional effects.

Term 2: pre-drilling mechanism baseline regularization. As shown in Equation (17), this term constrains the inverted parameters from deviating too far from the design baseline, preventing physically unreasonable results.

Term 3: offset-well statistical prior regularization. As shown in Equation (19), this term constrains the inverted parameters to conform to offset-well statistical laws, improving engineering rationality.

Term 4: temporal smoothness constraint. Parameters in adjacent windows should maintain continuity, suppressing high-frequency fluctuations induced by noise:

$${\mathcal{L}}_{\mathrm{smooth}}={\parallel {\mathbf{\theta }}_k-{\mathbf{\theta }}_{k-1}\parallel }_2^2,$$

(22)

where ${\mathbf{\theta }}_{k-1}$ is the inversion result from the previous window.

2.6.4. Solution Algorithm and Uncertainty Estimation

Minimization of the loss function is performed using the L-BFGS-B algorithm, which is suitable for small-to-medium-scale bounded optimization problems and exhibits fast convergence. Parameter upper and lower bounds are set according to physical meaning and engineering experience.

To quantify inversion uncertainty, a second-order Taylor expansion of the loss function is performed around the optimal solution, yielding an approximate posterior covariance matrix of the parameters:

$${\Sigma }_k\approx {\left({\nabla }^2\mathcal{L}\left({\mathbf{\theta }}_k^{*}\right)\right)}^{-1},$$

(23)

where ${\nabla }^2\mathcal{L}$ is the Hessian matrix of the loss function. This covariance matrix provides the uncertainty basis for subsequent probabilistic safety-constrained optimization. Note that the fitting loss in Equation (21) is built on weighted relative residuals, and the weighting coefficients $w_p$ and $w_t$ effectively encode the inverse of the observation noise variance. Accordingly, if the observation noise level ${\sigma }^2$ is known explicitly, the posterior covariance in Equation (23) should be scaled as ${\Sigma }_k\approx {\sigma }^2{\left({\nabla }^2\mathcal{L}\left({\mathbf{\theta }}_k^{*}\right)\right)}^{-1}$; otherwise the un-scaled form may under-estimate the parameter uncertainty by neglecting the noise-induced scale of the fitting residuals. In this study, the weights $w_p$ and $w_t$ are calibrated from the field observation signal-to-noise ratio, so the noise scale is implicitly absorbed, and Equation (23) provides a consistent relative uncertainty estimate.

2.7. FNO-Based Hydraulic Surrogate Model

2.7.1. Surrogate Model Positioning and Input–Output Definition

Numerical simulation of eccentric rotating-annular two-phase flow entails enormous computational cost and cannot support the thousands of online optimization iterations required under the sliding-window framework. This study adopts a one-dimensional Fourier neural operator (FNO) to construct a hydraulic surrogate model, learning the solution operator of the parameterized two-phase-flow governing equations and establishing a rapid mapping from rheological parameters and drilling operational parameters to the full-domain hydraulic response along the well depth.

The input vector of the surrogate model is defined as

$$X_{\mathrm{in}}=[\mathrm{PV},\mathrm{YP},C_s,e,\mathrm{WOB},\mathrm{RPM},Q,{\theta }_{\mathrm{dev}}].$$

(24)

The model outputs are one-dimensional discrete sequences along the well depth: the full-domain annular pressure-loss distribution, the equivalent circulating density $\mathrm{ECD}\left(z\right)$, the cuttings-bed height $h_b\left(z\right)$, and the bit hydraulic power sequence.

2.7.2. FNO Network Architecture

This study adopts a one-dimensional FNO architecture adapted to longitudinal well-depth field prediction. The network consists of an input lifting layer, multiple Fourier operator layers, nonlinear activation layers, and an output projection layer:

1.

Lifting layer: maps the low-dimensional parameter vector to the hidden channel dimension $d_{\mathrm{hidden}}$.

2.

Fourier operator core layers:

$${\left(\mathcal{F}u\right)}_k\left(z\right)=\sigma \left(Wu\left(z\right)+{\mathcal{F}}^{-1}\left(R_k\odot \mathcal{F}\left(u\right)\right)\left(z\right)\right),$$

(25)

where $\mathcal{F}$ and ${\mathcal{F}}^{-1}$ are the one-dimensional Fourier forward and inverse transforms, $R_k$ is the learnable kernel in Fourier space, ⊙ denotes the Hadamard product, W is the local linear transformation, and the activation function $\sigma$ is chosen as GELU. The model is configured with four Fourier operator layers, and the number of highest retained Fourier modes is truncated to balance accuracy and computational cost.

3.

Output projection layer: maps the high-dimensional hidden features to the hydraulic response field sequences.

2.7.3. Training Dataset Construction

All training samples are generated in batch by the eccentric annular two-phase hydraulic forward model described in Section 2, with parameters covering the engineering-reasonable ranges:

$$\begin{array}{cc}\hfill & \mathrm{PV}\in [10,80]\mathrm{mPa}·\mathrm{s},\mathrm{YP}\in [2,30]\mathrm{Pa},\hfill \\ \hfill & C_s\in [0,0.15],e\in [0,0.95],\hfill \\ \hfill & Q\in [15,40]\mathrm{L}/\mathrm{s},\mathrm{RPM}\in [40,180],\hfill \\ \hfill & {\theta }_{\mathrm{dev}}\in [0^{\circ },{90}^{\circ }].\hfill \end{array}$$

(26)

The dataset is split into 85% for training, 10% for validation, and 5% for testing; all input and output samples are uniformly subjected to min–max normalization to eliminate training bias caused by different physical dimensions.

2.7.4. Training Loss Function

The basic training adopts the field-averaged relative $L_2$ loss:

$${\mathcal{L}}_{\mathrm{FNO}}={\displaystyle \frac{1}{N_s}}{\displaystyle \sum _{i=1}^{N_s}}{\displaystyle \frac{{\parallel {\mathbf{Y}}_{\mathrm{pred},i}-{\mathbf{Y}}_{\mathrm{true},i}\parallel }_2}{{\parallel {\mathbf{Y}}_{\mathrm{true},i}\parallel }_2+\epsilon }},$$

(27)

where $N_s$ is the number of samples, ${\mathbf{Y}}_{\mathrm{pred},i}$ is the hydraulic field predicted by the surrogate model, ${\mathbf{Y}}_{\mathrm{true},i}$ is the numerical-simulation ground truth, and $\epsilon$ is a tiny constant to avoid division by zero.

An optional physics-informed regularization term can be introduced using the residual of the two-phase governing equations, forming a physics-enhanced FNO loss:

$${\mathcal{L}}_{\mathrm{phy}}={\lambda }_{\mathrm{data}}·{\mathcal{L}}_{\mathrm{FNO}}+{\lambda }_{\mathrm{pde}}·{\mathcal{L}}_{\mathrm{pde}},$$

(28)

where ${\mathcal{L}}_{\mathrm{pde}}$ is the residual loss of the continuity and momentum equations, effectively reducing prediction errors under extrapolation conditions beyond the training distribution. The main model in this study adopts a purely data-driven FNO; the effect of the physics-enhanced scheme can be compared in ablation experiments.

2.7.5. Training Strategy and Online Inference Workflow

The optimizer adopts AdamW with an initial learning rate of $1\times {10}^{-3}$, accompanied by cosine-annealing learning-rate decay. An early-stopping strategy is adopted: the validation loss is continuously monitored, and training is terminated after 30 consecutive epochs without improvement, saving the optimal weights.

The real-time inference workflow is as follows:

1.

Receive the real-time parameters ${\mathbf{\theta }}_k$ output from the inversion module described in Section 2.

2.

Concatenate with the current drilling operational parameters to form the complete input vector.

3.

Perform the FNO forward pass to output the full-domain ECD, cuttings-bed height, and annular pressure loss at the millisecond level.

4.

Feed the prediction results into the multi-objective optimization module to serve as the rapid evaluation kernel of the objective functions.

2.8. Hydraulic–Mechanical Coupled Multi-Objective Coordinated Optimization Framework

2.8.1. Multi-Objective Optimization Mathematical Model

The decision variables to be optimized are

$$\mathbf{x}={[\mathrm{WOB},\mathrm{RPM},Q]}^T.$$

(29)

Objective 1: maximize the rate of penetration ($\mathrm{ROP}$). The rate of penetration is the core indicator of drilling efficiency, and the drillability empirical model applicable to rotary drilling is adopted:

$$f_1\left(\mathbf{x}\right)=\mathrm{ROP}=K_d·{\mathrm{WOB}}^{a_1}·{\mathrm{RPM}}^{a_2}·e^{-\beta V_{\mathrm{risk}}},$$

(30)

where $K_d$ is the formation drillability coefficient, $a_1$ and $a_2$ are fitting exponents, and $\beta$ is the vibration-suppression coefficient. The optimization objective is to maximize $f_1\left(\mathbf{x}\right)$.

Objective 2: minimize the mechanical specific energy ($\mathrm{MSE}$). The mechanical specific energy reflects rock-breaking efficiency, defined as the energy required to break a unit volume of rock:

$$f_2\left(\mathbf{x}\right)=\mathrm{MSE}={\displaystyle \frac{\mathrm{WOB}}{A}}+{\displaystyle \frac{2\pi ·\mathrm{RPM}·T}{A·\mathrm{ROP}}},$$

(31)

where A is the bottom-hole area and T is the comprehensive bit torque. The optimization objective is to minimize $f_2\left(\mathbf{x}\right)$.

Objective 3: ECD safety-window deviation. The equivalent circulating density must be controlled within the safe window between formation pore pressure and fracture pressure. The safety deviation is defined as

$$f_3\left(\mathbf{x}\right)=max\left(0,{\mathrm{ECD}}_{max}-{\mathrm{ECD}}_{\mathrm{upper}}\right)+max\left(0,{\mathrm{ECD}}_{\mathrm{lower}}-{\mathrm{ECD}}_{min}\right),$$

(32)

where ${\mathrm{ECD}}_{\mathrm{upper}}$ and ${\mathrm{ECD}}_{\mathrm{lower}}$ are the upper and lower bounds of the safety window, respectively. The optimization objective is to minimize $f_3\left(\mathbf{x}\right)$, whose ideal value is zero.

Objective 4: minimize the cuttings-bed height. Excessive cuttings-bed height can cause stuck pipe and drag issues; the maximum cuttings-bed height along the full well section is taken:

$$f_4\left(\mathbf{x}\right)=\underset{z}{max}{h}_b(z,\mathbf{x},\mathbf{\theta }).$$

(33)

The optimization objective is to minimize $f_4\left(\mathbf{x}\right)$.

Objective 5: minimize the vibration risk. Excessive weight on bit and rotary speed can induce whirl and bit-bounce vibrations; an empirical vibration-risk model is adopted:

$$f_5\left(\mathbf{x}\right)=V_{\mathrm{risk}}(\mathrm{WOB},\mathrm{RPM}).$$

(34)

The optimization objective is to minimize $f_5\left(\mathbf{x}\right)$.

The constraints include the physical upper and lower bounds of each parameter and the equipment-capacity limits:

$${\mathrm{WOB}}_{min}\le \mathrm{WOB}\le {\mathrm{WOB}}_{max},{\mathrm{RPM}}_{min}\le \mathrm{RPM}\le {\mathrm{RPM}}_{max},Q_{min}\le Q\le Q_{max}.$$

(35)

The multi-objective optimization problem is uniformly converted into minimization form:

$$\underset{\mathbf{x}\in \Omega }{min}\mathbf{F}\left(\mathbf{x}\right)={[-f_1\left(\mathbf{x}\right),f_2\left(\mathbf{x}\right),f_3\left(\mathbf{x}\right),f_4\left(\mathbf{x}\right),f_5\left(\mathbf{x}\right)]}^T,$$

(36)

where $\Omega$ is the feasible domain.

2.8.2. Online Bayesian Multi-Objective Optimization Algorithm

Traditional multi-objective optimization algorithms such as NSGA-II require a large number of objective-function evaluations, incurring high computational cost and making them unsuitable for real-time applications. This study adopts a Pareto-frontier-based Bayesian multi-objective optimization algorithm, which intelligently selects evaluation points through an acquisition function, significantly reducing the number of function calls.

Surrogate model. An independent Gaussian-process (GP) surrogate model is constructed for each objective function, utilizing information from already-evaluated points to predict the objective values and uncertainties of unevaluated points.

Acquisition function. The Expected Hypervolume Improvement (EHVI) criterion is adopted, measuring the expected improvement in the hypervolume of the Pareto frontier contributed by a new evaluation point; the point with the maximum EHVI is selected for evaluation.

Mechanism-gradient acceleration. Leveraging the automatic-differentiation capability of the FNO surrogate model, the analytical gradient of the objective function with respect to the decision variables is computed, and the gradient information is incorporated into the acquisition function of Bayesian optimization to accelerate convergence.

2.8.3. Sliding-Window Rolling Update Mechanism

The optimization process proceeds synchronously with the inversion module, adopting a sliding-window rolling update strategy:

1.

After each inversion window is completed, the parameter estimate ${\mathbf{\theta }}_k$ and its confidence interval are updated.

2.

Based on the latest inversion results, the objective-function surface is recomputed.

3.

Online Bayesian optimization is executed to update the Pareto frontier and the optimal parameter recommendations.

4.

The optimization results are output to drilling operators or automatic control systems.

Because the inference speed of the FNO surrogate model is extremely fast, hundreds of objective-function evaluations can be completed within each window, ensuring that the optimization results are sufficiently converged.

2.8.4. Probabilistic Safety Constraints Based on Confidence Intervals

Inverted parameters possess uncertainty; direct use of mean values for optimization may lead to risk underestimation. This study constructs probabilistic safety constraints based on the confidence intervals of inverted parameters, ensuring that optimization results satisfy safety requirements at a certain confidence level. For the ECD safety constraint, the following must be satisfied:

$$\mathrm{P}\left(\mathrm{ECD}(\mathbf{x},\mathbf{\theta })\le {\mathrm{ECD}}_{\mathrm{upper}}\right)\ge 1-\alpha ,\mathrm{P}\left(\mathrm{ECD}(\mathbf{x},\mathbf{\theta })\ge {\mathrm{ECD}}_{\mathrm{lower}}\right)\ge 1-\alpha ,$$

(37)

where $\alpha$ is the risk level, typically taken as 0.05. Utilizing the posterior Gaussian-distribution assumption of the parameters, the distribution of ECD is approximated through a first-order Taylor expansion, converting the probabilistic constraints into deterministic constraints:

$$\mathrm{ECD}(\mathbf{x},{\mathbf{\mu }}_k)+z_{1-\alpha }·{\sigma }_{\mathrm{ECD}}(\mathbf{x},{\Sigma }_k)\le {\mathrm{ECD}}_{\mathrm{upper}},$$

(38)

where $z_{1-\alpha }$ is the quantile of the standard normal distribution and ${\sigma }_{\mathrm{ECD}}$ is the standard deviation of ECD, obtained through uncertainty propagation. Similarly, corresponding probabilistic constraints are constructed for other safety indicators such as the cuttings-bed height. The introduction of probabilistic safety constraints renders the optimization results robust, capable of absorbing the impact of inversion uncertainty.

2.8.5. Decision Scheme Selection

Multi-objective optimization yields a set of Pareto-optimal solutions; the final execution scheme must be selected according to actual field requirements. This study provides three decision modes:

1.

Safety-priority mode: selects the scheme with the largest ECD safety margin, applicable to high-risk conditions such as narrow density windows.

2.

Efficiency-priority mode: selects the scheme with the highest ROP, applicable to steady-inclination sections with wide safety windows.

3.

Balanced mode: adopts the TOPSIS method to normalize and weight each objective, selecting the scheme with the highest comprehensive score.

Operators can flexibly select the decision mode according to specific working conditions, or manually pick suitable parameter combinations directly from the Pareto frontier.

3. Results and Discussion

3.1. Example Setup

To validate the effectiveness of the proposed method, numerical simulation examples are constructed. All data are generated from the hydraulic PDE numerical simulation described in Section 2, without involving laboratory or field experiments. The base well parameters are listed in Table 1.

The dynamic variations of the four parameters during drilling are simulated: $\mathrm{PV}$ and $\mathrm{YP}$ slowly decrease with increasing well depth due to rising temperature, $C_s$ varies with ROP fluctuations, and e increases with increasing well deviation. Gaussian noise with an amplitude of 5% is added to the observation data to simulate actual measurement errors.

Three baseline methods are used for comparison: (i) traditional least-squares inversion (LS), which uses only pressure observations without prior constraints; (ii) a pure PINN prediction model, which directly learns the input–output mapping without an explicit inversion module; and (iii) offline NSGA-II optimization, which performs one-time optimization based on initial parameters without while-drilling updating.

Three groups of evaluation metrics are adopted: (i) inversion accuracy, quantified by the root-mean-square error (RMSE) and mean absolute percentage error (MAPE) of each parameter; (ii) surrogate-model accuracy, quantified by the full-field relative $L_2$ error and inference time; and (iii) optimization performance, quantified by the ROP improvement rate, MSE reduction rate, safety margin, and a comprehensive efficiency index.

3.2. Analysis of Dual-Prior Constrained Inversion Results

3.2.1. Inversion Accuracy Comparison

Table 2 compares the estimation accuracy of the four parameters obtained by different inversion methods.

As can be seen from Table 2, the traditional least-squares inversion yields the largest errors, with the MAPE of the four parameters ranging from 12.7% to 18.6%. After introducing a single prior constraint, the errors are reduced to some extent, with the offset-well prior performing slightly better than the baseline prior. The dual-prior constrained method proposed in this study achieves the optimal accuracy, with the MAPE of the four parameters reduced to 7.3%, 8.7%, 10.6%, and 8.2%, respectively, representing an average reduction of 42.7% compared with the least-squares method, which verifies the synergistic enhancement effect of the dual-prior constraints.

3.2.2. Temporal Stability Analysis

The least-squares inversion curve fluctuates severely, with multiple spikes, obviously disturbed by observation noise. In contrast, the dual-prior constrained inversion curve shows good agreement with the true values, and is temporally smooth without abnormal jumps. This benefits from the temporal smoothness constraint term and the regularization effect of the dual priors, enabling the inversion results to track the true variation trends while effectively suppressing noise interference.

Figure 1. Comparison of plastic-viscosity ($\mathrm{PV}$) inversion RMSE among least-squares inversion, single-prior (baseline/offset-well) inversion, and the proposed dual-prior-constrained method (unit: mPa·s). The dual-prior method achieves the lowest RMSE, demonstrating the synergistic effect of combining the pre-drill mechanistic baseline prior with the offset-well statistical prior.

Figure 1. Comparison of plastic-viscosity ($\mathrm{PV}$) inversion RMSE among least-squares inversion, single-prior (baseline/offset-well) inversion, and the proposed dual-prior-constrained method (unit: mPa·s). The dual-prior method achieves the lowest RMSE, demonstrating the synergistic effect of combining the pre-drill mechanistic baseline prior with the offset-well statistical prior.

Figure 2. Comparison of plastic-viscosity ($\mathrm{PV}$) inversion along well depth among different methods. The least-squares curve fluctuates sharply under observation noise, whereas the dual-prior-constrained curve closely follows the true variation without abnormal jumps, verifying the temporal stability of the proposed method.

Figure 2. Comparison of plastic-viscosity ($\mathrm{PV}$) inversion along well depth among different methods. The least-squares curve fluctuates sharply under observation noise, whereas the dual-prior-constrained curve closely follows the true variation without abnormal jumps, verifying the temporal stability of the proposed method.

3.2.3. Ablation Experiments

To verify the independent contribution of each loss term, ablation experiments are designed, removing each term from the loss function one by one to compare the changes in inversion performance, as summarized in Table 3.

The ablation results show that all three constraint terms make positive contributions to inversion accuracy. Among them, the offset-well prior contributes the most; after removal, the overall RMSE increases by 25.3%. The baseline prior is the second, with an increase of 17.2%. The smoothness constraint contributes relatively less, but still reaches 12.6%. When all constraints are removed leaving only the fitting term, the RMSE increases substantially by 74.7%, fully demonstrating the necessity of the multi-constraint mechanism.

3.3. Performance Analysis of FNO Surrogate Model

3.3.1. Prediction Accuracy

Table 4 compares the prediction accuracy of the FNO surrogate model and other methods on the test set.

The FNO model achieves the highest accuracy on all three outputs, with relative errors of ECD and pressure loss both below 2%, and the relative error of the cuttings-bed height at 3.24%, meeting engineering accuracy requirements. In contrast, PINN exhibits larger errors, especially in cuttings-bed prediction with an error approaching 9%, because PINN has inherent deficiencies in cross-condition generalization. CNN-Unet, as a traditional convolutional architecture, has intermediate accuracy. The response-surface method, as a classical surrogate model, yields the largest errors and struggles to handle high-dimensional nonlinear mappings.

3.3.2. Computational Efficiency

Table 5 compares the single-sample inference time of different methods.

The single-sample inference time of FNO is only 8 ms, achieving a speed-up of more than three orders of magnitude compared with traditional numerical simulation. Even compared with other deep-learning methods, FNO demonstrates a distinct speed advantage, providing sufficient computational margin for real-time optimization.

Figure 3. Comparison of single-sample inference time among conventional numerical simulation, PINN, CNN-Unet, and the proposed FNO surrogate model. The FNO model achieves millisecond-level inference, more than three orders of magnitude faster than conventional numerical simulation.

Figure 3. Comparison of single-sample inference time among conventional numerical simulation, PINN, CNN-Unet, and the proposed FNO surrogate model. The FNO model achieves millisecond-level inference, more than three orders of magnitude faster than conventional numerical simulation.

3.3.3. Generalization Capability Test

To verify the cross-condition generalization capability of the model, extreme conditions outside the training parameter range are selected for extrapolation testing. The results show that within a 10% parameter-extrapolation range, the prediction error of FNO remains within 5%, exhibiting good generalization performance. This benefits from the function-learning characteristics of neural operators, which possess stronger physical extrapolation capability compared with purely data-driven CNN methods.

3.4. Analysis of Multi-Objective Coordinated Optimization Results

3.4.1. Pareto Frontier Analysis

The Pareto frontier clearly demonstrates the trade-off relationship between drilling efficiency and energy consumption. The initial design parameters are located at the lower left of the Pareto frontier, indicating considerable optimization space. Through coordinated optimization, ROP can be improved by approximately 22% while maintaining MSE unchanged, or MSE can be reduced by approximately 18% while maintaining ROP unchanged.

Figure 4. Schematic of the Pareto front of the multi-objective optimization, illustrating the trade-off between drilling efficiency (ROP) and energy consumption (MSE). Coordinated optimization moves the initial operating point toward the Pareto front, where either ROP can be increased at constant MSE or MSE can be reduced at constant ROP.

Figure 4. Schematic of the Pareto front of the multi-objective optimization, illustrating the trade-off between drilling efficiency (ROP) and energy consumption (MSE). Coordinated optimization moves the initial operating point toward the Pareto front, where either ROP can be increased at constant MSE or MSE can be reduced at constant ROP.

3.4.2. Comparison of Different Optimization Strategies

Table 6 compares the comprehensive performance of different optimization strategies.

Single-parameter optimization yields limited improvement, with the comprehensive efficiency index increasing by only 3%–6%. Dual-parameter optimization shows significant improvement, with a 17% increase. The three-parameter coordinated optimization proposed in this study achieves the best results, with ROP improved by 24.3%, MSE reduced by 19.6%, and the comprehensive efficiency index increased by 28.5%. Offline NSGA-II performs inferior to the online coordinated optimization method because it cannot dynamically update with parameter variations, and its ECD safety margin is relatively small, posing certain risks.

3.4.3. Effect of Probabilistic Safety Constraints

To verify the robustness of the probabilistic safety constraints, the risk levels of deterministic optimization and probabilistically constrained optimization are compared. Under the condition of inversion-parameter uncertainty, the ECD exceedance probability of the deterministic optimization scheme is 12.7%, whereas that of the probabilistically constrained optimization scheme is only 3.2%, satisfying the 5% risk-design requirement. This demonstrates that confidence-interval-based probabilistic safety constraints can effectively absorb inversion uncertainty and guarantee the reliability of the optimization results.

3.5. Ablation Experiment Analysis

To verify the independent contribution of each core module one by one, four groups of ablation experiments are designed: basic inversion + offline optimization, single-prior constrained inversion + online optimization, dual-prior constrained inversion + PINN surrogate model, and the complete framework proposed in this study. The experimental results show that optimization improvements from a single module can only achieve local performance enhancement and cannot simultaneously account for inversion accuracy, inference efficiency, and dynamic optimization effectiveness.

To enable an intuitive comparison of the four ablation configurations, each dimension in Figure 5 is quantified by a normalized score. Four dimensions are defined: (i) inversion accuracy, measured by the average MAPE of the four inverted parameters; (ii) prediction speed, measured by the single-sample inference time of the surrogate; (iii) optimization gain, measured by the comprehensive efficiency index; and (iv) safety margin, measured by the ECD safety margin and the corresponding ECD violation probability. Because these dimensions have different physical units, each dimension is normalized independently: the best-performing configuration is anchored to a score of 100, and the remaining configurations are mapped to lower scores according to their relative gaps with respect to the best value. The resulting dimensionless scores range from 0 to 100 and are used only for relative comparison rather than as absolute accuracy metrics. The complete framework achieves the best raw indicators in the first three dimensions and is therefore scored 100 for them; in the safety dimension it is scored 96 because the probabilistically constrained scheme still retains a small residual ECD violation probability (3.2%), which satisfies the 5% risk-design requirement but does not completely eliminate the risk.

Based on these normalized scores, the independent contribution of each module is examined configuration by configuration. When only a single prior constraint is introduced, inversion robustness improves slightly, but the absence of offset-well statistical laws leads to large parameter-identification deviations under extreme working conditions. After replacing the surrogate model with PINN, the online-optimization iteration speed drops substantially, failing to meet real-time decision-making requirements. When the sliding-window update mechanism is removed, the optimization scheme becomes statically fixed, and dynamic adaptation capability is completely lost. Under the complete framework, dual-prior constraints guarantee high accuracy and stability of parameter inversion, the FNO surrogate model provides efficient forward-modeling support, and online Bayesian optimization with probabilistic constraints achieves dynamic, safe, and efficient coordinated parameter regulation. The synergistic coupling of the three modules jointly constitutes the core performance advantage of the proposed method, verifying the rationality and necessity of the three-layer framework design.

3.6. Discussion

3.6.1. Mechanistic Interpretation of the Inversion Improvement

The inversion accuracy gain originates from two coupled mechanisms. First, the dual observations resolve the ill-posedness inherent in single-channel inversion: standpipe pressure primarily constrains the frictional and hydraulic characteristics of the system, whereas rotary torque carries complementary information on the rheological and velocity-field state, so their joint fitting substantially improves parameter identifiability. Second, the two priors play complementary roles in stabilizing the solution. The pre-drill mechanistic baseline prior anchors the inverted parameters to the physically expected range, preventing convergence to non-physical solutions, while the offset-well statistical prior encodes block-scale correlations among parameters through the Mahalanobis distance, which suppresses the high-frequency fluctuations induced by observation noise. This explains the ablation results: removing the offset-well prior degrades the RMSE most severely (25.3%) because it carries the strongest statistical information, whereas the smoothness constraint contributes least (12.6%) because its role is partially redundant with the two priors under slowly varying conditions.

3.6.2. Mechanistic Interpretation of the Surrogate Efficiency and Generalization

The superiority of FNO over PINN originates from its operator-learning property. By parameterizing the integral kernel directly in Fourier space, FNO learns the mapping between function spaces rather than a discrete point-to-point mapping, endowing it with resolution invariance and strong cross-condition generalization. This is why FNO maintains a relative error below 5% even under a 10% parameter-extrapolation range, whereas PINN, which is essentially a condition-specific solver, degrades markedly when the operating condition deviates from the training distribution. The three-orders-of-magnitude inference speed-up further stems from the fact that an FNO forward pass is a single feed-forward computation, whereas numerical simulation requires iterative solution of the coupled PDE system.

3.6.3. Mechanistic Interpretation of the Coordinated Optimization Gain

The 28.5% comprehensive efficiency improvement of three-parameter optimization over single-parameter optimization reflects the exploitation of the intrinsic hydraulic–mechanical coupling among WOB, RPM, and flow rate. WOB and RPM jointly govern the rock-breaking rate and hence the cuttings generation rate, while the flow rate governs the cuttings-carrying capacity and the annular pressure loss; optimizing any one parameter in isolation shifts the operating point to a local, constrained optimum and frequently degrades the remaining indicators. The coordinated formulation explores the joint feasible region, and the probabilistic safety constraints additionally prevent the optimizer from recommending points whose true ECD may violate the safe window under inversion uncertainty, thereby achieving a robust balance among safety, efficiency, and energy consumption.

3.6.4. Limitations and Assumptions

The proposed method still has certain limitations. First, all experimental data in this study are generated by the in-house hydraulic PDE numerical simulation, without validation against field measured data; although the simulated conditions cover most conventional drilling scenarios, they cannot fully replicate complex geological disturbances, instrument anomalies, and sudden drilling-fluid performance changes, so the field adaptability of the method still requires further engineering verification. Second, the inverted parameters only focus on rheological parameters, cuttings concentration, and equivalent eccentricity, without considering special downhole complex conditions such as gas influx, wellbore collapse, and lost circulation, which limits the applicability under abnormal working conditions. Third, the FNO surrogate adopts a pure data-driven architecture; under extreme extrapolation conditions beyond the training parameter range, prediction accuracy exhibits slight decay, and the integration of physical constraints still has room for improvement. Fourth, the optimization decision mode relies on preset weighting and rule-based selection, without incorporating drilling-expert experience and real-time condition-aware intelligent decision-making, so the intelligence level can still be improved.

3.7. Future Works

To address the above limitations, future research can proceed along the following directions:

1.

Broader validation and engineering deployment. Collect field drilling data to fine-tune model parameters and constraint weights according to field working conditions, and complete engineering adaptation and deployment verification of the method across diverse geological basins.

2.

Extension of the inversion parameter dimensionality. Introduce variables such as gas-influx rate, lost-circulation volume, and wellbore-stability parameters to construct a multi-parameter intelligent inversion system for complex downhole conditions.

3.

Physics-enhanced surrogate modeling. Optimize the FNO architecture by integrating PDE physics residuals to construct a physics-enhanced FNO model, improving extrapolation prediction accuracy and physical rationality.

4.

Adaptive intelligent decision-making. Introduce reinforcement learning and expert systems to construct an adaptive intelligent decision-making mechanism, realizing fully automatic optimal parameter decisions under different formations and working conditions.

4. Conclusions

Aiming at the technical problems of low dynamic identification accuracy of hydraulic parameters, insufficient hydraulic-field prediction efficiency, isolated static optimization of drilling parameters, and the difficulty in balancing safety and efficiency during drilling in deep complex formations, this study proposes a WOB–RPM–flow-rate coordinated optimization method driven by dual-prior-constrained temporal inversion and an FNO surrogate model. A three-layer integrated while-drilling intelligent decision-making framework of “parameter inversion–hydraulic prediction–intelligent optimization” is constructed. The effectiveness of the method is verified through theoretical derivation and multiple groups of numerical simulation experiments. The main conclusions are summarized as follows:

1.

A dual-prior-constrained temporal inversion module is constructed, fusing the pre-drilling mechanism baseline prior and the offset-well statistical prior, combined with standpipe pressure–rotary torque dual-observation information and temporal smoothness constraints. A four-term weighted loss function is designed, effectively improving the ill-posedness of multi-parameter inversion. Compared with the traditional least-squares single-observation inversion method, the inversion RMSE is reduced by 42.7%, enabling accurate and stable identification of the dynamic variations of plastic viscosity, yield point, average annular cuttings concentration, and equivalent eccentricity, with excellent anti-noise performance and temporal stability.

2.

The Fourier neural operator is applied to eccentric rotating-annular two-phase-flow hydraulic prediction for the first time, and a one-dimensional FNO surrogate model adapted to the drilling-fluid hydraulics scenario is constructed. The model accurately maps the functional relationship from drilling operational parameters and rheological parameters to the full-field hydraulic response, with stable prediction accuracy maintained above 98.3%. The inference speed is improved by three orders of magnitude compared with traditional numerical simulation, providing efficient computational support for real-time online optimization.

3.

A hydraulic–mechanical coupled three-parameter multi-objective coordinated optimization framework is established, with drilling-efficiency maximization and energy-consumption and risk minimization as the core objectives. Combined with online Bayesian optimization, mechanism-gradient acceleration, and the sliding-window rolling update mechanism, dynamic coordinated optimization of WOB, RPM, and flow rate is realized. Probabilistic safety constraints based on the confidence intervals of inverted parameters effectively reduce optimization risks arising from parameter uncertainty. Compared with traditional single- and dual-parameter optimization schemes, the comprehensive drilling efficiency is improved by 28.5%, achieving a global optimal balance of drilling safety, efficiency, and energy consumption.

4.

Ablation experiments verify the independent value and coupling gain effect of each module in the three-layer framework: dual-prior constraints guarantee parameter-identification accuracy, the FNO surrogate model provides efficient forward-modeling support, and online probabilistic optimization achieves dynamic robust decision-making. The synergistic interaction of the three modules jointly constitutes the core advantage of the proposed method, which can provide a novel theoretical method and technical support for real-time intelligent regulation of safe and efficient drilling in deep and ultra-deep wells with narrow windows.

In the future, field drilling data can be combined to complete model engineering adaptation, the capability of identifying complex downhole abnormal working conditions can be expanded, model performance can be optimized by integrating physical constraints and intelligent algorithms, and the field deployment and application of drilling-fluid hydraulic intelligent inversion and parameter optimization technology can be further promoted.