Biological tissue viscoelastic optical coherence elastography method, system and equipment

Abstract

The invention discloses a biological tissue viscoelasticity optical coherence elastography method, system and equipment, and relates to the technical field of biological tissue viscoelasticity measurement. Non-contact biological tissue viscoelastic mechanical property detection is achieved through acoustic radiation force optical coherence elastography, the device can be applied to application scenes with biopsy requirements, distortion of a sound pressure field is reduced by applying the focal length-diameter ratio of the ultrasonic transducer and the frequency of an excitation signal, and the detection accuracy of the viscoelastic mechanical property of the biological tissue is improved. The capturing capability of tiny deformation of the biological tissue is improved, and the imaging quality and accuracy of the obtained displacement data of the biological tissue are improved, so that the accuracy of viscoelasticity calculation of the biological tissue is improved; meanwhile, viscoelasticity calculation is carried out on the phase velocity frequency dispersion curve of the biological tissue based on the surface acoustic wave frequency dispersion model, main energy distribution corresponding to elastic waves can be calculated, the detection and recognition effect is improved, and the accuracy of viscoelasticity calculation of the biological tissue is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of viscoelasticity measurement of biological tissue, and particularly relates to a biological tissue viscoelasticity optical coherence elastography method, system and device. BACKGROUND

[0002] At present, in the field of biological tissue disease diagnosis and treatment and bionic tissue preparation, the measurement of mechanical properties of biological tissue is of great significance.

[0003] In the prior art, traditional tensile testing can be used to evaluate the mechanical properties of the tissue, but it can only measure the overall elastic modulus of the ex vivo tissue, cannot perform biopsy and cannot test the layered structure of the tissue. In recent years, ultrasound elastography (UE), magnetic resonance imaging (MRI), laser speckle elastography (LSE) and computed tomography (CT) cannot perform accurate measurement due to slow imaging speed and low resolution.

[0004] ‌Optical coherence elastography (OCE) technology can distinguish tissues by measuring the microscopic elastic properties of the tissues, and the development of the detection method can quantitatively evaluate the elastic mechanical properties of the biological tissues.

[0005] However, the characteristics of the elastic waves generated under different excitation conditions will be different, and the existing optical coherence elastography technology has deficiencies in the aspects of elastic wave excitation, detection and identification, resulting in inaccurate viscoelasticity measurement of biological tissues. SUMMARY

[0006] Therefore, it is necessary to provide a biological tissue viscoelasticity optical coherence elastography method, system and device in view of the above technical problems.

[0007] The present application adopts the following technical solutions: The present application provides a biological tissue viscoelasticity optical coherence elastography method, comprising: obtaining displacement data of an elastic wave propagation process of a target biological tissue under acoustic radiation force; the acoustic radiation force is generated by an ultrasonic transducer with a sinusoidal signal of 0.1-1 kHz as an excitation signal, and the focal length-diameter ratio of the ultrasonic transducer is 0.8; performing two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue to determine the phase velocity dispersion curve of the target biological tissue; Based on the dispersion model of surface acoustic waves, the phase velocity dispersion curve of the target biological tissue is fitted by a fitting algorithm to obtain the elastic modulus and shear viscosity of the target biological tissue.

[0008] Optionally, the displacement data of the target biological tissue includes two-dimensional spatiotemporal slice data of the surface of the biological tissue, which is used to represent the displacement-time information of the surface of the target biological tissue. The step of performing a two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue to determine the phase velocity dispersion curve of the target biological tissue specifically includes: The displacement data of the target biological tissue is transformed from the spatiotemporal domain to the frequency-wavenumber domain by performing a two-dimensional discrete Fourier transform using the following formula: ; The maximum intensity value of the displacement data of the target biological tissue in the frequency-wavenumber domain is determined, and the signal below 10% of the maximum intensity value in the displacement data of the target biological tissue in the frequency-wavenumber domain is filtered out by threshold filtering to remove background noise. Based on the maximum surface acoustic wave number corresponding to each frequency in the displacement data of the target biological tissue in the filtered frequency-wavenumber domain, the phase velocity value at each frequency is determined by the following formula, thus obtaining the phase velocity dispersion curve of the target biological tissue: ; in, The energy distribution field of the target biological tissue in the frequency-wavenumber domain. For the time-space-displacement field of the target biological tissue surface, For the spatial coordinates of the elastic wave, For the time coordinate of the elastic wave, For wave number, For frequency, For spatial discrete sampling index, For time-discrete sampling index, The imaginary unit, The phase velocity value below, For frequency The corresponding maximum surface acoustic wave number.

[0009] Optionally, the dispersion model based on surface acoustic waves uses a fitting algorithm to fit the phase velocity dispersion curve of the target biological tissue to obtain the elastic modulus and shear viscosity of the target biological tissue, specifically including: Determine a fitting frequency interval according to the preset minimum frequency and the frequency interval length, fit the data of the phase velocity dispersion curve of the biological tissue in the fitting frequency interval to the acoustic surface wave dispersion model under viscoelastic, isotropic and uniform material, determine the initial range of the shear wave wavenumber according to the acoustic surface wave wavenumber, and find the optimal solution of the shear wave wavenumber satisfying the acoustic surface wave dispersion equation in the initial range of the shear wave wavenumber through traversal; According to the optimal solution of the shear wave wavenumber and the relationship between the shear wave wavenumber and the shear modulus and the shear viscosity, nonlinear least squares curve fitting is performed, the fitting correlation coefficient of the nonlinear least squares curve fitting result and the optimal solution of the shear wave wavenumber is better than the preset fitting correlation coefficient, and the elastic modulus and the shear viscosity of the target biological tissue are obtained.

[0010] Optionally, the target biological tissue is a dermis layer of skin, and the fitting frequency interval is 300-600 Hz.

[0011] The present application provides a kind of biological tissue viscoelastic optical coherence elastography system, comprising: Acoustic radiation force excitation module, comprising: ultrasonic transducer, for under the excitation of excitation signal, acoustic radiation force is generated, and the target biological tissue generates elastic wave;The focal length diameter ratio of the ultrasonic transducer is 0.8, and the sine signal of 0.1-1kHz is used as excitation signal; OCT acquisition module, for acquiring the displacement data of the elastic wave propagation process of the target biological tissue under the acoustic radiation force; Dispersion curve determination module, for carrying out two-dimensional discrete Fourier transform and threshold filtering to the displacement data of the target biological tissue, to determine the phase velocity dispersion curve of the target biological tissue; Inversion module, for fitting the phase velocity dispersion curve of the target biological tissue based on acoustic surface wave dispersion model through fitting algorithm, to obtain the elastic modulus and shear viscosity of the target biological tissue.

[0012] The present application provides a kind of computer readable storage medium, the storage medium stores computer program, the computer program is executed by processor and realizes the above-mentioned biological tissue viscoelastic optical coherence elastography method.

[0013] The present application provides a kind of computer equipment, including memory, processor and computer program stored on memory and can be run on processor, the processor realizes the above-mentioned biological tissue viscoelastic optical coherence elastography method when the program is executed.

[0014] The above-mentioned at least one technical scheme adopted by the present application can achieve the following beneficial effects: The application adopts acoustic radiation force optical coherence elastography to realize non-contact detection of the viscoelastic mechanical properties of biological tissues, can be applied to application scenarios with biopsy requirements, and by applying the focal length diameter ratio of the ultrasonic transducer and the frequency of the excitation signal, the distortion of the acoustic pressure field is reduced, the capture ability of the micro deformation of the biological tissue is improved, the imaging quality and accuracy of the obtained displacement data of the biological tissue are improved, and the accuracy of the viscoelastic calculation of the biological tissue is improved; at the same time, based on the viscoelastic calculation of the phase velocity dispersion curve of the biological tissue based on the acoustic surface wave dispersion model, the energy distribution corresponding to the elastic wave can be calculated, the detection and recognition effect is improved, and the accuracy of the viscoelastic calculation of the biological tissue is improved. BRIEF DESCRIPTION OF DRAWINGS

[0015] The drawings described herein are intended to provide further understanding of the present application, and form a part of the present application.

[0016] Figure 1 A biological tissue viscoelastic optical coherence elastography method provided by the present application is shown in the flowchart; Figure 2 The structure of an ARF-OCE system provided by the present application is shown in the schematic diagram; Figure 3 The process of establishing an ultrasonic transducer set model provided by the present application is shown in the schematic diagram; Figure 4 The surface of a concave transducer and its focusing area provided by the present application are shown in the schematic diagram; Figure 5 Different excitation signals provided by the present application are shown in the schematic diagram; Figure 6 A normalized acoustic pressure field provided by the present application is shown in the schematic diagram; Figure 7 A wave field simulation idea provided by the present application is shown in the schematic diagram; Figure 8 A three-layer model provided by the present application is shown in the schematic diagram; Figure 9 The acoustic radiation force field of a modulation signal and the force loading model of a single-layer acoustic field provided by the present application are shown in the schematic diagram; Figure 10 The propagation of transient elastic waves in a solid provided by the present application is shown in the schematic diagram; Figure 11 The grid division corresponding to different models provided by the present application is shown in the schematic diagram; Figure 12 A cloud image of an elastic wave provided by the present application is shown in the schematic diagram; Figure 13A model surface particle vibration characteristic schematic diagram provided by the present application; Figure 14 A different model displacement cloud map normalization processing schematic diagram provided by the present application; Figure 15 A elastic wave Fourier transform schematic diagram provided by the present application; Figure 16 A displacement-time gradient diagram of each model under different frequencies provided by the present application; Figure 17 A dispersion curve schematic diagram on the surface of different models provided by the present application; Figure 18 A phase velocity dispersion curve schematic diagram of an elastic tissue without viscosity provided by the present application; Figure 19 A phase velocity dispersion curve schematic diagram of an elastic tissue with viscosity provided by the present application; Figure 20 An experimental instrument specific parameter schematic diagram provided by the present application; Figure 21 An OCT acquisition mode schematic diagram provided by the present application; Figure 22 A dispersion curve iteration specific flow schematic diagram provided by the present application; Figure 23 An experimental result schematic diagram provided by the present application; Figure 24 An OCT system tracking elastic wave result schematic diagram provided by the present application; Figure 25 A single-layer phantom experimental result schematic diagram provided by the present application; Figure 26 A double-layer phantom experimental result schematic diagram provided by the present application; Figure 27 A three-layer phantom experimental result schematic diagram provided by the present application; Figure 28 A fitted stress-strain curve schematic diagram provided by the present application; Figure 29 An in-vitro pig biological tissue experiment schematic diagram provided by the present application; Figure 30 A biological tissue viscoelastic optical coherence elastography system schematic diagram provided by the present application; Figure 31 A computer device for realizing a biological tissue viscoelastic optical coherence elastography method schematic diagram provided by the present application. DETAILED DESCRIPTION

[0017] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions of the present application will be described clearly and completely below in combination with specific embodiments of the present application and corresponding drawings. Obviously, the described embodiments are only some of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative labor fall within the scope of protection of the present application.

[0018] At present, in the excitation of biological tissue elastic waves, the characterization of the biomechanical properties of biological tissues based on the mechanical loading of the OCE technology lays a foundation, but this loading method inevitably contacts with soft tissues, limiting the experimental test of fragile tissues, so it is necessary to develop a non-contact excitation method.

[0019] In the detection and identification of elastic waves, it is necessary to use appropriate wave propagation models and effective identification methods to evaluate the structure and mechanical properties of soft tissues using the results of elastic wave research. The characteristics of elastic waves generated under different excitation conditions may be different, and how to effectively detect elastic waves and determine the specific type of elastic waves is crucial for viscoelasticity measurement based on the elastic wave model, but few studies focus on the modeling and simulation of elastic waves in layered soft tissues, so it is necessary to develop an effective method for detecting and identifying elastic waves.

[0020] In the method of biological tissue viscoelasticity parameter inversion based on elastic waves, there are many different energy and non-overlapping frequency range wave modes in layered tissues, and the dispersion curve (relationship between propagation speed and frequency) depends on the thickness and elastic modulus of each layer. There are few studies from optical coherence tomography (OCT) raw data to viscoelasticity characterization, and many problems need to be solved in how to calculate the dispersion curve and select the frequency band to accurately obtain the viscoelasticity parameters.

[0021] The technical solutions provided by the embodiments of the present application will be described in detail below in combination with the drawings.

[0022] Figure 1 The flowchart of the biological tissue viscoelasticity optical coherence elastic imaging method in the present application specifically comprises the following steps: S101: Obtain displacement data of the propagation process of elastic waves of the target biological tissue under acoustic radiation force; the acoustic radiation force is generated by an ultrasonic transducer with a sinusoidal signal of 0.1-1 kHz as an excitation signal, and the focal length diameter ratio of the ultrasonic transducer is 0.8.

[0023] S102: Perform two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue to determine the phase velocity dispersion curve of the target biological tissue.

[0024] S103: Based on the dispersion model of the surface acoustic wave, the phase velocity dispersion curve of the target biological tissue is fitted by a fitting algorithm, and the elastic modulus and shear viscosity of the target biological tissue are obtained.

[0025] For the convenience of description, the following will be described taking the server as the execution subject. The server mentioned in the present application can be a server arranged in a service platform, or a device such as a desktop computer, a notebook computer, etc. capable of executing the scheme of the present application.

[0026] Generally, when performing viscoelastic optical coherence elastography on a target biological tissue, the focal length diameter ratio of the selected ultrasonic transducer and the excitation signal used will directly affect the characteristics of the acoustic pressure field. In one or more embodiments of the present application, an ultrasonic transducer using a 0.1-1 kHz sinusoidal signal as the excitation signal can be used to generate acoustic radiation force, and the focal length diameter ratio of the ultrasonic transducer can be 0.8, so as to obtain the displacement data of the elastic wave propagation process of the target biological tissue under the acoustic radiation force.

[0027] By applying the focal length diameter ratio of the ultrasonic transducer and the frequency of the excitation signal, the distortion of the acoustic pressure field can be reduced, the capture ability of the micro-deformation of the biological tissue can be improved, the imaging quality and accuracy of the obtained displacement data of the biological tissue can be improved, and thus the accuracy of the viscoelasticity calculation of the biological tissue can be improved.

[0028] Further, the viscoelasticity of the biological tissue can be calculated based on the dispersion model of the surface acoustic wave. The viscoelasticity calculation of the phase velocity dispersion curve of the biological tissue based on the dispersion model of the surface acoustic wave can calculate the energy distribution corresponding to the elastic wave, improve the detection and recognition effect, and improve the accuracy of the viscoelasticity calculation of the biological tissue.

[0029] Specifically, in one or more embodiments of the present application, the displacement data of the target biological tissue includes two-dimensional space-time slice data of the surface of the biological tissue, and the two-dimensional space-time slice data of the surface of the target biological tissue is used to represent the displacement-time information of the surface of the target biological tissue. When performing two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue to determine the phase velocity dispersion curve of the target biological tissue, the server can first perform two-dimensional discrete Fourier transform by the following formula to convert it from the space-time domain to the frequency-wave number domain: .

[0030] In the formula, is the energy distribution field of the target biological tissue in the frequency-wave number domain, is the time-space-displacement field of the surface of the target biological tissue, is the spatial coordinate of the elastic wave, is the time coordinate of the elastic wave, is the wave number, is a frequency, is a spatial discrete sampling index, is a time discrete sampling index, is an imaginary unit.

[0031] Then the maximum intensity value of the displacement data of the target biological tissue in the frequency-wave number domain is determined (the higher the intensity of a point (f_0, k_0), the higher the proportion of the component with the frequency f_0 and the wave number k_0), and the signals lower than 10% of the maximum intensity value in the displacement data of the target biological tissue in the frequency-wave number domain are filtered out by threshold filtering to remove background noise.

[0032] Thus, according to the maximum acoustic surface wave wave number corresponding to each frequency in the filtered displacement data of the target biological tissue in the frequency-wave number domain, the phase velocity value at each frequency is determined by the following formula to obtain the phase velocity dispersion curve of the biological tissue: . In the formula, is the phase velocity value at the frequency is the phase velocity value at the frequency is the phase velocity value at the frequency corresponding to the maximum acoustic surface wave wave number.

[0033] Then, the fitting frequency interval can be further determined according to the preset minimum frequency and the frequency interval length, the data of the phase velocity dispersion curve of the target biological tissue in the fitting frequency interval is fitted to the acoustic surface wave dispersion model under the viscoelastic, isotropic and uniform material, the initial range of the shear wave wave number is determined according to the acoustic surface wave wave number, and the optimal solution of the shear wave wave number satisfying the acoustic surface wave dispersion equation is found by iteration in the initial range of the shear wave wave number.

[0034] Finally, according to the optimal solution of the shear wave wave number and the relationship between the shear wave wave number and the shear modulus and the shear viscosity, a nonlinear least squares curve fitting is performed, the fitting correlation coefficient of the nonlinear least squares curve fitting result and the optimal solution of the shear wave wave number is better than the preset fitting correlation coefficient, and the elastic modulus and the shear viscosity of the target biological tissue are obtained.

[0035] The present application mainly optimizes the excitation parameters, determines the Rayleigh wave (i.e. acoustic surface wave) as the main elastic wave mode through elastic wave excitation and propagation simulation, integrates the ARF-OCE system based on this, and in actual application, the biological tissue viscoelastic optical coherent elastic imaging can be performed based on the ARF-OCE system.

[0036] The present application is further introduced from the following aspects, mainly including: Optimization of excitation parameters: Field II software is used to model the ultrasonic transducer with different geometric parameters, simulate the ultrasonic sound field under different excitation signals, and optimize the ultrasonic transducer with a focal length diameter ratio of about 0.8, and use a 0.1-1 kHz sinusoidal signal as the excitation signal to determine the optimal excitation control condition.

[0037] Elastic wave excitation and propagation simulation: based on COMSOL software, single-layer, double-layer, three-layer and viscoelastic biological tissue finite element models are constructed, the optimized excitation parameters are input, the propagation law of elastic wave in layered biological tissue structure is simulated, and through time-space-frequency multi-dimensional analysis method, the elastic wave is classified and identified, and the Rayleigh wave is determined as the main elastic wave mode.

[0038] Integration of experimental system: integrate the ARF-OCE system, including the OCT acquisition module and the acoustic radiation force excitation module, wherein the excitation module uses the optimized ultrasonic transducer and excitation signal, the OCT acquisition module uses the M-B mode to collect data, and realizes the imaging and elastic wave excitation detection of the layered microstructure of biological tissue. Figure 2 It is a structural schematic diagram of an ARF-OCE system in the application.

[0039] Dispersion curve calculation: two-dimensional discrete Fourier transform and threshold filtering are performed on the displacement data collected by OCT, and the phase velocity dispersion curve of Rayleigh wave is calculated.

[0040] Viscoelastic parameter inversion: based on the Rayleigh wave dispersion model, the fitting algorithm is used to fit the dispersion curve, and the elastic modulus and shear viscosity of the biological tissue are calculated.

[0041] (I) Finite element simulation and analysis of biological tissue: Simulation of acoustic radiation force field of ultrasonic transducer: 1. In Field II software, a set of ultrasonic transducer models is established, Figure 3 It is a process diagram for establishing a set of ultrasonic transducer models in the application.

[0042] 1) System initialization, mainly the basic settings of sampling frequency, ultrasonic transducer type, array element size and sound velocity used in the simulation process. Center frequency and focal length are user-defined parameters, the sampling frequency is set to 100MHz, the sound velocity is 1540m / s; then use the function xdc_concave to set the transducer as a concave focusing transducer, 1000 array elements are distributed, and the focusing area is set to a square area of 30mm*30mm, and 1000*1000 points are sampled. When the radius of the transducer is set to 20mm and the focal length is set to 45mm, the surface of the concave transducer and its focusing area can be obtained as Figure 4 , Figure 4Fig. 1 is a schematic diagram of a concave transducer surface and its focal region in the present application. The square region on the top is the focal region, and the bottom is the surface of the concave transducer.

[0043] 2) Design of excitation signal of transducer. The whole sensor is excited by the same excitation signal, which is a user-defined signal using the function xdc_excitation(). Taking the center frequency of 1 MHz as an example, Figure 5 Fig. 2 is a schematic diagram of different excitation signals in the present application, showing three different excitation signals, namely a Gaussian pulse signal, a modulated 1000 Hz continuous signal, and a sine windowed signal. In order to improve performance, the sine signal is windowed using a Hanning window.

[0044] 3) Calculation of sound pressure field. The distribution of sound field pressure at each position in the transducer space is calculated using the function clac_hp(). In order to more clearly display the sound pressure distribution under different parameters with images, the obtained sound pressure values are logarithmically compressed.

[0045] 2, simulation analysis of sound pressure field and selection of excitation parameters.

[0046] Three groups of focal length diameter ratios (F / D) under the condition of normalized sound pressure field (frequency of 1 MHz, Gaussian pulse) are simulated, and the results are shown in Figure 6 Fig. 3. Figure 6 Fig. 3 is a schematic diagram of simulation results of sound pressure field under different conditions in the present application. Figure 6 Figs. 3(a) to 3(c) correspond to the sound pressure field diagrams of F / D of 1.25, 1 and 0.8 respectively under the condition of Gaussian pulse excitation signal with a frequency of 1 MHz, Figure 6 Figs. 3(d) to 3(f) correspond to the sound pressure field diagrams obtained under different conditions when F / D=0.8, wherein Fig. 3(d) corresponds to the condition of Gaussian pulse excitation signal with a frequency of 0.5 MHz, Fig. 3(e) corresponds to the condition of sine wave excitation signal with a frequency of 1 MHz, and Fig. 3(f) corresponds to the condition of 1000 Hz sine wave excitation signal modulated by 1 MHz as a carrier. It can be seen from Figure 6 Figs. 3(a) to 3(c) that as F / D decreases, the focal point size gradually decreases. When F / D is about 0.8, the ideal imaging effect can be achieved. It can be seen from Figure 6 Figs. 3(c) and 3(d) that the lateral size of the focal point is about 2 mm and 1 mm respectively at 0.5 MHz and 1 MHz. This is consistent with the theoretical expectation that the higher the center frequency, the smaller the focal point. Considering the balance between resolution and penetration depth, a frequency of 1 MHz is sufficient to meet the requirements. It can be seen from Figure 6 Fig. 3(f) that under the condition of 1000 Hz modulated sine signal as excitation signal, the focal point size is more ideal, and the focal point shape is closer to an ellipse.

[0047] In the field of medical ultrasound imaging or material testing, high resolution and shallow focal depth are often required. Larger transducer diameter and shorter focal length help to produce smaller focal points, thus achieving a more focused beam. This configuration can improve spatial resolution, reduce the effects of tissue scattering and signal attenuation, and make the details of the target structure more clearly identifiable. Therefore, through this study, the ratio of F / D=0.8 is selected, and the carrier wave of 1MHz frequency is modulated into a sinusoidal wave of 1000Hz to generate the required acoustic radiation force field. Such parameter configuration aims to optimize the characteristics of the acoustic field to meet the needs of high-precision imaging, and provides a theoretical and experimental basis for the application of ultrasonic imaging technology.

[0048] In summary, the higher the center frequency, the smaller the acoustic field focal point, and the smaller the distortion of the sound pressure field when the focal length diameter ratio is near 0.8 and the sinusoidal modulation signal is used as the excitation signal. This lays a solid foundation for subsequent formulation of the excitation conditions of ARF-OCE experiments.

[0049] The basic process of wave field simulation based on COMSOL: In the present application, the wave field simulation in the solid is mainly realized by COMSOL Multiphysics simulation software. The solid mechanics module is based on the elastic wave equation to describe the propagation of acoustic waves in solid medium, as shown in the following formula, which describes the wavelength, frequency and propagation speed of acoustic waves in solid, etc. information: .

[0050] In the formula, the elastic dynamic equation uses tensor subscript notation: the displacement field component The comma in the subscript indicates the spatial partial derivative, the repeated subscript is summed according to Einstein's summation convention, and the free subscript represents the coordinate component independently; represents the gradient of the displacement divergence in the direction, represents the Laplacian of , represents the second-order time derivative of and represent the Lame constant, represents the external force component in the Cartesian coordinate system, represents the density of the medium.

[0051] In the simulation process, the velocity-strain formula is used to solve the control equation of general linear elastic material: ; ; ; Among them, is the velocity, is the density, is the stress tensor, is the strain tensor, is the elastic tensor (or stiffness tensor), is the body force, and the symbol “ ” represents the tensor product.

[0052] This equation is applicable for both isotropic and anisotropic materials, and dissipation can be added to the model in the form of Rayleigh damping, i.e., adding a mass and stiffness damping term to the right side of

[0053] In addition, it is also necessary to define the elastic parameters of the solid material, such as the elastic modulus, Poisson's ratio, etc. These parameters are used to calculate the wavelength of the sound wave in the material and the propagation speed, thereby affecting the propagation behavior of the sound wave; and define the boundary conditions to simulate the situation in the actual problem, which will affect the reflection, refraction and scattering of the sound wave in the solid; by adding a sound source to simulate the excitation of the sound wave. The solid mechanics module uses the finite element method to discretize and solve the elastic wave equation, in the process of discretization, the solid is divided into many small finite elements, and then the elastic wave equation is solved on these elements to obtain the propagation of the sound wave in the entire solid. In the post-processing process, the results can be post-processed, such as sound field distribution, wave speed map, frequency spectrum analysis, etc., which can be used to analyze the propagation characteristics of the sound wave in the solid.

[0054] Figure 7 is a schematic diagram of the simulation idea of a wave field in the present application, which mainly includes creating a geometric model, setting physical fields, parameters and boundary conditions, meshing, and post-processing.

[0055] Modeling and simulation of wave field in biological tissue: 1. Geometric model establishment In order to explore the influence of the layered characteristics of biological tissue on the propagation of elastic waves, single-layer, double-layer and three-layer tissue models are established.

[0056] Single-layer model: 15mm thickness, the elastic modulus is set to the value of the dermis layer.

[0057] Double-layer model: from top to bottom, the dermis layer with a thickness of 1.15mm and the subcutaneous tissue with a thickness of 13.85mm.

[0058] Three-layer model: from top to bottom, the epidermis layer with a thickness of 0.15mm, the dermis layer with a thickness of 1.15mm and the subcutaneous tissue with a thickness of 13.7mm.

[0059] For biological tissue, it is considered as an almost incompressible material. At this time ​, 0.499 in the simulation process. In order to analyze the elastic wave at different points over time, the generalized stretch coupling operator is defined at each layer, at the layered and in the vertical direction, which does not affect the subsequent mesh partitioning, and the required data can be exported for post-processing analysis after the simulation is completed. For the three-layer model defined by the generalized stretch coupling operator, the generalized stretch coupling operator is defined at equal to 15 mm, 14.925 mm, 14.85 mm, 14.275 mm, 13.7 mm, 13.2 mm between 6 mm and 9 mm, and in the vertical direction equal to 9 mm between 12 mm and 15 mm, which is used to analyze the time-domain characteristics of elastic waves at different positions. Figure 8 A three-layer model in the present application is shown in the figure.

[0060] 2. Determination of simulation parameters 1) Add physical field.

[0061] After the geometric model is established, the physical field of the simulation model is set to solid mechanics, and transient research is selected, so that the wave propagation can be analyzed in the time domain. In order to reduce the influence of boundary reflection, the left and right sides of the model are set to low reflection boundary conditions, so as to increase the absorption of the wave and reduce the influence on the collected elastic wave signal, and ensure the accuracy in the simulation process. According to the characteristics of biological tissues, the bottom side is set to a fixed constraint boundary, and the surface is set to a free boundary condition, and the interlayer is set to a solid-solid interface, so as to ensure the propagation of elastic waves.

[0062] 2) Apply excitation signal.

[0063] The sound field is calculated by Field II, and the sound field data is imported into COMSOL as an excitation source. The focused area is selected as the signal input source and point cloud data is used for input, which does not affect the construction of the grid. The sound pressure data is converted to a force field, and the most focused elliptical region of the force field is selected. The red elliptical region is the focus of attention, and the information in the elliptical region is extracted. First, find the maximum force and determine its coordinate value, find the most suitable focusing area, where the half major axis in the x-axis direction is 0.2 mm and the half major axis in the y-axis direction is 0.5 mm, and output the physical coordinates and force values of the points in the ellipse, which are loaded into COMSOL as input. In order to enable the acoustic radiation force to be loaded onto the surface and internal area of the tissue, the center horizontal coordinate of the extracted force field node information corresponds to the center position of the model, i.e. x = 10 mm; the maximum vertical coordinate of the force field corresponds to the upper surface of the model, i.e. the force field is tangent to the upper surface of the model. Figure 9The force loading model diagram of the acoustic radiation force field and the single-layer acoustic field of a modulated signal in the application.

[0064] 3) Dividing the grid size.

[0065] Compared with the quadrilateral grid, the triangular grid has stronger adaptability and can reduce the grid nodes, which can not only reduce the calculation cost and memory requirement, but also obtain better solution quality, so in the embodiment of the application, the triangular grid is adopted for division.

[0066] The relationship between the grid size and the grid density . .

[0067] The grid size of different biological tissue layers at different frequencies is different, and the relationship between the wavelength and the wave speed . .

[0068] For transient elastic waves in solids, at least 5 to 6 unit grids are required for each wavelength to analyze the wave propagation. Figure 10 The schematic diagram of analyzing wave propagation of transient elastic waves in solids in the application.

[0069] When the frequency is 1000 Hz, the Rayleigh wave speed is used for calculation when calculating the grid because the Rayleigh wave speed is smaller than the wave speed of other waves. The grid size required by different layers is different and is refined. Figure 11 The grid division schematic diagram corresponding to different models in the application, wherein (a) corresponds to a single-layer model; (b) corresponds to a double-layer model; and (c) corresponds to a three-layer model.

[0070] 4) Configuring the solver.

[0071] When the center frequency of the acoustic wave is set to , the period is , because the propagation of elastic waves needs time, in order to be able to capture the propagation trajectory of the wave, the total time T of the simulation process is set to 10 acoustic wave periods, .

[0072] When the propagation of elastic waves is described by using the displacement field or other physical fields, the smaller the time step is, the more time nodes are, and the more accurate the details of the elastic wave obtained are. In order to be able to observe the propagation of the elastic wave in detail, the time step is set to 1% of the acoustic wave period, .

[0073] Elastic wave simulation under different excitation models: 1. High-frequency excitation elastic wave simulation.

[0074] On the single-layer model, a 1MHz excitation signal, a Gaussian pulse signal, is used to load and get the wave propagation snapshots at different times.

[0075] Result 1) The energy of the elastic wave is mainly distributed on the surface, which is the Rayleigh wave, followed by the transverse wave, which has a slightly larger wave speed than the Rayleigh wave, and the longitudinal wave has the fastest speed. Figure 12 It is a cloud map of an elastic wave in the application, according to Figure 12 It can be seen that when , and , the longitudinal wave propagation, arrival at the bottom of the model and reflection are observed respectively.

[0076] Result 2) The Rayleigh wave mainly propagates on the surface, centered on the vibration source, while the transverse wave and the longitudinal wave propagate outward from the vibration source in a spherical shape. From the color change of the cloud map, it can be seen that the Rayleigh wave has the highest energy, followed by the transverse wave, and the longitudinal wave has the lowest energy.

[0077] In order to analyze the time-domain propagation characteristics of the wave, three particles on the model surface with an interval of 2mm are extracted, and the vibration characteristics of the particles are studied in detail. Figure 13 It is a schematic diagram of the vibration characteristics of a particle on the surface of a model in the application, and the results show that: as time goes on, the longitudinal wave separates from the mixed wave first due to its fastest speed, as shown in the red box; the transverse wave is second, shown in the yellow box; the Rayleigh wave has the largest amplitude, shown in the green box, and the longitudinal wave has the smallest amplitude, which confirms the significant difference in propagation characteristics of different types of elastic waves.

[0078] 2, acoustic radiation force field model simulation.

[0079] In order to show the characteristics of wave propagation at different frequencies and models, displacement cloud maps at different frequencies are selected for single-layer, double-layer and three-layer models. Specifically, when the frequency is 100Hz, the time point is 5.4ms; when the frequency is 1000Hz, the time point is 1.35ms; when the frequency is 5000Hz, the time point is 0.58ms. In order to more intuitively compare the propagation differences of the wave in different models, the displacement cloud map is normalized. Figure 14 It is a normalized displacement cloud map of different models in the application.

[0080] Analysis of the characteristics of elastic waves in different layered structures: 1, time domain analysis.

[0081] To further investigate the influence of interlayer interfaces on elastic wave propagation, this embodiment of the invention selected several particles along the depth direction in different models and conducted a detailed analysis of the time-domain signals of these particles. For models with different layers, this embodiment of the invention selected particles at different depths, 1 mm away from the center of the excitation point. In the single-layer model, three particles were selected at 1 mm intervals along the depth direction; in the layered model, corresponding particles were selected at the interlayer interfaces and the middle position of each layer.

[0082] When considering layered structures, at low frequencies, due to the longer wavelength, the displacement values ​​of two-layer and three-layer structures are generally greater than those of single-layer structures. At the specified frequencies, the displacements of the double-layer and triple-layer structures were approximately 63.86% and 64.33% greater than those of the single-layer structure, respectively. At high frequencies, the displacement values ​​of all three structures were generally smaller, with a maximum difference of 96.99% compared to the low-frequency case. Furthermore, in the presence of viscosity, the maximum displacement was 6.56% smaller than in the absence of viscosity. In the vertical direction, due to the distribution of the applied volume load along a certain length, the maximum displacement typically did not appear at the surface but rather in the interlayer locations. This phenomenon indicates that the displacement is significantly reduced at high frequencies and when the tissue has a certain degree of viscosity. Therefore, measurements should be performed at mid-to-low frequencies (0.1-1 kHz) for more effective data detection and analysis, which helps improve the accuracy of the experiment and the reliability of the data.

[0083] 2. Frequency domain analysis.

[0084] To analyze the frequency domain characteristics of elastic waves, Fourier transforms were performed on the time domain signals received at the same location points under different conditions, converting the signals from the time domain to the frequency domain (0-9000Hz). Figure 15 This is a schematic diagram of an elastic wave Fourier transform according to the present invention.

[0085] 3. Calculation of group velocity dispersion curve.

[0086] Statistical analysis was performed on the displacement data measured at various locations on the surface for different models, and the displacement-time (DST) values ​​for each model at different frequencies were obtained. Gradient plot, such as Figure 16 As shown, it provides quantitative information on the propagation speed of elastic waves, where the propagation speed can be determined by the slope of the wave's propagation trajectory. At this point, according to... The wave velocity determined by the slope of the displacement graph in the figure is the group velocity, which is equivalent to the wave velocity obtained by the time-of-flight method.

[0087] Based on the frequency sampling method in Table 1 below, dispersion curves were plotted on different model surfaces. Figure 17 This is a schematic diagram of dispersion curves on the surface of a different model in this invention.

[0088] Table 1 Frequency sampling method The analysis results show that the single-layer linear elastic model does not exhibit frequency dispersion within the test frequency range, and the average wave speed of all known points is calculated to be 2.8763 m / s, with a standard deviation of 0.0057 m / s. In contrast, the frequency dispersion curves of the double-layer and triple-layer models can be roughly divided into three parts: the data in region A exhibits distortion and lacks reliability; region B is a transition zone; and region C is a stable zone. In the double-layer model, the elastic moduli of the upper and lower layers are 27 kPa and 9 kPa, respectively, and the corresponding theoretical Rayleigh wave speeds are 2.8606 m / s and 1.6516 m / s, respectively. In the stable state C region (green box region), the wave speed is calculated to be 2.8849 ± 0.0029 m / s, with an error of 0.85%, which is suitable for characterizing the elastic properties of the tissue surface. In the transition phase of region B, the wavelength is long, and when the frequency is 1000 Hz, the wavelength is about 2.4 mm, which is about twice the thickness of the upper layer, and the calculation error reaches 16.7% at this time. When the frequency is as high as 3000 Hz, the wave speed begins to stabilize and approaches the theoretical Rayleigh surface wave speed. For the triple-layer model, the wave speed continues to increase with the increase of frequency, but it does not reach the thin layer at the top. It is expected that when the frequency is further increased, the wave speed will reach a stable state, and at this time the wave speed can be used to characterize the elastic properties of the model surface. Due to the influence of viscosity, the wave speed gradually changes with the frequency, although there is a stable trend, but using this wave speed to directly calculate the elastic modulus will introduce a large error, so it is necessary to accurately estimate the phase velocity of the viscoelastic material.

[0089] Simulation model viscoelasticity calculation: For elastic tissue without viscosity, take 5000 Hz as an example to analyze the wave speed, and the phase velocity dispersion curve obtained by the dispersion curve calculation method is as follows: Figure 18 FIG. 1 is a phase velocity dispersion curve of an elastic tissue without viscosity in the present application.

[0090] For tissue with viscosity, the previous analysis shows that the attenuation is large at a frequency of 5000 Hz, and the data has significant deviation. Therefore, data calculation is performed at a frequency of 1000 Hz to obtain more accurate results, Figure 19 FIG. 2 is a phase velocity dispersion curve of a tissue with viscosity in the present application.

[0091] Conclusion: 1) The change trend of phase velocity is consistent with the group velocity; 2) For the single-layer model, after discarding the data corresponding to the low-intensity frequency, the results show that the change of phase velocity is almost independent of frequency, and no frequency dispersion is exhibited; 3) For the double-layer model, the phase velocity gradually tends to be stable with the increase of frequency, and after reaching stability, the curves of the single-layer and double-layer models almost coincide, and the wave velocity is stable at about 2.8 m / s; 4) For the three-layer model, the frequency range of 300-600 Hz with higher curve fitting degree should be selected for analysis, which helps to improve the fitting accuracy of the model and ensure the reliability of the analysis results.

[0092] Taking the calculation of the single-layer elastic model as an example, it is found that when the frequency is in the range of 2000-2300 Hz, the fitting coefficient obtained by fitting the curve in this range is higher, the shear modulus is 8906.54±1.23 Pa, the shear viscosity is 0.06±0.01 Pa·s, and it is close to 0. The elastic modulus in the model is set to 27 kPa, and according to the formula, the shear modulus is about 9000 Pa. The relative error between the fitting result and the set parameters of the dermis layer is about 1.04%, and R2 is 0.997, indicating that the fitting quality is very high. It confirms the effectiveness of the measurement and calculation method, and also shows that the model and settings used can accurately reflect the viscoelastic properties of the single-layer biological tissue model.

[0093] The following Table 2 summarizes the shear modulus and shear viscosity estimated by the Rayleigh wave dispersion model for different models, as well as the shear modulus calculated by the wave velocity inversion after the dispersion curve is stable. As can be seen from the table, the shear modulus obtained by fitting and the shear modulus calculated by wave velocity without considering viscosity have good consistency, and the error compared with the actual set parameters is small and within an acceptable range. In contrast, the relative error of the inversion result without considering viscosity is smaller, about 0.71%, and compared with the double-layer model, the calculation result of the single-layer model is more accurate, but this method cannot obtain the viscosity value. In addition, it is found that the relative error of the calculation result of the three-layer model does not increase, which is considered to be due to the lack of stable frequency range at high frequency, and the result of selecting stable low-frequency data in the calculation. In the case where the wave velocity shows a stable trend, the fitting result in the stable stage is more accurate, and the error is relatively small. In the fitting process, the frequency range with smaller error should be selected for fitting according to the characteristics of the data. This can effectively improve the accuracy of the model estimation and ensure the reliability of the analysis results. By selecting the frequency range for fitting, the deviation can be minimized to obtain more accurate parameter estimation of the actual mechanical properties. In addition, the method of combining wave velocity inversion elastic modulus and fitting viscosity can be used for comprehensive analysis to obtain more accurate viscoelastic data.

[0094] Table 2 Inversion results of different models In summary, based on the proposed viscoelastic parameter inversion method, the wave field data in different models are calculated, and it is found that without considering the stratification and viscosity, the relative error of the phase velocity dispersion curve fitting is 1.04%, and the maximum relative error obtained by fitting in the calculation of the four models is not more than 9.74%, which proves the feasibility of the method for measuring the viscoelastic mechanical parameters of the dermis layer of biological tissues.

[0095] The ARF-OCE experimental technology for measuring the viscoelasticity of biological tissues is composed of an OCT acquisition module and an acoustic radiation force excitation module, and the experimental instrument is optimized according to the simulation results of the ultrasonic transducer, and the specific parameters are as shown in the table. Figure 20 The experimental instrument is optimized according to the simulation results of the ultrasonic transducer, as shown in the table.

[0096] Figure 21 The schematic diagram of one kind of OCT acquisition mode in the application is shown in the table. A-line is repeatedly acquired along time, and the acquisition points are composed of 1024 pixels along the longitudinal direction. The related information of the one-dimensional depth of the tissue is provided, and the effect of tracking the wave along the axial direction is provided. The time taken by the system to acquire two consecutive A-lines is called time resolution (Ts), and the A-line rate is 1 / Ts, which is the acquisition rate of the system. When the acquisition rate is 100 kHz, the time required for 1000 A-lines is 10 ms.

[0097] In order to analyze the process of the wave propagating along the axial direction, the M-Bmode acquisition mode can be used. The M-Bmode includes synchronous motion excitation and M-mode acquisition along each transverse position of the axial direction, so as to capture the wave propagation in two-dimensional space, as shown in part (b) of the table. Figure 21 In the system, the specific A-scan and B-scan number and voltage are set in Labview, and the motor is controlled to rotate by applying different voltages, so that different scanning lengths can be obtained. For a specific M-Bmode scan, the time taken by 1000 A-scans and n B-scans is 10n (ms).

[0098] In viscoelastic, isotropic and homogeneous materials (Kelvin-Voigt medium), the dispersion equation of surface acoustic wave (Rayleigh wave) can be expressed based on wave number as follows: ; , . ​​​​​

[0099] where, is the surface acoustic wave number, is the shear wave number, is the angular frequency, is the phase velocity of the surface acoustic wave, is the sample density, is the air density, and are the shear modulus and shear viscosity, respectively.

[0100] To estimate the viscoelastic parameters of the material, the obtained dispersion curve is iterated into the wave dispersion equation. The specific process is shown in Figure 22 .

[0101] The data acquisition is carried out by OCT. The phase stability of the experimental system during the detection process is the premise of determining whether the data is reliable or not. In order to reduce the interference of the phase information generated by the system itself, a correction arm is added as a reference. Before the experiment, the phase stability of the entire system still needs to be measured. The phase stability can be quantified in many ways. A common method is to calculate the phase noise or phase jitter, which can be represented by the standard deviation. The calculation formula of the standard deviation of the phase difference is as follows:

[0102] ; where, is the phase difference value of the th measurement, is the average value of all measured phase differences, is the number of measurements, which is equal to the number of A-Scans here, is the standard deviation of the phase difference, which represents the phase difference stability.

[0103] On the one hand, the experimental system was tested statically without applying excitation. The prepared biological phantom was placed without applying any excitation. The position of the objective lens was adjusted to focus the lens on the surface of the phantom. The length of the reference arm and the light intensity were adjusted to adjust the image position and improve the image quality. This experiment was carried out by single-point collection, i.e. M mode. The x and y voltages of the galvanometer were set to 0, and the A-scan was set to 2000. The results are shown in Figure 23 . Figure 23 (a) is the one-dimensional normalized amplitude along the depth direction; (b) is the OCT two-dimensional graph of M mode; (c) is the OCT phase difference graph; (d) is the phase difference curve of the 139th row in figure (c).

[0104] The ARF-OCE system is integrated to realize the detection and excitation of Rayleigh wave. Through the simulation results, the ultrasonic transducer with a focal length of 25 mm and a diameter of 30 mm is selected, and a low-frequency sinusoidal excitation signal of 600 Hz is designed, which contains the relevant parameters of other instruments. The system stability and synchronization of the OCT system are tested, and the phase stability of the experimental system is 0.2735 rad.

[0105] On the other hand, in order to test the synchronization of the imaging system and the excitation system under the loading state, dynamic testing is carried out on the experimental device. On the basis of static testing, the ultrasonic transducer is placed under the prepared biological phantom, and coupling agent is used as the coupling layer between the phantom and the transducer. The total thickness of the sample and the coupling layer is equal to the focal length of the transducer. The driving signal of the function generator is a continuous low-frequency sinusoidal signal modulated by a 1 MHz carrier, and the ultrasonic transducer is used as the excitation source to focus on the free surface-air interface of the phantom to generate elastic waves. The OCT system is used to track the elastic waves, and the results are shown in Figure 24 . Figure 24 Figures (a) and (c) are the phase difference images obtained after processing the continuous 600 Hz and 300 Hz sinusoidal signals, respectively; (b) and (d) are the original signals and filtered signals calculated from figures (a) and (c), respectively.

[0106] From Figure 24 (a) and (c), it can be reflected that different frequencies correspond to different wavelengths. Extract the data and perform filtering processing, and the results are shown in Figure 24 (b) and (d), the calculated filtered signals have frequencies of 600 Hz and 300 Hz, respectively, which are consistent with the input signals, so the system has good ability to track and identify the elastic waves generated in the phantom.

[0107] Phantom experiment verification: The results of the single-layer phantom experiment are shown in Figure 25 . For the single-layer model, the estimated elastic modulus and shear viscosity are 694.71 kPa and 0.78 Pa·s, respectively, by fitting the Rayleigh wave dispersion model.

[0108] The results of the double-layer phantom experiment are shown in Figure 26 . Fitting can obtain the elastic modulus and shear viscosity of 624.32 kPa and 0.41 Pa·s, respectively.

[0109] The results of the three-layer phantom experiment are shown in Figure 27 . Fitting can obtain the elastic modulus and shear viscosity of 746.82 kPa and 0.68 Pa·s, respectively.

[0110] Uniaxial tensile test verification: In order to determine the accuracy of the elastic modulus of the phantom obtained by the ARF-OCE experiment analysis, a micro tensile testing machine is used to carry out uniaxial tensile experiment on the prepared tensile phantom, the prepared tensile phantom is fixed at both ends of the tensile machine, the tensile rate is set to 0.01 m / s, first, pre-tension is carried out to make the test piece bear stress, the tensile length is set to 2 mm, the smaller tensile distance can assume that the behavior of the material at this time conforms to the linear elastic theory, the data after tension is calculated, and the result is shown in Figure 28 , Figure 28 It is a schematic diagram of fitting a stress-strain curve in the application. Table 3 is the tensile result of different degrees of silica gel in the application.

[0111] Table 3 Tensile result of silica gel with different degrees According to the experimental results of ARF-OCE, the relative error of the measurement of the elastic modulus can be obtained, the errors of the single-layer, double-layer and three-layer models based on the tensile result of 25-degree silica gel are 1.1%, 5.06% and 8.68% respectively, which is consistent with the conclusion of simulation, that is, the viscoelastic parameters of the dermis layer can be inverted, indicating the accuracy of the method.

[0112] Ex vivo pig biological tissue experiment: Figure 29 It is an ex vivo pig biological tissue experiment schematic diagram in the application, from the x-t diagram, the regular change can be accurately observed, which is consistent with the phenomenon of continuous excitation, the difference is that the phase velocity dispersion curve continuously increases in a certain frequency range, the elastic modulus and the shear viscosity are obtained by fitting, which are 2.09 MPa and 0.26 Pa·s respectively, and the tensile result of pigskin is 2.38 MPa, compared with the uniaxial tensile result, it is further proved that the ARF-OCE system and the viscoelastic inversion method designed in the research can be applied to the mechanical property measurement of biological tissue.

[0113] Based on the biological tissue viscoelastic optical coherence elastography method shown in Figure 1 , the following innovations are provided: Ultrasonic transducer: used for generating acoustic radiation force, the focal length diameter ratio is 0.8, the center frequency is 1 MHz, and a sinusoidal modulation signal of 0.1-1 kHz can be emitted.

[0114] OCT imaging module: including a swept-frequency laser source, an optical interference system, a signal acquisition and processing unit, used for acquiring phase information of biological tissue, and realizing imaging of layered microstructure of biological tissue.

[0115] Excitation control unit: used for controlling the excitation signal of the ultrasonic transducer, including a function generator and a signal amplifier.

[0116] Data processing unit: Used to process data acquired by OCT, including two-dimensional Fourier transform, threshold filtering, dispersion curve calculation and viscoelastic parameter inversion.

[0117] This invention conducts simulation studies on ARF (Advanced Radiofrequency) waves, using Field II to model transducers of different geometric dimensions and simulating the ultrasonic sound field under different excitation signals. The results show that when the focal length-to-diameter ratio is approximately 0.8, the results are ideal. Modulated low-frequency sinusoidal signals were used in subsequent simulations and experiments. Four biological tissue models were constructed using COMSOL, and time-domain and frequency-domain analyses of elastic waves were performed. The Rayleigh wave was identified as having the strongest energy on the model surface, and its wave velocity could be used to invert the elastic modulus of the tissue with an error of approximately 0.85%. The ARF-OCE system was optimized to achieve the excitation and detection of elastic waves, and the viscoelastic parameters of different phantoms were inverted. The calculated maximum error did not exceed 8.68%.

[0118] When applying the viscoelastic optical coherence elastography method for biological tissues provided by this invention, it is not necessary to consider... Figure 1 The steps shown are executed in sequence. The specific execution order of each step can be determined as needed, and this invention does not impose any restrictions on it.

[0119] The above describes a method for viscoelastic optical coherence elastography of biological tissues provided by one or more embodiments of the present invention. Based on the same idea, the present invention also provides a corresponding system for viscoelastic optical coherence elastography of biological tissues, such as... Figure 30 As shown.

[0120] Figure 30 A schematic diagram of a viscoelastic optical coherence elastography system for biological tissue provided by the present invention includes: The acoustic radiation force excitation module 201 includes: an ultrasonic transducer for generating acoustic radiation force under the excitation of an excitation signal, causing the target biological tissue to generate elastic waves; the ultrasonic transducer has a focal length-to-diameter ratio of 0.8 and uses a sine signal of 0.1-1kHz as the excitation signal. OCT acquisition module 202 is used to acquire displacement data of the target biological tissue during the elastic wave propagation process under the acoustic radiation force. The dispersion curve determination module 203 is used to perform two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue in order to determine the phase velocity dispersion curve of the target biological tissue. The inversion module 204 is used to fit the phase velocity dispersion curve of the target biological tissue based on the surface acoustic wave dispersion model and through a fitting algorithm to obtain the elastic modulus and shear viscosity of the target biological tissue.

[0121] The specific limitations of the bio-tissue viscoelastic optical coherence elastography system can refer to the limitations of the bio-tissue viscoelastic optical coherence elastography method described above, which will not be repeated here. Each module in the bio-tissue viscoelastic optical coherence elastography system described above can be realized by software, hardware, and a combination thereof, in whole or in part. The above-mentioned modules can be embedded in or independent of the processor in the computer device in hardware form, or can be stored in the memory in the computer device in software form, so as to be called and executed by the processor to perform the operations corresponding to each of the above modules.

[0122] The application also provides a computer readable storage medium storing a computer program, which can be used to execute the bio-tissue viscoelastic optical coherence elastography method described above. Figure 1 The bio-tissue viscoelastic optical coherence elastography method is provided.

[0123] The application also provides a computer device for implementing the bio-tissue viscoelastic optical coherence elastography method described above. Figure 31 The structure diagram of the computer device is shown in FIG. 1, which includes a processor, an internal bus, a network interface, a memory, and a non-volatile memory. Figure 31 At the hardware level, the computer device includes a processor, an internal bus, a network interface, a memory, and a non-volatile memory, and of course can also include other hardware required by the business. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs to implement the bio-tissue viscoelastic optical coherence elastography method described above. Figure 1 The bio-tissue viscoelastic optical coherence elastography method is provided.

[0124] A person of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer readable storage medium and can include the processes of the above-mentioned embodiments when executed. Any reference to a memory, storage, database, or other medium used in the embodiments provided by the application can include at least one of a non-volatile and volatile memory. The non-volatile memory can include a read-only memory (ROM), a tape, a floppy disk, a flash memory, or an optical memory. The volatile memory can include a random access memory (RAM) or an external cache memory. As an illustration but not limitation, the RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0125] The technical features of the above embodiments can be combined in any manner. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described, however, as long as the combinations of the technical features do not contradict each other, they should be considered to be within the scope of the present application.

Claims

1. A method of bio-tissue visco-elastic optical coherence elastography, characterized by, The method comprises the following steps: acquiring displacement data of an elastic wave propagation process of a target biological tissue under acoustic radiation force; the acoustic radiation force is generated by an ultrasonic transducer with a sine signal of 0.1-1 kHz as an excitation signal, and the focal length-diameter ratio of the ultrasonic transducer is 0.8; performing two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue to determine a phase velocity dispersion curve of the target biological tissue; based on a dispersion model of acoustic surface waves, fitting the phase velocity dispersion curve of the target biological tissue by a fitting algorithm to obtain the elastic modulus and shear viscosity of the target biological tissue.

2. The bio-tissue visco-elastic optical coherence elastography imaging method of claim 1, wherein, The displacement data of the target biological tissue comprises two-dimensional space-time slice data of a biological tissue surface, and the two-dimensional space-time slice data of the target biological tissue surface is used to represent displacement-time information of the target biological tissue surface; The two-dimensional discrete Fourier transform and threshold filtering are performed on the displacement data of the target biological tissue to determine the phase velocity dispersion curve of the target biological tissue, and specifically comprising: performing two-dimensional discrete Fourier transform on the displacement data of the target biological tissue by the following formula to convert it from the space-time domain to the frequency-wave number domain: ； determining the maximum intensity value of the displacement data of the target biological tissue in the frequency-wave number domain, and filtering out the signals lower than 10% of the maximum intensity value in the displacement data of the target biological tissue in the frequency-wave number domain by threshold filtering to remove background noise; according to the maximum acoustic surface wave wave number corresponding to each frequency in the filtered displacement data of the target biological tissue in the frequency-wave number domain, determining the phase velocity value at each frequency by the following formula to obtain the phase velocity dispersion curve of the target biological tissue: ； in, The energy distribution field of the target biological tissue in the frequency-wavenumber domain. For the time-space-displacement field of the target biological tissue surface, For the spatial coordinates of the elastic wave, For the time coordinate of the elastic wave, For wave number, For frequency, For spatial discrete sampling index, For time-discrete sampling index, The imaginary unit, For frequency The phase velocity value below, For frequency The corresponding maximum surface acoustic wave number.

3. The bio-tissue visco-elastic optical coherence elastography imaging method of claim 1, wherein, The fitting algorithm is used to fit the phase velocity dispersion curve of the target biological tissue based on the dispersion model of acoustic surface waves to obtain the elastic modulus and shear viscosity of the target biological tissue, and specifically comprising: determining a fitting frequency interval according to a preset minimum frequency and frequency interval length, fitting the data of the phase velocity dispersion curve of the target biological tissue in the fitting frequency interval to the dispersion model of acoustic surface waves under viscoelastic, isotropic and homogeneous materials, determining an initial range of shear wave wave number according to the acoustic surface wave wave number, and finding the optimal solution of the shear wave wave number satisfying the dispersion equation of the acoustic surface wave in the initial range of the shear wave wave number by iteration; according to the optimal solution of the shear wave wave number and the relationship between the shear wave wave number and the shear modulus and shear viscosity, performing nonlinear least squares curve fitting, making the fitting correlation coefficient of the nonlinear least squares curve fitting result and the optimal solution of the shear wave wave number better than a preset fitting correlation coefficient, and obtaining the elastic modulus and shear viscosity of the target biological tissue.

4. The bio-tissue visco-elastic optical coherence elastography imaging method of claim 3, wherein, The target biological tissue is the dermis of the skin, and the fitting frequency interval is 300-600 Hz.

5. A bio-tissue visco-elasticity optical coherence elastography system, characterized in that, The method comprises the following steps: an acoustic radiation force excitation module comprising: an ultrasonic transducer for generating acoustic radiation force under the excitation of an excitation signal to make the target biological tissue generate elastic waves; the focal length-diameter ratio of the ultrasonic transducer is 0.8, and a sine signal of 0.1-1 kHz is used as the excitation signal; an OCT acquisition module for acquiring displacement data of an elastic wave propagation process of a target biological tissue under acoustic radiation force; The frequency dispersion curve determination module is configured to perform two-dimensional discrete Fourier transform and threshold filtering on the displacement data of the target biological tissue to determine a phase velocity frequency dispersion curve of the target biological tissue. The inversion module is configured to perform fitting on the phase velocity frequency dispersion curve of the target biological tissue based on a surface acoustic wave frequency dispersion model by using a fitting algorithm to obtain the elastic modulus and shear viscosity of the target biological tissue.

6. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, and the computer program is executed by the processor to implement the method in any one of claims 1-4.

7. A computer device, comprising: The computer program is stored in the memory and executable on the processor, and the processor implements the method in any one of claims 1-4 when executing the program. The computer program is stored in the memory and executable on the processor, and the processor implements the method in any one of claims 1-4 when executing the program.

Images