Laminated plate impact event three-dimensional strain field inversion method based on binocular vision

Abstract

The invention discloses a laminated plate impact event three-dimensional strain field inversion method based on binocular vision, and the method comprises the steps: laying dynamic tracking mark points on the surface of an impacted carbon fiber laminated plate, synchronously collecting the image data of the whole impact process through a high-speed binocular camera, and generating a binocular vision image sequence with impact time sequence information; generating three-dimensional displacement field data of the impact area of the carbon fiber laminated plate based on the binocular vision image sequence; calling the anisotropic constitutive model of the carbon fiber laminated plate to generate a preliminary three-dimensional strain field inversion result; and performing coupling operation on the preliminary three-dimensional strain field inversion result, the impact load parameter and the laminated plate material characteristic parameter to generate a final carbon fiber laminated plate three-dimensional strain field inversion result consistent with the dynamic response of the impact event. According to the embodiment of the invention, accurate full-field representation of the dynamic strain of the impact event can be realized, and reliable data support is provided for impact resistance evaluation and structure optimization of the laminated plate.

Description

Technical Field

[0001] This invention belongs to the field of binocular vision technology, and in particular, it is a method for inverting the three-dimensional strain field of laminate impact events based on binocular vision. Background Technology

[0002] The dynamic strain response of carbon fiber laminates under impact loads is a key indicator for assessing their structural safety. However, traditional strain measurement methods have significant limitations: strain gauges can only acquire discrete point data and cannot reflect the overall strain distribution; monocular vision lacks spatial depth information, making it difficult to accurately invert three-dimensional strain; finite element simulations rely on simplified models, resulting in deviations from actual impact conditions. Especially under the influence of the anisotropic properties of laminates, the fiber layup direction leads to significant directional differences in strain distribution. Existing methods do not fully consider this characteristic, easily causing strain calculation errors and failing to accurately capture the local damage evolution during impact, thus hindering the impact-resistant design and optimization of laminate structures. Summary of the Invention

[0003] The purpose of this invention is to provide a three-dimensional strain field inversion method for laminate impact events based on binocular vision, so as to overcome the shortcomings of the prior art, realize accurate full-field characterization of dynamic strain of impact events, and provide reliable data support for the evaluation of impact resistance performance and structural optimization of laminates.

[0004] One embodiment of this application provides a method for inverting the three-dimensional strain field of a laminate impact event based on binocular vision, the method comprising: Dynamic tracking markers were placed on the surface of the impacted carbon fiber laminate. Image data of the entire impact process was acquired synchronously using a high-speed binocular camera. The magnitude of the impact load and the duration of the impact were recorded. The image coordinate system was unified through the binocular camera intrinsic parameter calibration algorithm to generate a binocular visual image sequence with impact time information. Based on the binocular vision image sequence, a hierarchical feature matching algorithm is used to extract the dynamic features of the marker points and the fiber texture of the laminate. The three-dimensional spatial coordinates of the feature points at each time sequence are calculated by stereo matching. The three-dimensional displacement vector of the feature points is solved by combining the difference between adjacent time sequence coordinates, and the three-dimensional displacement field data of the impact area of ​​the carbon fiber laminate is generated. The anisotropic constitutive model of carbon fiber laminate is invoked, the three-dimensional displacement field data is substituted into the strain tensor solution formula, and the spatial derivative is used to convert it into local strain components. Combined with the laminate layup direction to correct the strain calculation deviation, a preliminary three-dimensional strain field inversion result is generated. The preliminary three-dimensional strain field inversion results are coupled with the impact load parameters and the laminate material property parameters. Abnormal strain components are eliminated through strain-load correlation verification to generate the final three-dimensional strain field inversion results of carbon fiber laminate that are consistent with the dynamic response of the impact event.

[0005] Optionally, the process involves placing dynamic tracking markers on the surface of the impacted carbon fiber laminate, simultaneously acquiring image data of the entire impact process using a high-speed binocular camera, recording the magnitude and duration of the impact load, unifying the image coordinate system through a binocular camera intrinsic parameter calibration algorithm, and generating a binocular visual image sequence with impact timing information, including: The dynamic tracking markers are designed as concentric rings with alternating fluorescent and black colors. The rings are 2mm in diameter and 0.5mm apart. They are printed on a 50μm thick polyimide film using UV-curable adhesive to ensure that the markers have both high contrast and anti-motion blur characteristics under high-speed shooting. The output is the marker placement scheme. According to the marking point layout plan, the surface of the carbon fiber laminate is divided into an impact center area and an edge area. Marking points are pasted in the center area with a grid density of 5mm×5mm and in the edge area with a grid density of 10mm×10mm. The height difference of the marking points is controlled within 0.1mm using a laser thickness gauge. The laminate specimen with the marking points laid out is then output. The Zhang Zhengyou calibration method was used to perform intrinsic parameter calibration on the high-speed binocular camera. Twenty sets of images with different poses were acquired using a 12×9 checkerboard target. The focal length principal point coordinate distortion coefficients of the two cameras were calculated. The frame rate of the two cameras was locked at 1000fps and the trigger time difference was controlled within 1μs through the synchronous trigger module. The calibrated binocular camera system was then output. The laminate specimen with marked points is fixed on the impact test bench. The calibrated binocular camera system is started, and the drop hammer impact device is triggered synchronously. The magnitude and duration of the impact load are recorded by the force sensor. Binocular images of the entire impact process are acquired. The coordinate system of the two cameras is unified by the internal parameter data, and a binocular visual image sequence with impact time information is generated.

[0006] Optionally, based on the binocular vision image sequence, a hierarchical feature matching algorithm is used to extract the dynamic features of the marker points and the fiber texture of the laminate. The three-dimensional spatial coordinates of the feature points at each time step are calculated through stereo matching. The three-dimensional displacement vector of the feature points is then solved by combining the difference between adjacent time step coordinates, generating three-dimensional displacement field data of the carbon fiber laminate impact region, including: The binocular vision image sequence is preprocessed, adaptive bilateral filtering is used to remove high-speed imaging noise, HSV color space threshold segmentation is performed on the marked point region, Canny edge detection is performed on the fiber texture region, and the preprocessed image sequence is output. Based on the preprocessed image sequence, a hierarchical feature matching algorithm is adopted. The upper layer extracts ORB features from the marker points and performs brute-force matching, while the lower layer extracts histogram features of orientation gradients from the fiber texture and performs K-nearest neighbor matching to generate a sequence of matching pairs between marker points and texture features. The matching sequence is input into the stereo matching algorithm. The three-dimensional spatial coordinates of the feature points at each time step are calculated by triangulation in combination with the intrinsic parameters of the binocular camera. The coordinate accuracy is optimized by the bundle adjustment method, and the three-dimensional coordinate set of the feature points at each time step is output. For each time series feature point's three-dimensional coordinate set, the coordinate difference between adjacent time series is calculated to obtain the three-dimensional displacement vector. The discrete vector is extended into a continuous field using the inverse distance weighted interpolation method. Gaussian smoothing is used to eliminate interpolation errors, generating three-dimensional displacement field data of the carbon fiber laminate impact region.

[0007] Optionally, the step of calling the anisotropic constitutive model of the carbon fiber laminate, substituting the three-dimensional displacement field data into the strain tensor solution formula, converting it into local strain components through spatial derivative calculation, and correcting strain calculation errors by combining the laminate layup direction, generates preliminary three-dimensional strain field inversion results, including: Construct an anisotropic constitutive model of carbon fiber laminates, input ply parameters and material parameters, establish an elastic matrix based on the ply direction, and output the parameter set of the anisotropic constitutive model; The three-dimensional displacement field data was discretized into a 1mm×1mm×0.5mm grid. The spatial partial derivatives of the displacement components at each grid point were calculated using the second-order central difference method. The results were then substituted into the strain tensor solution formula to obtain 6 strain components, which were used as the initial strain component matrix. Based on the laminate ply direction, the local coordinate system rotation matrix corresponding to each grid point is calculated. The initial strain component matrix is ​​transformed into strain components in the local coordinate system through coordinate transformation. The calculation deviation caused by ply anisotropy is corrected, and the corrected strain component matrix is ​​output. The corrected strain component matrix is ​​mechanically consistent with the parameter set of the anisotropic constitutive model. Strain values ​​exceeding the elastic limit of the material are removed, and the removed region is completed by linear interpolation to generate preliminary three-dimensional strain field inversion results.

[0008] Optionally, the step of coupling the preliminary three-dimensional strain field inversion results with the impact load parameters and laminate material property parameters, and eliminating abnormal strain components through strain-load correlation verification to generate the final three-dimensional strain field inversion result of the carbon fiber laminate consistent with the dynamic response of the impact event includes: We organize the impact load parameters and the laminate material property parameters, construct a parameter association database, map parameters of different dimensions to the same dimension through normalization, and output a set of coupled operation parameters. The preliminary three-dimensional strain field inversion results are coupled with the coupling operation parameter set to establish a strain-load time domain correlation model. The cross-correlation coefficient R between each strain component and the load curve is calculated. R ≥ 0.8 is set as an effective correlation, and the strain-load correlation matrix is ​​output. Based on the strain-load correlation matrix, abnormal strain components with cross-correlation coefficient R < 0.8 or strain values ​​exceeding the preset range are marked. The rationality of the abnormal marking is verified by combining the stress concentration law in the impact area, and an abnormal strain marking map is output. Abnormal components in the abnormal strain marker map are removed, and weighted interpolation is used to complete the data in the removed area. Finally, the three-dimensional strain field of the carbon fiber laminate that is consistent with the dynamic response of the impact event is generated.

[0009] Another embodiment of this application provides a three-dimensional strain field inversion system for laminate impact events based on binocular vision, the system comprising: The acquisition module is used to set dynamic tracking markers on the surface of the impacted carbon fiber laminate. It uses a high-speed binocular camera to synchronously acquire image data of the entire impact process, and records the magnitude and duration of the impact load. The binocular camera intrinsic parameter calibration algorithm unifies the image coordinate system and generates a binocular visual image sequence with impact time information. The extraction module is used to extract the dynamic features of the marker points and the fiber texture of the laminate based on the binocular vision image sequence using a hierarchical feature matching algorithm, calculate the three-dimensional spatial coordinates of the feature points in each time sequence through stereo matching, and solve the three-dimensional displacement vector of the feature points by combining the difference between adjacent time sequence coordinates to generate three-dimensional displacement field data of the impact area of ​​the carbon fiber laminate. The calling module is used to call the anisotropic constitutive model of carbon fiber laminate, substitute the three-dimensional displacement field data into the strain tensor solution formula, convert it into local strain components through spatial derivative calculation, and combine the laminate layup direction to correct the strain calculation deviation and generate preliminary three-dimensional strain field inversion results. The generation module is used to couple the preliminary three-dimensional strain field inversion results with the impact load parameters and the laminate material property parameters, and eliminate abnormal strain components through strain-load correlation verification to generate the final three-dimensional strain field inversion result of carbon fiber laminate that is consistent with the dynamic response of the impact event.

[0010] Another embodiment of this application provides a storage medium storing a computer program, wherein the computer program is configured to execute the method described in any of the preceding claims when running.

[0011] Another embodiment of this application provides an electronic device including a memory and a processor, wherein the memory stores a computer program and the processor is configured to run the computer program to perform the method described in any of the preceding claims.

[0012] Compared with existing technologies, this invention provides a three-dimensional strain field inversion method for laminate impact events based on binocular vision. Dynamic tracking markers are placed on the surface of the impacted carbon fiber laminate, and image data of the entire impact process is simultaneously acquired using a high-speed binocular camera to generate a binocular vision image sequence with impact time-series information. Based on the binocular vision image sequence, three-dimensional displacement field data of the impact region of the carbon fiber laminate is generated. The anisotropic constitutive model of the carbon fiber laminate is called to generate preliminary three-dimensional strain field inversion results. The preliminary three-dimensional strain field inversion results are coupled with impact load parameters and laminate material property parameters to generate the final three-dimensional strain field inversion result of the carbon fiber laminate, consistent with the dynamic response of the impact event. This enables accurate full-field characterization of the dynamic strain of the impact event, providing reliable data support for the impact resistance performance evaluation and structural optimization of laminates. Attached Figure Description

[0013] Figure 1 A hardware structure block diagram of a computer terminal for a three-dimensional strain field inversion method for laminate impact events based on binocular vision, provided in an embodiment of the present invention; Figure 2 A flowchart illustrating a method for inverting the three-dimensional strain field of a laminate impact event based on binocular vision, provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the structure of a three-dimensional strain field inversion system for laminate impact events based on binocular vision, provided as an embodiment of the present invention. Detailed Implementation

[0014] The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0015] The present invention first provides a method for inverting the three-dimensional strain field of a laminate impact event based on binocular vision. This method can be applied to electronic devices, such as computer terminals, specifically ordinary computers.

[0016] The following detailed explanation uses a computer terminal as an example. Figure 1 This is a hardware structure block diagram of a computer terminal for a three-dimensional strain field inversion method for laminate impact events based on binocular vision, provided as an embodiment of the present invention. Figure 1 As shown, the computer device includes a processor, memory, and network interface connected via a system bus, wherein the memory may include non-volatile storage media and internal memory.

[0017] The non-volatile storage medium can store the operating system and computer program. The computer program includes program instructions that, when executed, cause the processor to perform any binocular vision-based three-dimensional strain field inversion method for laminate impact events.

[0018] The processor provides computing and control capabilities, supporting the operation of the entire computer device.

[0019] The internal memory provides an environment for the execution of computer programs in non-volatile storage media. When the computer program is executed by the processor, it enables the processor to execute any three-dimensional strain field inversion method for laminate impact events based on binocular vision.

[0020] This network interface is used for network communication, such as sending assigned tasks. Those skilled in the art will understand that... Figure 1 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0021] It should be understood that the processor can be a Central Processing Unit (CPU), but it can also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. Among these, a general-purpose processor can be a microprocessor or any conventional processor.

[0022] See Figure 2 The present invention provides a method for inverting the three-dimensional strain field of a laminate impact event based on binocular vision, which may include the following steps: S201: Dynamic tracking markers are placed on the surface of the impacted carbon fiber laminate. A high-speed binocular camera is used to synchronously acquire image data of the entire impact process. At the same time, the magnitude of the impact load and the duration of the impact are recorded. The image coordinate system is unified through the binocular camera intrinsic parameter calibration algorithm to generate a binocular visual image sequence with impact time information. Specifically, the dynamic tracking markers can be designed as concentric rings with alternating fluorescent and black colors, with a ring diameter of 2mm and a ring spacing of 0.5mm. They can be printed on a 50μm thick polyimide film using UV-curable adhesive to ensure that the markers have both high contrast and anti-motion blur characteristics under high-speed shooting, and output the marker placement scheme. The core of this step is designing marker points adapted for high-speed binocular vision acquisition. This requires addressing the issues of marker points becoming blurry and exhibiting low contrast during high-speed imaging, leading to feature extraction failures. Simultaneously, it's crucial to ensure that the marker points themselves do not affect the impact deformation of the laminate (e.g., thinness and good adhesion). Clearly defining the marker point structural parameters, material selection criteria, and performance verification standards is essential to developing an executable deployment plan.

[0023] I. Design and Principle of Marker Point Structure Parameters: Instead of a single dot, a concentric ring structure alternating between fluorescent and black was chosen. The core principle is to improve feature matching accuracy by utilizing the "gradient changes at the ring edges." The specific parameters are designed as follows: The ring diameter is 2mm: Combining the resolution of the high-speed binocular camera (e.g., 2048×2048 pixels, lens focal length 50mm, working distance 1m), the pixel size of the marker point in the image is calculated using the imaging formula: Imaging pixel size = (actual size of marker point × camera resolution) / (working distance × field of view corresponding to lens focal length). Substituting this into the formula, we get (2mm×2048 pixels) / (1000mm×40mm) ≈ 102 pixels (the field of view is calculated from a focal length of 50mm and a working distance of 1m, resulting in a horizontal field of view ≈ 45°, corresponding to a width of 40mm). The 102-pixel ring can still clearly distinguish the edges at 1000fps high-speed shooting, avoiding insufficient pixels due to its small size. Ring spacing 0.5mm: Fluorescent rings and black rings are arranged alternately with a ring spacing of 0.5mm to ensure that adjacent rings form a significant grayscale difference in the image (fluorescent ring grayscale value 200-255, black ring 30-50), with a contrast ratio ≥150, which meets the contrast requirements for feature extraction under high-speed imaging; if the spacing is too small (e.g., 0.2mm), the ring edges are easily blurred due to high-speed motion, and if the spacing is too large (e.g., 1mm), the overall size of the marker point is too large, occupying too many layers of plywood surface space and affecting the observation of fiber texture.

[0024] II. Marker Material Selection and Printing Process: The material needs to balance "high contrast, resistance to deformation, and no impact on the mechanical properties of the laminate". Specific selection criteria include: Substrate material: 50μm thick polyimide film, which has high toughness (elongation at break ≥70%) and low modulus (elastic modulus ≈2.5GPa, much lower than 150GPa of carbon fiber laminate). After being bonded to the surface of the laminate, it will not affect the impact deformation of the laminate due to its own stiffness; at the same time, it is resistant to high temperature (long-term operating temperature 260℃) and adaptable to the local temperature rise (≤50℃) that may occur in impact tests. Printing material: The fluorescent ring uses UV-curable fluorescent ink (wavelength 532nm, brightness ≥500cd / m² after excitation).2 The black ring is made of UV-curable black ink (light-shielding rate ≥95%), and both are printed on the film surface using UV-curable adhesive (model 3M DP460, shear strength ≥20MPa). The advantages of UV-curable adhesive are fast curing speed (curing in 30 seconds under UV light) and high bonding strength, ensuring that the marking points do not fall off when the laminate is impacted and deformed (maximum strain ≤5%). Printing process: High-precision screen printing (400 mesh screen) is used, with a printing accuracy of ±0.05mm, ensuring that the actual deviation of the ring diameter and spacing is ≤0.1mm, avoiding the impact of dimensional errors on the three-dimensional coordinate calculation.

[0025] III. Performance Verification and Deployment Scheme Output of Marker Points: The "anti-motion blur" and "high contrast" properties were verified through high-speed imaging tests: Anti-motion blur test: The marker point moving at a speed of 5m / s was photographed with a 1000fps camera (simulating the maximum surface velocity of the laminate under impact). The exposure time was set to 1μs (minimum exposure time of high-speed camera). The edge of the ring in the captured image was clear and there was no obvious ghosting (edge ​​blur width ≤ 2 pixels), which met the anti-blur requirements. Contrast test: Under different lighting conditions (natural light, backlight), the grayscale difference between the fluorescent ring and the black ring was measured. The difference was ≥150, which meets the contrast threshold (≥100) for feature extraction.

[0026] The final output marker placement scheme includes: "marker structure (2mm diameter concentric rings, 0.5mm ring spacing), material list (50μm polyimide film, UV-cured fluorescent / black ink, 3M DP460 adhesive), printing process (400 mesh screen printing, UV curing for 30s), and performance indicators (contrast ratio ≥150, anti-blurring edge width ≤2 pixels)".

[0027] According to the marking point layout plan, the surface of the carbon fiber laminate is divided into an impact center area and an edge area. Marking points are pasted in the center area with a grid density of 5mm×5mm and in the edge area with a grid density of 10mm×10mm. The height difference of the marking points is controlled within 0.1mm using a laser thickness gauge. The laminate specimen with the marking points laid out is then output. This step requires differentiating the placement of marker points based on the impact deformation characteristics of the laminate (larger deformation in the central area and smaller deformation in the edge area), while controlling the height difference of the marker points to avoid affecting the three-dimensional coordinate calculation, so as to ensure that the specimens after placement meet the test accuracy requirements.

[0028] I. Impact Zone Division and Mesh Density Design: The impact deformation of carbon fiber laminates is characterized by a large strain gradient in the central region (around the impact point) and a gentler strain gradient in the edge region. Therefore, a differentiated mesh density is adopted: Laminate specimen size: Select the common impact test size 300mm×300mm×3mm (carbon fiber T700 / epoxy resin, layup method [0° / 90°]4s), and set the impact center as the geometric center of the specimen (150mm×150mm). Region division: The impact center zone is a circular area with a radius of 50 mm centered on the impact center (area ≈ 7854 mm²). 2 The edge area is the remaining area outside the central area (area ≈ 300 × 300 - 7854 = 82146 mm). 2 ); Mesh density: 5mm×5mm mesh in the center area - 1 marker point is pasted on each mesh node, and the number of mesh nodes is (100mm / 5mm+1)×(100mm / 5mm+1)=21×21=441 (center area diameter 100mm). High density marker points can capture the details of severe deformation in the center area; 10mm×10mm mesh in the edge area - The number of nodes is (300mm / 10mm+1)×(300mm / 10mm+1)-441=31×31-441=961-441=520. Low density marker points can reflect edge deformation and avoid too many marker points from obscuring the fiber texture.

[0029] II. Marker point pasting process and height difference control: During the pasting process, ensure that the marked points are in contact with the surface of the laminate, with a height difference ≤ 0.1mm (to avoid deviations in 3D coordinate calculations due to height differences). Surface pretreatment: Wipe the surface of the laminate with alcohol to remove oil and dust, and then treat with plasma (500W power, 30s time) to improve the surface roughness (Ra from 0.05μm to 0.2μm) and enhance adhesion; Positioning and pasting: Using a CNC positioning platform (positioning accuracy ±0.1mm), mark the pasting position on the laminate surface according to the grid coordinates. Cut the polyimide film with the printed marks into small pieces of 5mm×5mm (corresponding to a single mark). Apply UV-curing adhesive to the back of the film (thickness ≤0.05mm), align it with the mark position and paste it. After pasting, apply pressure with a 200g weight for 10s to ensure adhesion. Height difference detection: A laser thickness gauge (model Keyence LK-G80, resolution 0.01mm, measurement range 0-10mm) is used to measure the height of the highest point of each marker relative to the surface of the laminate. For example, if the height of a marker in the center area is 0.08mm and the height of a marker in the edge area is 0.05mm, the maximum height difference is 0.03mm≤0.1mm, which meets the requirements. If the height of a marker is 0.12mm, it needs to be removed and re-attached to ensure that the height difference meets the standard.

[0030] III. Status output of the completed specimen setup: The output laminate specimen must include: "Specimen size 300mm×300mm×3mm, layup [0° / 90°] 4s, total number of markers 441+520=961, 5mm×5mm grid in the center area (441 markers), 10mm×10mm grid in the edge area (520 markers), height difference of markers ≤0.03mm, clean surface free of oil stains", and a test report from a laser thickness gauge (height data of each marker) to ensure traceability of subsequent tests.

[0031] The Zhang Zhengyou calibration method was used to perform intrinsic parameter calibration on the high-speed binocular camera. Twenty sets of images with different poses were acquired using a 12×9 checkerboard target. The focal length principal point coordinate distortion coefficients of the two cameras were calculated. The frame rate of the two cameras was locked at 1000fps and the trigger time difference was controlled within 1μs through the synchronous trigger module. The calibrated binocular camera system was then output. Binocular camera intrinsic parameter calibration is the core to ensure the accuracy of 3D coordinate calculation. Zhang Zhengyou's calibration method achieves high-precision calibration through a planar target, while the synchronous triggering module ensures that the timing of the images from the two cameras is consistent, avoiding stereo matching errors caused by asynchrony.

[0032] I. Implementation details of Zhang Zhengyou's calibration method: Zhang Zhengyou's calibration method is based on solving intrinsic parameters through "imaging of a planar target in different postures". Specific steps include: Target selection and parameters: Use a 12×9 checkerboard target, with each square measuring 20mm×20mm (dimensional accuracy ±0.01mm). The checkerboard is printed with alternating black and white (contrast ≥200). The target thickness is 5mm (to ensure flatness ≤0.05mm and avoid target bending affecting calibration). Image Acquisition: A high-speed binocular camera (Basler acA2040-90um, 2048×2048 pixel resolution, maximum frame rate 1000fps) was fixed on a tripod with a 500mm distance between the two cameras (baseline distance, optimized based on a 1m working distance to ensure a reasonable stereo matching parallax range). Twenty sets of target images in different poses were acquired—poses included: 0° (frontal), 15°, 30°, and 45° angles between the target plane and the camera's optical axis; the target was moved to the upper left, upper right, lower left, lower right, and center positions of the image. 3-5 frames were captured for each pose, for a total of 20 sets (ensuring coverage of the entire camera field of view to improve calibration accuracy). Intrinsic parameter calculation: Using the `calibrateCamera` function in OpenCV (based on Zhang Zhengyou's calibration algorithm), process 20 sets of images to calculate the intrinsic parameters of the two cameras. Left camera intrinsic parameters: focal length f_x1 = 1005 pixels, f_y1 = 1003 pixels (f_x and f_y are the focal lengths along the x and y axes, respectively; since the camera lens has no distortion, they are close); principal point coordinates (x01 = 1024 pixels, y01 = 1024 pixels) (the principal point is close to the image center, which is consistent with the camera design); radial distortion coefficients k1 = -0.012, k2 = 0.008; tangential distortion coefficients p1 = 0.001, p2 = -0.0005 (small absolute values ​​of distortion coefficients indicate slight lens distortion); Right camera intrinsic parameters: f_x2=1004 pixels, f_y2=1002 pixels, principal point (x02=1025 pixels, y02=1023 pixels), distortion coefficients k1=-0.011, k2=0.007, p1=0.0008, p2=-0.0006; After calibration, if the reprojection error is ≤0.5 pixels (reprojection error is the core indicator of calibration accuracy, and ≤1 pixel meets the requirements of high-speed vision testing), the calibration is deemed qualified.

[0033] II. Debugging and parameter locking of the synchronous trigger module: The two cameras need to acquire data synchronously to avoid the scene at the same moment being captured in different frames due to timing differences. The synchronization trigger module (National Instruments NI-9403) is implemented as follows: Frame rate locking: The trigger module sends a synchronization clock signal to the two cameras to lock the frame rate at 1000fps (frame interval 1ms, matching the dynamic response time of the laminate impact, the impact time is about 5ms, and 5 frames of data are needed to capture the impact process), with frame rate fluctuation ≤1fps, to ensure that the imaging rhythm of the two cameras is consistent. Trigger time difference control: Use an oscilloscope (Tektronix MDO3024, 200MHz bandwidth) to measure the trigger signals of the two cameras, adjust the delay parameters of the trigger module, and control the trigger time difference between the two cameras within 1μs (1μs is much smaller than the frame interval of 1ms and can be ignored, ensuring that the scene at the same moment is captured by both cameras simultaneously); if the time difference is 5μs, the trigger delay of the right camera needs to be adjusted to reduce the time difference to 0.8μs, which meets the requirements.

[0034] III. Output of the calibrated binocular camera system: The output system parameters include: "Left camera intrinsic parameters (f_x1=1005, f_y1=1003, x01=1024, y01=1024, k1=-0.012, etc.), right camera intrinsic parameters (f_x2=1004, f_y2=1002, etc.), baseline distance 500mm, frame rate 1000fps, synchronization time difference 0.8μs, reprojection error 0.4 pixels", and the intrinsic parameter file (XML format) is saved for subsequent use in image coordinate system one.

[0035] The laminate specimen with marked points is fixed on the impact test bench. The calibrated binocular camera system is started, and the drop hammer impact device is triggered synchronously. The magnitude and duration of the impact load are recorded by the force sensor. Binocular images of the entire impact process are acquired. The coordinate system of the two cameras is unified by the internal parameter data, and a binocular visual image sequence with impact time information is generated.

[0036] This step is the core of "hardware linkage and data acquisition". It requires the synchronous triggering of "camera-impact device-force sensor" and the use of an internal parameter unified coordinate system to ensure that the acquired image sequence has temporal and spatial consistency, laying the foundation for subsequent displacement field calculation.

[0037] I. Specimen Fixation and Impact Device Adjustment: Specimen Fixing: Fix the laminate specimen onto the clamp of the drop hammer impact test bench (model Instron CEAST 9350). The clamp adopts a vacuum adsorption type (adsorption force ≥500N) to ensure that the specimen does not loosen during the impact (loosening will cause additional displacement and affect strain inversion). The surface of the specimen is perpendicular to the optical axis of the camera (angle ≤1° to avoid perspective distortion in imaging). Impact device parameter settings: The drop hammer mass is set to 5kg (selected based on the laminate thickness of 3mm to ensure an impact energy of 5J, producing significant deformation but not damage), and the drop hammer height is 100mm (calculated using the energy formula E=mgh, E=5kg×9.8m / s). 2×0.1m=4.9J≈5J), impact head diameter 12mm (hemispherical, conforming to ASTM D7136 impact test standard).

[0038] II. Multi-device synchronous triggering and data acquisition: Synchronization of the camera, impact sensor, and force sensor is achieved through a trigger controller (integrated with a synchronization trigger module). Triggering logic: The trigger controller outputs three synchronous signals: the first signal starts the binocular camera to acquire data (starting 10ms in advance to capture the initial state before impact), the second signal releases the drop hammer (impact begins), and the third signal starts the force sensor to record data (synchronized with the release of the drop hammer). Force sensor data recording: A piezoelectric force sensor (model Kistler 9257B, range 0-20kN, sensitivity 50mV / kN, sampling frequency 10kHz) was installed on the top of the impact head to record the magnitude and duration of the impact load. Example data: maximum impact load 5.2kN, duration 4.8ms (load curve is a sine wave, rise time 1.2ms, peak time 0.5ms, fall time 3.1ms). The data was stored on a DAQ card (sampling frequency 10kHz), and the timestamp was synchronized with the camera image (error ≤1μs). Image acquisition: The camera starts acquiring images 10ms before impact (t=-10ms) and ends 20ms after impact (t=20ms), acquiring a total of 30,000 frames (1000fps×30ms). Each frame is 2048×2048 pixels in size and stored in RAW format (uncompressed, retaining the original grayscale information). Each frame is accompanied by a timestamp (accurate to μs, t=0 is the time when the impact begins).

[0039] III. Image Coordinate System and Sequence Generation: By applying camera intrinsics, the image from the right camera is transformed to the coordinate system of the left camera, thus eliminating the positional difference between the two cameras: Coordinate system method 1: The two camera extrinsic parameters (rotation matrix R, translation vector T) are calculated using Zhang Zhengyou's calibration method. The pixel coordinates (u_r, v_r) of the right camera image are converted into coordinates (u_l', v_l') in the left camera coordinate system. The conversion formula is based on the epipolar constraint of stereo vision and corrects the coordinate deviation caused by the different positions of the two cameras. Distortion correction: Apply the distortion coefficients from the intrinsic parameters to the left and right camera images respectively to remove radial and tangential distortion. For example, if a pixel in the left camera is offset by 2 pixels due to distortion, it will be corrected and returned to its correct position. Sequence generation: The corrected left and right camera images are sorted by timestamp, and each frame is labeled with "timestamp (e.g., t=-10ms, t=-9ms, ..., t=20ms), impact state (before impact / during impact / after impact), and load value (load size corresponding to the timestamp, e.g., load 0kN at t=0ms, load 3kN at t=1ms)", generating a binocular visual image sequence with impact timing information. For example, frame 1 (t=-10ms): left camera image (corrected), right camera image (corrected), load 0kN; frame 11 (t=0ms): left / right camera images, load 0kN (impact start); frame 16 (t=5ms): left / right camera images, load 0kN (impact end).

[0040] The final output binocular vision image sequence contains 30,000 corrected images from the left and right cameras. Each frame includes a timestamp, impact state, load value, and coordinate system 1 (with the left camera as the reference), providing high-quality data for subsequent feature extraction and displacement field calculation.

[0041] S202, based on the binocular vision image sequence, a hierarchical feature matching algorithm is used to extract the dynamic features of the marker points and the fiber texture of the laminate. The three-dimensional spatial coordinates of the feature points in each time sequence are calculated by stereo matching. The three-dimensional displacement vector of the feature points is solved by combining the difference between adjacent time sequence coordinates, and the three-dimensional displacement field data of the impact area of ​​the carbon fiber laminate is generated. Specifically, the binocular vision image sequence can be preprocessed, adaptive bilateral filtering can be used to remove high-speed imaging noise, HSV color space threshold segmentation can be performed on the marked point region, Canny edge detection can be performed on the fiber texture region, and the preprocessed image sequence can be output. During high-speed acquisition, binocular vision image sequences are susceptible to problems such as imaging noise (e.g., sensor thermal noise), insufficient contrast between marker points and background, and blurred fiber texture edges, leading to low accuracy in subsequent feature extraction. The core of preprocessing is "denoising while preserving edges, enhancing marker point features, and highlighting texture edges." Differentiated processing methods should be used for different regions (marker points, fiber textures) to ensure that the output of each step can directly support subsequent feature matching.

[0042] I. Adaptive bilateral filtering for noise reduction: High-speed camera (1000fps) sensors are prone to thermal noise (random fluctuations in grayscale values ​​of ±5-10) at high frame rates. Traditional Gaussian filtering blurs edges, while adaptive bilateral filtering balances noise reduction and edge preservation through "spatial domain weights + range domain weights." The specific implementation is as follows: Filter kernel parameter settings: Set the filter kernel size to 5×5 (to balance denoising effect and computational efficiency; a kernel that is too large will blur details, while a kernel that is too small will not denoise thoroughly), spatial domain standard deviation σ_d=5 (to control the influence range of the spatial neighborhood; the larger the value, the higher the weight of distant pixels), range domain standard deviation σ_r=20 (to control the influence of grayscale similarity; the larger the value, the higher the weight of pixels with large grayscale differences). Adaptive adjustment logic: For regions with large gray-level gradients in the image (such as marker edges and fiber textures), σ_r is automatically reduced to 10 (enhancing gray-level similarity filtering and preserving edges); for regions with smooth gray levels (such as laminate backgrounds), σ_r is increased to 30 (enhancing noise reduction). For example, if the gray-level gradient of a marker edge is 150 (the gray-level difference between the fluorescent ring and the black ring), after system recognition, σ_r is adjusted to 10, and the gray-level gradient of the edge remains at 145 after filtering, resulting in clear edges; if the gray-level gradient of the background region is 10, σ_r is adjusted to 30, and the gray-level fluctuation is reduced from ±8 to ±2 after filtering, significantly reducing noise. Denoising effect verification: The peak signal-to-noise ratio (PSNR) of the image was calculated and evaluated. Before preprocessing, PSNR = 28dB, and after preprocessing, PSNR = 35dB (PSNR ≥ 30dB meets the feature extraction requirements), so the denoising was deemed qualified.

[0043] II. HSV color space threshold segmentation of marked point regions: The markers have an alternating fluorescent-black structure, which is susceptible to reduced contrast in the RGB space due to lighting conditions. The HSV space (hue H, saturation S, lightness V) can more accurately distinguish color characteristics. Specific implementation details are as follows: HSV threshold range determined: Through sample calibration experiments, the HSV range of the fluorescent ring is: H=50-70 (yellow fluorescence, corresponding to a wavelength of 532nm), S=80-100 (high saturation, to avoid confusion with light-colored background areas), V=200-255 (high brightness, highlighting fluorescence characteristics); the HSV range of the black ring is: H=0-180 (no specific hue), S=0-20 (low saturation), V=0-50 (low brightness). Segmentation execution: For the image after adaptive bilateral filtering, first convert the RGB space to HSV space (conversion formula: H=arctan2 (GB,G+B)×(180 / π), S=1-3×min (R,G,B) / (R+G+B), V=(R+G+B) / 3, where R, G, and B are gray values ​​of 0-255), and then segment the fluorescent ring and black ring regions according to the above threshold range to generate a binary image (the marked area is white and the background is black). Segmentation effect verification: Statistical analysis of 100 frames of images showed that the segmentation accuracy of the marked point region was ≥98% (i.e., 98% of the marked point pixels were correctly identified), and the misidentification rate was ≤1% (the proportion of background pixels that were misidentified as marked points), which meets the region localization requirements for subsequent ORB feature extraction.

[0044] III. Improved Canny edge detection for fiber texture regions: The carbon fiber laminate has a parallel, elongated fiber texture. Traditional Canny edge detection is prone to missed detections due to discontinuous texture edges. The improved solution enhances detection performance through "multi-scale gradient fusion + adaptive thresholding," specifically implemented as follows: Multi-scale gradient calculation: Image gradients are calculated using Sobel operators at two scales, 3×3 and 5×5 (3×3 operator detects fine textures, 5×5 operator detects coarse textures), and the gradient values ​​at the two scales are fused according to weights (3×3 gradient weight 0.6, 5×5 gradient weight 0.4) to enhance the continuity of texture edges; Adaptive threshold setting: Traditional Canny uses fixed high and low thresholds (e.g., 200, 100). The improved version calculates the optimal threshold for the image using the Otsu algorithm, and then determines the high and low thresholds proportionally (high threshold = Otsu threshold × 1.2, low threshold = Otsu threshold × 0.5). For example, for a frame of fiber texture image, the Otsu threshold is 150, the high threshold is 180, and the low threshold is 75, which can suppress false edges caused by noise while preserving weak edges of the texture. Edge thinning: For the detected edge image, a morphological thinning algorithm (skeleton extraction) is used to thin wide edges (3-5 pixels) into single-pixel edges, which facilitates subsequent HOG feature extraction. For example, the original texture edge width is 4 pixels, and after thinning it is 1 pixel, improving the edge localization accuracy to ±0.5 pixels.

[0045] The final output preprocessed image sequence contains "denoised grayscale image, marker point binary segmentation image, and fiber texture edge image" in each frame, providing clear and accurate feature regions for subsequent hierarchical feature matching.

[0046] Based on the preprocessed image sequence, a hierarchical feature matching algorithm is adopted. The upper layer extracts ORB features from the marker points and performs brute-force matching, while the lower layer extracts histogram features of orientation gradients from the fiber texture and performs K-nearest neighbor matching to generate a sequence of matching pairs between marker points and texture features. Layered feature matching takes into account the different characteristics of "marker points (high contrast, discrete distribution)" and "fiber texture (low contrast, continuous distribution)", and adopts differentiated feature extraction and matching algorithms to avoid missed detections or false matches caused by a single algorithm, and ensure that feature points can still be stably matched throughout the entire impact process (when deformation is severe).

[0047] I. Upper Layer: ORB Feature Extraction and Brute-Force Matching of Marker Points ORB (Oriented FAST and Rotated BRIEF) features combine the speed of FAST corner detection with the robustness of BRIEF descriptors, making them suitable for marker matching in high-speed, dynamic scenes. Implementation details: ORB Feature Extraction: Corner detection: The FAST-9 algorithm is used (to detect corners where the gray-level difference between 9 consecutive neighboring pixels is greater than a threshold). The corner response threshold is set to 20 (to control the number of corners and avoid too many redundant corners). In the binary segmentation map of the marker points, corners are extracted only from the white areas (marked points). 8-12 corners are extracted from each marker point (to ensure feature richness). For example, 10 corners are extracted from a fluorescent ring area and distributed on the edge of the ring. Descriptor generation: For each corner point, calculate its principal direction (based on the gray centroid of the corner point's neighborhood), and then generate a 256-bit BRIEF descriptor (generated by randomly selecting 128 pairs of pixels for gray-level comparison). The descriptor has rotation invariance (due to principal direction alignment) and scale invariance (achieved through an image pyramid with 4 layers and a scale factor of 1.2). Brute-Force Matching: Matching logic: For the left and right images of the stereo camera, calculate the Hamming distance (measures the difference between binary descriptors, the smaller the value, the higher the matching degree) between the ORB descriptor of the left image and all ORB descriptors of the right image. Match filtering: A Hamming distance threshold of 50 is set (through experimental calibration, matching pairs with a Hamming distance ≤ 50 are considered valid matches). Cross-validation is used (only matches where the best result of the left image matches the best result of the right image, and vice versa), to eliminate one-way mismatches. For example, the Hamming distance between the descriptor of marker A in the left image and marker B in the right image is 42, and the cross-validation shows they match, so they are retained as a matching pair; the Hamming distance between the descriptor of marker C in the left image and marker D in the right image is 65, exceeding the threshold, so they are discarded. Matching effect verification: The number of matching pairs of marker points in each frame of the image is ≥90% (that is, 90% of the marker point corner points are successfully matched), and the matching accuracy is ≥95% (verified by reprojection error, the error is ≤1 pixel), which meets the requirements of stereo matching.

[0048] II. Lower Layer: HOG Feature Extraction and K-Nearest Neighbor Matching of Fiber Texture HOG (Histogram of Oriented Gradients) features can effectively describe the directional distribution of textures and are suitable for matching continuous fiber textures. Specific implementation details are as follows: HOG Feature Extraction: Image segmentation: Divide the fiber texture edge map into “cell units (8×8 pixels) - blocks (2×2 cell units)”, with a block step size of 4 pixels (adjacent blocks overlap by 50% to enhance feature continuity). For example, a 2048×2048 image is divided into 256×256 cell units and 253×253 blocks. Gradient calculation: For each cell unit, calculate the gradient direction (0-180°, divided into 18 intervals of 10°) and gradient magnitude of the pixel to generate an 18-dimensional gradient histogram; for each block, concatenate the histograms of 4 cell units to generate a 72-dimensional block feature; for the entire image, concatenate all block features to generate 253×253×72≈4.6×10 6 1×10 HOG feature vectors (in practical applications, only texture regions are extracted, reducing the dimension to 1×10). 6 ); K-Nearest Neighbor (K-NN) matching: K value setting: K=2 (select the 2 candidate matching pairs with the highest matching degree and filter by distance ratio) to avoid false matching when K=1; Matching and filtering: Calculate the Euclidean distance between the HOG features of the left and right images. Let the distances between the first two candidate matching pairs be d1 and d2. If d1 / d2 ≤ 0.8 (distance ratio threshold, determined experimentally), then retain the matching pair corresponding to d1; otherwise, discard it. For example, the candidate matching pairs for texture block P in the left image are block Q (d1=10) and block R (d2=13) in the right image. If d1 / d2 ≈ 0.77 ≤ 0.8, retain the PQ matching pair; if d2 = 12.6, d1 / d2 ≈ 0.79 ≤ 0.8, still retain it; if d2 = 12.7, d1 / d2 ≈ 0.787 ≤ 0.8, retain it, ensuring matching robustness. Matching effect verification: The density of matching pairs of fiber textures is ≥50 / mm. 2 (That is, at least 50 matching pairs per square millimeter), and the orientation consistency of the matching pairs is ≥90% (the proportion of matching pairs with the same texture orientation), to ensure that subsequent stereo matching can cover the entire texture area.

[0049] III. Generation of matching sequences: The ORB matching pairs of marker points and the HOG matching pairs of fiber textures are organized according to "time sequence + left and right camera coordinates". A matching pair set is generated for each frame of image. For example, "time sequence t=0ms (impact start): 150 marker point matching pairs (left coordinate (u1,v1), right coordinate (u2,v2)), 20,000 texture matching pairs (left coordinate (u3,v3), right coordinate (u4,v4))". All frames are sorted by timestamp to form a matching pair sequence from t=-10ms to t=20ms, which provides a continuous and complete feature correspondence for subsequent stereo matching.

[0050] The matching sequence is input into the stereo matching algorithm. The three-dimensional spatial coordinates of the feature points at each time step are calculated by triangulation in combination with the intrinsic parameters of the binocular camera. The coordinate accuracy is optimized by the bundle adjustment method, and the three-dimensional coordinate set of the feature points at each time step is output. The core of stereo matching is to use the parallax information of a binocular camera to convert two-dimensional image matching pairs into three-dimensional spatial coordinates. Triangulation is a classic method to achieve this conversion. Bundle adjustment is used to eliminate intrinsic parameter errors and matching errors, improve coordinate accuracy, and ensure that the three-dimensional coordinates can accurately reflect the actual position of the laminate.

[0051] I. Preparation of internal parameters for stereo matching: From the calibrated binocular camera system, extract key intrinsic and extrinsic parameters: left camera focal length f_x1=1005 pixels, f_y1=1003 pixels, principal point coordinates (x01=1024, y01=1024); right camera focal length f_x2=1004 pixels, f_y2=1002 pixels, principal point coordinates (x02=1025, y02=1023); baseline distance between the two cameras b=500mm (horizontal distance between the optical centers of the left and right cameras); extrinsic rotation matrix R is an identity matrix (the optical axes of the two cameras are parallel); translation vector T=(-b,0,0) (the right camera is at position b to the left of the left camera). These parameters need to be stored in advance in the parameter library of the stereo matching algorithm for coordinate calculation.

[0052] II. Calculation of three-dimensional spatial coordinates using triangulation: Triangulation is based on a "left and right camera projection model." It calculates the three-dimensional coordinates (X, Y, Z) of a point in space by matching the disparity (the difference in pixel coordinates between the left and right images). Specific implementation: Pixel coordinate normalization: Convert the matching pixel coordinates of the left and right images into normalized image coordinates (eliminating the influence of focal length and principal point). The normalized coordinates of the left image are (u1', v1') = ( (u1-x01) / f_x1, (v1-y01) / f_y1 ), and the normalized coordinates of the right image are (u2', v2') = ( (u2-x02) / f_x2, (v2-y02) / f_y2 ). Parallax calculation and depth determination: Since the optical axes of the two cameras are parallel, the parallax d = u1' - u2' (only horizontal parallax is considered, vertical parallax has been eliminated due to calibration). The depth Z of the spatial point (along the direction of the camera optical axis, i.e. the distance from the laminate to the camera) is calculated by the formula Z = f_x1 × b / d. 3D coordinate calculation: The formulas for calculating the X (horizontal) and Y (vertical) coordinates of a spatial point in the left camera coordinate system are X = Z × u1', Y = Z × v1'. For example, the left pixel coordinates (u1 = 1050, v1 = 1024) and right pixel coordinates (u2 = 1000, v2 = 1024) of a marker point matching pair are substituted into the intrinsic parameters for calculation: Normalized coordinates: u1'=(1050-1024) / 1005≈0.0259, v1'=(1024-1024) / 1003=0; u2'=(1000-1025) / 1004≈-0.0249, v2'=0; The parallax d = 0.0259 - (-0.0249) = 0.0508; Depth Z = 1005 × 500 mm / 0.0508 ≈ 1005 × 500 / 0.0508 ≈ 9.9 × 10 6 mm = 9900mm = 9.9m (consistent with the camera's working distance of 10m, with an error of ≤1%). 3D coordinates: X=9900mm×0.0259≈256.4mm, Y=9900mm×0=0mm, Z=9900mm, that is, the coordinates of the marker point in the left camera coordinate system are (256.4, 0, 9900) mm.

[0053] III. Optimizing coordinate accuracy using bundle adjustment method: The coordinate accuracy of triangulation is affected by intrinsic parameter errors and matching errors. Bundle adjustment optimizes the 3D coordinates and camera parameters by minimizing reprojection errors. Specifically: Reprojection error is defined as follows: The calculated 3D coordinates (X,Y,Z) are reprojected onto the left and right camera images to obtain the projected pixel coordinates (u1_proj,v1_proj) and (u2_proj,v2_proj). The reprojection error is the Euclidean distance between the actual matching pair pixel coordinates and the projected coordinates, i.e., e=√[(u1-u1_proj)]. 2 +(v1-v1_proj) 2 ] + √[(u2-u2_proj) 2 +(v2-v2_proj) 2 ]; The objective function is minΣe, which aims to minimize the sum of reprojection errors of all matching pairs. Optimization solution: The Levenberg-Marquardt algorithm is used for iterative solution (the number of iterations is ≤50, and the convergence threshold is 1e-6). For example, the reprojection error of a certain marker point was 0.8 pixels before optimization, and it was reduced to 0.1 pixels after optimization. The 3D coordinate error was reduced from ±5mm to ±0.5mm. Optimization effect verification: After optimization, the average reprojection error of all feature points is ≤0.3 pixels, and the three-dimensional coordinate accuracy is ≤±1mm (meets the requirements for laminate impact displacement measurement, and the impact displacement is usually ≥1mm).

[0054] IV. Output of the 3D coordinate set of each temporal feature point: The optimized 3D coordinates are organized according to "time series + feature point type". A coordinate set is generated for each frame. For example, "time series t=0ms: 150 3D coordinates of marker points (such as marker point A (256.4,0,9900) mm, marker point B (260.1,5,9900) mm), 20,000 3D coordinates of texture feature points (such as texture point P (258.2,2,9900) mm)". All time series coordinate sets are sorted by timestamp to form continuous coordinate data from t=-10ms to t=20ms, which provides a basis for subsequent displacement field calculations.

[0055] For each time series feature point's three-dimensional coordinate set, the coordinate difference between adjacent time series is calculated to obtain the three-dimensional displacement vector. The discrete vector is extended into a continuous field using the inverse distance weighted interpolation method. Gaussian smoothing is used to eliminate interpolation errors, generating three-dimensional displacement field data of the carbon fiber laminate impact region.

[0056] The three-dimensional displacement field is the direct input for strain field inversion. It is necessary to convert the displacement vectors (differences between adjacent time-series coordinates) of discrete feature points into continuous field data, while eliminating interpolation noise to ensure that the displacement field can accurately reflect the spatial distribution law of impact deformation of laminated plates.

[0057] I. Calculation of three-dimensional displacement vectors: The displacement vector reflects the positional changes of a feature point in adjacent time sequences. The calculation object is "the three-dimensional coordinates of the same feature point in different time sequences." Specific implementation: Temporal matching: For each feature point (e.g., marker point A), find its corresponding three-dimensional coordinates in the coordinate set of adjacent temporal sequences (e.g., t=0ms and t=1ms) to ensure that they are the same physical point (determined by the consistency of the ORB descriptor or the positional continuity of the texture features). Displacement component calculation: The three-dimensional displacement vector contains three components: X (horizontal), Y (vertical), and Z (depth). The calculation formulas are ΔX = X(t + Δt) - X(t), ΔY = Y(t + Δt) - Y(t), and ΔZ = Z(t + Δt) - Z(t), where Δt = 1ms (the interval between adjacent time steps, due to the camera frame rate of 1000fps). For example, the coordinates of marker point A at t = 0ms are (256.4, 0, 9900) mm, and the coordinates at t = 1ms are (257.0, 0.2, 9899.8) mm. The displacement vectors are ΔX = 0.6mm, ΔY = 0.2mm, and ΔZ = -0.2mm, indicating that the point moves 0.6mm in the positive X direction, 0.2mm in the positive Y direction, and 0.2mm in the negative Z direction (closer to the camera). Displacement vector filtering: Set displacement thresholds (e.g., the absolute values ​​of ΔX, ΔY, and ΔZ are ≤10mm, and the maximum displacement during impact is usually ≤5mm) and remove abnormal vectors that exceed the threshold (e.g., displacement of 100mm caused by matching errors). For example, the displacement of a certain texture point ΔX=15mm exceeds the threshold and is judged as abnormal, so the vector is removed.

[0058] II. Constructing a continuous displacement field using inverse distance weighted interpolation: The displacement vectors of discrete feature points cannot cover all areas of the laminate, so they need to be extended to a continuous field through interpolation. Inverse distance weighted interpolation (IDW) can assign different weights according to the distance of feature points, ensuring that the interpolation result is close to the actual deformation. Specific implementation: Interpolation grid setting: Based on the laminate size (300mm×300mm), the interpolation grid is set to 1mm×1mm (horizontal X and Y step size 1mm), with a total of 300×300=90000 grid points. For each grid point, the three displacement components ΔX, ΔY, and ΔZ need to be calculated. Weight calculation: For each grid point (x, y), search for all discrete feature points within a 50mm radius around it (ensure that the interpolation region has enough feature points to support it), calculate the Euclidean distance d_i from the grid point to each feature point, and the weight w_i = 1 / (d_i^p), where p = 2 (weight decay coefficient, the larger p is, the higher the weight of the nearest point), and the weight after normalization is w_i' = w_i / Σw_i; Displacement component interpolation: The displacement component of a grid point is the weighted average of the displacement components of surrounding feature points, i.e., ΔX_grid=Σ(w_i'×ΔX_i), and ΔY_grid and ΔZ_grid are similarly calculated. For example, there are 4 feature points around the grid point (257,1) mm, with distances of d1=1mm, d2=2mm, d3=3mm, and d4=4mm, and weights w1=1 / 1=1, w2=1 / 4=0.25, w3=1 / 9≈0.111, and w4=1 / 16≈0.0625. After normalization, w1'=1 / (1+0.25+0.111+0.0625)≈0.706, w2'≈0.176, w3'≈0.078, and w4'≈0.044; the ΔX_grid of the feature point... The values ​​are 0.6, 0.5, 0.4, and 0.3 mm respectively. After interpolation, ΔX_grid = 0.706×0.6 + 0.176×0.5 + 0.078×0.4 + 0.044×0.3 ≈ 0.424 + 0.088 + 0.031 + 0.013 ≈ 0.556 mm. Interpolation effect verification: Comparing the displacement of the interpolated grid points with the actual discrete feature point displacements, the error is ≤0.05mm, indicating that the interpolation result can accurately reflect the deformation trend.

[0059] III. Gaussian smoothing to eliminate interpolation errors: Inverse distance weighted interpolation may introduce local noise due to uneven distribution of feature points (e.g., displacement fluctuation of ±0.2mm at a certain grid point). Gaussian smoothing eliminates noise through neighborhood weighted averaging. Specific implementation details are as follows: Smoothing kernel parameter settings: Smoothing kernel size 3×3, Gaussian standard deviation σ=1.0 (controls the degree of smoothing; the larger σ is, the more obvious the smoothing), the weight of pixels within the kernel is calculated by a Gaussian function, for example, the weight matrix of a 3×3 kernel is: [0.075, 0.124, 0.075; 0.124, 0.204, 0.124; [0.075, 0.124, 0.075] (The weighted sum is 1 to ensure that the total displacement remains unchanged); Smoothing: For the interpolated displacement field, the displacement component of each grid point is equal to the sum of the products of the displacement components of all grid points in its 3×3 neighborhood and their corresponding weights. For example, the ΔX of the grid point (257,1) mm is 0.556 mm, and the ΔX of the grid points in the neighborhood are 0.55, 0.56, 0.55, 0.57, 0.556, 0.54, 0.56, 0.55, and 0.57 mm respectively. After smoothing, ΔX = 0.075×0.55 + 0.124×0.56 + ... + 0.075×0.57 ≈ 0.555 mm, and the fluctuation is reduced from ±0.2 mm to ±0.05 mm. Final displacement field output: The smoothed displacement field data contains ΔX, ΔY, and ΔZ components of 90,000 grid points, organized in time sequence (t=-10ms to t=20ms). For example, "t=1ms impact region displacement field: grid points (250-300, 0-50) mm (impact center region) ΔX=0.5-1.0mm, ΔY=0.1-0.3mm, ΔZ=-0.1-0.2mm; edge region ΔX=0.1-0.3mm", clearly reflecting the pattern of large deformation in the impact center region and small deformation in the edge region, providing continuous and smooth displacement data for subsequent strain field inversion.

[0060] S203, call the anisotropic constitutive model of carbon fiber laminate, substitute the three-dimensional displacement field data into the strain tensor solution formula, convert it into local strain components through spatial derivative calculation, and combine the laminate layup direction to correct the strain calculation deviation, and generate preliminary three-dimensional strain field inversion results. Specifically, an anisotropic constitutive model of carbon fiber laminates can be constructed by inputting ply parameters and material parameters, establishing an elastic matrix based on the ply direction, and outputting the parameter set of the anisotropic constitutive model. Carbon fiber laminates exhibit significant anisotropy due to the different fiber layup directions (the mechanical properties along the fiber direction differ greatly from those perpendicular to the fiber direction). Traditional isotropic constitutive models cannot accurately describe their mechanical response; therefore, anisotropic constitutive models need to be constructed. The core of this step is to "transform the layup structure and material properties into mathematical matrices," providing a mechanical foundation for subsequent strain calculations. This requires clearly defining the layup parameters, the basis for selecting material parameters, and the method for constructing the elastic matrix to ensure that the model can accurately reflect the mechanical behavior of the laminate.

[0061] I. Definition and selection of layup parameters and material parameters: Ply parameters: Based on common carbon fiber laminate impact test specimens, the layup configuration was set to [0° / 90°] 4s, where: “0°” indicates that the fibers are laid along the length of the laminate (X-axis of the global coordinate system), and “90°” indicates that they are laid along the width (Y-axis of the global coordinate system). “ / ” indicates the direction separation between adjacent plies, “4” indicates 4 cycles of alternating 0° and 90° (8 layers in total), and “s” indicates symmetrical plies (i.e., the ply sequence from bottom to top is 0°, 90°, 0°, 90°, 90°, 0°, 90°, 0°) to ensure that the laminate has no coupled bending deformation. The thickness of a single layer is t0 = 0.375 mm, and the total thickness is t = 8 × 0.375 = 3 mm. The distribution of each layer in the thickness direction (Z axis) is as follows: layer 1 (Z = 0~0.375 mm), layer 2 (Z = 0.375~0.75 mm), ..., layer 8 (Z = 2.625~3 mm).

[0062] Material parameters: T700 carbon fiber / epoxy resin composite material (commonly used in aerospace) was selected. Its anisotropic material parameters were measured through standard tests (ASTM D3039, ASTM D3518), and the specific values ​​are as follows: Longitudinal elastic modulus E 11 =150GPa (maximum stiffness along the fiber direction); Transverse elastic modulus E 22 =10GPa (perpendicular to the fiber direction, the stiffness is only 1 / 15 of that in the longitudinal direction). In-plane shear modulus G 12 =5GPa (shear stiffness in the XY plane); Longitudinal Poisson's ratio ν 12 =0.2 (the proportion of transverse contraction during longitudinal stretching); Transverse Poisson's ratio ν 21 =ν 12 ×(E 22 / E 11 =0.2×(10 / 150)≈0.013 (because E) 22 Much smaller than E 11 (The Poisson's ratio in the horizontal direction is much smaller than that in the vertical direction). Thickness-related parameters (simplified, impact problems are mainly characterized by in-plane strain): E 33 =8GPa, G 13 =G 23 =4GPa, ν 13 =ν 23 =0.15.

[0063] II. Construction of the elasticity matrix (based on the ply direction): The core of the anisotropic constitutive model is the elastic matrix [Q] (also known as the stiffness matrix), which describes the linear relationship between stress and strain ([σ]=[Q][ε]). It needs to be constructed according to the ply direction (local coordinate system) and is divided into the in-plane elastic matrix (describing the mechanical behavior in the XY plane, which is the main focus of impact problems) and the three-dimensional elastic matrix (including the thickness direction). Here, we focus on constructing the in-plane elastic matrix [Q]: Derivation of the formula for the in-plane elasticity matrix: For orthogonal anisotropic materials (fiber direction orthogonal to the perpendicular direction), the in-plane elastic matrix [Q] is: [Q] = [ [E 11 / (1-ν 12 ν 21 ), ν 12 E 22 / (1-ν 12 ν 21 ), 0], [ν 12 E 22 / (1-ν 12 ν 21 ), E 22 / (1-ν 12 ν 21 ), 0], [0, 0, G 12 ] ]; Wherein, the denominator (1-ν) 12 ν 21 ) is the Poisson ratio coupling term, because ν 12 ν 21 ≈0.2×0.013≈0.0026, denominator≈0.9974, can be approximated as 1, but should be retained for precise calculation.

[0064] Substitute the parameters to calculate the elasticity matrix: Substitute the material parameters of T700 into the formula: First row, first column (Q) 11 ):E 11 / (1-ν 12 ν 21 )=150GPa / (1-0.0026)≈150.39GPa; First row, second column (Q) 12 ):ν 12 E 22 / (1-ν 12 ν 21 )=0.2×10GPa / 0.9974≈2.005GPa; Second row, second column (Q)22 ):E 22 / (1-ν 12 ν 21 = 10 GPa / 0.9974 ≈ 10.026 GPa; Third row, third column (Q) 66 ): G 12 =5GPa (In-plane shear stiffness, corresponding to shear stress τ) 12 With shear strain γ 12 (relationship) All other off-diagonal elements are 0, and the final in-plane elasticity matrix [Q] is: [Q] ≈ [ [150.39, 2.005, 0], [2.005, 10.026, 0], [0, 0, 5] ]; Three-dimensional elasticity matrix extension: Considering the thickness direction (Z-axis), the three-dimensional elastic matrix [Q] 33 For a 6×6 matrix, add a stiffness term (Q) in the thickness direction. 33 =E 33 / (1-ν 13 ν 31 -ν 23 ν 32 -2ν 12 ν 23 ν 31 )), shear stiffness term (Q) 44 =G 13 Q 55 =G 23 ), etc., Q is calculated here. 33 ≈8.1 GPa, Q 44 =4GPa, Q 55 =4GPa, the remaining coupling terms are negligible due to their extremely small values, the three-dimensional elastic matrix needs to be distributed according to the layup of the laminate thickness, and constructed separately for each layer (the [Q] matrix of the 0° layer and the 90° layer only needs to be swapped). 11 With E 22 ν 12 With ν 21 That is, since the Y-axis of the local coordinate system of the 90° layer is the fiber direction.

[0065] III. Output of the parameter set for the anisotropic constitutive model: The parameter set includes a "layout parameter table + material parameter table + elasticity matrix table", as shown in the example below: "Parameter set for the anisotropic constitutive model of carbon fiber laminate (T700 / epoxy resin):" Ply parameters: Ply method [0° / 90°] 4s, single layer thickness 0.375mm, total thickness 3mm, number of layers 8, ply sequence (Z direction): 0° (0-0.375mm), 90° (0.375-0.75mm), 0° (0.75-1.125mm), 90° (1.125-1.5mm), 90° (1.5-1.875mm), 0° (1.875-2.25mm), 90° (2.25-2.625mm), 0° (2.625-3mm); Material parameters: E 11 =150GPa, E 22 =10GPa, E 33 =8GPa, G 12 =5GPa, G 13 =G 23 =4GPa, ν 12 =0.2, ν 21 ≈0.013, ν 13 =ν 23 =0.15; Elasticity matrix: 0° layer 3D elasticity matrix [Q] 0° (6×6, unit GPa): First row [150.39,2.005,1.21,0,0,0], Second row [2.005,10.026,0.81,0,0,0], Third row [1.21,0.81,8.1,0,0,0], Fourth row [0,0,0,4,0,0], Fifth row [0,0,0,0,4,0], Sixth row [0,0,0,0,0,5]; 90° layer 3D elasticity matrix [Q] 90 ° (Exchange E) 11 With E 22 The first row is [10.026, 2.005, 0.81, 0, 0, 0], the second row is [2.005, 150.39, 1.21, 0, 0, 0], and the rest is the same as layer 0°. The three-dimensional displacement field data was discretized into a 1mm×1mm×0.5mm grid. The spatial partial derivatives of the displacement components at each grid point were calculated using the second-order central difference method. The results were then substituted into the strain tensor solution formula to obtain 6 strain components, which were used as the initial strain component matrix. The three-dimensional displacement field is continuous field data (such as the displacement of 90,000 in-plane grid points output in step three), but strain calculation requires "spatial partial derivatives." Therefore, it needs to be discretized to a smaller grid first, and then differentiated using the numerical difference method to finally obtain the six components of the strain tensor (three normal strains and three shear strains). This step needs to solve three core problems: "how to discretize the continuous displacement field," "how to accurately calculate the partial derivatives," and "how to correspond the strain components to the tensor," to ensure the accuracy of the initial strain calculation.

[0066] I. Mesh discretization of the three-dimensional displacement field: Grid size setting basis: The discretization mesh needs to balance computational accuracy and efficiency. The mesh step size in the in-plane directions (X, Y) is Δx=Δy=1mm (consistent with the interpolation mesh in step three to avoid errors caused by re-interpolation), and the step size in the thickness direction (Z) is Δz=0.5mm (since the total thickness of the laminate is 3mm, Δz=0.5mm can cover 6 thickness layers, which can reflect the strain differences between layers and avoid excessive computation due to too many meshes). The final discretized mesh is 1mm(X)×1mm(Y)×0.5mm(Z), covering the X range of 0~300mm, Y range of 0~300mm, and Z range of 0~3mm of the laminate, with a total of 300×300×6=540000 mesh points.

[0067] Discrete assignment of displacement components: The three-dimensional displacement field contains three displacement components for each spatial point: u (displacement in the X direction, along the fiber / length direction), v (displacement in the Y direction, along the width direction), and w (displacement in the Z direction, along the thickness direction). During discretization, the u, v, and w values ​​of the continuous displacement field are assigned to the grid point coordinates (x_i, y_j, z_k), where x_i = i × Δx (i = 0, 1, ..., 299), y_j = j × Δy (j = 0, 1, ..., 299), and z_k = k × Δz (k = 0, 1, ..., 5). For example, the displacement components of the grid point (10mm, 20mm, 1mm) (i = 10, j = 20, k = 2) are assigned u = 0.6mm, v = 0.2mm, and w = -0.1mm (taken from the displacement field data in step three, with the w value in the thickness direction supplemented by linear interpolation).

[0068] II. Calculation of spatial partial derivatives using the second-order central difference method: Strain is the spatial rate of change of displacement. It requires calculating the first-order partial derivatives (normal strain) and mixed partial derivatives (shear strain) of u, v, and w with respect to X, Y, and Z. The second-order central difference method is more accurate than the first-order method and can effectively reduce discretization errors. Specific formulas and calculation examples are as follows: First-order partial derivatives (related to normal strain): The partial derivative of u with respect to X ( The formula describes the rate of change of displacement along the X-direction. ≈ [u (x+Δx,y,z) - 2u (x,y,z) + u (x-Δx,y,z)] / (Δx 2 The partial derivative of normal strain is the first derivative, and the formula for the first derivative of the second-order central difference is: ≈ [u (x+Δx,y,z) - u (x-Δx,y,z)] / (2Δx), where u (x+Δx,y,z) is the u value of the adjacent grid point on the right, u (x-Δx,y,z) is the u value of the adjacent grid point on the left, and 2Δx is the difference step size.

[0069] Example: For grid point (10mm, 20mm, 1mm), u = 0.6mm; for the left point (9mm, 20mm, 1mm), u = 0.5mm; for the right point (11mm, 20mm, 1mm), u = 0.7mm; Δx = 1mm. Substituting these values, we get... =(0.7-0.5) / (2×1)=0.2mm / mm=0.2 (dimensionless, unit of normal strain); Similarly, the partial derivative of v with respect to Y ( )≈ [v (x,y+Δy,z) - v (x,y-Δy,z)] / (2Δy), where in the example v=0.2mm, v=0.3mm at the upper point (10mm,21mm,1mm), v=0.1mm at the lower point (10mm,19mm,1mm), and Δy=1mm. =(0.3-0.1) / (2×1)=0.1; The partial derivative of w with respect to Z ( )≈ [w(x,y,z+Δz) - w(x,y,z-Δz)] / (2Δz), where in the example w=-0.1mm, at the upper point (10mm,20mm,1.5mm) w=-0.15mm, and at the lower point (10mm,20mm,0.5mm) w=-0.05mm, Δz=0.5mm, thus obtaining =(-0.15 - (-0.05)) / (2×0.5)=(-0.1) / 1=-0.1.

[0070] Mixed partial derivatives (related to shear strain): Shear strain is the sum of the cross rates of change of displacement in two directions, and the formula is: γ_xy= + (Shear strain in the XY plane), where ≈[u (x,y+Δy,z)-u(x,y-Δy,z)] / (2Δy), ≈[v (x+Δx,y,z)-v (x-Δx,y,z)] / (2Δx); Example: Partial derivatives of u with respect to Y: u(x,y+Δy,z)=0.65mm, u(x,y-Δy,z)=0.55mm, Therefore... =(0.65-0.55) / (2×1)=0.05; Partial derivatives of v with respect to X: v(x+Δx,y,z)=0.25mm, v(x-Δx,y,z)=0.15mm, thus... =(0.25-0.15) / (2×1)=0.05; therefore γ_xy=0.05+0.05=0.1; Similarly, γ_yz= + (YZ plane shear strain), in the example =0.02, =0.03, therefore γ_yz=0.05; γ_zx= + (ZX plane shear strain), in the example =0.01, =0.02, so γ_zx=0.03.

[0071] III. Generation of the initial strain component matrix: The strain tensor is a 6-dimensional vector, arranged in the order of "normal strain → shear strain", i.e., [ε] = [ε_xx,ε_yy, ε_zz, γ_xy, γ_yz, γ_zx]^T, where: ε_xx= (X-direction normal strain), ε_yy= (Normal strain in the Y direction), ε_zz= (Z-direction normal strain); γ_xy, γ_yz, and γ_zx are shear strains (engineering shear strains, which are twice the tensor shear strains).

[0072] The six strain components of each grid point are organized according to the correspondence between "grid point coordinates - strain components" to form an initial strain component matrix. The matrix dimension is 540000×6 (540000 grid points, each with 6 strain components). For example, the row vector of the initial strain component matrix for grid point (10mm, 20mm, 1mm) is [0.2, 0.1, -0.1, 0.1, 0.05, 0.03], where ε_xx=0.2 (2000με, micro-strain, 1με=10). -6 ), ε_yy=0.1 (1000με), ε_zz=-0.1 (-1000με, compressive strain), γ_xy=0.1 (1000με), etc.

[0073] The initial strain component matrix must be labeled with "grid point coordinates, the difference step size used in the calculation, and the source of the displacement field data" to ensure that subsequent corrections and verifications are traceable.

[0074] Based on the laminate ply direction, the local coordinate system rotation matrix corresponding to each grid point is calculated. The initial strain component matrix is ​​transformed into strain components in the local coordinate system through coordinate transformation. The calculation deviation caused by ply anisotropy is corrected, and the corrected strain component matrix is ​​output. The initial strain components are calculated based on a "global coordinate system" (X-axis along the laminate length, Y-axis along the width). However, the anisotropy of the carbon fiber laminate is based on the "fiber orientation" (local coordinate system). The local coordinate systems for different layup orientations (0°, 90°) have angles with the global coordinate system. Directly using the global strain will lead to deviations in the mechanical response calculation (e.g., the fiber orientation of the 90° layer is the global Y-axis, and the global ε_xx is actually the transverse strain of that layer, not the longitudinal strain). Therefore, a coordinate transformation is needed to convert the global strain into local strain to correct for the anisotropy deviation.

[0075] I. Determination of Local Coordinate System and Rotation Angle: Local coordinate system definition: For each ply (0° ply or 90° ply), define a local coordinate system (1-2-3): 1st axis: along the fiber direction (the 1st axis of the 0° layer is consistent with the global X-axis, and the 1st axis of the 90° layer is consistent with the global Y-axis); 2-axis: Perpendicular to the fiber direction and in the plane of the laminate (the 2-axis of the 0° layer is consistent with the global Y-axis, and the 2-axis of the 90° layer is consistent with the global X-axis); 3rd axis: along the thickness direction of the laminate (consistent with the global Z-axis, no rotation).

[0076] Setting the rotation angle θ: The rotation angle θ is the angle between the global coordinate system's X-axis and the local coordinate system's X-axis, with counterclockwise being positive. 0° layer: θ=0° (axis 1 = X-axis, axis 2 = Y-axis); 90° layer: θ = 90° (axis 1 = Y-axis, axis 2 = -X-axis); Since the thickness direction (Z-axis) of the laminate does not rotate, θ only affects the in-plane strain (ε_xx, ε_yy, γ_xy), and the thickness direction strain (ε_zz, γ_yz, γ_zx) does not need to be transformed.

[0077] II. Calculation of the rotation matrix: The coordinate transformation of strain components needs to be achieved through the "strain transformation matrix [T]", with the formula [ε]_local = [T]^T [ε]_global [T], where [T] is the coordinate rotation matrix (3×3, in-plane part). For in-plane strains (ε_xx, ε_yy, γ_xy), the in-plane part (2×2) of the transformation matrix [T] is: [T] = [ [cosθ, -sinθ], [sinθ, cosθ] ]; The corresponding strain transformation coefficient matrix [C] (used for engineering strain vector transformation, a 6×6 matrix, with only the in-plane portion being non-zero) is as follows: C 11 =cos 2 θ, C 12 =sin 2 θ, C 16 =sinθcosθ; C 21 =sin 2 θ, C 22 =cos 2 θ, C 26 =-sinθcosθ; C 61 =2sinθcosθ,C 62 =-2sinθcosθ,C 66 =cos 2 θ - sin 2 θ; The remaining components (C) 13 C 14 All of these values ​​are 0 (no rotation in the thickness direction).

[0078] Calculation of the rotation matrix for the 0° layer (θ=0°): cos0°=1, sin0°=0, substituting these values, we get: C 11 =1 2 =1, C 12 =0 2 =0, C 16 =0; C 21 =0 2 =0, C 22 =1 2 =1, C 26 =0; C 61 =0, C 62 =0, C 66 =1 2 -0 2 =1; Therefore, [C] is the identity matrix, and [ε]_local = [ε]_global, meaning that the local strain of the 0° layer is consistent with the global strain and does not require correction.

[0079] Calculation of the 90° layer rotation matrix (θ=90°): cos90°=0, sin90°=1, substituting them, we get: C 11 =0 2 =0, C 12 =1 2 =1, C 16 =1×0=0; C 21 =1 2 =1, C 22 =0 2 =0, C 26 =-1×0=0; C 61 =2×1×0=0,C 62 =-2×1×0=0,C 66 =0 2 -1 2 =-1; Thickness component C 33 =C 44 =C 55 =1, and the rest are 0, therefore the strain transformation formula for the 90° layer is: ε 11 (Local longitudinal normal strain) = C 11 ε_xx + C 12 ε_yy + C 16 γ_xy=0×ε_xx +1×ε_yy +0×γ_xy=ε_yy; ε 22(Local transverse normal strain) = C 21 ε_xx + C 22 ε_yy + C 26 γ_xy=1×ε_xx +0×ε_yy +0×γ_xy=ε_xx; γ 12 (Local in-plane shear strain) = C 61 ε_xx + C 62 ε_yy + C 66 γ_xy=0+0+(-1)×γ_xy=-γ_xy; Thickness strain ε 33 =ε_zz,γ 13 =γ_yz,γ 23 =γ_zx (no change).

[0080] III. Calculation and Deviation Correction Example of Local Strain Components: Taking a grid point (10mm, 20mm, 0.5mm) in a 90° layer (Z=0.5mm, belonging to the 2nd layer, 90° layup) as an example, its global initial strain components are [ε_xx=0.1, ε_yy=0.2, ε_zz=-0.1, γ_xy=0.1, γ_yz=0.05, γ_zx=0.03], which can be converted into local strain: ε 11 =ε_yy=0.2 (local longitudinal strain, along the fiber direction, original global strain in the Y direction); ε 22 =ε_xx=0.1 (Local transverse strain, perpendicular to the fiber direction, original global strain in the X direction); γ 12 =-γ_xy=-0.1 (Local shear strain, sign reversed due to 90° rotation of the coordinate system); The strain in the thickness direction remains constant: ε 33 =-0.1, γ 13 =0.05, γ 23 =0.03; The corrected local strain components are [0.2, 0.1, -0.1, -0.1, 0.05, 0.03].

[0081] If not corrected, directly using the global strain ε_xx=0.1 as the longitudinal strain of the 90° layer (which should actually be the transverse strain) will lead to deviations in subsequent stress calculations (longitudinal stress σ). 11 =Q 11 ε_xx+Q 12ε_yy≈150.39×0.1+2.005×0.2≈15.44GPa, while the corrected σ 11 =Q 11 ε 11 +Q 12 ε 22 ≈150.39×0.2+2.005×0.1≈30.28GPa, with a deviation of 50%, therefore the correction step is crucial.

[0082] IV. Output of the corrected strain component matrix: The global strain of all grid points is converted into local strain according to the angle θ of their respective ply, forming a corrected strain component matrix. The matrix dimension remains 540000×6, but the strain component of each grid point now corresponds to the local coordinate system (fiber direction) of its ply. For example, the strain components of the grid points in the 0° layer are consistent with the initial matrix, the strain components of the grid points in the 90° layer have been converted to θ=90°, and the strain components of the grid points in the thickness direction (such as Z=1mm, 4th layer 90°) are also converted according to the corresponding θ.

[0083] When outputting, you need to mark "the ply assignment (0° / 90°), rotation angle θ, and transformation matrix [T] of each grid point" to ensure that the elastic matrix of the corresponding ply can be called in subsequent mechanical verification.

[0084] The corrected strain component matrix is ​​mechanically consistent with the parameter set of the anisotropic constitutive model. Strain values ​​exceeding the elastic limit of the material are removed, and the removed region is completed by linear interpolation to generate preliminary three-dimensional strain field inversion results.

[0085] The corrected strain components must satisfy the "material mechanical property constraints," meaning the strain value cannot exceed the material's elastic limit (if it does, the material enters the plastic stage, and the constitutive model is no longer applicable). Simultaneously, it must be consistent with the stiffness characteristics reflected by the elastic matrix (e.g., longitudinal strain should be less than transverse strain because longitudinal stiffness is greater, resulting in smaller deformation under the same load). This step involves eliminating outliers through mechanical consistency verification and then supplementing the data to ensure the physical rationality of the initial strain field.

[0086] I. Determination of the elastic limit strain of a material: Based on the standard tensile test of T700 carbon fiber / epoxy resin composites, the elastic limit strain (proportional limit strain, beyond which stress-strain is no longer linear) is: Longitudinal elastic limit strain ε 11 ^lim=1.5%=15000με (along the fiber direction, high stiffness and high elastic limit); Transverse elastic limit strain ε 22^lim=0.5%=5000με (perpendicular to the fiber direction, low stiffness, low elastic limit); In-plane shear elastic limit strain γ 12 ^lim=1.0%=10000με; Elastic limit strain ε in the thickness direction 33 ^lim=0.8%=8000με, shear elastic limit strain γ 13 ^lim=γ 23 ^lim=0.9%=9000με.

[0087] These limit values ​​serve as thresholds for verification; exceeding them indicates abnormal strain (which may be caused by displacement field interpolation errors or differential calculation errors).

[0088] II. Execution of Mechanical Consistency Verification: Single-component limit check: Traverse the corrected strain component matrix, and for each of the six strain components at each grid point, determine whether the corresponding elastic limit is exceeded: If ε 11 >15000με or ε 11 If the value is less than -15000με (the compression limit and the tensile limit are the same, so this is a simplified treatment), it is considered abnormal. If ε 22 >5000με or ε 22 <-5000με is considered abnormal; The same applies to shear strain, such as γ. 12 >10000με or γ 12 <-10000με is considered abnormal.

[0089] Example: Local strain ε at a grid point in a 90° layer 22 =6000με (exceeding 5000με) is identified as an abnormal strain component and marked as "to be removed".

[0090] Stiffness consistency verification: By combining the elasticity matrix of the anisotropic constitutive model, the "matching between strain components and stiffness" is verified—the direction with larger stiffness (such as longitudinal E). 11 =150GPa), the strain should be less than in the direction of lower stiffness (such as the transverse direction E). 22 =10GPa), that is, for the 0° layer, |ε 11 |≤|ε 22 |;For the 90° layer, |ε 11 |(original ε_yy)≤|ε 22 | (original ε_xx), if this rule is violated, it is judged as abnormal.

[0091] Example: Local strain ε at a grid point in layer 0° 11 =8000με,ε 22 =5000με,|ε 11 |>|ε 22 |(Longitudinal strain is greater than transverse strain, and E) 11 >E 22 If a statement is found to be contradictory, it is considered abnormal and must be removed even if it does not exceed the limit.

[0092] III. Removal of Abnormal Strain Components and Data Completion: Removal rules: For strain components marked as anomalous, they are directly removed from the matrix, and the location is marked with "NaN" (Not Numerical), for example, the ε of a certain grid point. 22 =6000με (abnormal), this component is removed and marked as NaN, while the other components are retained.

[0093] Linear interpolation completion: The strain values ​​at NaN locations are completed using "3×3×3 neighborhood linear interpolation". Specifically, for an anomalous grid point (x_i, y_j, z_k), 3×3×3 normal grid points (26 neighbors) are selected around it, and the completed value is calculated using a distance-weighted average. The formula is: ε_comp = Σ(ε_nei × w_nei) / Σw_nei. Where ε_nei is the normal strain value of the neighboring grid point, and w_nei = 1 / d_nei (d_nei is the Euclidean distance between the anomalous point and its neighbors; the closer the distance, the greater the weight).

[0094] Example: ε of the outlier (x=10, y=20, z=1) 22 =NaN, ε of the three surrounding neighboring points 22 Given 4500με (d=1mm, w=1), 4800με (d=√2mm, w≈0.707), and 4600με (d=√3mm, w≈0.577), then ε_comp=(4500×1 + 4800×0.707 + 4600×0.577) / (1+0.707+0.577)≈(4500+3393.6+2654.2) / 2.284≈10547.8 / 2.284≈4618με (within the 5000με limit, which is reasonable).

[0095] IV. Generation of preliminary three-dimensional strain field inversion results: The completed strain component matrix is ​​the preliminary three-dimensional strain field inversion result, which includes "spatial distribution + temporal variation" - according to the X (0~300mm), Y (0~300mm), Z (0~3mm) spatial coordinates of the laminate, and the temporal sequence of the impact event (t=-10ms~t=20ms), with each spatiotemporal point corresponding to 6 local strain components.

[0096] Example result description: "At t=1ms (impact load rise stage), the local longitudinal strain ε at the 0° layer grid point in the impact center region (X=140~160mm, Y=140~160mm) is..." 11 =8000~12000με (not exceeding 15000με), transverse strain ε 22 =3000~5000με (close to the 5000με limit), in-plane shear strain γ 12 =4000~6000με; ε of the region corresponding to the 90° layer 11 =7000~11000με, ε 22 =3500~4800με; the strain in the edge region (X=0~50mm, Y=0~50mm) is generally less than 1000με, which conforms to the principle of stress concentration at the impact center and small deformation at the edge. The results need to be stored in a structured form of "spatiotemporal grid - strain components" (such as HDF5 format), including strain field data for each time series and verification logs (number of outliers removed, completion method), to provide a foundation for the generation of the final inversion results.

[0097] S204, the preliminary three-dimensional strain field inversion results are coupled with the impact load parameters and the laminate material property parameters. Abnormal strain components are eliminated through strain-load correlation verification to generate the final three-dimensional strain field inversion result of carbon fiber laminate that is consistent with the dynamic response of the impact event.

[0098] Specifically, impact load parameters and laminate material property parameters can be organized to build a parameter association database. Through normalization, parameters of different dimensions can be mapped to the same dimension, and a set of coupled operation parameters can be output. This step is fundamental to the coupled computation, and its core purpose is to address the problem that "the dimensions of impact loads (force, time) and material properties (stiffness, ultimate strain) differ greatly, making direct correlation impossible." It is necessary to first clarify the specific content and source of both types of parameters, and then eliminate the influence of dimensions through normalization to form a structured parameter library, providing a unified input standard for subsequent strain-load correlation.

[0099] I. Parameter preparation and database construction: Extraction and definition of impact load parameters: The impact load parameters are derived from the data recorded by the force sensor in the above steps (piezoelectric force sensor, sampling frequency 10kHz). Key indicators reflecting the dynamic characteristics of the impact need to be extracted, and each parameter must be labeled with its physical meaning and the basis for its value. Maximum impact load F_max: The maximum force during the impact process, reflecting the impact intensity. The example value is 5.2kN (drop hammer mass 5kg, height 100mm, impact energy ≈5J, matching the 3mm thickness of the laminate). Total impact time t_total: The total time from when the load rises to 0 until it falls back to 0. Example value is 4.8ms (the complete time from when the impact head contacts and separates from the laminate). Load rise time t_rise: The time it takes for the load to rise from 0 to F_max, reflecting the impact loading rate. The example value is 1.2ms (rapid impact, high loading rate, and easy local stress concentration in the laminate). Load fall time t_fall: The time from F_max to 0, the example value is 3.1ms (unloading process, laminate elastic recovery stage); Peak load time t_peak: The time when the load reaches F_max, the example value is 1.7ms (the critical moment of the impact dynamic response, corresponding to the peak strain of the laminate).

[0100] A summary of the material properties of laminates: The material property parameters are derived from the anisotropic constitutive model constructed in the above steps. Parameters that are strongly correlated with strain verification need to be selected to ensure that subsequent coupling calculations can reflect the "constraint of material properties on the strain-load relationship". Longitudinal elastic modulus E 11 Stiffness along the fiber direction determines the sensitivity of the longitudinal strain response to load; the example value is 150 GPa (standard value for T700 carbon fiber / epoxy resin). Transverse elastic modulus E 22 Stiffness perpendicular to the fiber direction, example value is 10 GPa; In-plane shear modulus G 12 XY plane shear stiffness, example value is 5GPa; Longitudinal elastic limit strain ε 11 ^lim: The proportional limit along the fiber direction. If exceeded, the material enters the plastic stage. Example value is 15000με (1.5%). Transverse elastic limit strain ε 22 ^lim: The proportional limit perpendicular to the fiber direction, with an example value of 5000με (0.5%).

[0101] Parameter association database construction: Organized logically according to "parameter category - physical meaning - value - dimension - source", a textual database is formed, for example: Parameters are associated with the database: Impact load parameter class: F_max: 5.2kN (dimension: force), source: force sensor record (t=1.7ms); t_total: 4.8ms (unit: time), source: total duration of the load curve from 0 to 0; t_rise: 1.2ms (unit: time), source: duration of load from 0 to 5.2kN; t_peak: 1.7ms (unit: time), source: the moment corresponding to the peak value of the load curve; Material property parameters: E 11 150 GPa (dimension: stress), source: T700 / Epoxy Resin Standard Tensile Test; ε 11 ^lim: 15000με (dimensionless), source: material proportional limit test; ε 22 ^lim: 5000με (dimensionless), source: same as above. II. Normalization Process (Min-Max Standardization): The dimensions of different parameters differ greatly (e.g., F_max = 5.2kN vs ε). 11 Since ^lim=15000με), all parameters need to be mapped to the [0,1] interval using Min-Max normalization, as shown in the formula: x_norm = (x - x_min) / (x_max - x_min): Where x is the original value of the parameter, x_min is the theoretical minimum value of the parameter (usually 0, such as load or time) or the minimum threshold of the material / equipment, and x_max is the reasonable maximum value of the parameter (the upper limit of the equipment range or common material limit), ensuring that the normalized value has physical meaning (non-negative and not exceeding 1).

[0102] Example of normalization of impact load parameters: F_max normalization: x_max is taken as the maximum load range of the impact test bench, 10kN (the upper limit of the equipment capacity), x_min=0kN, substituting them, we get: F_max_norm = (5.2 - 0) / (10 - 0) = 0.52; t_total normalization: x_max takes the upper limit of common impact time as 10ms (laminate impact is usually ≤10ms), x_min=0ms, substituting them gives: t_total_norm = (4.8 - 0) / (10 - 0) = 0.48; Normalized t_peak: x_max = 10ms, x_min = 0ms, substituting, we get: t_peak_norm = (1.7 - 0) / (10 - 0) = 0.17.

[0103] Example of normalization of material property parameters: E 11 Normalization: x_max is taken as the upper limit of the common longitudinal modulus of carbon fiber composites, 200 GPa, and x_min = 50 GPa (lower limit of low modulus composites). Substituting these values, we get: E 11 _norm = (150 - 50) / (200 - 50) = 100 / 150 ≈ 0.67; ε 11 Normalization using ^lim: x_max = 20000με (a common upper limit for the elastic limit of materials), x_min = 0με, substituting these values ​​gives: ε 11 ^lim_norm = (15000 - 0) / (20000 - 0) = 0.75; ε 22 Normalization using ^lim: x_max = 10000με, x_min = 0με, substituting these values ​​gives: ε 22 ^lim_norm = (5000 - 0) / (10000 - 0) = 0.5.

[0104] Output of the coupled operation parameter set: All normalized parameters are categorized and organized according to "impact load - material properties" to form a parameter set that can be directly used for coupled calculations, for example: "Coupled operation parameter set (normalized, all dimensions are [0,1]):" Impact load normalization parameters: F_max_norm=0.52, t_total_norm=0.48, t_rise_norm=0.24 (x_max=5ms), t_peak_norm=0.17; Material property normalization parameter: E 11 _norm=0.67, E 22 _norm=0.1 (E 22 =10GPa, x_max=100GPa), G 12 _norm=0.17 (G 12 =5GPa, x_max=30GPa), ε 11 ^lim_norm=0.75, ε 22 ^lim_norm=0.5. The preliminary three-dimensional strain field inversion results are coupled with the coupling operation parameter set to establish a strain-load time domain correlation model. The cross-correlation coefficient R between each strain component and the load curve is calculated. R ≥ 0.8 is set as an effective correlation, and the strain-load correlation matrix is ​​output. The core of this step is to verify the dynamic consistency between the strain field and the impact load through "time-domain correlation"—during the impact, the rise / fall trend of strain should be synchronized with the load curve (e.g., the load peak corresponds to the strain peak), and the cross-correlation coefficient R is the key indicator for quantifying this synchronicity. It is necessary to clarify the construction logic of the correlation model, the calculation method of R, and the criteria for determining effective correlation, ensuring that only strain components matching the dynamic response of the load are retained.

[0105] I. Establishment of the strain-load time-domain correlation model: The time-domain correlation model uses "time" as a link to establish a point-to-point correlation between the "strain time series curve" and the "impact load time series curve" of each grid point in the initial strain field. The model assumes that during the elastic deformation stage, the time-domain change of strain should be linearly correlated with the time-domain change of load (which conforms to the dynamic extension of Hooke's Law).

[0106] Extraction of time series data: Strain time series ε(t): Extracting the variation of a strain component of a specific grid point with time from the preliminary three-dimensional strain field inversion results. The time range is consistent with the load record (t=-10ms~t=20ms, 1ms interval, a total of 31 time points). Example: Grid point (150mm, 150mm, 1mm) in the impact center area (0° layer, longitudinal strain ε 11 Timing of ε from t=-10ms to t=0ms 11=0 (no deformation before impact); 8000με at t=1ms; 12000με at t=1.7ms (peak value, consistent with the peak load time); 8000με at t=3ms; 0 at t=5ms (elastic recovery after impact); ε at t=6ms~t=20ms 11 =0 (no residual strain).

[0107] Load timing F(t): Extracted from force sensor records, with 31 time points at 1ms intervals. Example: F=0 at t=-10ms~t=0ms; 3kN at t=1ms; 5.2kN (peak) at t=1.7ms; 3kN at t=3ms; 0 at t=5ms; F=0 at t=6ms~t=20ms.

[0108] Mathematical expression of the association model: The model is constructed using linear regression and residual analysis, assuming ε(t) = a×F(t) + b + e(t), where: 'a' is a proportionality coefficient (reflecting the sensitivity of the load to the strain, inversely proportional to the material stiffness, E) 11 The larger the value, the smaller the value of a). b is a constant term (usually 0, since both strain and load are 0 before impact); e(t) is the residual (model error; the smaller the value, the higher the correlation).

[0109] In the example, linear regression yields a = 2307.7 με / kN (12000 με / 5.2kN ≈ 2307.7), b = 0, and the maximum value of the residual e(t) is 500 με (≤ 5% × 12000 με, the error is acceptable), indicating that the strain component has a good linear correlation with the load.

[0110] II. Calculation of the cross-correlation coefficient R (Pearson correlation coefficient): The cross-correlation coefficient R is used to quantify the degree of linear correlation between ε(t) and F(t), with a value ranging from -1 to 1. The closer R is to 1, the stronger the correlation; when R < 0.8, the correlation is considered weak (strain changes may be caused by noise or calculation errors, not by load). The calculation formula is: R = COV(ε(t), F(t)) / √[VAR(ε(t)) × VAR(F(t))]; in: COV (ε(t), F (t)) is the covariance of ε(t) and F (t), reflecting the consistency of their overall changing trends; VAR(ε(t)) and VAR(F(t)) are the variances of ε(t) and F(t), respectively, reflecting their respective degrees of dispersion.

[0111] Example of calculation steps (based on the time series data above): Selecting valid data from t=0ms to t=5ms (a total of 6 time points: t0=0, t1=1, t2=1.7, t3=3, t4=5), the calculations are as follows: Calculate the mean: μ_ε = (0 + 8000 + 12000 + 8000 + 0) / 5 = 5600με; μ_F = (0 + 3 + 5.2 + 3 + 0) / 5 = 2.24kN; Calculate the covariance (COV): COV = [Σ(ε_i - μ_ε)(F_i - μ_F)] / (n - 1) = [(0-5600)(0-2.24) + (8000-5600)(3-2.24) + (12000-5600)(5.2-2.24) + (8000-5600)(3-2.24) + (0-5600)(0-2.24)] / (5-1) = [12544 + 1824 + 18944 + 1824 + 12544] / 4 = 47680 / 4 = 11920 (με·kN); Calculate the variance VAR: VAR(ε) = [Σ(ε_i - μ_ε) 2 ] / (n - 1) = [(0-5600) 2 + (8000-5600) 2 +(12000-5600) 2 + (8000-5600) 2 + (0-5600) 2 ] / 4 = (31360000 + 5760000 + 40960000 + 5760000 + 31360000) / 4 = 115200000 / 4 = 28800000 (με 2 ); VAR(F) = [Σ(F_i - μ_F) 2 ] / (n - 1) = [(0-2.24)2 + (3-2.24) 2 + (5.2-2.24) 2 + (3-2.24) 2 + (0-2.24) 2 ] / 4 = (5.0176 + 0.5776 + 8.7616 + 0.5776 +5.0176) / 4 = 19.952 / 4 = 4.988 (kN 2 ); Calculate R: R = 11920 / √(28800000 × 4.988) = 11920 / √143654400 ≈ 11920 / 11985.6 ≈ 0.994 (close to 1, extremely high correlation).

[0112] Examples of R calculation results for different strain components: Calculate R for other strain components at the same grid point: ε 22 (Transverse strain): Time-series peak 4800 με (t=1.7 ms), R=0.98 (good correlation); γ 12 (In-plane shear strain): Time-series peak value 6000με (t=1.7ms), R=0.92 (good correlation); ε 33 (Thickness direction strain): Time-series peak value -1200με (compression, t=1.7ms), R=0.75 (<0.8, weak correlation, because the thickness direction deformation is constrained by the laminate and has low correlation with in-plane load). γ 13 (ZX plane shear strain): time-series peak 800με, R=0.78 (<0.8, weak correlation).

[0113] III. Output of the strain-load correlation matrix: The R values ​​of the six strain components of all grid points are organized according to "grid point coordinates - strain component - R value - associated state" to form a strain-load correlation matrix (textualized form, no table). Example: Strain-load correlation matrix (t=-10ms~t=20ms, mesh size 1mm×1mm×0.5mm): Impact center zone (X=140~160mm, Y=140~160mm, Z=0~3mm): Grid point (150,150,1) (0° layer): ε11 _R=0.994 (valid), ε 22 _R=0.98 (valid), ε 33 _R=0.75 (invalid), γ 12 _R=0.92 (valid), γ 13 _R=0.78 (invalid), γ 23 _R=0.76 (invalid); Grid point (150, 150, 1.5) (90° layer): ε 11 _R=0.97 (Valid, ε) 11 Peak value 11000με), ε 22 _R=0.96 (valid), ε 33 _R=0.77 (invalid), γ 12 _R=0.91 (valid), γ 13 _R=0.79 (invalid), γ 23 _R=0.77 (invalid); Edge zone (X=20~30mm, Y=20~30mm, Z=0~3mm): Grid point (25,25,1) (0° layer): ε 11 _R=0.6 (invalid, ε) 11 Peak value 6000με, exceeding the pattern of small edge strain), ε 22 _R=0.58 (invalid), ε 33 _R=0.65 (invalid), γ 12 _R=0.62 (invalid), γ 13 _R=0.61 (invalid), γ 23 _R=0.63 (invalid); Association state definition: R≥0.8 is "valid", R<0.8 is "invalid".

[0114] Based on the strain-load correlation matrix, abnormal strain components with cross-correlation coefficient R < 0.8 or strain values ​​exceeding the preset range are marked. The rationality of the abnormal marking is verified by combining the stress concentration law in the impact area, and an abnormal strain marking map is output. This step requires double screening of abnormal strains: first, "dynamic correlation anomalies (R<0.8)" and second, "static numerical anomalies (exceeding the preset range)". At the same time, the rationality of the anomalies is verified by combining the basic laws of impact mechanics (central stress concentration and small edge stress) to avoid misjudging "reasonable weak correlation components" (such as thickness direction strain) as anomalies.

[0115] I. Marking rules for abnormal strain components: Type I anomalies: Dynamic correlation anomalies with R < 0.8: Directly label all components in the strain-load correlation matrix where R < 0.8, example: ε of the grid points in the central area 33 γ 13 γ 23 (R=0.75~0.79); All strain components (R=0.58~0.65) at grid points in the edge region.

[0116] Category 2 anomalies: Static numerical anomalies where strain values ​​exceed the preset range. The preset range is based on the "stress concentration law in the impact zone + material elastic limit" setting, and the range varies for different regions: Impact center zone (X=140~160mm, Y=140~160mm): ε 11 ≤15000με (elastic limit of the material), ε 22 ≤5000με,γ 12 ≤10000με; Transition zone (X=120~140mm, Y=120~140mm; X=160~180mm, Y=160~180mm): ε 11 ≤8000με, ε 22 ≤3000με, γ 12 ≤6000με; Edge zone (X<120mm or X>180mm; Y<120mm or Y>180mm): ε 11 ≤2000με, ε 22 ≤1500με, γ 12 ≤2000με; Thickness direction and shear strain (ε) 33 γ 13 γ 23 There are no strict numerical limits, but they need to match the stress level of the region (the central area can be slightly larger, and the edge area needs to be extremely small).

[0117] In the example, the ε of the edge region grid point (25,25,1) 11 =6000με (exceeding the ≤2000με range in the edge region), marked as a numerical anomaly; ε at grid point (150,150,1) in the central region. 11 =12000με (≤15000με), no numerical anomalies.

[0118] II. Verification of the rationality of abnormal markings (in conjunction with the law of impact stress concentration): Under impact loading, the stress-strain distribution of the laminate follows a pattern of "concentration at the center and attenuation at the edges"—the strain is greatest at the impact center due to direct load, and decreases exponentially with distance from the center. Anomaly markers are verified based on this pattern. Verification of the rationality of weakly correlated components in the central region: ε in the central area 33 (R=0.75), γ 13 (R=0.78) Although R<0.8, the strain in the thickness direction (ε) 33 Constrained by the interlaminar structure of the laminate, the deformation is small (-1200με), and the load mainly causes in-plane deformation (ε). 11 ε 22 The correlation between the thickness direction and the load should be relatively weak, so it is marked as "reasonably weak correlation" and is not considered an anomaly to be removed.

[0119] Verification of the rationality of anomalous components in the edge region: ε of the edge region grid point (25,25,1) 11 =6000με (exceeding ≤2000με) and R=0.6 (<0.8), which contradicts the rule of "small edge strain" - the edge area is far from the impact center, and the load has been greatly attenuated by this point, so it is impossible to generate a large strain of 6000με. Therefore, it is marked as "unreasonable anomaly" and needs to be removed later.

[0120] Verification of the reasonableness of the transitional differentiation quantity: ε of the transition zone grid point (130,130,1) 11 =7500με (≤8000με), R=0.82 (≥0.8), which conforms to the rule that "the strain in the transition zone is between the center and the edge", and is marked as "normal and effective".

[0121] III. Output of Abnormal Strain Marker Diagram: Anomaly strain marker maps are based on "spatial location + strain components + anomaly type + rationality," using text to describe the distribution of anomalies in different areas. Example: Abnormal strain marking diagram (corresponding to the peak load time at t=1.7ms, laminate X=0~300mm, Y=0~300mm, Z=0~3mm): Impact center zone (140~160mm, 140~160mm): ε of all grid points 33 γ 13 γ 23 Marked as 'Weak correlation anomaly (R=0.75~0.79), Reasonableness: Reasonable (weak correlation between thickness / shear direction and in-plane load)'; No numerical outliers were found. Transition zone (120~140mm, 120~140mm; 160~180mm, 160~180mm): γ of a few grid points 13 γ 23 Marked as 'weak association anomaly (R=0.78~0.79), reasonableness: reasonable'; No numerical outliers were found. Edge zone (X<120mm or X>180mm; Y<120mm or Y>180mm): ε of all grid points 11 ε 22 γ 12 Marked as 'numerical + correlation dual anomaly (ε)' 11 =4000~6000με>2000με, R=0.58~0.65<0.8), Reasonableness: Unreasonable; ε of some grid points 33 γ 13 γ 23 Marked as 'Association Anomaly (R=0.61~0.65<0.8), Reasonableness: Unreasonable (no load transfer in the edge area, strain should be close to 0)'; Anomaly type definitions: 'Weak association anomaly' (only R < 0.8, reasonable value), 'Numerical anomaly' (only out of range, R ≥ 0.8), 'Double anomaly' (R < 0.8 and out of range).

[0122] Abnormal components in the abnormal strain marker map are removed, and weighted interpolation is used to complete the data in the removed area. Finally, the three-dimensional strain field of the carbon fiber laminate that is consistent with the dynamic response of the impact event is generated.

[0123] This step is the final stage of the inversion process. It is necessary to remove "unreasonable abnormal components" and fill in the data gaps through interpolation to ensure that the final strain field meets the requirements of "dynamic synchronization with load, static conformity to stress law, and spatial continuity without gaps", so as to provide reliable data for subsequent impact damage analysis of laminated plates.

[0124] I. Rules for eliminating abnormal strain components: Only components marked "Reasonable: Unreasonable" in the abnormal strain marker diagram are removed, while components marked "Reasonable Weak Correlation" (such as ε in the central region) are retained. 33 To avoid excessive filtering that could lead to the loss of useful information: Removed components: all strain components in the edge region (double anomalies), and components in the transition region with R < 0.8 and values ​​outside the range; Rejection operation: Mark the strain value of the rejected object as "NaN" (non-numerical), indicating that the data at that location is invalid and needs to be completed later; Reserved object: ε in the central area 11 ε 22 γ 12 (Effective association), ε in the central region 33 γ 13 γ 23 (Reasonable weak association), effective association components in the transition zone.

[0125] Example: ε of the edge grid point (25,25,1) 11 ε 22 γ 12 ε 33 γ 13 γ 23 All are labeled as NaN; the ε of the grid point (150,150,1) in the central region. 33 γ 13 γ 23 Retain the original values ​​(-1200με, 800με, 750με).

[0126] II. Weighted Interpolation Completion (Inverse Distance Weighted Interpolation, IDW): For regions marked as NaN, inverse distance weighted interpolation is used for completion—using the strain values ​​of surrounding normal grid points, the completion value is calculated according to the principle of "the closer the distance, the greater the weight", to ensure that the strain distribution after completion is continuous and conforms to the stress attenuation law.

[0127] Interpolation parameter settings: Search neighborhood: For each NaN grid point, search for normal grid points within a 50mm radius around it (to ensure sufficient data support and avoid interpolation bias). Weight formula: weight w_i = 1 / (d_i^p), where d_i is the Euclidean distance (unit: mm) between the NaN point and the i-th normal point, p=2 (weight decay coefficient, the larger p is, the higher the weight of the nearest point, and the closer the interpolation is to the local trend); Normalized weights: w_i' = w_i / Σw_i (to ensure that all weights sum to 1 and avoid interpolation values ​​from exceeding a reasonable range).

[0128] Completion example (ε of edge grid point (25,25,1)) 11 ): Search for normal points within the neighborhood (the boundary between the transition zone and the edge zone): Point A (120, 120, 1): ε 11=2000με,d_A=√((120-25) 2 +(120-25) 2 )=√(95 2 +95 2 =√18050≈134.3mm; Point B (120, 100, 1): ε 11 =1800με,d_B=√((120-25) 2 +(100-25) 2 )=√(95 2 +75 2 ) = √14800 ≈ 121.7 mm; Point C (100, 120, 1): ε 11 =1700με,d_C=√((100-25) 2 +(120-25) 2 )=√(75 2 +95 2 ) = √14800 ≈ 121.7 mm; Calculate the weights: w_A=1 / (134.3 2 )≈1 / 18036≈0.000055; w_B=1 / (121.7 2 )≈1 / 14811≈0.0000675; w_C=1 / (121.7 2 )≈0.0000675; Σw=0.000055+0.0000675+0.0000675≈0.00019; w_A'=0.000055 / 0.00019≈0.289; w_B'=0.0000675 / 0.00019≈0.355; w_C'=0.0000675 / 0.00019≈0.355; Calculation of complete value: ε 11 The completion is calculated as follows: 2000×0.289 + 1800×0.355 + 1700×0.355 ≈ 578 + 639 + 603.5 ≈ 1820.5με (≈1820με, which is consistent with the pattern within the ≤2000με range of the edge region).

[0129] Full area completion effect: After completion, the strain values ​​in the edge region are all in the range of 0~2000με, the transition region is in the range of 2000~8000με, and the central region is in the range of 8000~12000με. The strain decreases continuously from the center to the edge, with no obvious data gaps.

[0130] III. Generation of the final three-dimensional strain field inversion results: After integrating and completing the strain data, the final results are output according to three dimensions: spatial distribution, temporal consistency, and material consistency, ensuring that the results fully match the dynamic response, material properties, and mechanical laws of the impact event. Final inversion results of the three-dimensional strain field of carbon fiber laminate (T700 / epoxy resin, layup [0° / 90°] 4s, impact energy 5J): Spatial distribution characteristics (t=1.7ms, peak load time): Impact center zone (140~160mm, 140~160mm, Z=0~3mm): 0° layer: ε 11 =10000~12000με (R=0.98~0.994), ε 22 =4500~5000με (R=0.96~0.98), γ 12 =5500~6000με (R=0.91~0.92); 90° layer: ε 11 =9000~11000με (R=0.95~0.97),ε 22 =4000~4800με (R=0.94~0.96), γ 12 =5000~5500με (R=0.90~0.91); Thickness direction: ε 33 =-1000~-1200με (R=0.75~0.77, weak correlation), γ 13 =750~800με (R=0.78~0.79), γ 23 =700~750με (R=0.76~0.77); Transition zone (120~140mm, 120~140mm; 160~180mm, 160~180mm): 0° layer: ε 11 =5000~8000με (R=0.82~0.90), ε 22 =2500~3000με (R=0.81~0.88); Edge zone (X<120mm or X>180mm; Y<120mm or Y>180mm): All layers: ε 11 =500~1820με (R=0.80~0.85, after completion), ε 22 =400~1500με (R=0.80~0.84); Timing consistency verification: The peak time for all strain components is t=1.7ms (consistent with the peak time of the load). The strain rise time (0~1.7ms) is matched with the load rise time (1.2ms), and the fall time (1.7~5ms) is matched with the load fall time (3.1ms). After impact (t>5ms), the strain returns to 0 with no residual strain, which is consistent with the elastic deformation characteristics. Material consistency verification: All strain values ​​did not exceed the material's elastic limit (ε). 11 ≤15000με, ε 22 ≤5000με); Longitudinal strain (ε) 11 ) greater than the transverse strain (ε) 22 ), conforming to E 11 >E 22 Stiffness characteristics; Applications of the results: This can be used to analyze stress concentration areas (central region) and interlaminar deformation patterns (ε in the thickness direction) under impact on laminates. 33 This provides strain input for impact damage prediction.

[0131] Another embodiment of the present invention provides a three-dimensional strain field inversion system for laminate impact events based on binocular vision, see [link to documentation]. Figure 3 The system may include: The acquisition module 301 is used to set dynamic tracking markers on the surface of the impacted carbon fiber laminate, and to use a high-speed binocular camera to synchronously acquire image data of the entire impact process. At the same time, it records the magnitude of the impact load and the duration of the impact. The binocular camera intrinsic parameter calibration algorithm unifies the image coordinate system and generates a binocular visual image sequence with impact time information. Extraction module 302 is used to extract the dynamic features of the marker points and the fiber texture of the laminate based on the binocular vision image sequence using a hierarchical feature matching algorithm, calculate the three-dimensional spatial coordinates of the feature points in each time sequence through stereo matching, and solve the three-dimensional displacement vector of the feature points by combining the difference between adjacent time sequence coordinates to generate three-dimensional displacement field data of the impact area of ​​carbon fiber laminate. Module 304 is used to call the anisotropic constitutive model of carbon fiber laminate, substitute the three-dimensional displacement field data into the strain tensor solution formula, convert it into local strain components through spatial derivative calculation, and combine the laminate layup direction to correct the strain calculation deviation and generate preliminary three-dimensional strain field inversion results. The generation module 305 is used to couple the preliminary three-dimensional strain field inversion results with the impact load parameters and the laminate material property parameters, and eliminate abnormal strain components through strain-load correlation verification to generate the final three-dimensional strain field inversion result of carbon fiber laminate that is consistent with the dynamic response of the impact event.

[0132] This invention also provides a storage medium storing a computer program, wherein the computer program is configured to execute the steps in any of the above method embodiments when running.

[0133] This invention also provides an electronic device, including a memory and a processor, wherein the memory stores a computer program, and the processor is configured to run the computer program to perform the steps in any of the above method embodiments.

[0134] Specifically, the aforementioned electronic device may further include a transmission device and an input / output device, wherein the transmission device is connected to the aforementioned processor, and the input / output device is connected to the aforementioned processor.

[0135] The above description, based on the embodiments shown in the figures, details the structure, features, and effects of the present invention. The above description is only a preferred embodiment of the present invention, but the present invention is not limited to the scope of implementation shown in the figures. Any changes made in accordance with the concept of the present invention, or equivalent embodiments modified to have equivalent changes, that do not exceed the spirit covered by the specification and figures, should be within the protection scope of the present invention.

Claims

1. A method for inverting the three-dimensional strain field of a laminated plate impact event based on binocular vision, characterized in that, The method includes: Dynamic tracking markers were placed on the surface of the impacted carbon fiber laminate. Image data of the entire impact process was acquired synchronously using a high-speed binocular camera. The magnitude of the impact load and the duration of the impact were recorded. The image coordinate system was unified through the binocular camera intrinsic parameter calibration algorithm to generate a binocular visual image sequence with impact time information. Based on the binocular vision image sequence, a hierarchical feature matching algorithm is used to extract the dynamic features of the marker points and the fiber texture of the laminate. The three-dimensional spatial coordinates of the feature points at each time sequence are calculated by stereo matching. The three-dimensional displacement vector of the feature points is solved by combining the difference between adjacent time sequence coordinates, and the three-dimensional displacement field data of the impact area of ​​the carbon fiber laminate is generated. The anisotropic constitutive model of carbon fiber laminate is invoked, the three-dimensional displacement field data is substituted into the strain tensor solution formula, and the spatial derivative is used to convert it into local strain components. Combined with the laminate layup direction to correct the strain calculation deviation, a preliminary three-dimensional strain field inversion result is generated. The preliminary three-dimensional strain field inversion results are coupled with the impact load parameters and the laminate material property parameters. Abnormal strain components are eliminated through strain-load correlation verification to generate the final three-dimensional strain field inversion results of carbon fiber laminate that are consistent with the dynamic response of the impact event.

2. The method according to claim 1, characterized in that, The process involves placing dynamic tracking markers on the surface of the impacted carbon fiber laminate, simultaneously acquiring image data of the entire impact process using a high-speed binocular camera, recording the magnitude and duration of the impact load, and unifying the image coordinate system through a binocular camera intrinsic parameter calibration algorithm to generate a binocular visual image sequence with impact time-series information, including: The dynamic tracking markers are designed as concentric rings with alternating fluorescent and black colors. The rings are 2mm in diameter and 0.5mm apart. They are printed on a 50μm thick polyimide film using UV-curable adhesive to ensure that the markers have both high contrast and anti-motion blur characteristics under high-speed shooting. The output is the marker placement scheme. According to the marking point layout plan, the surface of the carbon fiber laminate is divided into an impact center area and an edge area. Marking points are pasted in the center area with a grid density of 5mm×5mm and in the edge area with a grid density of 10mm×10mm. The height difference of the marking points is controlled within 0.1mm using a laser thickness gauge. The laminate specimen with the marking points laid out is then output. The Zhang Zhengyou calibration method was used to perform intrinsic parameter calibration on the high-speed binocular camera. Twenty sets of images with different poses were acquired using a 12×9 checkerboard target. The focal length principal point coordinate distortion coefficients of the two cameras were calculated. The frame rate of the two cameras was locked at 1000fps and the trigger time difference was controlled within 1μs through the synchronous trigger module. The calibrated binocular camera system was then output. The laminate specimen with marked points is fixed on the impact test bench. The calibrated binocular camera system is started, and the drop hammer impact device is triggered synchronously. The magnitude and duration of the impact load are recorded by the force sensor. Binocular images of the entire impact process are acquired. The coordinate system of the two cameras is unified by the internal parameter data, and a binocular visual image sequence with impact time information is generated.

3. The method according to claim 2, characterized in that, Based on the binocular vision image sequence, a hierarchical feature matching algorithm is used to extract the dynamic features of the marker points and the fiber texture of the laminate. The three-dimensional spatial coordinates of the feature points at each time step are calculated through stereo matching. The three-dimensional displacement vectors of the feature points are then solved by combining the differences between adjacent time step coordinates, generating three-dimensional displacement field data of the impact region of the carbon fiber laminate, including: The binocular vision image sequence is preprocessed, adaptive bilateral filtering is used to remove high-speed imaging noise, HSV color space threshold segmentation is performed on the marked point region, Canny edge detection is performed on the fiber texture region, and the preprocessed image sequence is output. Based on the preprocessed image sequence, a hierarchical feature matching algorithm is adopted. The upper layer extracts ORB features from the marker points and performs brute-force matching, while the lower layer extracts histogram features of orientation gradients from the fiber texture and performs K-nearest neighbor matching to generate a sequence of matching pairs between marker points and texture features. The matching sequence is input into the stereo matching algorithm. The three-dimensional spatial coordinates of the feature points at each time step are calculated by triangulation in combination with the intrinsic parameters of the binocular camera. The coordinate accuracy is optimized by the bundle adjustment method, and the three-dimensional coordinate set of the feature points at each time step is output. For each time series feature point's three-dimensional coordinate set, the coordinate difference between adjacent time series is calculated to obtain the three-dimensional displacement vector. The discrete vector is extended into a continuous field using the inverse distance weighted interpolation method. Gaussian smoothing is used to eliminate interpolation errors, generating three-dimensional displacement field data of the carbon fiber laminate impact region.

4. The method according to claim 3, characterized in that, The process involves calling the anisotropic constitutive model of the carbon fiber laminate, substituting the three-dimensional displacement field data into the strain tensor solution formula, converting it into local strain components through spatial derivative calculations, and correcting strain calculation errors based on the laminate layup direction to generate preliminary three-dimensional strain field inversion results, including: Construct an anisotropic constitutive model of carbon fiber laminates, input ply parameters and material parameters, establish an elastic matrix based on the ply direction, and output the parameter set of the anisotropic constitutive model; The three-dimensional displacement field data was discretized into a 1mm×1mm×0.5mm grid. The spatial partial derivatives of the displacement components at each grid point were calculated using the second-order central difference method. The results were then substituted into the strain tensor solution formula to obtain 6 strain components, which were used as the initial strain component matrix. Based on the laminate ply direction, the local coordinate system rotation matrix corresponding to each grid point is calculated. The initial strain component matrix is ​​transformed into strain components in the local coordinate system through coordinate transformation. The calculation deviation caused by ply anisotropy is corrected, and the corrected strain component matrix is ​​output. The corrected strain component matrix is ​​mechanically consistent with the parameter set of the anisotropic constitutive model. Strain values ​​exceeding the elastic limit of the material are removed, and the removed region is completed by linear interpolation to generate preliminary three-dimensional strain field inversion results.

5. The method according to claim 4, characterized in that, The process involves coupling the preliminary three-dimensional strain field inversion results with impact load parameters and laminate material property parameters, eliminating anomalous strain components through strain-load correlation verification, and generating the final three-dimensional strain field inversion result of the carbon fiber laminate consistent with the dynamic response of the impact event. This includes: We organize the impact load parameters and the laminate material property parameters, construct a parameter association database, map parameters of different dimensions to the same dimension through normalization, and output a set of coupled operation parameters. The preliminary three-dimensional strain field inversion results are coupled with the coupling operation parameter set to establish a strain-load time domain correlation model. The cross-correlation coefficient R between each strain component and the load curve is calculated. R ≥ 0.8 is set as an effective correlation, and the strain-load correlation matrix is ​​output. Based on the strain-load correlation matrix, abnormal strain components with cross-correlation coefficient R < 0.8 or strain values ​​exceeding the preset range are marked. The rationality of the abnormal marking is verified by combining the stress concentration law in the impact area, and an abnormal strain marking map is output. Abnormal components in the abnormal strain marker map are removed, and weighted interpolation is used to complete the data in the removed area. Finally, the three-dimensional strain field of the carbon fiber laminate that is consistent with the dynamic response of the impact event is generated.

6. A three-dimensional strain field inversion system for laminated plate impact events based on binocular vision, characterized in that, The system includes: The acquisition module is used to set dynamic tracking markers on the surface of the impacted carbon fiber laminate. It uses a high-speed binocular camera to synchronously acquire image data of the entire impact process, and records the magnitude and duration of the impact load. The binocular camera intrinsic parameter calibration algorithm unifies the image coordinate system and generates a binocular visual image sequence with impact time information. The extraction module is used to extract the dynamic features of the marker points and the fiber texture of the laminate based on the binocular vision image sequence using a hierarchical feature matching algorithm, calculate the three-dimensional spatial coordinates of the feature points in each time sequence through stereo matching, and solve the three-dimensional displacement vector of the feature points by combining the difference between adjacent time sequence coordinates to generate three-dimensional displacement field data of the impact area of ​​the carbon fiber laminate. The calling module is used to call the anisotropic constitutive model of carbon fiber laminate, substitute the three-dimensional displacement field data into the strain tensor solution formula, convert it into local strain components through spatial derivative calculation, and combine the laminate layup direction to correct the strain calculation deviation and generate preliminary three-dimensional strain field inversion results. The generation module is used to couple the preliminary three-dimensional strain field inversion results with the impact load parameters and the laminate material property parameters, and eliminate abnormal strain components through strain-load correlation verification to generate the final three-dimensional strain field inversion result of carbon fiber laminate that is consistent with the dynamic response of the impact event.

7. The system according to claim 6, characterized in that, The acquisition module is specifically used for: The dynamic tracking markers are designed as concentric rings with alternating fluorescent and black colors. The rings are 2mm in diameter and 0.5mm apart. They are printed on a 50μm thick polyimide film using UV-curable adhesive to ensure that the markers have both high contrast and anti-motion blur characteristics under high-speed shooting. The output is the marker placement scheme. According to the marking point layout plan, the surface of the carbon fiber laminate is divided into an impact center area and an edge area. Marking points are pasted in the center area with a grid density of 5mm×5mm and in the edge area with a grid density of 10mm×10mm. The height difference of the marking points is controlled within 0.1mm using a laser thickness gauge. The laminate specimen with the marking points laid out is then output. The Zhang Zhengyou calibration method was used to perform intrinsic parameter calibration on the high-speed binocular camera. Twenty sets of images with different poses were acquired using a 12×9 checkerboard target. The focal length principal point coordinate distortion coefficients of the two cameras were calculated. The frame rate of the two cameras was locked at 1000fps and the trigger time difference was controlled within 1μs through the synchronous trigger module. The calibrated binocular camera system was then output. The laminate specimen with marked points is fixed on the impact test bench. The calibrated binocular camera system is started, and the drop hammer impact device is triggered synchronously. The magnitude and duration of the impact load are recorded by the force sensor. Binocular images of the entire impact process are acquired. The coordinate system of the two cameras is unified by the internal parameter data, and a binocular visual image sequence with impact time information is generated.

8. The system according to claim 7, characterized in that, The extraction module is specifically used for: The binocular vision image sequence is preprocessed, adaptive bilateral filtering is used to remove high-speed imaging noise, HSV color space threshold segmentation is performed on the marked point region, Canny edge detection is performed on the fiber texture region, and the preprocessed image sequence is output. Based on the preprocessed image sequence, a hierarchical feature matching algorithm is adopted. The upper layer extracts ORB features from the marker points and performs brute-force matching, while the lower layer extracts histogram features of orientation gradients from the fiber texture and performs K-nearest neighbor matching to generate a sequence of matching pairs between marker points and texture features. The matching sequence is input into the stereo matching algorithm. The three-dimensional spatial coordinates of the feature points at each time step are calculated by triangulation in combination with the intrinsic parameters of the binocular camera. The coordinate accuracy is optimized by the bundle adjustment method, and the three-dimensional coordinate set of the feature points at each time step is output. For each time series feature point's three-dimensional coordinate set, the coordinate difference between adjacent time series is calculated to obtain the three-dimensional displacement vector. The discrete vector is extended into a continuous field using the inverse distance weighted interpolation method. Gaussian smoothing is used to eliminate interpolation errors, generating three-dimensional displacement field data of the carbon fiber laminate impact region.

9. A storage medium, characterized in that, The storage medium stores a computer program, wherein the computer program is configured to execute the method of any one of claims 1-5 when it is run.

10. An electronic device comprising a memory and a processor, characterized in that, The memory stores a computer program, and the processor is configured to run the computer program to perform the method of any one of claims 1-5.

Images