A strain adaptive calculation method for optical coherence elastic imaging

Abstract

The present invention discloses a strain adaptive calculation method for optical coherence elastic imaging, which reconstructs the strain field distribution inside a composite material with high precision through an adaptive rotation model of a differential phase combined with a light vector method. The present invention adopts an algorithm combining an adaptive rotation model with a vector method, and does not have the problem of being unable to improve the strain calculation signal-to-noise ratio under multiplicative noise when using a least squares method. In addition, the present invention finds the optimal rotation angle by constructing an adaptive rotation model of a differential phase, and then combines the light vector method to reduce its own lateral inhomogeneity to a certain extent, thereby solving the problem that the light vector method cannot accurately construct the strain field distribution inside a composite material with a low signal-to-noise ratio under strong noise and complex deformation conditions, especially when the lateral changes are abnormally non-uniform.

Description

Technical Field

[0001] The invention relates to the field of optical coherence elastic imaging, and in particular to a strain adaptive calculation method for optical coherence elastic imaging. Background Art

[0002] There are two existing optical coherence elastography strain calculation schemes: the least squares method and the light vector method.

[0003] Among them, the least square method is essentially equivalent to smoothing the noisy signal, which can only be used to suppress additive noise, but cannot improve the strain calculation signal-to-noise ratio under multiplicative noise. This shortcoming makes the least square method very limited in the internal tomography measurement of materials. Compared with the least square method, the light vector method cleverly uses the characteristics of vector superposition and has the ability to suppress additive and multiplicative noise through weighting and phase smoothing; however, when the light vector method faces complex deformation of composite materials, the differential phase of optical coherence elastic imaging will show lateral inhomogeneity; at this time, the phase information will be lost in the horizontal direction of the light vector method, resulting in low resolution and poor accuracy of strain imaging in the horizontal direction.

[0004] In summary, in the field of strain calculation of optical coherence elasticity, there is still a lack of a strain calculation method that can perform well under strong noise and complex deformation conditions. Summary of the invention

[0005] The purpose of the present invention is to provide a strain adaptive calculation method for optical coherence elastic imaging, which is used to accurately calculate the strain field distribution under the conditions of low signal-to-noise ratio and non-uniform variation of the differential wrapped phase field.

[0006] In order to achieve the above tasks, the present invention adopts the following technical solutions:

[0007] A strain adaptive calculation method for optical coherence elastic imaging comprises the following steps:

[0008] Step 1, constructing tomographic interference light field signals before and after deformation for the composite material respectively; constructing a differential interference light field signal based on the tomographic interference light field signal, and determining a differential wrapped phase field corresponding to the deformation field;

[0009] Step 2, rotating the differential wrapped phase field counterclockwise around its center according to the set step size;

[0010] Step 3, constructing a sliding window for the differential wrapped phase field corresponding to each rotation angle, averaging the differential interference light field signal in the selected sliding window in the horizontal direction, and obtaining the differential interference light field signal after horizontal averaging;

[0011] Step 4: for the differential interference light field signals averaged in the horizontal direction, based on the multiplication rule of complex numbers, the differential wrapped phases of the preset pixel intervals are multiplied in pairs to construct the second-order differential interference light field signals in the depth direction;

[0012] Step 5, averaging and normalizing the second-order differential interference light field signal in the depth direction to obtain a normalized second-order differential interference light field signal after averaging in the depth direction;

[0013] Step 6, calculating the strain field distribution in the depth direction on the normalized second-order differential interference light field signal averaged in the depth direction;

[0014] Step 7, using the same method as steps 1 to 6, wherein the horizontal average in step 3 is replaced by the depth average, and the depth average in step 5 is replaced by the horizontal average, so as to construct the strain field distribution in the horizontal direction;

[0015] The strain field distributions in the horizontal direction and the depth direction are rotated clockwise again to obtain the strain field distributions in the horizontal direction and the depth direction that are the same as the differential wrapped phase field position before the rotation;

[0016] Step 8, integrating the strain field distribution in the horizontal direction and the depth direction after rotation in the horizontal direction and the depth direction respectively, and obtaining the corresponding unwrapped differential phase field in the horizontal direction and the depth direction;

[0017] Step 9, performing phase wrapping on the unwrapped differential phase field to obtain differential wrapped phase fields in the horizontal direction and the depth direction; and then obtaining a restored differential wrapped phase field based on the multiplication rule of complex numbers;

[0018] Step 10, solving the correlation coefficient between the differential wrapped phase field and the differential wrapped phase field;

[0019] Step 11, according to the above steps 2 to 9, the restored complete differential wrapped phase field in step 9 is obtained at different rotation angles, and then through step 10, a corresponding set of correlation coefficients can be obtained;

[0020] The least squares fitting method was used to obtain the polynomial function between the rotation angle and the correlation coefficient;

[0021] Step 12, solving the polynomial function to obtain an optimal solution, and then substituting the solution into the polynomial function to obtain an optimal rotation angle;

[0022] Step 13, according to the optimal rotation angle obtained in step 12, obtain the strain field distribution in the depth direction of the composite material through steps 3 to 6, and then rotate clockwise to restore the position of the strain field distribution in the depth direction of the material to before rotation. At this time, the strain field distribution under complex deformation of the material can be obtained.

[0023] Furthermore, the tomographic interference light field signals before and after the deformation are expressed as:

[0024] a1(m,j)=A1(m,j)×exp(iφ1(m,j))

[0025] a2(m,j)=A2(m,j)×exp(iφ2(m,j))

[0026] Wherein, a1(m, j) represents the tomographic interference light field signal before deformation, a2(m, j) represents the tomographic interference light field signal after deformation, (m, j) represents the pixel coordinates of the interference light field spatial image, i is an imaginary unit, exp represents exponential operation, A1(m, j) represents the interference spectrum amplitude before deformation, φ1(m, j) represents the phase field of the composite material before deformation, A2(m, j) represents the interference spectrum amplitude after deformation, and φ2(m, j) represents the phase field of the composite material after deformation;

[0027] The complex multiplication criterion is used to construct the differential interference light field signal:

[0028] a2(m,j)a1(m,j) ★ ≡b(m,j)=B(m,j)exp[i·Φ(m,j)]

[0029] In the above formula, ★ represents the conjugate operation, b(m, j) represents the differential interference light field signal, B(m, j) = A2(m, j)A1(m, j) represents the amplitude information, and Φ(m, j) represents the phase field of the differential interference light field signal;

[0030] The differential wrapped phase field corresponding to the deformation field is: Δφ(m, j)=Arg{exp[i·Φ(m, j)]}, where Arg represents the phase of a complex number.

[0031] Furthermore, the rotation angle is θ, θ=θ0+Δθ, and the rotation stops when the rotation angle reaches 360°; wherein the range of θ is controlled between 0° and 360°, θ0 is the starting angle, and Δθ is the interval between continuous rotation angles θ.

[0032] Furthermore, the horizontally averaged differential interference light field signal is expressed as:

[0033]

[0034] Where N X 、N z represents the width of the sliding window in the horizontal and depth directions, B(j) represents the amplitude of the differential interference light field signal after horizontal averaging, represents the amplitude and phase of the differential interference light field signal after being averaged in the horizontal direction, and m and j represent the horizontal and vertical coordinates of the pixels in the spatial image of the differential interference light field.

[0035] Furthermore, the second-order differential interference light field signal in the depth direction is expressed as:

[0036]

[0037] in, B(j+1), Graph represents the horizontally averaged differential interference light field signal at j+1, the horizontally averaged differential interference light field signal amplitude, and the horizontally averaged differential interference light field signal amplitude phase.

[0038] Furthermore, the strain field distribution in the depth direction is expressed as:

[0039]

[0040] In the formula, Represents the normalized second-order differential interference light field signal after averaging in the depth direction The phase, d represents the distance between pixels, and the correlation coefficient of the ratio Where λ0 represents the central wavelength of the laser and n represents the refractive index of the composite material.

[0041] Furthermore, the recovered differential wrapped phase field The atmosphere is shown as:

[0042]

[0043]

[0044] In the formula, ★ represents the conjugate operation, and Arg{} represents the phase of the complex number.

[0045] Compared with the prior art, the present invention has the following technical features:

[0046] The solution of the present invention can reconstruct the strain field distribution inside the composite material with high precision by combining the adaptive rotation model of the differential phase with the light vector method. Since the algorithm used in the present invention that combines the adaptive rotation model with the vector method is developed based on the vector method, there is no problem that the signal-to-noise ratio of strain calculation under multiplicative noise cannot be improved when the least squares method is used. In addition, the adaptive rotation model of the differential phase is constructed, the optimal rotation angle is found, and then combined with the light vector method, the lateral inhomogeneity itself is reduced to a certain extent, and the problem that the light vector method cannot accurately construct the strain field distribution inside the low signal-to-noise ratio composite material under strong noise and complex deformation conditions, especially the abnormally uneven lateral changes. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 It is a schematic diagram of the process of the present invention;

[0048] Figure 2 A phase diagram showing an abnormally non-uniform lateral change in stiffness of a composite material in one embodiment of the present invention;

[0049] Figure 3 The strain field diagram of the composite material calculated by the light vector method;

[0050] Figure 4 The strain field component map of the composite material is reconstructed by the method of the present invention. DETAILED DESCRIPTION

[0051] The present invention aims to provide a strain adaptive calculation method for optical coherence elastic imaging. By constructing an adaptive rotation model of differential phase combined with a light vector method, the strain field distribution inside a low signal-to-noise ratio composite material can be reconstructed with high precision under complex material deformation conditions, especially when the lateral changes are abnormally uneven.

[0052] The present invention provides a strain adaptive calculation method for optical coherence elastic imaging, such as Figure 1 As shown, the following steps are included:

[0053] Step 1: construct tomographic interference light field signals before and after deformation for a composite material with complex deformation; construct differential interference light field signals based on the tomographic interference light field signals, and determine the differential wrapping phase field corresponding to the deformation field. The deformation refers to the deformation of the composite material caused by some factors (such as pressure, temperature, etc.).

[0054] a1(m,j)=A1(m,j)×exp(iφ1(m,j))

[0055] a2(m,j)=A2(m,j)×exp(iφ2(m,j))

[0056] Wherein, a1(m, j) represents the tomographic interference light field signal before deformation, a2(m, j) represents the tomographic interference light field signal after deformation, (m, j) represents the pixel coordinates of the interference light field spatial image, i is an imaginary unit, exp represents exponential operation, A1(m, j) represents the interference spectrum amplitude before deformation, φ1(m, j) represents the phase field of the composite material before deformation, A2(m, j) represents the interference spectrum amplitude after deformation, and φ2(m, j) represents the phase field of the composite material after deformation;

[0057] Use the complex multiplication criterion to construct the differential interference light field signal:

[0058] a2(m,j)a1(m,j) ★≡b(m,j)=B(m,j)exp[i·Φ(m,j)] (1)

[0059] In the above formula, ★ represents the conjugate operation, b(m, j) represents the differential interference light field signal, B(m, j) = A2(m, j)A1(m, j) represents the amplitude information, and Φ(m, j) represents the phase field of the differential interference light field signal.

[0060] The differential wrapped phase field corresponding to the deformation field is: Δφ(m, j)=Arg{exp[i·Φ(m, j)]}, where Arg represents the phase of a complex number.

[0061] Step 2, rotate the differential wrapped phase field counterclockwise around its center, and the rotation angle is θ (θ = θ0 + Δθ, stop when the rotation angle reaches 360°; where the range of θ is controlled between 0° and 360°, θ0 is the starting angle (set to 0°), and Δθ is the interval between consecutive rotation angles θ (set to 1°)).

[0062] Step 3: To suppress the phase and amplitude effects of speckle noise, a sliding window of size N is constructed for the differential wrapped phase field with a rotation angle of θ. X ×N z , the differential interference light field signal b(m, j) is averaged horizontally in the selected sliding window to obtain the differential interference light field signal after horizontal averaging:

[0063]

[0064] Where N X 、N z represents the width of the sliding window in the horizontal and depth directions, B(j) represents the amplitude of the differential interference light field signal after horizontal averaging, represents the amplitude and phase of the differential interference light field signal after being averaged in the horizontal direction, and m and j represent the horizontal and vertical coordinates of the pixels in the spatial image of the differential interference light field.

[0065] Step 4: for the differential interference light field signal after horizontal averaging Based on the multiplication rule of complex numbers, the differential wrapped phase with a pixel interval of 1 is multiplied two by two to construct the second-order differential interference light field signal in the depth direction:

[0066]

[0067] in, B(j+1), It represents the horizontally averaged differential interference light field signal at j+1, the horizontally averaged differential interference light field signal amplitude, and the horizontally averaged differential interference light field signal amplitude phase.

[0068] Step 5: In view of the influence of the noise in the depth direction, the second-order differential interference light field signal c(j) in the depth direction is averaged and normalized along the depth direction to obtain a normalized second-order differential interference light field signal after averaging in the depth direction:

[0069]

[0070] In the above formula, N z represents the depth width of the selected sliding window, and |c(j)| represents the amplitude of the second-order differential interference light field signal in the depth direction.

[0071] Step 6: According to the gradient distribution of the strain along the depth direction, based on the linear mapping relationship between the differential wrapped phase field and the displacement, the strain field distribution in the depth direction is calculated on the normalized second-order differential interference light field signal averaged in the depth direction:

[0072]

[0073] In the formula, Represents the normalized second-order differential interference light field signal after averaging in the depth direction The phase, d represents the distance between pixels, and the correlation coefficient of the ratio Where λ0 represents the central wavelength of the laser and n represents the refractive index of the composite material.

[0074] Step 7, using the same method as steps 1 to 6, wherein the horizontal average in step 3 is replaced by the depth average, and the depth average in step 5 is replaced by the horizontal average, so as to construct the strain field distribution in the horizontal direction;

[0075] Let ε1 and ε2 be the obtained strain field distribution in the horizontal direction and depth direction respectively. Rotate ε1 and ε2 clockwise again with a rotation angle of θ to obtain the strain field distribution in the horizontal direction and depth direction with the same position of the differential wrapped phase field before rotation.

[0076] Step 8: Based on the nature of strain: the gradient distribution of displacement along the depth direction, Integrate in the horizontal and depth directions respectively to obtain the corresponding unwrapped differential phase fields in the horizontal and depth directions

[0077]

[0078] In the formula, dm and dj represent Integration operator in the horizontal and depth directions.

[0079] Step 9: Based on the relationship between the wrapped phase and the unwrapped phase, the unwrapped differential phase field in step 8 is Perform phase wrapping to obtain differential wrapped phase fields in the horizontal and depth directions Then, based on the multiplication rule of complex numbers, the recovered differential wrapped phase field is obtained

[0080]

[0081]

[0082] In the formula, ★ represents the conjugate operation, and Arg{} represents the phase of the complex number.

[0083] Step 10: Solve the differential wrapped phase field in step 9 by using the existing method based on the Pearson correlation coefficient principle. The correlation coefficient T between the differential wrapped phase field Δφ(m, j) in step 1 corr .

[0084] Step 11, according to the above steps 2 to 9, at different rotation angles θ (0° to 360°), obtain the restored complete differential wrapped phase field in step 9 The two correspond one to one, and then after step 10, a corresponding set of correlation coefficients T can be obtained. corrj (j = 0, 1, 2, ..., 360), the least squares fitting method is used to obtain the rotation angle θ and the correlation coefficient T corr The polynomial function f between θ (T corr ):

[0085] f θ (T corr )=k0T corr0 +k1T corr1 +k2T corr2 2 +…+k n T corrn n , (9)

[0086] Where k0, k1, ..., k n represents the function expansion coefficient obtained using the least squares fitting method, and n represents the order of the polynomial function.

[0087] Step 12: Use the quasi-Newton method to solve the polynomial function f θ (T corr ), and get the optimal solution T corr * , and then substitute it into the polynomial function f θ (T corr), and get the optimal rotation angle θ * :

[0088] T corr * =argmin θ f θ (T corr ) (10)

[0089] θ * =f θ (T corr * ) (11)

[0090] In the formula, argmin θ Represents solving the polynomial function f θ (T corr ) is the optimization operator for θ in .

[0091] Step 13: According to the optimal rotation angle θ obtained in step 12 * , the strain field distribution in the depth direction of the composite material is obtained through steps 3 to 6, and then the strain field distribution in the depth direction of the composite material is rotated clockwise by θ * , so that the strain field distribution position in the depth direction of the material is restored to the position before rotation, at this time the strain field distribution under complex deformation of the material can be obtained.

[0092] The technical solution of the present invention is described below in conjunction with a specific embodiment of the change in the stiffness of a composite material.

[0093] like Figure 2 As shown, it is the phase field distribution when the stiffness of the composite material changes abnormally and unevenly in the lateral direction. In this embodiment, the central wavelength of the laser is 8.4×10 -7 m, the light source bandwidth is 5×10 -8 m, the material refractive index is n = 1, and the selected sliding window size is 10 × 10 pixels.

[0094] Calculation using the light vector method Figure 2 , the reconstructed strain field distribution can be obtained, such as Figure 3 As shown in the figure, it can be seen that the strain field distribution obtained by the light vector method is very poor; this also confirms the problem of the light vector method: under the condition of strong noise and complex deformation, it is impossible to reconstruct the strain field distribution inside the composite material with low signal-to-noise ratio with high precision.

[0095] The strain field distribution of the composite material is calculated by the method proposed in the present invention. The results are as follows: Figure 4 As shown in the figure, it can be clearly seen that, compared with the light vector method, the method proposed in the present invention can reconstruct the strain field distribution of the composite material with high precision.

[0096] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. These modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application, and should all be included in the protection scope of the present application.

Claims

1. A strain adaptive calculation method for optical coherence elastic imaging, characterized in that: The following steps are involved: Step 1, constructing tomographic interference light field signals before and after deformation for the composite material respectively; constructing a differential interference light field signal based on the tomographic interference light field signal, and determining a differential wrapped phase field corresponding to the deformation field; Step 2, rotating the differential wrapped phase field counterclockwise around its center according to the set step size; Step 3, constructing a sliding window for the differential wrapped phase field corresponding to each rotation angle, averaging the differential interference light field signal in the selected sliding window in the horizontal direction, and obtaining the differential interference light field signal after horizontal averaging; Step 4: for the differential interference light field signals averaged in the horizontal direction, based on the multiplication rule of complex numbers, the differential wrapped phases of the preset pixel intervals are multiplied in pairs to construct the second-order differential interference light field signals in the depth direction; Step 5, averaging and normalizing the second-order differential interference light field signal in the depth direction to obtain a normalized second-order differential interference light field signal after averaging in the depth direction; Step 6, calculating the strain field distribution in the depth direction on the normalized second-order differential interference light field signal averaged in the depth direction; Step 7, using the same method as steps 1 to 6, wherein the horizontal average in step 3 is replaced by the depth average, and the depth average in step 5 is replaced by the horizontal average, so as to construct the strain field distribution in the horizontal direction; The strain field distributions in the horizontal direction and the depth direction are rotated clockwise again to obtain the strain field distributions in the horizontal direction and the depth direction that are the same as the differential wrapped phase field position before the rotation; Step 8, integrating the strain field distribution in the horizontal direction and the depth direction after rotation in the horizontal direction and the depth direction respectively, and obtaining the corresponding unwrapped differential phase field in the horizontal direction and the depth direction; Step 9, performing phase wrapping on the unwrapped differential phase field to obtain differential wrapped phase fields in the horizontal direction and the depth direction; and then obtaining a restored differential wrapped phase field based on the multiplication rule of complex numbers; Step 10, solving the correlation coefficient between the differential wrapped phase field and the differential wrapped phase field; Step 11, according to the above steps 2 to 9, the restored complete differential wrapped phase field in step 9 is obtained at different rotation angles, and then through step 10, a corresponding set of correlation coefficients can be obtained; The least squares fitting method was used to obtain the polynomial function between the rotation angle and the correlation coefficient; Step 12, solving the polynomial function to obtain an optimal solution, and then substituting the solution into the polynomial function to obtain an optimal rotation angle; Step 13, according to the optimal rotation angle obtained in step 12, obtain the strain field distribution in the depth direction of the composite material through steps 3 to 6, and then rotate clockwise to restore the position of the strain field distribution in the depth direction of the material to before rotation. At this time, the strain field distribution under complex deformation of the material can be obtained.

2. The strain adaptive calculation method for optical coherence elastic imaging according to claim 1, characterized in that: The tomographic interference light field signals before and after the deformation are expressed as: a1(m,j)=A1(m,j)×exp(iφ1(m,j)) a2(m,j)=A2(m,j)×exp(iφ2(m,j)) Among them, a1(m,j) represents the tomographic interference light field signal before deformation, a2(m,j) represents the tomographic interference light field signal after deformation, (m,j) represents the pixel coordinates of the interference light field space image, i is an imaginary unit, exp represents the exponential operation, A1(m,j) represents the interference spectrum amplitude before deformation, φ1(m,j) represents the phase field of the composite material before deformation, A2(m,j) represents the interference spectrum amplitude after deformation, and φ2(m,j) represents the phase field of the composite material after deformation; The complex multiplication criterion is used to construct the differential interference light field signal: a2(m,j)a1(m,j)*≡b(m,j)=B(m,j)exp[i·Φ(m,j)] In the above formula, * represents the conjugate operation, b(m,j) represents the differential interference light field signal, B(m,j)=A2(m,j)A1(m,j) represents the amplitude information, and Φ(m,j) represents the phase field of the differential interference light field signal; The differential wrapped phase field corresponding to the deformation field is: Δφ(m,j)=Arg{exp[i·Φ(m,j)]}, where Arg represents the phase of a complex number.

3. The strain adaptive calculation method for optical coherence elastic imaging according to claim 1, characterized in that: The rotation angle is θ, θ=θ0+Δθ, and the rotation stops when the rotation angle reaches 360°; wherein the range of θ is controlled between 0° and 360°, θ0 is the starting angle, and Δθ is the interval between continuous rotation angles θ.

4. The strain adaptive calculation method for optical coherence elastic imaging according to claim 2, characterized in that: The horizontally averaged differential interference light field signal is expressed as: Where N X 、N z represents the width of the sliding window in the horizontal and depth directions, B(j) represents the amplitude of the differential interference light field signal after horizontal averaging, represents the amplitude and phase of the differential interference light field signal after being averaged in the horizontal direction, and m, j represents the horizontal and vertical coordinates of the pixel in the spatial image of the differential interference light field.

5. The strain adaptive calculation method for optical coherence elastic imaging according to claim 4, characterized in that: The second-order differential interference light field signal in the depth direction is expressed as: in, B(j+1), It represents the horizontally averaged differential interference light field signal at j+1, the horizontally averaged differential interference light field signal amplitude, and the horizontally averaged differential interference light field signal amplitude phase.

6. The strain adaptive calculation method for optical coherence elastic imaging according to claim 1, characterized in that: The strain field distribution in the depth direction is expressed as: In the formula, Represents the normalized second-order differential interference light field signal after averaging in the depth direction The phase, d represents the distance between pixels, and the correlation coefficient of the ratio Where λ0 represents the central wavelength of the laser and n represents the refractive index of the composite material.

7. The strain adaptive calculation method for optical coherence elastic imaging according to claim 5, characterized in that: The recovered differential wrapped phase field It is expressed as: Wherein, * represents the conjugate operation, and Arg{} represents the phase of the complex number.

Images