A fracture morphology inversion method and device combining full time steps and physical models

Abstract

This invention discloses a method and apparatus for crack morphology inversion combining full-time step and physical model, comprising: inputting crack height and crack offset in the crack height direction into a physical model; first, processing the crack height through a crack propagation sub-model to output crack geometric parameters for the full-time step; then, processing the crack geometric parameters and crack offset in the crack height direction through a fiber strain response sub-model to output axial strain rate at simulated measurement points for the full-time step; adjusting the axial strain rate at simulated measurement points based on positive and negative strain rate scaling factors; calculating the loss value between the adjusted simulated axial strain rate and the measured axial strain rate; if the loss value is greater than a preset loss threshold, adjusting the corresponding optimization control parameters until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached; and outputting the target crack geometric parameters for the full-time step based on the finally adjusted optimization control parameters. This invention can improve the accuracy and efficiency of crack morphology inversion.

Description

Technical Field

[0001] This invention relates to the field of oil and gas field development technology, and in particular to a method and apparatus for fracture morphology inversion that combines full time step and physical model. Background Technology

[0002] In hydraulic fracturing engineering in oil and gas fields, accurate inversion of fracture morphology (fracture geometry / fracture geometric parameters) is the core basis for objectively evaluating the fracturing effect, optimizing construction parameters, predicting production capacity, and guiding subsequent well network deployment. It is of great significance for improving oil and gas recovery rate and development efficiency.

[0003] However, existing crack morphology inversion methods suffer from the prominent problem of "difficulty in balancing high accuracy and high efficiency": in pursuit of inversion accuracy, a large number of spatial and temporal inversion optimizations are required for parameters such as crack length and crack aperture, resulting in high computational complexity and high resource consumption, which cannot meet the real-time requirements of engineering; in order to improve inversion efficiency, single-time-step data is mainly used for crack geometric parameter inversion, but single-time-step inversion will disrupt the continuity of crack propagation, making it difficult to guarantee the physical rationality of the inversion results, ultimately leading to insufficient inversion accuracy and large deviation from the actual crack morphology.

[0004] There is currently no effective solution to the above problems. Summary of the Invention

[0005] This specification provides a method and apparatus for crack morphology inversion that combines full-time stepping with a physical model, in order to solve the problem that existing crack morphology inversion techniques cannot balance inversion accuracy and computational efficiency.

[0006] Firstly, embodiments of this specification provide a crack morphology inversion method combining full-time step and physical model, including: Obtain optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors; The crack height and the crack offset in the crack height direction are input into the physical model, which includes a crack propagation sub-model and an optical fiber strain response sub-model. First, the crack height is processed by the crack propagation sub-model to output the crack geometric parameters at the full time step. Then, the crack geometric parameters and the crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate at the simulated measurement point at the full time step. The crack geometric parameters include crack length and crack aperture. The axial strain rate of the simulated measuring point is adjusted based on the positive and negative strain rate scaling factors, and the loss value between the adjusted axial strain rate of the simulated measuring point and the measured axial strain rate of the measuring point is calculated. If the loss value is greater than the preset loss threshold, the corresponding optimization control parameters are adjusted until the loss value is less than or equal to the preset loss threshold, or the preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters for the entire time step are output.

[0007] In some embodiments, processing the crack height through a crack propagation sub-model to output crack geometry parameters at the full time step includes: The crack height is processed according to the following formula to output the crack length at the full time step:

[0008] The crack height is processed according to the following formula to output the crack aperture at the full time step:

[0009] Where L(t) is the crack length at time t; G is the shear modulus; q is the flow rate across the crack cross section; and v is the Poisson's ratio of the rock. denoted as η, where η is the dynamic viscosity coefficient of the fracturing fluid; h is the fracture height; w(t) is the fracture aperture at time t; and a and b are constants.

[0010] In some embodiments, the step of processing the crack geometry parameters and the crack offset in the crack height direction using a fiber optic strain response sub-model to output the axial strain rate at the simulated measurement point across the entire time step includes: Based on the crack length and the crack offset in the crack height direction, the influence coefficients in the fiber strain response sub-model are corrected. The influence coefficients depend on the relative positional distance between the crack element and the fiber measuring point. The axial displacement of the fiber optic measuring point at the full time step is determined based on the corrected influence coefficient and the crack opening. The fiber gauge length is set based on the crack length, and the axial strain of the fiber measuring point is determined at the full time step according to the fiber gauge length and the axial displacement of the fiber measuring point. Based on the axial strain at the fiber optic measuring point and the simulation calculation time interval, the axial strain rate of the simulated measuring point at the full time step is determined and output.

[0011] In some embodiments, determining the axial displacement of the fiber optic measuring point across the entire time step based on the corrected influence coefficient and the crack aperture includes: The axial displacement of the fiber optic measuring point throughout the entire time step is determined using the following formula:

[0012] The step of determining the axial strain of the fiber optic measuring point across the entire time step based on the fiber gauge length and the axial displacement of the fiber optic measuring point includes: The axial strain at the fiber optic measuring point throughout the entire time step is determined using the following formula:

[0013] The step of determining and outputting the axial strain rate of the simulated measuring point across the entire time step based on the axial strain of the fiber optic measuring point and the simulation calculation time interval includes: The axial strain rate at the simulated measurement points throughout the entire time step is determined using the following formula:

[0014] Among them, u f (t) represents the axial displacement of the fiber optic measuring point at time t; M represents the total number of crack elements; v represents the Poisson's ratio of the rock; I1 and I2 are the corrected influence coefficients; w i (t) represents the crack opening of the i-th crack element at time t; L represents the axial strain at the fiber optic measuring point at time t. g z is the fiber gauge length; z is the center position of the fiber measuring point; u f (t)(z+L g / 2)- u f (t)(zL g / 2) represents the relative displacement difference between the two ends of the fiber gauge length; Let be the axial strain rate at the simulated measuring point at time t; for Axial strain at the fiber optic measuring point at time t; This is for simulating the calculation time interval.

[0015] In some embodiments, the positive and negative strain rate scaling factors include a positive strain rate scaling factor and a negative strain rate scaling factor; correspondingly, adjusting the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factors includes: If the axial strain rate of the simulated measuring point is greater than or equal to zero, the axial strain rate of the simulated measuring point is adjusted based on the positive strain rate scaling factor, and the product of the positive strain rate scaling factor and the axial strain rate of the simulated measuring point is used as the adjusted axial strain rate of the simulated measuring point. If the axial strain rate of the simulated measuring point is less than zero, the axial strain rate of the simulated measuring point is adjusted based on the negative strain rate scaling factor. The product of the negative strain rate scaling factor and the axial strain rate of the simulated measuring point is used as the adjusted axial strain rate of the simulated measuring point.

[0016] In some embodiments, the fracture propagation sub-model also outputs the fracturing fluid pressure in the injection hole; correspondingly, if the loss value is greater than a preset loss threshold, the corresponding optimization control parameters are adjusted, including: If the loss value is greater than the preset loss threshold, the error type is determined by combining the error characteristics of the injection hole fracturing fluid pressure and the adjusted simulated axial strain rate and the measured axial strain rate, and the corresponding optimized control parameters are adjusted according to the error type.

[0017] In some embodiments, determining the error type by combining the error characteristics of the fracturing fluid pressure at the injection hole and the adjusted simulated axial strain rate at the measuring point with the measured axial strain rate at the measuring point includes: If the fracturing fluid pressure in the injection hole is within the preset engineering reasonable range, and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are characterized by an overall amplitude deviation, then it is determined to be an amplitude deviation. If the fracturing fluid pressure in the injection hole exceeds the preset reasonable engineering range, and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are manifested as deviations in waveform width, signal arrival time, or spatial distribution, then it is determined to be a spatiotemporal morphological deviation.

[0018] In some embodiments, adjusting the corresponding optimization control parameters according to the error type includes: If the error type is amplitude deviation, then adjust the positive and negative strain rate scaling factors; If the error type is spatiotemporal morphological deviation, then adjust the crack height and the crack offset in the crack height direction.

[0019] Secondly, embodiments of this specification also provide a crack morphology inversion device that combines full-time stepping with a physical model, comprising: The acquisition module is used to acquire optimized control parameters, which include crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. The output module is used to input the crack height and the crack offset in the crack height direction into the physical model. The physical model includes a crack propagation sub-model and an optical fiber strain response sub-model. First, the crack height is processed by the crack propagation sub-model to output the crack geometric parameters at the full time step. Then, the crack geometric parameters and the crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate of the simulated measurement point at the full time step. The crack geometric parameters include crack length and crack aperture. The loss value calculation module is used to adjust the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factor, and calculate the loss value between the adjusted axial strain rate of the simulated measuring point and the measured axial strain rate of the actual measuring point. The adjustment iteration module is used to adjust the corresponding optimization control parameters if the loss value is greater than a preset loss threshold, until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters for the entire time step are output.

[0020] Thirdly, embodiments of this specification also provide a computer-readable storage medium storing computer program instructions that, when executed by a processor, implement the steps of the aforementioned crack morphology inversion method combining full time step and physical model.

[0021] This specification provides a method and apparatus for crack morphology inversion combining full-time step and physical model. First, optimized control parameters are obtained, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. Then, the crack height and crack offset in the crack height direction are input into the physical model, which includes a crack propagation sub-model and an optical fiber strain response sub-model. The crack height is first processed by the crack propagation sub-model to output the crack geometric parameters for the full-time step. Then, the crack geometric parameters and crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate of the simulated measuring point for the full-time step. The crack geometric parameters include crack length and crack aperture. Finally, the axial strain rate of the simulated measuring point is adjusted based on the positive and negative strain rate scaling factors, and the loss value between the adjusted simulated axial strain rate and the measured axial strain rate is calculated. Finally, if the loss value is greater than a preset loss threshold, the corresponding optimization control parameters are adjusted until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters for the entire time step are output. In this embodiment, only four key physical parameters need to be obtained: crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. This can effectively reduce computational complexity and lay the foundation for improving inversion efficiency. Through the physical model, a quantitative mapping of "optimized control parameters → crack geometry parameters → simulated strain rate" is established, which can ensure the physical rationality of the inversion results. Moreover, the inversion process is calculated for the entire time step, which can avoid the fragmentation problem of "independent fitting in a single time step," fully capture the continuous evolution law of crack geometry parameters, and eliminate the need to optimize the crack geometry parameters for each time step separately, greatly reducing the amount of computation and balancing efficiency and accuracy. By adjusting the axial strain rate of the simulated measuring point using positive and negative strain rate scaling factors, the loss value between the simulated and measured axial strain rates is calculated. When the loss value exceeds a preset loss threshold, the corresponding optimized control parameters are adjusted, and the data is re-input into the physical model for quantitative mapping of "optimized control parameters → crack geometry parameters → simulated strain rate," which effectively improves the inversion accuracy. Finally, the target crack geometry parameters output across the entire time step provide accurate and comprehensive data support for engineering applications. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings: Figure 1A flowchart illustrating a crack morphology inversion method combining full time step and physical model provided in the embodiments of this specification; Figure 2 A schematic diagram of the crack propagation and fiber strain response sub-model provided in the embodiments of this specification; Figure 3 Single-time-step strain rate comparison images provided in the embodiments of this specification; Figure 4 A comparison image of the crack length between the real crack and the inverted crack provided in the embodiments of this specification; Figure 5 A comparison image of the crack aperture of a real crack and an inverted crack provided for the embodiments of this specification; Figure 6 This is a schematic diagram of the structural composition of a crack morphology inversion device that combines full time step and physical model, as provided in the embodiments of this specification. Figure 7 This is a schematic diagram of the structural composition of the electronic device provided in the embodiments of this specification. Detailed Implementation

[0023] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this specification, and not all embodiments. Based on the embodiments in this specification, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this specification.

[0024] As mentioned earlier, accurate inversion of fracture morphology (fracture geometry / fracture geometry parameters) is the core basis for objectively evaluating the fracturing effect, optimizing construction parameters, predicting production capacity, and guiding subsequent well network deployment. It is of great significance for improving oil and gas recovery and development efficiency.

[0025] Crack morphology inversion methods based on low-frequency distributed acoustic sensing (LF-DAS) have become an important tool for evaluating fracturing effectiveness. Existing crack morphology inversion methods mainly rely on fitting fiber strain response characteristics. They invert crack geometric parameters by comparing LF-DAS strain rate waterfall plots generated from the fiber strain response with field-measured data and iteratively optimizing the results. However, existing methods primarily rely on single-time-step data, making it difficult to accurately reflect the dynamic strain rate and crack morphology evolution during crack propagation. Furthermore, traditional methods require extensive spatial and temporal optimization of parameters such as crack width and length, which consumes significant computational resources.

[0026] In summary, existing crack morphology inversion methods suffer from the prominent problem of "difficulty in balancing high accuracy and high efficiency".

[0027] To address the aforementioned issues, this specification provides a method and apparatus for crack morphology inversion that combines full-time step data with a physical model. By simultaneously considering data from all time steps and integrating it with a physical model (crack propagation sub-model and fiber strain response sub-model), the method effectively solves the crack morphology inversion bias problem in traditional methods and significantly reduces computational costs, thereby achieving high-precision and high-efficiency crack morphology inversion.

[0028] It is understood that the methods described in the embodiments of this specification can be applied to electronic devices, which can refer to electronic devices with data computing, processing, and storage capabilities. These electronic devices can be terminals such as PCs (Personal Computers), tablets, smartphones, wearable devices, and intelligent robots; they can also be servers. A server can be an independent physical server, a server cluster or distributed system composed of multiple physical servers, or a cloud server providing cloud computing services.

[0029] See Figure 1 As shown in the embodiments of this specification, a crack propagation sub-model and an optical fiber strain response sub-model method are provided. In specific implementations, this method may include the following: S101: Obtain optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors.

[0030] In some embodiments, the aforementioned optimization control parameters can be the core input for iterative optimization or the core control variables for crack geometry parameters. These parameters may include crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. Crack height can be a core physical parameter of the crack propagation sub-model (PKN model), directly determining the dynamic evolution of crack length and aperture, eliminating the need for separate optimization of these parameters. Crack offset in the crack height direction can be used to correct the influence coefficients in the fiber strain response sub-model (i.e., to correct the relative positional distance between crack elements and fiber measuring points), thereby ensuring the simulation accuracy of the axial strain rate at the measuring points. Positive and negative strain rate scaling factors can be used to further adjust the axial strain rate at the simulated measuring points, avoiding mismatches between simulated and measured data due to simplified geological parameters (such as assuming homogeneous rock).

[0031] The initial values ​​of the above-mentioned optimized control parameters can be set based on engineering experience. For example, the initial value of the crack height can be referenced to the formation thickness (e.g., 10-15m when the formation thickness is 20m), the initial value of the offset can be 0m (when there is no clear offset), and the initial value of the scaling factor can be 1.0 (without initial amplitude correction). Subsequently, the initial values ​​of the crack height and the initial value of the crack offset in the crack height direction can be input into the physical model first, and then the axial strain rate of the simulated measuring point can be adjusted based on the initial positive and negative strain rate scaling factors. If the error between the simulated axial strain rate and the measured axial strain rate is too large, the above initial values ​​should be adjusted.

[0032] By obtaining the crack height, the crack offset in the crack height direction, and the positive and negative strain rate scaling factors, the approach of "optimizing massive spatiotemporal parameters such as crack length and aperture one by one" in existing technologies is abandoned. This can reduce the computational complexity from the source and lay the foundation for improving the inversion efficiency.

[0033] S102: Input the crack height and the crack offset in the crack height direction into the physical model. The physical model includes a crack propagation sub-model and an optical fiber strain response sub-model. First, process the crack height through the crack propagation sub-model and output the crack geometric parameters at the full time step. Then, process the crack geometric parameters and the crack offset in the crack height direction through the optical fiber strain response sub-model and output the axial strain rate of the simulated measurement point at the full time step. The crack geometric parameters include crack length and crack aperture.

[0034] In some embodiments, the fracture propagation sub-model described above can adopt the PKN model, assuming a constant fracture height, an elastic rock medium, and a Newtonian fluid. Through continuity equations (fracturing fluid continuity equations within the fracture), fracture propagation control equations (hydraulic fracturing fracture propagation control equations), and boundary conditions, a quantitative mapping of "fracture height → fracture length / aperture" is established to ensure that the output parameters conform to the physical laws of hydraulic fracturing. Given parameters such as fracture height, the fracture length, fracture aperture, and intra-fracture pressure (or injection hole fracturing fluid pressure) can be calculated at each moment.

[0035] The continuity equation (the fracturing fluid continuity equation within the fracture) can be as follows: (1) Where q(x,t) represents the flow rate through the cross-section of the crack; q l (x,t) represents the fluid loss rate per unit length of the crack; A(x,t) represents the cross-sectional area of ​​the crack.

[0036] This equation is the mass conservation equation, which describes the mass balance relationship when fracturing fluid flows within a fracture.

[0037] The governing equations for fracture propagation (the governing equations for fracture propagation in hydraulic fracturing) can be as follows: (2) (3) Among them, C e Let G be the elastic modulus, v be the Poisson's ratio of the rock, and h be the crack height; and w be the crack opening. denoted as ν, where ν is the dynamic viscosity coefficient of the fracturing fluid; G is the shear modulus of the rock.

[0038] This equation is a dynamic coupling equation, which can describe the mechanical process of fluid pressure driving crack opening.

[0039] The above boundary conditions can be as follows: ① Injection end (x=0) boundary: Given flow rate or given pressure: (4) (5) ② Crack tip (x=L(t)) boundary: (6) (7) This equation defines the physical boundary constraints for solving the system of equations, and can be used to constrain the solutions to ensure physical correctness.

[0040] Formulas (1)-(3) under boundary conditions (4)-(7) can yield the fracture aperture distribution w(x,t), pressure distribution characteristics p(x,t), and flow distribution q(x,t) at any given time. Based on this, the fracture length, representative aperture, and fracturing fluid pressure in the injection hole are further extracted and given by (8)-(10): The crack length L(t) can be: (8) This equation can be used to calculate the total length of the crack in the direction of extension.

[0041] The crack aperture w(t) can be: (9) This equation can be used to calculate the crack opening.

[0042] The fracturing fluid pressure p(t) in the injection hole can be: (10) This equation can be used to calculate the pressure of the fracturing fluid at the injection point.

[0043] Among them, the crack propagation sub-model (PKN model) mentioned above can include the above formulas (8)-(10), and their basic formulas are the above formulas (1)-(7).

[0044] In some embodiments, the process of processing the crack height through the crack propagation sub-model and outputting the crack geometry parameters at the full time step in S102 above may, in specific implementation, include: The crack height is processed according to the following formula to output the crack length at the full time step:

[0045] The crack height is processed according to the following formula to output the crack aperture at the full time step:

[0046] Where L(t) is the crack length at time t (i.e., the crack length over the entire time step); G is the shear modulus; q is the flow rate across the crack cross section; and v is the Poisson's ratio of the rock. denoted as , where h is the dynamic viscosity coefficient of the fracturing fluid; h is the fracture height; w(t) is the fracture aperture at time t (i.e., the fracture aperture over the entire time step); a and b are constants. Based on engineering application experience with the PKN model, the value of a can range from 0.65 to 0.70 (preferably 0.68), and the value of b can range from 2.4 to 2.6 (preferably 2.5), which can ensure the matching degree between the analytical solution and the actual engineering data.

[0047] Specifically, the crack height h can be processed using the above formula (8) to output the crack length L(t) for the entire time step. This crack length L(t) is related to time. Proportional to the crack height h 4 The crack height is inversely proportional to the crack height, which conforms to the physical law that "the greater the crack height, the greater the fluid flow resistance, and the slower the propagation." The crack height h can be processed using the above formula (9) to output the crack aperture w(t) for the entire time step. This crack aperture w(t) is inversely proportional to time. It is directly proportional to the fracture height h and inversely proportional to the fracture height h, which can reflect the characteristic that "the greater the fracture height, the stronger the mechanical constraint, and the greater the difficulty of opening". The fracture height h can also be processed by the above formula (10) to output the injection hole fracturing fluid pressure p(t). The injection hole fracturing fluid pressure p(t) can be used to help judge the error type so as to adjust the corresponding optimization control parameters according to the error type. This part will be explained separately later, and will not be repeated here.

[0048] By providing quantitative calculation methods for crack length and aperture, the output of the crack propagation sub-model can have a clear mathematical basis, avoiding the subjectivity of parameter calculation and providing accurate input for subsequent strain rate conversion.

[0049] In some embodiments, the aforementioned fiber strain response sub-model can be constructed based on the displacement discontinuity method (DDM). Its core is to transform "crack geometric parameters" into "fiber-monitorable strain rates," achieving the transformation from "physical parameters to monitoring data." The output of the crack propagation sub-model (such as crack length, crack aperture, or crack width) can serve as the input source for the fiber strain response sub-model. The offset in the crack height direction directly corrects the latter's spatial correlation parameters (influence coefficients), forming a complete link of "crack geometric parameters → spatial correction → strain rate output." Specifically, the crack length and crack aperture output from the aforementioned crack propagation sub-model can be used as source terms to calculate the axial displacement (fiber measuring point axial displacement) u at each measuring point location in the fiber according to the elasticity equation. f (t), and then differentiating with respect to displacement, the simulated strain rate (simulated axial strain rate of the measuring point) at each measuring point of the optical fiber at each time moment is derived, denoted as .

[0050] Among them, the axial displacement u of the fiber optic measuring point f (t) can be: (11) Among them, w i (t) represents the crack aperture of the i-th crack element at time t, calculated from the previous PKN model; I1 and I2 are the corrected influence coefficients (pure geometric terms, depending on the relative position distance between the crack element and the fiber optic measuring point, where I1 can be called the first influence coefficient and I2 can be called the second influence coefficient); M is the total number of crack elements; v is the Poisson's ratio of the rock. This equation can be used to convert "crack opening width" into "rock displacement".

[0051] Axial strain at the above fiber optic measuring points It can be: (12) Among them, L g The fiber gauge length (i.e., the spatial resolution when measuring strain with fiber optic cables); u f (t)(z+L g / 2)- u f (t)(zL g / 2) represents the relative displacement difference between the two ends of the fiber gauge length; This equation can be used to convert "rock displacement" into "fiber axial strain".

[0052] Axial strain rate at the above simulated measurement points It can be: (13) in, , These are the current time t and the next time, respectively. Calculated strain value (axial strain at fiber optic measuring point). The time interval for simulation calculation (i.e., the simulation calculation time interval); This equation can be used to convert "static strain" into "dynamic strain rate".

[0053] The fiber strain response sub-model mentioned above can include the above formulas (11)-(13). The three formulas form a progressive relationship of "crack opening → axial displacement → axial strain → axial strain rate". The output of each step is the input of the next step, and all parameters are related to the crack length, opening, and optimization control parameters (offset) output by the PKN model, forming a complete parameter transfer link.

[0054] In some embodiments, the process of processing the crack geometry parameters and the crack offset in the crack height direction using the fiber optic strain response sub-model in S102 above, and outputting the axial strain rate of the simulated measurement point throughout the time step, may include, in specific implementation: S1: Based on the crack length and the crack offset in the crack height direction, correct the influence coefficient in the fiber strain response sub-model. The influence coefficient depends on the relative positional distance between the crack element and the fiber measuring point. S2: Determine the axial displacement of the fiber optic measuring point throughout the time step based on the corrected influence coefficient and the crack opening. S3: Set the fiber gauge length based on the crack length, and determine the axial strain of the fiber measuring point at the full time step according to the fiber gauge length and the axial displacement of the fiber measuring point. S4: Based on the axial strain of the fiber optic measuring point and the simulation calculation time interval, determine and output the axial strain rate of the simulated measuring point for the entire time step.

[0055] In some embodiments, determining the axial displacement of the fiber optic measuring point across the entire time step based on the corrected influence coefficient and the crack aperture in S2 above may, in specific implementation, include: The axial displacement of the fiber optic measuring point throughout the entire time step is determined using the following formula:

[0056] In the above-mentioned S3, determining the axial strain of the fiber optic measuring point across the entire time step based on the fiber gauge length and the axial displacement of the fiber optic measuring point can, in specific implementation, include: The axial strain at the fiber optic measuring point throughout the entire time step is determined using the following formula:

[0057] In S4 above, determining and outputting the axial strain rate of the simulated measuring point based on the axial strain of the fiber optic measuring point and the simulation calculation time interval can, in specific implementation, include: The axial strain rate at the simulated measurement points throughout the entire time step is determined using the following formula:

[0058] Among them, u f (t) represents the axial displacement of the fiber optic measuring point at time t (i.e., the axial displacement of the fiber optic measuring point throughout the entire time step); M represents the total number of crack elements; v represents the Poisson's ratio of the rock; I1 and I2 are the corrected influence coefficients; w i (t) represents the crack aperture of the i-th crack element at time t (i.e., the crack aperture of the i-th crack element at the full time step). L represents the axial strain at the fiber optic measuring point at time t (i.e., the axial strain at the fiber optic measuring point throughout the entire time step); g z is the fiber gauge length; z is the center position of the fiber measuring point; u f (t)(z+L g / 2)- u f (t)(zL g / 2) represents the relative displacement difference between the two ends of the fiber gauge length; Let be the axial strain rate of the simulated measuring point at time t (i.e., the axial strain rate of the simulated measuring point throughout the entire time step). for Axial strain at the fiber optic measuring point at time t; This is for simulating the calculation time interval.

[0059] Specifically, the aforementioned influence coefficients I1 and I2 (referred to as influence coefficients before correction and modified influence coefficients after correction) are purely geometric terms. Their essence is the "displacement contribution weight of the crack element to the fiber optic measuring point," which is entirely determined by the relative spatial position between the crack element and the fiber optic measuring point. They are the core geometric parameters of the fiber optic strain response sub-model. For example, each crack element deforms the surrounding rock during propagation. Crack elements closer to the fiber optic measuring point contribute more to the displacement of the measuring point, and their corresponding influence coefficients are larger; conversely, those farther away contribute less. This can be based on the crack length L(t) and the crack offset v in the crack height direction. h Make corrections, for example: (1) The crack length L(t) can define the range of influence: only crack elements within 0≤x≤L(t) (crack extension direction) contribute displacement to the measuring point, and the influence coefficient of elements outside the range can be 0; (2) Offset v in the seam height direction hThe relative position can be corrected: when the crack shifts upwards / downwards, the vertical distance between the measuring point and the crack element changes, and I1 and I2 are recalculated through geometric integration (if the shift increases, the vertical distance increases, and the influence coefficient decreases).

[0060] The aforementioned fiber gauge length L g It is the spatial resolution of fiber optic strain measurement, which can meet L g ≤0.1 L min (t)(L min (t) represents the minimum crack length across the entire time step, and its purpose is: (1) Ensure that the rock displacement distribution within the gauge length is approximately linear to avoid distortion of strain calculation due to excessive gauge length; (2) Ensure that the gauge length can capture the local strain changes of crack propagation and improve the spatial resolution of the simulation data.

[0061] For example, if the minimum crack length L across the entire time step min If the gauge length is 8m, then the gauge length can be set to 0.5~0.8m.

[0062] After correcting the influence coefficients in the fiber strain response sub-model based on the crack length and the crack offset in the crack height direction, the corrected influence coefficients I1, I2 and the crack aperture w of the i-th crack element at time t can be processed by the above formula (11). i (t), outputting the axial displacement u of the fiber optic measuring point throughout the full time step. f (t). After setting the fiber gauge length based on the crack length, the fiber gauge length L can be processed by the above formula (12). g and axial displacement u of the fiber optic measuring point f (t), outputting the axial strain at the fiber optic measuring point throughout the full time step. Then, you can set the simulation calculation time interval. With the data acquisition interval consistent with that of the optical fiber (e.g., 5s), the axial strain rate of the simulated measuring point can be calculated by the "strain difference between adjacent time steps", that is, by the above formula (13), and finally the axial strain rate of the simulated measuring point for the entire time step is output. .

[0063] By correcting the influence coefficient and setting a reasonable gauge length, the spatial and temporal resolution of the simulated strain rate can be ensured to match the characteristics of the measured data, thus improving the effectiveness of error comparison. By clarifying the quantitative calculation methods for each physical quantity, the implementation of the fiber optic strain response sub-model can be made operable and repeatable, avoiding simulation data deviations caused by ambiguity in the calculation process.

[0064] The physical models described above (crack propagation sub-model and fiber strain response sub-model) allow for the simultaneous calculation of parameters across all time steps, rather than independent calculations at each time step. This ensures the continuity of crack propagation (e.g., t). n+1 The crack length at time t is based on n (Evolution of extended results at time).

[0065] S103: Adjust the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factors, and calculate the loss value between the adjusted axial strain rate of the simulated measuring point and the measured axial strain rate of the actual measuring point.

[0066] In some embodiments, prior to S103 above, a distributed fiber optic sensing monitoring system may be deployed in the monitoring well. Raw DAS data (DAS data refers to full-band acoustic or vibration signals acquired by Distributed Acoustic Sensing (DAS) technology) can be acquired, and low-LF-DAS data (LF-DAS data refers to low-frequency distributed acoustic sensing (DAS) data, a subset of DAS data that focuses on acquiring and analyzing acoustic signals with lower frequencies (typically from a few hertz to tens of hertz)) can be extracted from it. Filtering and denoising algorithms are then used to remove downhole noise interference, retaining the effective low-frequency signals during fracture propagation, generating a clear strain rate waterfall plot, from which the measured strain rate data is obtained and denoted as... .

[0067] In some embodiments, the positive and negative strain rate scaling factors in S103 above may include a positive strain rate scaling factor and a negative strain rate scaling factor; correspondingly, adjusting the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factors in S103 above may, in specific implementation, include: If the axial strain rate of the simulated measuring point is greater than or equal to zero, the axial strain rate of the simulated measuring point is adjusted based on the positive strain rate scaling factor, and the product of the positive strain rate scaling factor and the axial strain rate of the simulated measuring point is used as the adjusted axial strain rate of the simulated measuring point. If the axial strain rate of the simulated measuring point is less than zero, the axial strain rate of the simulated measuring point is adjusted based on the negative strain rate scaling factor. The product of the negative strain rate scaling factor and the axial strain rate of the simulated measuring point is used as the adjusted axial strain rate of the simulated measuring point.

[0068] Specifically, positive strain rate corresponds to the strain increase caused by crack tensile propagation, while negative strain rate corresponds to the strain decrease caused by crack slowdown or local contraction. The two have different physical causes and their amplitude deviation characteristics may differ (e.g., the amplitude of positive strain rate is too small and the amplitude of negative strain rate is too large). More accurate amplitude correction can be achieved by segmented scaling.

[0069] Due to the idealization of geological models, the amplitude of directly calculated simulated data may not match the actual data. Scaling factors can be used to correct the simulated data. As a new Substitute the loss function. The specific correction formula can be as follows: (14) in, This is the positive strain rate scaling factor; It is a negative strain rate scaling factor; The adjusted axial strain rate at the simulated measuring point.

[0070] In some embodiments, a loss function based on the Frobenius norm can be constructed to quantify the overall error across all time steps and all measurement points. Its formula can be as follows: (15) in, The Frobenius norm is used to take the square root of the sum of squares of the "difference between simulated and measured strain rates" across all time steps and all measurement points. This can be used to quantify the overall error (avoiding local interference from single time step / single measurement point errors); T is the total number of time steps (e.g., if the fracturing monitoring duration is 50s and the time interval is 5s, then T=10); N is the total number of fiber optic measurement points (e.g., if 20 measurement points are arranged in the monitoring well, then N=20). The loss value between the adjusted simulated axial strain rate and the measured axial strain rate can be calculated using the above formula (15), thus avoiding interference from the global optimal solution by the error of a single time step or measuring point.

[0071] S104: If the loss value is greater than the preset loss threshold, adjust the corresponding optimization control parameters until the loss value is less than or equal to the preset loss threshold, or until the preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, output the target crack geometry parameters for the entire time step.

[0072] In some embodiments, the fracture propagation sub-model described above can also output the fracturing fluid pressure in the injection hole; correspondingly, the adjustment of the corresponding optimization control parameters in S104 if the loss value is greater than a preset loss threshold can, in specific implementations, include: If the loss value is greater than the preset loss threshold, the error type is determined by combining the error characteristics of the injection hole fracturing fluid pressure and the adjusted simulated axial strain rate and the measured axial strain rate, and the corresponding optimized control parameters are adjusted according to the error type.

[0073] In some embodiments, the above-mentioned determination of the error type based on the error characteristics of the injection hole fracturing fluid pressure and the adjusted simulated axial strain rate at the measuring point compared with the measured axial strain rate at the measuring point may, in specific implementations, include: If the fracturing fluid pressure in the injection hole is within the preset engineering reasonable range, and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are characterized by an overall amplitude deviation, then it is determined to be an amplitude deviation. If the fracturing fluid pressure in the injection hole exceeds the preset reasonable engineering range, and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are manifested as deviations in waveform width, signal arrival time, or spatial distribution, then it is determined to be a spatiotemporal morphological deviation.

[0074] Specifically, the fracture height h can be processed using the above formula (10) to output the injection hole fracturing fluid pressure p(t) at time t. Its value can directly reflect the rationality of the fracture height (e.g., the larger the fracture height, the greater the flow resistance inside the fracture, and the higher the injection hole pressure). The loss value calculated by the above formula (15) can be compared with the preset loss threshold. If the loss value is greater than the preset loss threshold, the injection hole fracturing fluid pressure can be compared with the engineering reasonable range of the injection hole fracturing fluid pressure or the field measured wellhead pressure. Based on the comparison result and combined with the error characteristics between the adjusted simulated measuring point axial strain rate and the measured measuring point axial strain rate (e.g., the overall amplitude deviation or the waveform width, signal arrival time or spatial distribution deviation), the error type between the simulated measuring point axial strain rate and the measured measuring point axial strain rate (e.g., amplitude deviation or spatiotemporal morphology deviation) can be determined, and the corresponding optimized control parameters can be adjusted according to the error type.

[0075] The reasonable engineering range for the fracturing fluid pressure at the injection hole can be ±10% of the measured wellhead pressure. For example, if the measured wellhead pressure is 12 MPa, the reasonable range is 10.8~13.2 MPa. This range can be determined based on the allowable pressure fluctuation value during engineering construction to ensure the rationality of the constraints. If the fracturing fluid pressure at the injection hole is within the preset reasonable engineering range, and the error characteristic between the adjusted simulated axial strain rate and the measured axial strain rate is an overall amplitude deviation, then it is determined to be an amplitude deviation. If the fracturing fluid pressure at the injection hole exceeds the preset reasonable engineering range, and the error characteristic between the adjusted simulated axial strain rate and the measured axial strain rate is a deviation in waveform width, signal arrival time, or spatial distribution, then it is determined to be a spatiotemporal morphological deviation.

[0076] When the fracturing fluid pressure in the injection hole is within the preset reasonable engineering range, it indicates that the current fracture height setting conforms to the flow resistance and mechanical constraints of hydraulic fracturing, meaning the fracture geometry is physically reasonable. Under the premise of physically reasonable geometry, the simulated and measured strain rates only show "consistent spatiotemporal distribution with only overall amplitude deviation," indicating that the error does not stem from incorrect geometry but from idealized assumptions in the physical model (such as assuming homogeneous rock and ignoring differences in fracturing fluid loss). These assumptions lead to a systematic deviation in the amplitude between simulated and measured data. In this case, simply correcting the amplitude of the simulated data using positive and negative strain rate scaling factors is sufficient to match the simulated and measured data, without changing the physically reasonable geometry.

[0077] When the injection hole pressure exceeds the reasonable range for engineering applications, it indicates that the current fracture height setting violates the physical laws of hydraulic fracturing. For example: High pressure: Larger fracture height → Excessive flow resistance within the fracture → Slower fracture propagation rate. Low pressure: Smaller fracture height → Insufficient flow resistance within the fracture → Faster fracture propagation rate. Errors in geometric parameters directly lead to discrepancies between the fracture propagation rate and its actual spatial location, resulting in a mismatch between the simulated strain rate's "waveform width, signal arrival time, and spatial distribution" and the measured data. For example: Larger fracture height → Slower fracture propagation → Later arrival time of the strain rate peak and a narrower range of response measurement points. Incorrect fracture height offset → Deviated fracture spatial location → Strain rate peak measurement point does not match the measured value. In this case, simply adjusting the amplitude cannot resolve the contradiction at the physical law level. It is necessary to adjust the fracture height and offset to bring the geometric parameters back to the reasonable physical range in order to fundamentally correct the spatiotemporal deviation.

[0078] In some embodiments, the above-mentioned adjustment of the corresponding optimization control parameters according to the error type may, in specific implementation, include: If the error type is amplitude deviation, then adjust the positive and negative strain rate scaling factors; If the error type is spatiotemporal morphological deviation, then adjust the crack height and the crack offset in the crack height direction.

[0079] Specifically, if the loss value is greater than the preset loss threshold If the error type is amplitude deviation (the error is mainly manifested as "overall amplitude is too large / too small" (which can be judged by the energy ratio or peak ratio of the positive and negative parts)), then adjust the normal strain rate scaling factor according to the following formula. With negative strain rate scaling factor : (16) in, / This is the positive / negative strain rate scaling factor after the (k+1)th iteration update; / The positive / negative strain rate scaling factor for the k-th iteration can be used to correct for amplitude deviations in the simulated strain rate. / This is the direction coefficient for the positive / negative strain rate scaling factor (takes a value of +1 or -1, which controls the adjustment direction of the positive / negative strain rate scaling factor). If the simulated peak positive strain rate is smaller than the measured value, (Increase the normal strain rate scaling factor) If it is greater than the measured value, (Reduce the negative strain rate scaling factor), the same applies to the negative strain rate scaling factor; The adjustment step size for scaling the system can range from 0.1 to 0.3, and can be used to control the adjustment range and avoid overcorrection.

[0080] If the loss value is greater than the preset loss threshold If the error type is a spatiotemporal morphological deviation (the error mainly manifests as geometric mismatches such as "waveform too narrow / too wide", arrival time offset, and incorrect spatial position), then adjust the crack height h and the crack offset v in the crack height direction according to the following formula. h : (17) in, / This represents the offset of the crack height / crack height direction after the (k+1)th iteration update; / is the offset of crack height / crack height direction in the kth iteration, which is the core parameter affecting crack geometry; / The direction coefficient (value is +1 or -1) represents the offset in the height / crack height direction. If the injection hole pressure is too high, it indicates that the crack height is too large. (Reduce the height) If the pressure is too low (indicating that the crack height is too small). (Increase the height). If the simulated peak measurement point number is greater than the actual measured peak (peak value deviates to the right),... (Reduce the offset) If it is less than the actual measurement (peak value is off to the left). (Increase the offset); This is the adjustment step size for the geometric parameters, and its value can range from 0.2. 0.5m can be used to control and adjust the shape, avoiding excessive parameter fluctuations.

[0081] By constructing a dual judgment system of injection hole fracturing fluid pressure physical constraint + strain rate error feature classification, and combining it with the targeted adjustment logic of error type-optimized control parameters, an iterative optimization closed loop that is accurate, efficient and in line with engineering physics laws is formed. This can fundamentally solve the technical pain points of traditional fracture morphology inversion methods, such as blind iteration, difficulty in balancing accuracy and efficiency, and inversion results that violate the physical laws of hydraulic fracturing. Overall, it achieves multiple improvements in iterative efficiency, inversion accuracy, engineering practicality and result reliability.

[0082] In some embodiments, the loss value in the k-th iteration can be... When, or when the preset iteration stopping condition is reached (such as when the maximum number of iterations is reached). When the iteration stops, the target fracture geometry parameters (i.e., the real fracture geometry parameters obtained by inversion) can be output based on the final adjusted optimized control parameters for the entire time step. At the same time, the dynamic curves of the fracture geometry parameters (fracture length, fracture height, fracture width) evolving over time can also be output, along with the optimized full-time step fracturing fluid pressure in the injection hole, the adjusted axial strain rate at the simulated measuring points, and the convergence curve of the loss value.

[0083] The output fracture geometry parameters can be used to evaluate fracturing effectiveness, such as quantifying the fracture stimulation range and whether the aperture meets the requirements for proppant embedding and oil and gas flow, thus determining whether the fracturing operation has achieved the design goals. They can also be used to optimize fracturing parameters, such as adjusting subsequent fracturing injection flow rate and fracturing fluid viscosity based on fracture propagation dynamics (e.g., length growth rate, aperture variation patterns) to improve stimulation effects. Furthermore, they can be used to formulate subsequent development plans, such as guiding well network deployment (e.g., avoiding areas with water channeling risk) and predicting production capacity (calculating oil and gas flow capacity based on fracture parameters), providing core data support for optimizing development plans.

[0084] The output fracturing fluid pressure in the injection hole can be used to: compare with the measured wellhead pressure in the field to further verify the engineering rationality of the inversion results; if the pressure data match, it can be confirmed that the fracture parameters are closer to reality. The output adjusted simulated strain rate can be used for: quantitative evaluation of inversion accuracy (such as calculating peak ratio and energy ratio), and can also serve as benchmark data for subsequent model iteration optimization. The output loss value convergence curve can be used for: verifying the effectiveness of the technical solution, providing data support for the method's promotion, and can also serve as the basis for subsequent adjustment iteration stopping conditions and parameter adjustment step sizes.

[0085] Through the above-mentioned scheme, the present invention can significantly improve the accuracy of fracture geometry inversion, reduce computational costs, and is effectively applicable to complex downhole environments, providing key technical support for the accurate evaluation and optimization of fracturing effects.

[0086] The various embodiments in this specification are described in a progressive manner. Similar or identical parts between embodiments can be referred to interchangeably. Each embodiment focuses on describing the differences from other embodiments. For details, please refer to the foregoing descriptions of the relevant processing embodiments; they will not be repeated here.

[0087] The foregoing description of this method is for illustrative purposes only and describes specific embodiments. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps described in the claims may be performed in a different order than those shown in the embodiments and still achieve the desired results. Furthermore, the processes depicted in the drawings do not necessarily require the specific or sequential order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are possible or may be advantageous.

[0088] In a specific implementation scenario, refer to Figure 2 As shown, it intuitively demonstrates the physical basis and signal transmission link of the inversion method of this invention, which is the basic support for understanding the technical logic of "physical model forward modeling - strain rate error matching - crack parameter inversion" of this invention. Figure 2 A three-dimensional coordinate system was established with the fractured well as the origin: the X-axis represents the fracture extension direction (corresponding to fracture length), the Z-axis represents the vertical direction (corresponding to fracture height), and the Y-axis is horizontal and perpendicular to the fracture extension direction (corresponding to the well distance between the fractured well and the monitoring well). The grid region represents the three-dimensional fracture volume generated by fracturing, and its length and other geometric parameters are the core inversion objects of this invention.

[0089] Distributed optical fibers are laid along the inner wall of the monitoring well casing. When the crack expands, the surrounding rock elements undergo three-dimensional deformation (as shown in the enlarged box in the upper left corner). / / ) and rotation ( / / The deformation is transmitted to the casing through the cement ring, causing axial strain in the casing (as shown in the upper right magnified box). The optical fiber is bound to the casing and synchronously senses this strain, forming a signal transmission link of "crack propagation → rock deformation → casing strain → optical fiber sensing", which can provide a data source of "measured axial strain rate at the actual measurement point" for inversion.

[0090] Specific implementation examples are as follows: a. Generate synthetic fiber data: Fiber data is generated by coupling the crack propagation sub-model (PKN model) and the fiber strain response sub-model, which can be used to simulate distributed acoustic wave transmission inductive variable rate signals.

[0091] b. Constructing a forward model (i.e., the physical model mentioned above): Initial values ​​are set for the crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. The crack geometric parameters are obtained by simulating the crack propagation process using the PKN model, and the fiber strain rate data is calculated using the fiber strain response sub-model. A fiber strain waterfall plot is then plotted. The forward model is used to simulate and obtain the spatiotemporal evolution characteristics of displacement, strain, and strain rate during crack propagation.

[0092] c. Parameter Optimization and Iteration: An iterative optimization algorithm is used to input the strain rate data synthesized from the optical fiber and the simulated strain rate data calculated using the PKN model and the optical fiber strain response sub-model into the loss function to check if the accuracy requirements are met. If the accuracy requirements are not met, the model parameters are adjusted, and steps b and c are repeated to minimize the loss function until the calculated simulated strain rate meets the accuracy requirements. Specifically, the inversion algorithm automatically iterates to find the optimal solution: (1) Initialization: Based on the actual situation, set the parameter set for crack height h and crack offset v in the crack height direction. h With positive and negative strain rate scaling factors and initial value and step size With stop threshold .

[0093] Kth iteration: (2) Forward calculation: The input is fed into the physical model, and the final output is the simulated strain rate (simulated axial strain rate at the measuring point). .

[0094] (3) Amplitude correction: Adjusted based on positive and negative strain rate scaling factors That is, scaling according to the positive and negative signs respectively. .

[0095] (4) Calculate the loss: Compared with actual measurement By comparing and calculating the losses of the two, we can obtain... .

[0096] (5) When the amplitude is mismatched, adjust and :like If the error is mainly manifested as "overall amplitude being too large / too small" (which can be judged by the energy ratio or peak ratio of the positive and negative parts), then update it first:

[0097] in: It is determined by whether the simulation is "larger or smaller than the actual measurement".

[0098] After updating, return to steps (3)-(4) and recalculate. .

[0099] (6) When the shapes do not match, adjust h and v. h If in the current Below, after several times Even after the update If the error is mainly manifested as geometric mismatch such as "waveform too narrow / too wide, arrival time offset, incorrect spatial position", then update:

[0100] After updating, return to step (2) and repeat the forward operation, then repeat steps (3)-(6).

[0101] (7) Stopping condition: when or the maximum number of iterations is reached. The iteration stops when the parameter values ​​at this point are the true parameters obtained through inversion.

[0102] d. Output the inversion results: Output the crack geometry parameters (crack length, crack width) and their dynamic curves over time. The inversion algorithm then terminates. Specifically, after parameter optimization and iteration stop, the optimal crack height determined in step c is used. The PKN model is run again, and the final output is the fracture geometry parameter curve (fracture length and fracture width) and its dynamic curve of evolution over time, realizing the visualization of fracture morphology and providing data support for fracturing effect evaluation and optimization design.

[0103] e. Further model validation: Inversion was performed using field data to extract LF-DAS data, and filtering and denoising algorithms were used to remove downhole noise interference to complete data processing. The data inversion results were compared to further validate the accuracy of the inversion model.

[0104] See Figure 3 As shown, the horizontal axis represents the time step (corresponding to the time series of fracturing operations), and the vertical axis represents the axial strain rate (unit: ×10). -7 The curves show that the overall trend, peak position, and waveform width of the simulated axial strain rate and the measured axial strain rate are highly consistent. The amplitude deviation of the positive and negative strain rates is controlled within 10%, which verifies the amplitude correction effect of the positive and negative strain rate scaling coefficients in this invention and the calculation accuracy of the full-time step physical model, indicating that the simulated data has achieved high-precision matching with the measured data.

[0105] See Figure 4As shown, the horizontal axis represents the time step, and the vertical axis represents the fracture half-length (i.e., the fracture length mentioned above, in meters). The curve represents the monitored value of the actual fracture half-length, and the dots represent the fracture half-length values ​​obtained by the method of this invention. The curve shows that the fracture half-length obtained by the inversion exhibits a linear growth trend with the actual value, and the numerical deviation is less than 3%. This indicates that the logic of calculating the fracture length by correlating fracture height in this invention is effective, and the inversion results can accurately reflect the dynamic process of horizontal fracture propagation, providing a reliable basis for assessing the reservoir stimulation range.

[0106] See Figure 5 As shown, the horizontal axis represents the time step, and the vertical axis represents the fracture aperture (unit: mm). The curve represents the monitored value of the actual fracture aperture, and the dots represent the fracture aperture values ​​obtained by the method of this invention. The curve shows that the trend of the fracture aperture obtained by the inversion is completely consistent with the actual value. In the early stage, the aperture increases rapidly with the injection of fracturing fluid, and then tends to stabilize in the later stage, with an overall deviation of less than 5%. This result verifies the rationality of calculating the fracture aperture of the entire time step through the fracture propagation sub-model of this invention, and the accuracy of the fracture geometric parameters after iterative optimization.

[0107] This invention addresses key issues in existing distributed fiber optic sensing fracture inversion methods, such as difficulty in capturing strain rate, dynamic changes in fracture morphology, and high computational costs, by proposing a novel inversion method. The core of this method lies in fusing the classical PKN hydraulic fracturing physical model with full-time-step data fitting, constructing a new inversion approach. Instead of fitting data from a single time step in isolation, this method applies physical constraints to the entire dynamic process of fracture propagation through the PKN model. By simplifying the inversion parameters from a massive number of spatiotemporal variables to four key physical parameters (fracture height, location offset, and positive and negative strain rate scaling factors), this invention effectively improves computational efficiency without sacrificing accuracy. Experiments show that, whether using synthetic data or field data (HFTS), the fracture morphology and strain response inverted by this method highly match the actual situation, significantly outperforming traditional single-time-step inversion methods. Therefore, this invention provides a more efficient, accurate, and reliable technical tool for real-time evaluation and optimization design of fracturing effects in oil and gas fields, possessing significant engineering application value.

[0108] By using the PKN model as a constraint, a full-time-step fitting method was applied for the first time to the inversion of fracture geometry in hydraulic fracturing. This method overcomes the core problems of insufficient accuracy and high computational cost of existing inversion techniques, and provides an efficient, accurate and reliable method for inverting fracture geometry. It provides a solid data foundation for the accurate evaluation of hydraulic fracturing effects and subsequent optimization design.

[0109] Although this specification provides the following examples or appendices Figure 6The method or apparatus structure shown may include more or fewer combined operation steps or module units based on conventional or non-creative labor. In steps or structures where there is no logically necessary causal relationship, the execution order of these steps or the module structure of the apparatus is not limited to the execution order or module structure shown in the embodiments or drawings of this specification. When the method or module structure is applied in actual devices, servers, or terminal products, it can be executed sequentially or in parallel according to the method or module structure shown in the embodiments or drawings (e.g., in parallel processor or multi-threaded processing environments, or even distributed processing or server cluster implementation environments). Based on the above-described crack morphology inversion method combining full time step and physical model, this specification also proposes an embodiment of a crack morphology inversion apparatus combining full time step and physical model. Figure 6 As shown, the device may specifically include the following modules: The acquisition module 601 can be used to acquire optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. The output module 602 can be used to input the crack height and the crack offset in the crack height direction into the physical model. The physical model includes a crack propagation sub-model and an optical fiber strain response sub-model. First, the crack height is processed by the crack propagation sub-model to output the crack geometric parameters of the entire time step. Then, the crack geometric parameters and the crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate of the simulated measuring point of the entire time step. The crack geometric parameters include crack length and crack aperture. The loss value calculation module 603 can be used to adjust the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factor, and calculate the loss value between the adjusted axial strain rate of the simulated measuring point and the measured axial strain rate of the actual measuring point. The adjustment iteration module 604 can be used to adjust the corresponding optimization control parameters if the loss value is greater than the preset loss threshold, until the loss value is less than or equal to the preset loss threshold, or the preset iteration stop condition is reached, and output the target crack geometry parameters for the entire time step based on the finally adjusted optimization control parameters.

[0110] In some embodiments, the output module 602 described above can be specifically used to process the crack height according to the following formula and output the crack length at the full time step:

[0111] The crack height is processed according to the following formula to output the crack aperture at the full time step:

[0112] Where L(t) is the crack length at time t; G is the shear modulus; q is the flow rate across the crack cross section; and v is the Poisson's ratio of the rock. denoted as η, where η is the dynamic viscosity coefficient of the fracturing fluid; h is the fracture height; w(t) is the fracture aperture at time t; and a and b are constants.

[0113] In some embodiments, the output module 602 can also be used to correct the influence coefficients in the fiber strain response sub-model based on the crack length and the crack offset in the crack height direction, wherein the influence coefficients depend on the relative positional distance between the crack element and the fiber measuring point; determine the axial displacement of the fiber measuring point throughout the time step based on the corrected influence coefficients and the crack opening; set the fiber gauge length based on the crack length; determine the axial strain of the fiber measuring point throughout the time step based on the fiber gauge length and the axial displacement of the fiber measuring point; and determine and output the axial strain rate of the simulated measuring point throughout the time step based on the axial strain of the fiber measuring point and the simulation calculation time interval.

[0114] In some embodiments, the output module 602 described above can also be used to determine the axial displacement of the fiber optic measuring point throughout the time step according to the following formula:

[0115] The step of determining the axial strain of the fiber optic measuring point across the entire time step based on the fiber gauge length and the axial displacement of the fiber optic measuring point includes: The axial strain at the fiber optic measuring point throughout the entire time step is determined using the following formula:

[0116] The step of determining and outputting the axial strain rate of the simulated measuring point across the entire time step based on the axial strain of the fiber optic measuring point and the simulation calculation time interval includes: The axial strain rate at the simulated measurement points throughout the entire time step is determined using the following formula:

[0117] Among them, u f (t) represents the axial displacement of the fiber optic measuring point at time t; M represents the total number of crack elements; v represents the Poisson's ratio of the rock; I1 and I2 are the corrected influence coefficients; w i (t) represents the crack opening of the i-th crack element at time t; L represents the axial strain at the fiber optic measuring point at time t. g z is the fiber gauge length; z is the center position of the fiber measuring point; u f (t)(z+L g / 2)- u f (t)(zL g / 2) represents the relative displacement difference between the two ends of the fiber gauge length; Let be the axial strain rate at the simulated measuring point at time t; for Axial strain at the fiber optic measuring point at time t; This is for simulating the calculation time interval.

[0118] In some embodiments, the positive and negative strain rate scaling factors in the loss value calculation module 603 may include a positive strain rate scaling factor and a negative strain rate scaling factor. Accordingly, the loss value calculation module 603 may specifically be used to adjust the axial strain rate of the simulated measuring point based on the positive strain rate scaling factor if the axial strain rate of the simulated measuring point is greater than or equal to zero, and take the product of the positive strain rate scaling factor and the axial strain rate of the simulated measuring point as the adjusted axial strain rate of the simulated measuring point; if the axial strain rate of the simulated measuring point is less than zero, adjust the axial strain rate of the simulated measuring point based on the negative strain rate scaling factor, and take the product of the negative strain rate scaling factor and the axial strain rate of the simulated measuring point as the adjusted axial strain rate of the simulated measuring point.

[0119] In some embodiments, the fracture propagation sub-model also outputs the fracturing fluid pressure in the injection hole; correspondingly, the adjustment iteration module 604 can be specifically used to determine the error type and adjust the corresponding optimization control parameters according to the error type if the loss value is greater than a preset loss threshold, by combining the error characteristics of the fracturing fluid pressure in the injection hole and the axial strain rate of the adjusted simulated measuring point with the axial strain rate of the measured measuring point.

[0120] In some embodiments, the above-mentioned adjustment iteration module 604 can also be used to determine an amplitude deviation if the fracturing fluid pressure at the injection hole is within a preset engineering reasonable range and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are characterized by an overall amplitude deviation; and to determine a spatiotemporal morphology deviation if the fracturing fluid pressure at the injection hole exceeds the preset engineering reasonable range and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are characterized by deviations in waveform width, signal arrival time, or spatial distribution.

[0121] In some embodiments, the above-mentioned adjustment iteration module 604 can also be used to adjust the positive and negative strain rate scaling coefficients if the error type is amplitude deviation; and to adjust the crack height and the crack offset in the crack height direction if the error type is spatiotemporal morphology deviation.

[0122] As can be seen from the above, the crack morphology inversion device combining full time step and physical model provided in the embodiments of this specification can improve the accuracy of crack geometry inversion, significantly reduce computational costs and improve inversion efficiency.

[0123] This specification also provides an electronic device based on the above-described crack morphology inversion method combining full-time step and physical model, including a processor and a memory for storing processor-executable programs / instructions. Specifically, the processor can execute the following steps according to the program / instructions: obtaining optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors; inputting the crack height and crack offset in the crack height direction into a physical model, which includes a crack propagation sub-model and an optical fiber strain response sub-model; first processing the crack height through the crack propagation sub-model, and outputting the crack morphology inversion method at the full-time step. The parameters are then processed using a fiber optic strain response sub-model to determine the crack geometry parameters and the crack offset in the crack height direction, outputting the axial strain rate of the simulated measuring point at the full time step. The crack geometry parameters include crack length and crack aperture. The axial strain rate of the simulated measuring point is adjusted based on the positive and negative strain rate scaling factors, and the loss value between the adjusted simulated axial strain rate and the measured axial strain rate is calculated. If the loss value is greater than a preset loss threshold, the corresponding optimization control parameters are adjusted until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters at the full time step are output.

[0124] To execute the above instructions more accurately, please refer to... Figure 7 As shown in the embodiments of this specification, another specific electronic device is also provided, wherein the electronic device includes a network communication port 701, a processor 702, and a memory 703. The above structures are connected by internal cables so that the various structures can perform specific data interaction.

[0125] Specifically, the network communication port 701 can be used to acquire optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. The processor 702 is specifically used to input the crack height and the crack offset in the crack height direction into a physical model. The physical model includes a crack propagation sub-model and an optical fiber strain response sub-model. First, the crack height is processed by the crack propagation sub-model to output the crack geometric parameters at the full time step. Then, the crack geometric parameters and the crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate of the simulated measuring point at the full time step. The crack geometric parameters include crack length and crack aperture. The axial strain rate of the simulated measuring point is adjusted based on the positive and negative strain rate scaling factors. The loss value between the adjusted simulated axial strain rate and the measured axial strain rate is calculated. If the loss value is greater than a preset loss threshold, the corresponding optimization control parameters are adjusted until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometric parameters at the full time step are output.

[0126] The memory 703 can be used to store the corresponding instruction program.

[0127] In this embodiment, the network communication port 701 can be a virtual port bound to different communication protocols, thereby enabling the sending or receiving of different data. For example, the network communication port can be a port responsible for web data communication, a port responsible for FTP data communication, or a port responsible for email data communication. Furthermore, the network communication port can also be a physical communication interface or communication chip. For example, it can be a wireless mobile network communication chip, such as GSM or CDMA; it can also be a Wi-Fi chip; or it can be a Bluetooth chip.

[0128] In this embodiment, the processor 702 can be implemented in any suitable manner. For example, the processor can take the form of a microprocessor or processor and a computer-readable medium storing computer-readable program code (e.g., software or firmware) executable by the (micro)processor, logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers, and embedded microcontrollers, etc. This specification is not limiting.

[0129] In this embodiment, the memory 703 may include multiple layers. In a digital system, anything that can store binary data can be a memory. In an integrated circuit, a circuit with storage function but no physical form is also called a memory, such as RAM, FIFO, etc. In a system, a storage device with a physical form is also called a memory, such as a memory stick, TF card, etc.

[0130] This specification also provides a computer storage medium based on the above-described crack morphology inversion method combining full-time step and physical model. The computer storage medium stores a computer program / instruction that, when executed, performs the following: acquiring optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors; inputting the crack height and crack offset in the crack height direction into a physical model, which includes a crack propagation sub-model and an optical fiber strain response sub-model; first processing the crack height through the crack propagation sub-model to output the crack geometric parameters at the full-time step; and then... The fiber strain response sub-model processes the crack geometry parameters and the crack offset in the crack height direction, outputting the axial strain rate of the simulated measuring point at the full time step. The crack geometry parameters include crack length and crack aperture. The axial strain rate of the simulated measuring point is adjusted based on the positive and negative strain rate scaling factors, and the loss value between the adjusted simulated axial strain rate and the measured axial strain rate is calculated. If the loss value is greater than a preset loss threshold, the corresponding optimization control parameters are adjusted until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters at the full time step are output.

[0131] In this embodiment, the storage medium includes, but is not limited to, Random Access Memory (RAM), Read-Only Memory (ROM), Cache, Hard Disk Drive (HDD), or Memory Card. The memory can be used to store computer program instructions. The network communication unit can be an interface configured according to standards specified in the communication protocol for network connection communication.

[0132] In this embodiment, the specific functions and effects implemented by the program instructions stored in the computer storage medium can be explained in comparison with other implementation methods, and will not be repeated here.

[0133] While this specification provides the steps of operation for the methods described in the embodiments or flowcharts, more or fewer steps may be included based on conventional or non-inventive means. The order of steps listed in the embodiments is merely one possible order of execution among many steps and does not represent the only possible order. In actual device or client product execution, the methods shown in the embodiments or drawings may be executed sequentially or in parallel (e.g., in a parallel processor or multi-threaded processing environment, or even a distributed data processing environment). The terms "comprising," "including," or any other variations thereof are intended to cover a non-exclusive inclusion, such that a process, method, product, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, product, or apparatus. Without further limitations, the presence of other identical or equivalent elements in a process, method, product, or apparatus that includes said elements is not excluded. The terms "first," "second," etc., are used to denote names and do not indicate any particular order.

[0134] Those skilled in the art will also know that, besides implementing the controller using purely computer-readable program code, the same functions can be achieved by logically programming the method steps, making the controller function as logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers (PLCs), and embedded microcontrollers. Therefore, such a controller can be considered a hardware component, and the devices within it used to implement various functions can also be considered structures within that hardware component. Alternatively, the devices used to implement various functions can be considered as both software modules implementing the method and structures within a hardware component.

[0135] This specification can be described in the general context of computer-executable instructions that are executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, classes, etc., that perform a specific task or implement a specific abstract data type. This specification can also be practiced in distributed computing environments, where tasks are performed by remote processing devices connected via a communication network. In distributed computing environments, program modules can reside in local and remote computer storage media, including storage devices.

[0136] As can be seen from the above description of the embodiments, those skilled in the art can clearly understand that this specification can be implemented by means of software plus necessary general-purpose hardware platforms. Based on this understanding, the technical solutions of this specification can essentially be embodied in the form of a software product. This computer software product can be stored in a storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, mobile terminal, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments of this specification.

[0137] The various embodiments in this specification are described in a progressive manner. Similar or identical parts between embodiments can be referred to interchangeably. Each embodiment focuses on its differences from other embodiments. This specification can be used in numerous general-purpose or special-purpose computer system environments or configurations. Examples include: personal computers, server computers, handheld or portable devices, tablet devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable electronic devices, network PCs, minicomputers, mainframe computers, and distributed computing environments including any of the above systems or devices, etc.

[0138] Although this specification has been described by way of examples, those skilled in the art will recognize that many variations of this specification are possible without departing from its spirit, and it is intended that the appended claims cover such variations without departing from the spirit of this specification.

Claims

1. A crack morphology inversion method combining full-time step and physical model, characterized in that, include: Obtain optimized control parameters, including crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors; The crack height and the crack offset in the crack height direction are input into the physical model, which includes a crack propagation sub-model and an optical fiber strain response sub-model. First, the crack height is processed by the crack propagation sub-model to output the crack geometric parameters at the full time step. Then, the crack geometric parameters and the crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate at the simulated measurement point at the full time step. The crack geometric parameters include crack length and crack aperture. The axial strain rate of the simulated measuring point is adjusted based on the positive and negative strain rate scaling factors, and the loss value between the adjusted axial strain rate of the simulated measuring point and the measured axial strain rate of the measuring point is calculated. If the loss value is greater than the preset loss threshold, the corresponding optimization control parameters are adjusted until the loss value is less than or equal to the preset loss threshold, or the preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters for the entire time step are output.

2. The method according to claim 1, characterized in that, The process of processing the crack height through a crack propagation sub-model to output the crack geometry parameters at the full time step includes: The crack height is processed according to the following formula to output the crack length at the full time step: The crack height is processed according to the following formula to output the crack aperture at the full time step: Where L(t) is the crack length at time t; G is the shear modulus; q is the flow rate across the crack cross section; and v is the Poisson's ratio of the rock. denoted as η, where η is the dynamic viscosity coefficient of the fracturing fluid; h is the fracture height; w(t) is the fracture aperture at time t; and a and b are constants.

3. The method according to claim 1, characterized in that, The process of processing the crack geometry parameters and the crack offset in the crack height direction using a fiber optic strain response sub-model to output the axial strain rate at the simulated measurement points across the entire time step includes: Based on the crack length and the crack offset in the crack height direction, the influence coefficients in the fiber strain response sub-model are corrected. The influence coefficients depend on the relative positional distance between the crack element and the fiber measuring point. The axial displacement of the fiber optic measuring point at the full time step is determined based on the corrected influence coefficient and the crack opening. The fiber gauge length is set based on the crack length, and the axial strain of the fiber measuring point is determined at the full time step according to the fiber gauge length and the axial displacement of the fiber measuring point. Based on the axial strain at the fiber optic measuring point and the simulation calculation time interval, the axial strain rate of the simulated measuring point at the full time step is determined and output.

4. The method according to claim 3, characterized in that, The determination of the axial displacement of the fiber optic measuring point across the entire time step based on the corrected influence coefficient and the crack aperture includes: The axial displacement of the fiber optic measuring point throughout the entire time step is determined using the following formula: The step of determining the axial strain of the fiber optic measuring point across the entire time step based on the fiber gauge length and the axial displacement of the fiber optic measuring point includes: The axial strain at the fiber optic measuring point throughout the entire time step is determined using the following formula: The step of determining and outputting the axial strain rate of the simulated measuring point across the entire time step based on the axial strain of the fiber optic measuring point and the simulation calculation time interval includes: The axial strain rate at the simulated measurement points throughout the entire time step is determined using the following formula: Among them, u f (t) represents the axial displacement of the fiber optic measuring point at time t; M represents the total number of crack elements; v represents the Poisson's ratio of the rock; I1 and I2 are the corrected influence coefficients; w i (t) represents the crack opening of the i-th crack element at time t; L represents the axial strain at the fiber optic measuring point at time t. g z is the fiber gauge length; z is the center position of the fiber measuring point; u f (t)(z+L g / 2)- u f (t)(zL g / 2) represents the relative displacement difference between the two ends of the fiber gauge length; Let be the axial strain rate at the simulated measuring point at time t; for Axial strain at the fiber optic measuring point at time t; This is for simulating the calculation time interval.

5. The method according to claim 1, characterized in that, The positive and negative strain rate scaling factors include a positive strain rate scaling factor and a negative strain rate scaling factor; correspondingly, adjusting the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factors includes: If the axial strain rate of the simulated measuring point is greater than or equal to zero, the axial strain rate of the simulated measuring point is adjusted based on the positive strain rate scaling factor, and the product of the positive strain rate scaling factor and the axial strain rate of the simulated measuring point is used as the adjusted axial strain rate of the simulated measuring point. If the axial strain rate of the simulated measuring point is less than zero, the axial strain rate of the simulated measuring point is adjusted based on the negative strain rate scaling factor. The product of the negative strain rate scaling factor and the axial strain rate of the simulated measuring point is used as the adjusted axial strain rate of the simulated measuring point.

6. The method according to claim 1, characterized in that, The fracture propagation sub-model also outputs the fracturing fluid pressure in the injection hole; correspondingly, if the loss value is greater than a preset loss threshold, the corresponding optimization control parameters are adjusted, including: If the loss value is greater than the preset loss threshold, the error type is determined by combining the error characteristics of the injection hole fracturing fluid pressure and the adjusted simulated axial strain rate and the measured axial strain rate, and the corresponding optimized control parameters are adjusted according to the error type.

7. The method according to claim 6, characterized in that, The error type is determined by combining the error characteristics of the fracturing fluid pressure at the injection hole and the adjusted axial strain rate at the simulated measuring point with the measured axial strain rate at the actual measuring point, including: If the fracturing fluid pressure in the injection hole is within the preset engineering reasonable range, and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are characterized by an overall amplitude deviation, then it is determined to be an amplitude deviation. If the fracturing fluid pressure in the injection hole exceeds the preset reasonable engineering range, and the error characteristics between the adjusted simulated axial strain rate and the measured axial strain rate are manifested as deviations in waveform width, signal arrival time, or spatial distribution, then it is determined to be a spatiotemporal morphological deviation.

8. The method according to claim 7, characterized in that, The step of adjusting the corresponding optimization control parameters according to the error type includes: If the error type is amplitude deviation, then adjust the positive and negative strain rate scaling factors; If the error type is spatiotemporal morphological deviation, then adjust the crack height and the crack offset in the crack height direction.

9. A crack morphology inversion device combining full-time step and physical model, characterized in that, include: The acquisition module is used to acquire optimized control parameters, which include crack height, crack offset in the crack height direction, and positive and negative strain rate scaling factors. The output module is used to input the crack height and the crack offset in the crack height direction into the physical model. The physical model includes a crack propagation sub-model and an optical fiber strain response sub-model. First, the crack height is processed by the crack propagation sub-model to output the crack geometric parameters at the full time step. Then, the crack geometric parameters and the crack offset in the crack height direction are processed by the optical fiber strain response sub-model to output the axial strain rate of the simulated measurement point at the full time step. The crack geometric parameters include crack length and crack aperture. The loss value calculation module is used to adjust the axial strain rate of the simulated measuring point based on the positive and negative strain rate scaling factor, and calculate the loss value between the adjusted axial strain rate of the simulated measuring point and the measured axial strain rate of the actual measuring point. The adjustment iteration module is used to adjust the corresponding optimization control parameters if the loss value is greater than a preset loss threshold, until the loss value is less than or equal to the preset loss threshold, or a preset iteration stop condition is reached. Based on the finally adjusted optimization control parameters, the target crack geometry parameters for the entire time step are output.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by a processor, they implement the steps of the method according to any one of claims 1 to 8.

Images