基于最小二乘的量化评估页岩气压裂改造效果方法

Methods to Quantitatively Evaluate the Effect of Shale Gas Fracturing Stimulation Based on Least Squares

摘要

在页岩气压裂改造过程中，水平井分段多簇压裂技术是实现页岩气高效开发的关键技术之一，但这项技术缺乏低成本、可量化的评价方法。页岩气压裂改造效果的好坏通常通过井下射孔的有效开启个数来判断，有效孔眼开启得多，则压裂改造效果好，反之则差。针对目前页岩气压裂改造技术评估的瓶颈，本文设计一种基于非线性最小二乘的量化评估模型，该模型结合井下各项摩阻理论与孔眼扩展方程，实现了有效开启射孔数的高效计算。首先，根据压裂过程中阶梯降排量测试中的总摩阻（井筒摩阻、孔眼摩阻、近井摩阻之和）与压裂液排量的对应关系，构造了页岩气摩阻拟合的非线性最小二乘目标函数，从而拟合出各项摩阻系数，尤其是孔眼摩阻系数；然后，结合孔眼摩阻方程和孔眼扩展方程，创建了射孔有效开启孔数和磨蚀后射孔平均直径的计算模型；最后，将该模型应用于工业生产，根据有效开启射孔数和磨蚀后射孔平均直径的计算结果，得到该井段水平井分段多簇压裂改造的实施效果评估依据。在摩阻拟合方面，相比原有模型，本文模型具有可区分所有拟合参数的优点，并且添加了符合工程实际的约束条件，使得拟合结果符合工程需求。所提出的评估方法具有数学理论完备的优点，且实施成本低、效率高，可作为评价压裂改造的新方法。

Abstract

Objective Shale gas fracturing constitutes the cornerstone of contemporary energy extraction, with horizontal well-staged multi-cluster fracturing technology emerging as a pivotal technique for achieving efficient shale gas development. Despite its critical role in meeting global energy demands, the industry faces a persistent challenge: the lack of cost-effective, quantifiable methodologies to assess the effectiveness of fracturing. Conventional approaches rely predominantly on post-fracturing active perforation counts as a proxy for stimulation effectiveness. Although empirical evidence indicates a positive correlation between the number of active perforations and production enhancement, this oversimplified metric fails to capture the inherent complexity of hydraulic fracturing dynamics. The process involves complex interactions among geological formations, fluid rheology, wellbore configurations, and operational parameters. Practical limitations at field sites, where only total friction and flow rate are measurable, further compel engineers to neglect quantitative analysis of individual friction components (wellbore friction vs. perforation friction). Hence, traditional methods rely heavily on subjective experiential judgment, resulting in compromised accuracy in active perforation assessment and suboptimal fracturing design. This study develops a comprehensive, physics-based quantitative model to accurately evaluate the effectiveness of fracturing stimulation, enabling data-driven optimization of shale gas extraction processes. Methods This study established a novel quantitative evaluation model based on principles of fluid mechanics and mathematical optimization theory. The methodology utilized nonlinear least squares optimization and proceeded through two integrated computational phases: 1) Friction coefficient fitting: A dedicated nonlinear least squares objective function was constructed to resolve friction components during staged multi-cluster fracturing. Using data obtained from step-down discharge tests, the model analyzed the relationship between total friction (comprising wellbore friction, perforation friction, and near-wellbore friction) and the fracturing fluid flow rate. The optimization targeted the perforation friction coefficient k_perf as the primary unknown variable. Engineering-informed constraints, such as realistic friction ranges and fluid behavior boundaries, were incorporated to ensure physically meaningful solutions. Advanced global optimization algorithms, including the Levenberg-Marquardt method, were applied to effectively address this non-convex problem. 2) Active perforation assessment: A computational method was developed based on the fitted perforation friction coefficient and established perforation erosion equations. This method identified the precise number of active perforations and calculated their average diameter after erosion during the fracturing process. Advanced optimization algorithms were employed to efficiently address both the fitting and calculation tasks. The model integrated comprehensive fluid mechanics theories related to downhole friction with existing perforation erosion models. Results and Discussions The model was deployed in shale gas wells across the Sichuan Basin in China. During the step-down tests, a high level of agreement was observed between the predicted and measured friction pressure curves, confirming the model's robustness under complex field conditions. It delivered accurate quantitative outputs for both the number of active perforations and their average diameter after abrasion, overcoming the subjectivity and inaccuracy associated with traditional methods. The proposed model demonstrated several advantages over conventional approaches: 1) Enhanced accuracy and relevance: The model ensured highly accurate and practically applicable results by carefully defining fitting parameters and incorporating engineering constraints. 2) Robust theoretical foundation: It was grounded in established mathematical theory and principles of fluid mechanics. 3) Practicality and efficiency: The model featured low implementation costs and high computational efficiency, presenting a viable and promising alternative for field evaluations. Its significance to the industry lay in addressing the major challenge of the absence of a cost-effective, quantifiable assessment method. It offered detailed insights into the effectiveness of fracturing (in terms of the number and quality of perforations), enabling engineers to improve the fracturing stimulation process for improved production results. As the demand for shale gas continued to increase, this innovative approach proved critical to improving industry efficiency and sustainability. The model established a foundation for future developments, supporting more efficient and environmentally sustainable shale gas extraction practices by enabling a deeper understanding of the fracturing process and its results.

Graphical abstract

关键词

页岩气开采 / 水平井分段多簇压裂 / 摩阻拟合 / 非线性最小二乘 / 有效开启孔数

Key words

shale gas extraction / horizontal well staged multi-cluster fracturing / friction fitting / nonlinear least squares / active perforations

本刊网刊

中国四川盆地富含页岩气资源，其天然气产量正在快速增长，对中国石油和天然气行业发挥着重要作用^[1]。 然而，页岩储层通常极为致密，这些储层的物理性质差、渗透性低且流动性差，在短期内降低开发成本仍然是一个巨大的挑战。围绕页岩油气生产过程中各个生产环节的基于不同目标的量化建模^[2‒3]与数值计算^[4‒5]具有重要意义，特别地，压裂阶段的优化建模结果可以指导生产并降低成本。水平井分段多簇压裂技术是实现页岩气高效开发的关键技术之一^[6‒7]，它可以利用水力喷射、桥塞射孔联作、裸眼滑套等方法，在相对较短的时间内完成多个储层的压裂作业，并最大限度地减少对储层的伤害，实现多层合采，形成复杂的裂缝网格，从而提高储层改造效果和单井产量，最大限度地提高地质储层的可利用性。

在水力压裂过程中，存在一些常见问题，例如：井筒中的套管变形、工具入井困难、无法实施分段多簇射孔和下入桥塞等，这些问题极大降低了压裂施工的效率^[8‒9]。同时，受射孔孔眼进液不均匀、多个裂缝相互干扰、地应力条件及储层非均质等条件的影响^[10]，压裂施工面临一次性改造射孔簇开启不充分、各簇裂缝同步和均匀扩展困难、水平井段改造不完全等关键工程问题^[11]。为了有效提高分段多簇压裂改造效果并提高油气产量，近年来发展了压裂暂堵转向技术，该技术利用泵注暂堵剂进入井底，以阻塞液体主要流通通道并增加井底压力，从而促进新裂缝的形成和扩展。作业结束后，暂堵剂完全溶解，解除封堵^[12‒14]。这项技术可以实现对原有裂缝的改造和新裂缝的生成，从而有效提高射孔孔眼的开启效率，并增加储层的有效渗透率和产能^[15]。

尽管暂堵转向技术可以有效提升页岩改造效果，包括最大化储层改造体积、裂缝网格复杂性和改造均匀性，但目前暂堵转向施工主要依据现场施工经验，对于暂堵前后射孔开启情况缺乏定量认识，缺少系统、量化的设计和评估方法，这限制了暂堵转向施工的设计和压裂改造效果^[16]。为识别有效开启的射孔孔眼，目前已经发展了井下电视^[17]、钻井取心观察^[18]等方法，但实施成本高，施工难度大。因此，如何以低成本获得井下射孔孔眼开启识别方法是当前研究的一个重要方向。考虑射孔摩阻与有效开启孔数的关系，以及有效开启孔数对水平井分段多簇压裂施工效果的影响，有必要建立页岩气摩阻拟合模型，以更快、更好地计算磨蚀后的射孔数及平均射孔直径。

为此，本文建立一个新模型计算有效开启孔数及磨蚀后的平均射孔直径。本文的主要贡献是：1）基于总摩阻的构成及各项摩阻的基础理论，建立页岩气摩阻拟合的非线性最小二乘模型；2）结合孔眼摩阻的达西公式及射孔随时间磨蚀的趋势（孔眼扩展方程），利用第1）步得到的孔眼摩阻系数，建立孔眼有效开启孔数及磨蚀后平均孔眼直径的计算模型；3）利用射孔孔眼有效开启孔数的计算结果判断当前的压裂效果，以优化和调整非常规油气储层的开发和生产，为下一步暂堵转向技术的实施提供依据，从而增加油气产量和产能。将本文的方法应用于四川盆地页岩气开采，结果显示，效果明显优于原有方案，可协助工程师在施工现场迅速评估下一步压裂方案的优劣，从而优化并调整后续页岩气的开采过程。

1 水平井压裂摩阻构成与拟合模型

1.1 分段压裂水平井中的各项摩阻

水平井中总摩阻的构成为^[19‒21]：

P_{\mathrm{f}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=P_{\mathrm{s}}-P_{\mathrm{0}}=P_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}+P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}+P_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}（1）

式中，P_{\mathrm{s}}为施工总压力，P_{\mathrm{0}}为停泵压力，P_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}为井筒摩阻，P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}为孔眼摩阻，P_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}为近井摩阻，单位均为MPa。

井筒摩阻的计算公式来源于流体力学中的达西公式^[22]：

P_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=\lambda \frac{l}{\mathrm{2}d}\rho {\left(\frac{\mathrm{4}q}{\mathrm{\pi }{d}^{\mathrm{2}}}\right)}^{\mathrm{2}}=k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}^{\mathrm{2}}（2）

式中：\lambda定量地描述了流体流动状态（雷诺数）和管道内壁相对粗糙程度共同作用对摩擦阻力的影响程度，它把流体力学中复杂的流动特性引入到了工程计算中；l为管道长度；d为井筒直径；\rho为压裂液密度；q为压裂液流量（即后文简称的排量）；k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}为系数。式（2）中第2个等式k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}^{\mathrm{2}}将达西公式高度简化，极大地便利了现场工程计算，k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}的物理根源即为包含了\lambda \mathrm{、}\rho \mathrm{、}l\mathrm{、}d的复杂组合。在实际生产中，受压力、温度等情况的影响，实际的井筒摩阻情况较为复杂，不易定量计算^[23]。对于工程上使用的五寸半套管，参考行业中清水的千米摩阻与排量的对应关系，见表1。

表1表示当清水排量确定时，五寸半套管每千米所具有的摩阻。工程中，通常控制降阻剂的用量使得液体的摩阻满足需求，降阻率一般为70%～80%^[24]，因此，实际施工中通过降阻剂后的摩阻是清水的20%～30%。将裂缝井段的深度范围乘以降阻后的千米摩阻，即可得到井筒摩阻的大致范围。

射孔孔眼摩阻的计算公式为^[25]：

P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\frac{\alpha \rho q^{\mathrm{2}}}{n_{\mathrm{p}}^{\mathrm{2}}{d}_{\mathrm{p}}^{\mathrm{4}}{C}_{\mathrm{d}}^{\mathrm{2}}}（3）

式中，\alpha为经验参数，n_{\mathrm{p}}为孔眼开启数，d_{\mathrm{p}}为磨蚀后的平均孔眼直径，C_{\mathrm{d}}为孔眼流动系数^[26‒27]。令系数k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\alpha \rho /\left(n_{\mathrm{p}}^{\mathrm{2}}{d}_{\mathrm{p}}^{\mathrm{4}}{C}_{\mathrm{d}}^{\mathrm{2}}\right)，则孔眼摩阻的计算公式为：

P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}^{\mathrm{2}}（4）

对于近井摩阻的计算，被广泛应用的摩阻与排量的幂律关系有^[28]

P_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}^{\beta }（5）

式中：\beta通常取0.5；系数k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}定量地表征了特定井在当前状态下（射孔孔眼数量、尺寸、效率、清洁度、近井地层特性等）产生的近井流动阻力强度，无法仅靠理论精确计算，而必须通过现场压裂测试来拟合确定。

1.2 工程上常用的摩阻拟合模型

结合方程（1）、（2）、（4）、（5）与多次降排量测试得到的排量与其对应总摩阻的值，可以得到一个关于3个摩阻系数（k_well、k_perf、k_nearwell）的最小二乘模型（\beta =\mathrm{0.5}）^[29]：

\mathrm{m}\mathrm{i}\mathrm{n}\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^m}{f_i}^{\mathrm{2}},f_i=k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\mathrm{2}}+k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}_i^{\mathrm{2}}+k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\mathrm{0.5}}-P_i（6）

式中，m为阶梯降排量测试的台阶数，q_i为阶梯排量，P_i为对应的摩阻。

观察优化式（6），结合线性最小二乘算法，需拟合的参数k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}与k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}的系数均为q_i^{\mathrm{2}}，故该优化模型只能求出使式（6）最小的k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}+k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}，而不能将二者区分；即使根据表1和井深范围可以得到井筒摩阻系数k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}的一个大致范围，但对应的孔眼摩阻系数也是一个范围，在这个范围内的解都是式（6）的最优解。因此，需要对其做进一步改进。

2 改进的摩阻拟合的非线性最小二乘模型

首先，在页岩气开采过程中，井下情况十分复杂，射孔孔眼摩阻的流动指数事实上可能并不是准确的值（即2），因此，本文决定将孔眼摩阻流动指数作松弛处理，即将其也作为需要拟合的参数\gamma，由此，方程（4）变为：

P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}^{\gamma },\gamma \in [\mathrm{1.5,2.5}]（7）

其次，对近井筒裂缝几何形状不同假设的理论推导表明，实际的\beta取值在0.25～1.00之间^[28]。

事实上，对于在固定宽度平行板间的流动，压降P（摩阻）与排量q之间的关系为^[30]：

P=\frac{{\mathrm{\pi }}^{\mathrm{2}}\mu }{{\omega }^{\mathrm{3}}}q（8）

式中，\mu为压裂液黏度，\omega为平行板模型的固定宽度。由式（8）可知，此时q的指数为1，对应于式（5）有\beta =\mathrm{1}。

在裂缝宽度取决于流体压力情况下的流动中，PKN（Perkins‒Kern‒Nordgren）模型给出了另外一个公式^[31‒33]：

q=\frac{PH\mathrm{16}{L}^{\mathrm{2}}{\left(p_{\mathrm{f}}+\frac{P}{\mathrm{2}}-{\sigma }_{\mathrm{h}\mathrm{m}\mathrm{i}\mathrm{n}}\right)}^{\mathrm{3}}}{\mathrm{3}\mu {E^{\text{'}}}^{\mathrm{3}}}（9）

式中：H和L分别为裂缝高度和长度；p_f为裂缝内的流体压力；{\sigma }_{\mathrm{h}\mathrm{m}\mathrm{i}\mathrm{n}} 为最小原地应力；E^{\text{'}}=E{\left(\mathrm{1}-\nu \right)}^{\mathrm{2}}，E和\nu分别为杨氏模量和泊松比。在这种情况下，排量与摩阻的关系为q\propto P^{\mathrm{4}}+P^{\mathrm{3}}+P^{\mathrm{2}}+P， \beta \in (\mathrm{0.25,1.00})。

如果近井筒裂缝的几何形状是固定宽度的曲线而不是平行板，那么有^[26]

P=\frac{{\mathrm{\pi }}^{\mathrm{2}}\mu }{{\omega }^{\mathrm{3}}}\times R（10）

式中，R=\theta \sqrt[]{\frac{E^{\mathrm{3}}\mu q}{H}}{\left(\frac{\mathrm{1}}{{\sigma }_{\mathrm{h}\mathrm{m}\mathrm{i}\mathrm{n}}\left(\frac{{\sigma }_{\mathrm{H}\mathrm{m}\mathrm{a}\mathrm{x}}}{{\sigma }_{\mathrm{h}\mathrm{m}\mathrm{i}\mathrm{n}}}-\mathrm{1}\right)}\right)}^{\mathrm{2}}，\theta为试验系数，{\sigma }_{\mathrm{H}\mathrm{m}\mathrm{a}\mathrm{x}}为最大水平应力。此时，对应于式（5），有\beta =\mathrm{0.5}。

在实际工程中，近井筒压裂的几何形状可能更加复杂，并且具有时间依赖性，恒定的某个\beta不能很好地表示流动类型或近井筒的几何形状。故在本文模型中，将\beta也作为未知参数进行拟合，范围为\left[\mathrm{0.25,1.00}\right]。

由于每个未知参数表示一定的物理意义，并与实际情况紧密相关，故本文决定对每个拟合参数进行相应的约束。

首先，根据第1.1节，将实际施工中的井深范围与表1对比可得到井筒摩阻的大致范围，再由井筒摩阻公式反推得到井筒摩阻系数的约束区间。例如：该施工井段的井深范围为a～b（单位\mathrm{k}\mathrm{m}），稳定的施工排量为q_{\mathrm{1}}（单位{\mathrm{m}}^{\mathrm{3}}\cdot \mathrm{m}\mathrm{i}{\mathrm{n}}^{-\mathrm{1}}），查阅表1可得对应的五寸半套管的千米摩阻为p_{\mathrm{1}}（单位\mathrm{M}\mathrm{P}\mathrm{a}）；根据实际施工过程中的降阻率，可以得到对应的井筒摩阻范围大致为\mathrm{0.2}a{p}_{\mathrm{1}}～\mathrm{0.3}b{p}_{\mathrm{1}}（单位\mathrm{M}\mathrm{P}\mathrm{a}）；再由式（2）可以反推得到井筒摩阻系数的约束区间为：

k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}\in \left[\frac{\mathrm{0.2}a{p}_{\mathrm{1}}}{q_{\mathrm{1}}^{\mathrm{2}}},\frac{\mathrm{0.3}b{p}_{\mathrm{1}}}{q_{\mathrm{1}}^{\mathrm{2}}}\right]（11）

对于孔眼摩阻系数和近井摩阻系数，并没有明确的范围限制，首先保证是正实数，其次在工程上主要通过控制井筒摩阻的范围，结合施工中总摩阻的大小，得到井筒摩阻大于孔眼摩阻与近井摩阻之和的关系，可用数学公式表示为：

k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\mathrm{2}}\ge k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}_i^{\gamma }+k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\beta },i=\mathrm{1,2},\cdots ,m（12）

因此，最终改进的优化问题为：

\begin{array}{c}\mathrm{m}\mathrm{i}\mathrm{n}\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^m}{f_i}^{\mathrm{2}},\\ f_i=k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\mathrm{2}}+k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}_i^{\gamma }+k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\beta }-P_i,\\ \mathrm{s}.\mathrm{t}.\frac{\mathrm{0.2}a{p}_{\mathrm{1}}}{q_{\mathrm{1}}^{\mathrm{2}}}\le k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}\le \frac{\mathrm{0.3}b{p}_{\mathrm{1}}}{q_{\mathrm{1}}^{\mathrm{2}}},\\ \mathrm{1.5}\le \gamma \le \mathrm{2.5},\\ \mathrm{0.25}\le \beta \le \mathrm{1},\\ k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}-k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}_i^{\gamma -\mathrm{2}}-k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\beta -\mathrm{2}}\ge \mathrm{0}\end{array}（13）

对于式（13），定义\mathit{x}=[k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}},k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}},\gamma ,k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{,\beta ]}^{\mathrm{T}}，且有函数 f_i(\mathit{x})=k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\mathrm{2}}+k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{q}_i^{\gamma }+k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{q}_i^{\beta }-P_i,i=\mathrm{1,2},\cdots ,m，则可将式（13）转化为矩阵形式，即为求解使式（14）最小的\mathit{x}：

F(\mathit{x})=\frac{\mathrm{1}}{\mathrm{2}}\mathit{f}{(\mathit{x})}^{\mathrm{T}}\mathit{f}(\mathit{x})（14）

式中，\mathit{f}(\mathit{x})=[f_{\mathrm{1}}(\mathit{x}),f_{\mathrm{2}}(\mathit{x}),\cdots ,f_m{(\mathit{x})]}^{\mathrm{T}}。

根据非线性目标函数的特性，采用梯度下降法进行迭代求解。故需要找到使目标函数F(\mathit{x})下降的方向和步长。考虑在\mathit{x}处对\mathit{f}进行泰勒展开：

\begin{array}{c}\mathit{f}(\mathit{x}+\mathit{h})=\mathit{f}(\mathit{x})+\mathit{J}(\mathit{x})\mathit{h}+O({‖\mathit{h}‖}^{\mathrm{2}})\approx \\ \mathit{f}(\mathit{x})+\mathit{J}(\mathit{x})\mathit{h}=\mathit{l}(\mathit{h})\end{array}（15）

式中，\mathit{J}(\mathit{x})为\mathit{f}(\mathit{x})的Jacobi矩阵，\mathit{J}\in {\mathbb{R}}^{m\times m}，它的第i行、第j列的分量表示为\mathit{J}{(\mathit{x})}_{ij}=\frac{\partial f_i(\mathit{x})}{\partial x_j}，其中，x_j为\mathit{x}的第j个分量。\mathit{J}(\mathit{x})具体表示为\mathit{J}=[{\mathit{J}}_{\mathrm{1}},{\mathit{J}}_{\mathrm{2}},{\mathit{J}}_{\mathrm{3}},{\mathit{J}}_{\mathrm{4}},{\mathit{J}}_{\mathrm{5}}]，其中：

\begin{array}{c}{\mathit{J}}_{\mathrm{1}}={\mathit{Q}}^{\mathrm{2}},{\mathit{J}}_{\mathrm{2}}={\mathit{Q}}^{\gamma },{\mathit{J}}_{\mathrm{3}}=k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{\mathit{Q}}^{\gamma }\mathrm{l}\mathrm{n}\mathit{Q},{\mathit{J}}_{\mathrm{4}}={\mathit{Q}}^{\beta },\\ {\mathit{J}}_{\mathrm{5}}=k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{\mathit{Q}}^{\beta }\mathrm{l}\mathrm{n}\mathit{Q}\end{array}（16）

式中，\mathit{Q}=[q_i]为列向量，分别表示降排量测试中的总摩阻及对应的排量。

从而，对目标函数，有

\begin{array}{l}F(\mathit{x}+\mathit{h})\approx L(\mathit{h})=\frac{\mathrm{1}}{\mathrm{2}}\mathit{l}{(\mathit{h})}^{\mathrm{T}}\mathit{l}(\mathit{h})=\\ \frac{\mathrm{1}}{\mathrm{2}}{\mathit{f}}^{\mathrm{T}}\mathit{f}+{\mathit{h}}^{\mathrm{T}}{\mathit{J}}^{\mathrm{T}}\mathit{f}+\frac{\mathrm{1}}{\mathrm{2}}{\mathit{h}}^{\mathrm{T}}{\mathit{J}}^{\mathrm{T}}\mathit{J}\mathit{h}=\\ F(\mathit{x})+{\mathit{h}}^{\mathrm{T}}{\mathit{J}}^{\mathrm{T}}\mathit{f}+\frac{\mathrm{1}}{\mathrm{2}}{\mathit{h}}^{\mathrm{T}}{\mathit{J}}^{\mathrm{T}}\mathit{J}\mathit{h}\end{array}（17）

高斯牛顿法的步长{\mathit{h}}_{\mathbf{g}\mathbf{n}}表示为{\mathit{h}}_{\mathbf{g}\mathbf{n}}=\underset{\mathit{h}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}L(\mathit{h})，则对L(\mathit{h})求导，有\mathit{L}\mathrm{\text{'}}(\mathit{h})={\mathit{J}}^{\mathrm{T}}\mathit{f}+{\mathit{J}}^{\mathrm{T}}\mathit{J}\mathit{h},\mathit{L}″(\mathit{h})={\mathit{J}}^{\mathrm{T}}\mathit{J}。

令\mathit{L}\mathrm{\text{'}}(\mathit{h})=\mathrm{0}，有

{\mathit{J}}^{\mathrm{T}}\mathit{J}{\mathit{h}}_{\mathbf{g}\mathbf{n}}=-{\mathit{J}}^{\mathrm{T}}\mathit{f}（18）

由于{\mathit{J}}^{\mathrm{T}}\mathit{J}只能保证半正定，若{\mathit{J}}^{\mathrm{T}}\mathit{J}为奇异矩阵，则式（18）无解，此时高斯牛顿法失效，故在高斯牛顿法的基础上增加了一项阻尼项，称为阻尼‒高斯牛顿法 (damping‒Gauss‒Newton method)，中文称为列文伯格‒马夸尔特方法（Levenberg‒Marquardt method，LM）^[34]：

\left({\mathit{J}}^{\mathrm{T}}\mathit{J}+\kappa \mathit{I}\right){\mathit{h}}_{\mathbf{l}\mathbf{m}}=-{\mathit{J}}^{\mathrm{T}}\mathit{f},\mu \ge \mathrm{0}（19）

由式（19）可知，若LM方法的步长κ足够大，则{\mathit{J}}^{\mathrm{T}}\mathit{J}+\kappa \mathit{I}必为正定矩阵，式（19）一定有解。

如果目标函数F(\mathit{x})直接对x_j求导，则有\frac{\partial F(\mathit{x})}{\partial x_j}={\displaystyle \sum _{i=\mathrm{1}}^m}{f}_j(\mathit{x})\frac{\partial f_i(\mathit{x})}{\partial x_j}，由此，F(\mathit{x})的梯度可以表示为：

\mathit{g}=\nabla F(\mathit{x})=\mathit{J}{(\mathit{x})}^{\mathrm{T}}\mathit{f}(\mathit{x})（20）

由以上方法得到的\mathit{h}只能保证L(\mathit{h})下降，而L(\mathit{h})只是F(\mathit{x}+\mathit{h})的近似，并不能保证F(\mathit{x}+\mathit{h})下降，因此引入一个参数\xi来衡量下降的质量：

\xi =\frac{F(\mathit{x})-F(\mathit{x}+\mathit{h})}{L(\mathrm{0})-L(\mathit{h})}（21）

此方法迭代优化的停止条件如下。

1） 若{\mathit{x}}^{\mathrm{*}}为局部最小值点，则\nabla F(\mathit{x})=\mathit{g}({\mathit{x}}^{\mathrm{*}})=\mathrm{0}。由此可以采用如下停止条件：

{‖\mathit{g}‖}_{\mathrm{\infty }}\le {ϵ}_{\mathrm{1}}（22）

式中，{ϵ}_{\mathrm{1}}为充分小的正数，使得存在\mathit{x}使‖\mathit{g}‖趋近于0。

2） 若迭代前后两次的\mathit{x}变化很小，也可认为达到了局部最小值点，即：

‖{\mathit{x}}_{\mathrm{n}\mathrm{e}\mathrm{w}}-\mathit{x}‖\le {ϵ}_{\mathrm{2}}(‖\mathit{x}‖+{ϵ}_{\mathrm{2}})（23）

式中，{ϵ}_{\mathrm{2}}为充分小的正数， x_new为更新后的值。

3） 设置最大迭代次数：

（24）

3 有效开启孔数与磨蚀后孔眼平均直径的计算模型

实际工程中孔眼摩阻系数的计算公式为^[35]：

{\tilde{k}}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\frac{\mathrm{0.236}\mathrm{9}\tilde{\rho }}{n_{\mathrm{p}}^{\mathrm{2}}{\tilde{d}}_{\mathrm{p}}^{\mathrm{4}}{C}_{\mathrm{d}}^{\mathrm{2}}}（25）

式中，{\tilde{k}}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}表示孔眼摩阻系数（\mathrm{p}\mathrm{s}\mathrm{i}\cdot \mathrm{m}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}/\mathrm{g}\mathrm{a}\mathrm{l})，\tilde{\rho }为压裂液密度（\mathrm{l}\mathrm{b}/\mathrm{g}\mathrm{a}\mathrm{l}），{\tilde{d}}_{\mathrm{p}}为平均孔眼直径（\mathrm{i}\mathrm{n}），n_{\mathrm{p}}为有效开启孔数，C_{\mathrm{d}}为孔眼流量系数。

将式（25）中的所有参数换为本文所用单位，有^[27]

k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\frac{\mathrm{2.232}\mathrm{6}\times {\mathrm{10}}^{-\mathrm{7}}\rho }{n_{\mathrm{p}}^{\mathrm{2}}{d}_{\mathrm{p}}^{\mathrm{4}}{C}_{\mathrm{d}}^{\mathrm{2}}}（26）

式中，各项参数的单位分别为k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}（\mathrm{M}\mathrm{P}\mathrm{a}\cdot \mathrm{m}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}/{\mathrm{m}}^{\mathrm{6}}）、\rho（\mathrm{g}/\mathrm{c}{\mathrm{m}}^{\mathrm{3}}）、d_{\mathrm{p}}（\mathrm{m}）。

式（25）、（26）中的C_{\mathrm{d}}会随着磨蚀的加深而增大，其变化规律为^[36]：

\frac{\partial C_{\mathrm{d}}}{\partial t}=\zeta c{\left(\frac{\mathrm{4}q}{n_{\mathrm{p}}\mathrm{\pi }{d}_{\mathrm{p}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\left(\mathrm{1}-\frac{C_{\mathrm{d}}}{C_{\mathrm{d}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}\right)（27）

式中：\zeta为根据实际情况由经验获取的系数；c为支撑剂浓度，\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}；C_{\mathrm{d}}的范围为\mathrm{0.6}0～\mathrm{0.95}，C_{\mathrm{d}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{0.95}；t为时间。通常在压裂后，认为有效开启孔数的射孔均被完全磨蚀，即直接取C_{\mathrm{d}}=C_{\mathrm{d}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{0.95}。

而孔眼扩展方程为^[37]：

\frac{\partial d_{\mathrm{p}}}{\partial t}=\eta c{\left(\frac{\mathrm{4}q}{n_{\mathrm{p}}\mathrm{\pi }{d}_{\mathrm{p}}^{\mathrm{2}}}\right)}^{\mathrm{2}}（28）

式中，\eta与式（27）中的\zeta一样，为经验系数。

式（28）同时对时间积分，利用初始条件d_{\mathrm{p}}{|}_{t=\mathrm{0}}=d_{\mathrm{0}}（d₀为已知的平均初始直径），有

d_{\mathrm{p}}^{\mathrm{5}}=\eta \frac{\mathrm{80}c{q}^{\mathrm{2}}t}{n_{\mathrm{p}}^{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}}+d_{\mathrm{0}}^{\mathrm{5}}（29）

结合式（26）、（29），有

n_{\mathrm{p}}^{\mathrm{2}}=\frac{\alpha \rho }{k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{d}_{\mathrm{p}}^{\mathrm{4}}{C}_{\mathrm{d}}^{\mathrm{2}}}=\frac{\mathrm{80}\eta c{q}^{\mathrm{2}}t}{(d_{\mathrm{p}}^{\mathrm{5}}-d_{\mathrm{0}}^{\mathrm{5}}){\mathrm{\pi }}^{\mathrm{2}}}（30）

利用式（30）的后两式，得到关于d_{\mathrm{p}}的方程：

d_{\mathrm{p}}^{\mathrm{5}}-\frac{\mathrm{80}\eta c{q}^{\mathrm{2}}t{k}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{C}_{\mathrm{d}}^{\mathrm{2}}}{\alpha \rho {\mathrm{\pi }}^{\mathrm{2}}}{d}_{\mathrm{p}}^{\mathrm{4}}-d_{\mathrm{0}}^{\mathrm{5}}=\mathrm{0}（31）

求解这个关于d_{\mathrm{p}}的5次方程，即可得到磨蚀后的平均孔眼直径；将d_{\mathrm{p}}的值代入式（30）中即可求出有效开启孔数n_{\mathrm{p}}。

由伽罗瓦理论可知，一般的5次方程没有求根公式，故采用牛顿迭代法来求解式（31）。牛顿迭代法的最大优点在于它在方程f(x)=\mathrm{0}的单根附近具有平方收敛，而且该法还可以用来求解方程的重根、复根，此时线性收敛，且可以通过一些方法变成超线性收敛。

定　理　Newton‒Raphson法的局部收敛性定理^[38]。若目标函数f，2阶连续可微且其Hessian矩阵{\nabla }^{\mathrm{2}}f(x)在满足充分条件的最优值点x^{\mathrm{*}}的邻域内Lipschitz连续，考虑迭代x_{k+\mathrm{1}}=x_k+p_k，其中p_k=-{\nabla }^{\mathrm{2}}f(x_k{)}_{}^{-\mathrm{1}}\nabla f(x_k)。则

1） 如果迭代初值x_{\mathrm{0}}足够接近x^{\mathrm{*}}，则牛顿法产生的迭代点列\left\{x_i\right\}收敛到x^{\mathrm{*}}；

2） \{x_k\} 2次收敛到x^{\mathrm{*}}；

3） \{‖\nabla f(x_k)‖\} 2次收敛到0。

式（31）中，d_{\mathrm{p}}>d_{\mathrm{0}}，因此，迭代初值取x_{\mathrm{0}}=d_{\mathrm{0}}。令f(x)=x^{\mathrm{5}}-\frac{\mathrm{80}\eta c{q}^{\mathrm{2}}t{k}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{C}_{\mathrm{d}}^{\mathrm{2}}}{\alpha \rho {\mathrm{\pi }}^{\mathrm{2}}}{x}^{\mathrm{4}}-d_{\mathrm{0}}^{\mathrm{5}}，则f\mathrm{\text{'}}(x)=\mathrm{5}{x}^{\mathrm{4}}-\frac{\mathrm{320}\eta c{q}^{\mathrm{2}}t{k}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}{C}_{\mathrm{d}}^{\mathrm{2}}}{\alpha \rho {\mathrm{\pi }}^{\mathrm{2}}}x^{\mathrm{3}}，f(x)在x_n处的1阶泰勒展开为f(x)\approx f(x_n)+f\text{'}(x_n)(x-x_n),令右端等于0，有x=x_n-\frac{f(x_n)}{f\mathrm{\text{'}}(x_n)}，写为迭代格式为x_{n+\mathrm{1}}=x_n-\frac{f(x_n)}{f\mathrm{\text{'}}(x_n)}。将上述算法的结果作为d_{\mathrm{p}}的近似d_{\mathrm{p}}^{\mathrm{*}}，代入式（30），有

n_p=\sqrt[]{\frac{\alpha \rho }{k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}}}{d}_{\mathrm{p}}^{\mathrm{*}\mathrm{2}}{C}_{\mathrm{d}}（32）

因此，评价压裂效果的算法如下所示。

算　法　基于最小二乘的量化评估页岩气压裂改造效果的算法

输入：随机初值采样次数N，最大迭代次数M，\mathit{Q}=\left[q_i\right]，\mathit{P}=\left[P_i\right]，\left[a,b\right]，c，t，\rho =\mathrm{1}，C_{\mathrm{d}}=\mathrm{0.95}，\alpha =\mathrm{2.232}\mathrm{6}\times{\mathrm{10}}^{-\mathrm{7}}， \eta =\mathrm{1.01}\times {\mathrm{10}}^{-\mathrm{13}}，{ϵ}_{\mathrm{1}}={ϵ}_{\mathrm{2}}={\mathrm{10}}^{-\mathrm{10}}，i≔\mathrm{0}，k≔\mathrm{0}，\nu ≔\mathrm{2}，n≔\mathrm{0}，x_{\mathrm{0}}=d_{\mathrm{0}}。

输出：k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}，k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，\gamma，\beta，P_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}，P_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，n_{\mathrm{p}}， d_{\mathrm{p}}。

1. while do

2. i=i+\mathrm{1}

3. 由表1和式（11）计算井筒摩阻系数的约束区间

4. \begin{array}{l}\forall k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{0}}\in \left[\frac{\mathrm{0.2}a{p}_{\mathrm{1}}}{q_{\mathrm{1}}^{\mathrm{2}}},\frac{\mathrm{0.3}b{p}_{\mathrm{1}}}{q_{\mathrm{1}}^{\mathrm{2}}}\right],{\gamma }^{\mathrm{0}}\in \left[\mathrm{1.5,2.5}\right],\\ {\beta }^{\mathrm{0}}\in \left[\mathrm{0.25,1.00}\right],\end{array}

5. \mathit{x}≔\left[k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{0}},k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}^{\mathrm{0}},{\gamma }^{\mathrm{0}},k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{0}},{\beta }^{\mathrm{0}}\right],

\mathit{A}=\mathit{J}{(\mathit{x})}^{\mathrm{T}}\mathit{J}(\mathit{x}),\mathit{g}≔\mathit{J}{(\mathit{x})}^{\mathrm{T}}\mathit{f}(\mathit{x})

6. found≔\left({‖\mathit{g}‖}_{\mathrm{\infty }}\le {ϵ}_{\mathrm{1}}\right),\mu ≔{\mathrm{10}}^{-\mathrm{8}}\cdot \mathrm{m}\mathrm{a}\mathrm{x}{a}_{jj}

7. While notfound and do

8. k=k+\mathrm{1},\mathrm{解方}\mathrm{程}(\mathit{A}+\mu \mathit{I})\mathit{h}=-\mathit{g}

9. if ‖\mathit{h}‖\le {ϵ}_{\mathrm{2}}(‖\mathit{x}‖+{ϵ}_{\mathrm{2}}) then

10. found≔true

11. else

12. {\mathit{x}}_{\mathbf{n}\mathbf{e}\mathbf{w}}≔\mathit{x}+\mathit{h} and {\mathit{x}}_{\mathbf{n}\mathbf{e}\mathbf{w}}满足式（12）

13. \theta ≔\frac{F(\mathit{x})-F({\mathit{x}}_{\mathbf{n}\mathbf{e}\mathbf{w}})}{L(\mathrm{0})-L(\mathit{h})}

14. if \theta >\mathrm{0} then

15. \mathit{x}≔{\mathit{x}}_{\mathbf{n}\mathbf{e}\mathbf{w}}，\mathit{A}≔\mathit{J}{(\mathit{x})}^{\mathrm{T}}\mathit{J}(\mathit{x}),

\begin{array}{l}\mathit{g}≔\mathit{J}{(\mathit{x})}^{\mathrm{T}}\mathit{f}(\mathit{x}),found≔\left({‖\mathit{g}‖}_{\mathrm{\infty }}\le {ϵ}_{\mathrm{1}}\right),\\ \mu ≔\mu \cdot \mathrm{m}\mathrm{a}\mathrm{x}\{\frac{\mathrm{1}}{\mathrm{3}},\mathrm{1}-{\left(\mathrm{2}\theta -\mathrm{1}\right)}^{\mathrm{3}}\},\nu ≔\mathrm{2}\end{array}

16. else

17. \mu ≔\mu \cdot \nu ,\nu ≔\mathrm{2}\cdot \nu

18. end if

19. end if

20. {\mathit{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}}≔\mathit{x},loss=F({\mathit{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}})

21. end while

22. Loss=\mathrm{m}\mathrm{i}\mathrm{n}\{loss\},\mathit{x}\mathrm{为}Loss\mathrm{对应}\mathrm{的}{\mathit{x}}_{\mathbf{r}\mathbf{e}\mathbf{s}}

从式（2）、（5）、（7）中分别计算P_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}、 P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}、 P_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}

23. 计算f\mathrm{\text{'}}(x_{\mathrm{0}})

24. error=\mathrm{0.1}

25. while error>{\mathrm{10}}^{-\mathrm{5}} do

26. n=n+\mathrm{1}

27. x_{n+\mathrm{1}}=x_n-\frac{f(x_n)}{f\mathrm{\text{'}}(x_n)}

28. error=|x_{n+\mathrm{1}}-x_n|

29. end while

30. d_{\mathrm{p}}=x_{n+\mathrm{1}}

31. 从式（32）计算n_{\mathrm{p}}

4 现场试验效果

为了说明模型的有效性，将工程上常用的摩阻拟合模型与本文模型的计算结果进行对比。以四川盆地某段页岩施工实例的加砂压裂施工曲线为例，如图1所示。

根据工程中的记录和图1，第1次阶梯降排量过程从7点19分开始至7点27分结束，共有5个阶梯，对应的压裂液排量为\mathit{Q}=[\mathrm{11.00,9.00,7.00,5.00,3.04}]，相应的施工压力为{\mathit{P}}_{\mathbf{s}}=[\mathrm{44.0,41.3,38.4,36.2,34.5}]，停泵压力P_{\mathrm{0}}=\mathrm{33}\mathrm{M}\mathrm{P}\mathrm{a}。

工程上常用的摩阻拟合模型（式（6））直接使用已有的曲线拟合函数，拟合的参数范围均设置为(\mathrm{0},\mathrm{\infty })，运行多次，均得到不同结果，见表2。

由表2可知，每次运行得到的结果都不同，并且发现3次结果均满足k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}+k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}为定值。常用模型的算法特性决定了拟合结果不唯一，并且均为算法的最优解。同时，原有模型并不遵循施工实际的客观规律，即井筒摩阻大于孔眼摩阻与近井摩阻之和。而本文提出的模型优化了已有模型的不足。

输入井深范围[a,b]=[\mathrm{3.6,3.7}] km，由表1和式（11）可以得到井筒摩阻系数的范围。多次选择初值后，通过比较得到的最优结果为：k_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=\mathrm{0.075}\mathrm{02}，k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\mathrm{0.011}\mathrm{12}，k_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=\mathrm{0.506}\mathrm{8}，\gamma =\mathrm{1.657}\mathrm{3}，\beta =\mathrm{0.500}\mathrm{7}。

此时的目标函数（损失函数）的值为Loss=F(\mathit{x})=\mathrm{0.116}\mathrm{8}。当q=\mathrm{11} m³/min时，拟合的井筒摩阻、孔眼摩阻及近井摩阻为P_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=\mathrm{9.077}\mathrm{1}\mathrm{M}\mathrm{P}\mathrm{a}，P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\mathrm{0.591}\mathrm{3}\mathrm{M}\mathrm{P}\mathrm{a}，P_{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=\mathrm{1.683}\mathrm{4}\mathrm{M}\mathrm{P}\mathrm{a}。

相应的井筒摩阻、孔眼摩阻、近井摩阻、总摩阻的拟合曲线与实际总摩阻的散点图如图2所示。

由图2可知，本文方法的计算结果并不会随着运行次数的增多而改变，并且遵循了实际施工的客观规律。

输入 c=\mathrm{105.1}，t=\mathrm{5}\mathrm{749}，d_{\mathrm{0}}=\mathrm{11.7}，根据之前拟合得到的孔眼摩阻系数k_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}}=\mathrm{0.011}\mathrm{12}，计算可得有效开启孔数为n_{\mathrm{p}}=\mathrm{27}，磨蚀后的平均孔眼直径为d_{\mathrm{p}}=\mathrm{11.932}\mathrm{4}\mathrm{m}。

5 结　论

根据页岩气施工过程中总摩阻的构成原理，建立了页岩气中井筒摩阻系数、孔眼摩阻系数、近井摩阻系数与孔眼摩阻流动指数、近井摩阻流动指数的非线性最小二乘模型。该模型符合施工实际情况，解决了常用摩阻拟合模型最优解不唯一、不符合施工情况等问题；通过得到的孔眼摩阻系数与孔眼摩阻流动指数计算出拟合的孔眼摩阻大小，联立孔眼摩阻系数的计算公式与孔眼扩展方程，可利用非线性方程迭代法求得有效开启孔数与磨蚀后的孔眼平均直径。该模型从数学上量化评估了页岩气压裂改造效果，并为下一步的暂堵优化技术的实施提供了依据，且该方法实施成本低、效率高。

未来可考虑根据工程总结的经验进一步细化模型的参数约束条件以降低模型误差；同时，可根据工程经验选取更恰当的模型迭代初值以提高计算效率、改善计算结果。

参考文献

[1] Zhang Liehui, He Xiao, Li Xiaogang,et al.Shale gas exploration and development in the Sichuan basin:Progress,challenge and countermeasures[J].Natural Gas Industry B,2022,9(2):176‒186. doi:10.1016/j.ngib.2021.08.024

[2] Hu Chaolang, Lu Jing, He Xiaoming.Numerical solutions of a hypersingular integral equation with application to productivity formulae of horizontal wells at constant wellbore pressure[J].International Journal of Numerical Analysis and Modeling,2014,5(3):269‒288.

[3] Hu Chaolang, Lu Jing, He Xiaoming.Productivity formulae of an infinite-conductivity hydraulically fractured well producing at constant wellbore pressure based on numerical solutions of a weakly singular integral equation of the first kind[J].Mathematical Problems in Engineering,2012,2012(1):428596. doi:10.1155/2012/428596

[4] Hu Chaolang, He Xiaoming, Tao Lü.Euler-Maclaurin expansions and approximations of hypersingular integrals[J].Discrete and Continuous Dynamical Systems‒Series B,2015,20(5):1355‒1375. doi:10.3934/dcdsb.2015.20.1355

[5] Hu Chaolang, Lu Tao.Approximations of hypersingular integrals for negative fractional exponent[J].Journal of Computational Mathematics,2018,36(5):627‒643. doi:10.4208/jcm.1703-m2016-0544

[6] Wang Yonghui, Lu Yongjun, Li Yongping,et al.Progress and application of fracturing reconstruction technology in unconventional reservoirs[J].Acta Petrolei Sinica,2012,33(Supp1):149‒158.

[7] 王永辉,卢拥军,李永平,等.非常规储层压裂改造技术进展及应用[J].石油学报,2012,33():149‒158.

[8] Zhang Linsen, Zhao Zhongcheng, Li Xiaowen.Application of refracturing technology in the old well of oilfield[J].China Foreign Energy,2006,11(4):42‒45.

[9] 张林森,赵忠诚,李效文.油田老井重复压裂技术研究应用[J].中外能源,2006,11(4):42‒45.

[10] Cao Xuejun, Wang Minggui, Kang Jie,et al.Fracturing technologies of deep shale gas horizontal wells in the Weirong Block,southern Sichuan Basin[J].Natural Gas Industry B,2020,7(1):64‒70. doi:10.1016/j.ngib.2019.07.003

[11] Liang Xing, Zhu Juhui, Shi Xiaozhi,et al.Staged fracturing of horizontal shale gas wells with temporary plugging by sand filling[J].Natural Gas Industry B,2017,4(2):134‒140. doi:10.1016/j.ngib.2017.07.017

[12] Wu Kan, Olson J E.Mechanisms of simultaneous hydraulic-fracture propagation from multiple perforation clusters in horizontal wells[J].SPE Journal,2016,21(3):1000‒1008. doi:10.2118/178925-pa

[13] Miller C, Waters G, Rylander E.Evaluation of production log data from horizontal wells drilled in organic shales[C]//Proceedings of the North American Unconventional Gas Conference and Exhibition.Texas:SPE,2011:SPE-144326-MS. doi:10.2118/144326-ms

[14] Wang Bo.Study on plugging and steering law of temporary plugging fracturing fractures[D].Beijing:China University of Petroleum(Beijing),2019.

[15] 王博.暂堵压裂裂缝封堵与转向规律研究[D].北京:中国石油大学(北京),2019.

[16] Zhou Fujian, Su Hang, Liang Xingyuan,et al.Integrated hydraulic fracturing techniques to enhance oil recovery from tight rocks[J].Petroleum Exploration and Development,2019,46(5):1065‒1072. doi:10.1016/s1876-3804(19)60263-6

[17] Xiong Chunming, Shi Yang, Zhou Fujian,et al.High efficiency reservoir stimulation based on temporary plugging and diverting for deep reservoirs[J].Petroleum Exploration and Development,2018,45(5):948‒954. doi:10.1016/s1876-3804(18)30098-3

[18] Wang Bo, Zhou Fujian, Yang Chen,et al.Experimental st-udy on injection pressure response and fracture geometry during temporary plugging and diverting fracturing[J].SPE Journal,2020,25(2):573‒586. doi:10.2118/199893-pa

[19] Guo Jianchun, Zhao Feng, Zhan Li,et al.Recent advances and development suggestions of temporary plugging and diverting fracturing technology for shale gas reservoirs in the Sichuan Basin[J].Petroleum Drilling Techniques,2023,51(4):170‒183.

[20] 郭建春,赵峰,詹立,等.四川盆地页岩气储层暂堵转向压裂技术进展及发展建议[J].石油钻探技术,2023,51(4):170‒183.

[21] 李海军,刘恒,李力民,等.套损检测技术在储气库的应用[C]//2017年全国天然气学术年会论文集, 2017.

[22] Zhang Chuanju, Xu Hongliang.Influence and correction of drilling coring[J].China Petroleum and Chemical Standard and Quality,2023,43(3):137‒139.

[23] 张传举,徐宏亮.钻井取心的影响与校正[J].中国石油和化工标准与质量,2023,43(3):137‒139.

[24] Cramer D D.The application of limited-entry techniques in massive hydraulic fracturing treatments[C]//Proceedings of the SPE Production Operations Symposium.Oklahoma:SPE,1987:SPE-16189-MS. doi:10.2118/16189-ms

[25] Lord D L, Shah S N, Rein R G,et al.Study of perforation friction pressure employing a large-scale fracturing flow simulator[C]//Proceedings of the SPE Annual Technical Conference and Exhibition.Louisiana:SPE,1994:SPE-28508-MS. doi:10.2118/28508-ms

[26] Kim G H, Wang J Y.Interpretation of hydraulic fracturing pressure in low-permeability gas formations[C]//Proceedings of the SPE Production and Operations Symposium.Oklahoma:SPE,2011:SPE-141525-MS. doi:10.2118/141525-ms

[27] Howell G W, Weathers T W.Aerospace fluid component designers' handbook[M].Air Force Rocket Propulsion Laboratory, Research and Technology Division,Air Force Systems Command,1967. doi:10.21236/ad0817250

[28] 周宁杰,孙涛.计算不同管壁粗糙度产生的井筒摩阻压降[J].内蒙古石油化工,2021,47(10):34‒37.

[29] Xiong Ying, Liu Youquan, Mei Zhihong,et al.Slick water technology of high salinity resistance for shale gas development in Sichuan[J].Chemical Engineering of Oil & Gas,2019,48(3):62‒65. doi:10.3969/j.issn.1007-3426.2019.03.011

[30] 熊颖,刘友权,梅志宏,等.四川页岩气开发用耐高矿化度滑溜水技术研究[J].石油与天然气化工,2019,48(3):62‒65. doi:10.3969/j.issn.1007-3426.2019.03.011

[31] Crump J B, Conway M W.Effects of perforation-entry friction on bottomhole treating analysis[J].Journal of Petroleum Technology,1988,40(8):1041‒1048. doi:10.2118/15474-pa

[32] Romero J, Mack M G, Elbel J L.Theoretical model and numerical investigation of near-wellbore effects in hydraulic fracturing[J].SPE Production & Facilities,2000,15(2):76‒82. doi:10.2118/63009-pa

[33] El‒Rabba A M, Shah S N, Lord D L.New perforation pressure-loss correlations for limited-entry fracturing treatments[J].SPE Production & Facilities,1999,14(1):63‒71. doi:10.2118/54533-pa

[34] Weijers L, Wright C A, Sugiyama H,et al.Simultaneous propagation of multiple hydraulic fractures—Evidence,impact and modeling implications[C]//Proceedings of the Internatio-nal Oil and Gas Conference and Exhibition in China.Beijing:SPE,2000:SPE-64772-MS. doi:10.2118/64772-ms

[35] Mondal S, Ugueto G, Huckabee P,et al.Uncertainties in step-down test interpretation for evaluating completions effectiveness and near wellbore complexities[C]//Proceedings of the 7th Unconventional Resources Technology Conference.Denver:American Association of Petroleum Geologists,2019. doi:10.15530/urtec-2019-1141

[36] Mondal S, Zhang Min, Huckabee P,et al.Advancements in step down tests to guide perforation cluster design and limited entry pressure intensities—Learnings from field tests in multiple basins[C]//Proceedings of the SPE Hydraulic Fracturing Technology Conference and Exhibition.Virtual:SPE,2021:D021S006R002. doi:10.2118/204147-ms

[37] Perkins T K, Kern L R.Widths of hydraulic fractures[J].Journal of Petroleum Technology,1961,13(9):937‒949. doi:10.2118/89-pa

[38] Nordgren R P.Propagation of a vertical hydraulic fracture[J].Society of Petroleum Engineers Journal,1972,12(4):306‒314. doi:10.2118/3009-pa

[39] Massaras L V, Dragomir A, Chiriac D.Enhanced fracture entry friction analysis of the rate step-down test[C]//Proceedings of the SPE Hydraulic Fracturing Technology Conference.Texas:SPE,2007:SPE-106058-MS. doi:10.2118/106058-ms

[40] Madsen K, Nielsen H B, Tingleff O.Methods for non-linear least squares problems[J].Informatics and Mathematical Modelling Technical University of Denmark,2004.

[41] Willingham J D, Tan H C, Norman L R.Perforation friction pressure of fracturing fluid slurries[C]//Proceedings of the Low Permeability Reservoirs Symposium.Denver:SPE,1993:SPE-25891-MS. doi:10.2523/25891-ms

[42] Long Gongbo, Liu Songxia, Xu Guanshui,et al.Modeling of perforation erosion for hydraulic fracturing applications[C]//Proceedings of the SPE Annual Technical Conference and Exhibition.Houston:SPE,2015:D021S015R008. doi:10.2118/174959-ms

[43] Vincent M C, Miller H B, Milton‒Tayler D,et al.Erosion by proppant:A comparison of the erosivity of sand and ceramic proppants during slurry injection and flowback of proppant[C]//Proceedings of the SPE Annual Technical Conference and Exhibition.Houston:SPE,2004:SPE-90604-MS. doi:10.2523/90604-ms

[44] Jorge N, Stephen J W.Numerical optimization[M].Berlin:Spinger,2006. doi:10.1007/978-3-540-35447-5

基金资助

国家自然科学基金项目(11971337)