A farmland road surface spectrum detection method based on a small spherical mobile platform

Abstract

The application provides a farmland road surface spectrum detection method based on a small spherical mobile platform, and belongs to the technical field of farmland road surface spectrum detection.The application synchronously collects triaxial acceleration, angular velocity and pressure data when driving on the farmland road surface, calculates centrifugal acceleration and tangential acceleration generated by rotational motion by using an angular velocity vector measured by a gyroscope and a sensor position vector, subtracts these kinematic interference components from original acceleration to realize kinematic compensation, fuses the compensated acceleration with an acceleration estimation true value provided by a pressure sensor by using a Kalman filter to realize noise reduction and spectrum analysis to obtain a smooth displacement power spectrum density curve, finally determines a road surface unevenness grade by dynamic time warping matching and calculates a road surface unevenness coefficient, and solves the technical problem that kinematic interference components generated by sensor rotation pollute acceleration measurement values, so that road surface vibration characteristics cannot be accurately extracted.

Description

Technical Field

[0001] This invention belongs to the field of farmland pavement spectrum detection technology, and more specifically, relates to a farmland pavement spectrum detection method based on a small spherical mobile platform. Background Technology

[0002] Farmland pavement spectrum detection is a crucial technique for assessing the smoothness of roads used by agricultural machinery. Traditional methods employ laser profilometers, inertial measurement units, or vehicle-mounted accelerometers to collect vibration data on farmland surfaces and obtain the pavement power spectral density curve through spectrum analysis to determine the pavement unevenness level. However, in actual testing, when the testing platform travels on complex farmland surfaces, it generates rotational motions such as pitch, roll, and yaw. These rotational motions introduce additional kinematic interference components, such as centrifugal and tangential acceleration, at the accelerometer measurement locations. These interference components are superimposed on the actual vibration acceleration caused by pavement excitation. In other words, existing technologies suffer from the technical problem of kinematic interference components generated by sensor rotational motion contaminating acceleration measurements and hindering the accurate extraction of pavement vibration characteristics. Summary of the Invention

[0003] In view of this, the present invention provides a method for detecting the pavement spectrum of farmland based on a small spherical mobile platform, which can solve the technical problem in the prior art that the acceleration measurement value is mixed with kinematic interference components due to the rotational motion of the sensor during the detection of the pavement spectrum of farmland, thus making it impossible to accurately extract the pure vibration acceleration caused by the pavement excitation.

[0004] This invention is implemented as follows: It provides a method for detecting the pavement spectrum of farmland using a small spherical mobile platform. An accelerometer, tilt sensor, gyroscope, pressure sensor, GNSS / BeiDou positioning module, remote control module, temperature sensor, and humidity sensor are installed on the spherical mobile platform. The remote control module controls the spherical mobile platform to move along the farmland pavement to be tested and simultaneously collect data. Kinematic compensation is performed on the raw longitudinal axis acceleration collected by the accelerometer to obtain the longitudinal axis linear acceleration. The longitudinal axis linear acceleration is input into a Kalman filter for noise reduction to obtain the optimal longitudinal axis acceleration estimate sequence. The optimal longitudinal axis acceleration estimate sequence is then integrated twice and a high-pass filter is applied. The filter obtains the vibration displacement sequence; the vibration displacement sequence is segmented and weighted, and the Welch method is used to calculate the smooth displacement power spectral density curve; the smooth displacement power spectral density curve is dynamically time-warped and matched with the standard road surface power spectral density curve to obtain the road surface roughness level; the frequency index and road surface roughness coefficient are obtained by linear fitting based on the double logarithmic coordinate relationship between the smooth displacement power spectral density curve and the spatial frequency; the pose map is constructed using the spatial position coordinate sequence recorded by the GNSS Beidou positioning module to obtain the globally consistent trajectory estimation; the road surface roughness level, road surface roughness coefficient, frequency index, smooth displacement power spectral density curve and spatial position coordinates are correlated to draw the road surface spectrum image.

[0005] The accelerometer, tilt sensor, and gyroscope are installed on the inner platform of the spherical mobile platform, the pressure sensor is installed at the bottom of the spherical mobile platform, and the GNSS Beidou positioning module, remote control module, temperature sensor, and humidity sensor are installed on the outside of the spherical mobile platform.

[0006] The synchronously acquired data includes the GNSS Beidou positioning module recording the spatial coordinates of the spherical mobile platform in real time, the pressure sensor collecting the pressure value along the longitudinal axis of the bottom of the spherical mobile platform, the accelerometer collecting the three-axis raw acceleration and quaternion attitude data of the inner platform, the gyroscope measuring the three-axis angular velocity of the spherical mobile platform, the tilt sensor measuring the pitch and roll angles of the spherical mobile platform, and the temperature and humidity sensors recording the ambient temperature and humidity values.

[0007] Before the kinematic compensation, the motion state of the spherical moving platform is determined based on the absolute value of the time derivative of the pressure value collected by the pressure sensor, and the sensor data stream is segmented. When the absolute value of the time derivative of the pressure value is less than the preset motion state threshold, it is determined to be a uniform motion state. Zero speed correction is performed at the beginning and end of the uniform motion state segment to force the cumulative speed deviation and cumulative displacement deviation to zero.

[0008] Wherein, the preset motion state threshold ∈ [1, 10] N / s, and the sensor data stream consists of the triaxial raw acceleration collected by the accelerometer, the triaxial angular velocity measured by the gyroscope, and the pitch and roll angles measured by the tilt sensor.

[0009] Specifically, the kinematic compensation involves calculating the centrifugal acceleration component and the tangential acceleration component using the three-axis angular velocity vector measured by the gyroscope and the position vector of the accelerometer relative to the center of mass of the spherical moving platform. The longitudinal axis linear acceleration is obtained by subtracting the centrifugal acceleration component, the tangential acceleration component, and the gravitational acceleration component from the original longitudinal axis acceleration.

[0010] Specifically, the calculation of the centrifugal acceleration component involves performing a cross product operation between the triaxial angular velocity vector and the position vector to obtain a first intermediate vector, and then performing a cross product operation between the triaxial angular velocity vector and the first intermediate vector. The negative value of the result is the centrifugal acceleration component vector.

[0011] Specifically, the calculation of the tangential acceleration component involves taking the time derivative of the three-axis angular velocity vector to obtain the angular acceleration vector, and then performing a cross product operation between the angular acceleration vector and the position vector to obtain the tangential acceleration component vector.

[0012] Specifically, the Kalman filter performs noise reduction processing by dividing the pressure value collected by the pressure sensor by the mass of the spherical moving platform plus the gravitational acceleration component as the true value of the longitudinal axis acceleration estimation. The Kalman filter then fuses the true value of the longitudinal axis acceleration estimation with the longitudinal axis linear acceleration to output the optimal longitudinal axis acceleration estimation value.

[0013] The two cumulative integrations are performed using the trapezoidal integration rule, and the high-pass filter is used to remove the low-frequency drift trend of the vertical axis displacement sequence. The cutoff frequency of the high-pass filter is ∈ [0.01, 0.1] Hz.

[0014] Specifically, the segmented weighting involves dividing the vibration displacement sequence into multiple vibration displacement data segments according to a preset segment length, applying a Hanning window function to each vibration displacement data segment to obtain a weighted vibration displacement data segment, with a preset segment length ∈ [256, 2048] data points and an overlap rate between adjacent vibration displacement data segments ∈ [50%, 75%].

[0015] Specifically, the Welch method involves performing a fast Fourier transform on the weighted vibration displacement data segments and calculating the single-segment displacement power spectral density. The single-segment displacement power spectral density results of multiple vibration displacement data segments are then averaged to obtain a smooth displacement power spectral density curve.

[0016] Specifically, the dynamic time warping matching involves dynamically warping the smooth displacement power spectral density curve with the A to H grade standard pavement power spectral density curves in the standard pavement roughness classification table, calculating the dynamic time warping distance between the smooth displacement power spectral density curve and each standard pavement power spectral density curve, and selecting the pavement roughness level of the standard pavement power spectral density curve corresponding to the minimum dynamic time warping distance as the pavement roughness level of the farmland pavement to be tested.

[0017] Specifically, the frequency index is calculated by measuring the slope of the fitted straight line, and the road surface roughness coefficient is calculated by measuring the reference spatial frequency 0.1. The power spectral density value corresponding to the fitted straight line is used to obtain the road surface roughness coefficient.

[0018] Specifically, the construction of the pose graph involves using the spatial coordinates of the spherical mobile platform at key moments as pose graph nodes, and the relative motion transformations between adjacent pose graph nodes as pose graph edge constraints. When the spherical mobile platform is detected to be returning to a visited spatial position, a closure constraint edge is added between the starting pose graph node and the returning pose graph node. The globally consistent trajectory estimate is obtained by minimizing the weighted sum of squared residuals of all pose graph edge constraints using the Levenberg-Marquardt optimization algorithm.

[0019] The road surface spectrum image displays ambient temperature and humidity values, as well as power spectral density curves for standard E-grade and standard F-grade road surfaces.

[0020] This invention utilizes a spherical mobile platform carrying a multi-sensor fusion system. It employs gyroscopes to measure triaxial angular velocities in real time and calculates angular acceleration. Combining this with the position vector of the accelerometer relative to the center of mass, it accurately calculates the centrifugal and tangential acceleration components generated by rotational motion at the sensor position based on rigid body kinematics. These additional kinematic interference components and gravitational acceleration components are subtracted from the original longitudinal axis acceleration, thereby extracting the longitudinal axis linear acceleration purely caused by road surface excitation. Furthermore, this invention employs a Kalman filter to fuse the kinematically compensated acceleration with the true acceleration estimate provided by the pressure sensor, eliminating residual measurement noise and outputting the optimal acceleration estimate. Accurate vibration displacement sequences are obtained through trapezoidal integration and high-pass filtering. Finally, the Welch method is used to calculate the power spectral density curve reflecting the actual road surface roughness characteristics. In summary, this invention solves the technical problem mentioned in the background art where kinematic interference components generated by sensor rotational motion contaminate acceleration measurements, leading to inaccurate extraction of road surface vibration characteristics. Attached Figure Description

[0021] Figure 1 This is a time-domain waveform diagram of the vibration displacement of a spherical moving platform.

[0022] Figure 2 Comparison chart before and after global trajectory optimization.

[0023] Figure 3 This is a schematic diagram of the internal structure of a spherical mobile platform.

[0024] Figure 4 This is a schematic diagram of the chip and sensing layer structure of a spherical mobile platform. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0026] like Figure 1 The diagram shown is a flowchart of a farmland road surface spectrum detection method based on a small spherical mobile platform provided by the present invention. This method includes the following steps:

[0027] S01. Install an acceleration sensor, tilt sensor and gyroscope on the inner platform of the spherical mobile platform, install a pressure sensor on the bottom of the spherical mobile platform, and install a GNSS Beidou positioning module, remote control module, temperature sensor and humidity sensor on the outside of the spherical mobile platform.

[0028] S02. Start the spherical mobile platform and control it to move on the farmland road surface to be tested through the remote control module. The GNSS Beidou positioning module records the spatial position coordinates of the spherical mobile platform in real time, and the pressure sensor collects the pressure value in the longitudinal direction of the bottom of the spherical mobile platform simultaneously.

[0029] S03. The accelerometer collects the three-axis raw acceleration and quaternion attitude data of the inner platform, the gyroscope synchronously measures the three-axis angular velocity of the spherical moving platform, the tilt sensor measures the pitch and roll angles of the spherical moving platform, and the temperature and humidity sensors record the ambient temperature and humidity values.

[0030] S04. Determine the motion state of the spherical moving platform based on the absolute value of the time derivative of the pressure value collected by the pressure sensor, and segment the sensor data stream composed of the three-axis raw acceleration collected by the accelerometer, the three-axis angular velocity measured by the gyroscope, and the pitch angle and roll angle measured by the tilt sensor. When the absolute value of the time derivative of the pressure value is less than the preset motion state threshold, it is determined to be a uniform motion state. At the beginning and end of the uniform motion state segment, zero speed correction is performed to force the cumulative velocity deviation and cumulative displacement deviation to zero.

[0031] S05. Perform kinematic compensation on the original longitudinal axis acceleration collected by the accelerometer. Calculate the centrifugal acceleration component and tangential acceleration component using the three-axis angular velocity vector measured by the gyroscope and the position vector of the accelerometer relative to the center of mass of the spherical moving platform. Subtract the centrifugal acceleration component, tangential acceleration component, and gravitational acceleration component from the original longitudinal axis acceleration to obtain the longitudinal axis linear acceleration.

[0032] S06. Input the kinematically compensated longitudinal axis linear acceleration into the Kalman filter for noise reduction. Divide the pressure value collected by the pressure sensor by the mass of the spherical moving platform plus the gravitational acceleration component as the true value of the longitudinal axis acceleration estimation. The Kalman filter fuses the true value of the longitudinal axis acceleration estimation with the longitudinal axis linear acceleration to output the optimal longitudinal axis acceleration estimation value.

[0033] S07. The optimal longitudinal acceleration estimate sequence output by the Kalman filter is converted into a longitudinal displacement sequence by two cumulative integrations using the trapezoidal integral rule. A high-pass filter is then applied to the longitudinal displacement sequence to remove the low-frequency drift trend, resulting in the vibration displacement sequence of the longitudinal axis of the spherical moving platform relative to the average position of the longitudinal axis.

[0034] S08. Divide the vibration displacement sequence into multiple vibration displacement data segments according to the preset segment length. Apply the Hanning window function to each vibration displacement data segment to obtain a weighted vibration displacement data segment. Use the Welch method to perform a fast Fourier transform on the weighted vibration displacement data segment and calculate the single-segment displacement power spectral density. Average the single-segment displacement power spectral density results of multiple vibration displacement data segments to obtain a smooth displacement power spectral density curve.

[0035] S09. Perform dynamic time warping matching between the smooth displacement power spectral density curve and the standard pavement power spectral density curves of grades A to H in the standard pavement roughness classification table. Calculate the dynamic time warping distance between the smooth displacement power spectral density curve and each standard pavement power spectral density curve. Select the pavement roughness grade of the standard pavement power spectral density curve corresponding to the minimum dynamic time warping distance as the pavement roughness grade of the farmland pavement to be tested.

[0036] S10. Based on the logarithmic coordinate relationship between the smooth displacement power spectral density curve and the spatial frequency, a straight line is fitted to obtain the fitted line. The slope of the fitted line is calculated to obtain the frequency index, and the reference spatial frequency 0.1 is calculated. The power spectral density value corresponding to the fitted straight line is used to obtain the road surface roughness coefficient;

[0037] S11. Construct a pose graph using the spatial position coordinate sequence recorded by the GNSS Beidou positioning module. Use the spatial position coordinates of the spherical mobile platform at key moments as pose graph nodes. Use the relative motion transformation between adjacent pose graph nodes as pose graph edge constraints. When the spherical mobile platform is detected to return to the visited spatial position, add a closure constraint edge between the starting pose graph node and the returning pose graph node. Minimize the weighted residual sum of squares of all pose graph edge constraints using the Levenberg-Marquardt optimization algorithm to obtain a globally consistent trajectory estimate.

[0038] S12. Associate the road surface roughness grade, road surface roughness coefficient, frequency index, smooth displacement power spectral density curve with the spatial location coordinates in the global consistent trajectory estimation, draw the road surface spectrum image, and display the ambient temperature value, ambient humidity value, standard E-grade road surface power spectral density curve and standard F-grade road surface power spectral density curve next to the road surface spectrum image.

[0039] The preset motion state threshold is used to determine whether the spherical moving platform is in a uniform motion state. When the absolute value of the time derivative of the pressure value is less than the preset motion state threshold, it indicates that the force change at the bottom of the spherical moving platform is slow and is determined to be in a uniform motion state. The preset motion state threshold ranges from 1 N / s to 10 N / s.

[0040] The zero-speed correction refers to the technical means of eliminating integral drift by forcibly returning the cumulative speed deviation and cumulative displacement deviation to zero when the spherical moving platform is in a state of rest or uniform motion, in which theoretically the speed change should be zero. The cumulative speed deviation refers to the cumulative speed value calculated by integration from the last zero-speed correction time to the current time, and the cumulative displacement deviation refers to the cumulative displacement value calculated by double integration from the last zero-speed correction time to the current time.

[0041] The kinematic compensation refers to the process of measuring the three-axis angular velocity vector and calculating the angular acceleration vector of the spherical moving platform, using rigid body kinematics to calculate the additional acceleration generated by the rotational motion at the acceleration sensor position, and subtracting the additional acceleration from the original longitudinal axis acceleration to extract the longitudinal axis linear acceleration caused purely by road surface excitation.

[0042] The centrifugal acceleration component refers to the acceleration directed towards the center of rotation generated by the accelerometer due to the positional offset of the accelerometer from the center of rotation when the spherical moving platform rotates around a certain axis. The calculation formula for the centrifugal acceleration component is as follows: perform a cross product operation between the three-axis angular velocity vector and the position vector to obtain the first intermediate vector, then perform a cross product operation between the three-axis angular velocity vector and the first intermediate vector, and take the negative value of the result to obtain the centrifugal acceleration component vector.

[0043] The tangential acceleration component refers to the acceleration generated by the accelerometer due to positional shift when the angular velocity of the spherical moving platform changes, i.e., when an angular acceleration vector exists. This acceleration is perpendicular to the plane containing the position vector and the three-axis angular velocity vector. The calculation formula for the tangential acceleration component is as follows: the angular acceleration vector is obtained by taking the time derivative of the three-axis angular velocity vector, and the tangential acceleration component vector is obtained by performing a cross product operation between the angular acceleration vector and the position vector.

[0044] The position vector refers to the spatial vector of the accelerometer installation position relative to the center of mass of the spherical moving platform. The three components of the position vector represent the distances of the accelerometer relative to the center of mass in the longitudinal, transverse, and vertical directions of the spherical moving platform, respectively. The position vector is obtained by measurement through mechanical design drawings.

[0045] The Kalman filter described is a recursive state estimation algorithm. It establishes a state-space model of the system, predicts the state at the next time step using the longitudinal linear acceleration and the system dynamics model, then calculates the Kalman gain based on the prediction error covariance matrix and the measurement error covariance matrix. Finally, it fuses the predicted value and the true longitudinal acceleration estimate to obtain the optimal longitudinal acceleration estimate, and updates the error covariance matrix for the recursive calculation at the next time step. The state-space model includes state equations and observation equations. The state equations describe the evolution of the system state over time, while the observation equations describe the relationship between the measured values ​​and the system state.

[0046] The estimated true value of the longitudinal axis acceleration refers to the longitudinal axis support force acceleration obtained by dividing the pressure value collected by the pressure sensor by the mass of the spherical moving platform, plus the gravitational acceleration component. This value should be equal to the theoretical true value of the original longitudinal axis acceleration when the bottom of the spherical moving platform contacts the ground, and therefore serves as a reference benchmark for longitudinal axis acceleration measurement.

[0047] The trapezoidal integral rule is a numerical integration method that divides a continuous integration interval into several smaller intervals. Within each smaller interval, the integral value is approximated by the area of ​​a trapezoid, that is, by multiplying the average of the optimal longitudinal acceleration estimates at two adjacent times by the time interval. The total integral value is obtained by accumulating the trapezoidal areas between all the smaller intervals. The first integration yields the velocity, and the second integration yields the longitudinal displacement sequence.

[0048] The high-pass filter is a filter that allows high-frequency signals to pass through while attenuating low-frequency signals. In this method, it is used to remove the slow drift trend caused by the accumulation of integral errors in the longitudinal displacement sequence and retain the rapidly changing vibration components caused by road surface excitation. The cutoff frequency of the high-pass filter is set in the range of 0.01Hz to 0.1Hz.

[0049] The preset segment length refers to the number of data points contained in each segment when the vibration displacement sequence is divided into multiple vibration displacement data segments. The preset segment length ranges from 256 to 2048 data points, and the overlap rate between adjacent vibration displacement data segments is set from 50% to 75%.

[0050] The single-segment displacement power spectral density refers to the square of the frequency domain amplitude obtained by performing a fast Fourier transform on a single weighted vibration displacement data segment, divided by the data segment length and the sampling frequency, representing the energy distribution of the displacement signal at different spatial frequencies.

[0051] The standard pavement roughness classification table refers to the pavement classification standard specified in the international standard ISO 8608. Based on the pavement roughness coefficient, the pavement is divided into 8 levels from A to H, and each level corresponds to a standard pavement power spectral density curve.

[0052] The dynamic time warping matching is an algorithm for measuring the similarity between two time series. It allows the time axis to be non-linearly stretched to find the best alignment. By constructing a cumulative distance matrix, a dynamic programming method is used to search for the path that minimizes the sum of the distances between the corresponding points of the smooth displacement power spectral density curve and the standard road surface power spectral density curve.

[0053] The dynamic time warping distance refers to the minimum cumulative distance obtained by the dynamic time warping matching algorithm. The smaller the dynamic time warping distance, the more similar the smooth displacement power spectral density curve is to the standard road surface power spectral density curve.

[0054] The aforementioned double logarithmic coordinate relationship refers to plotting the horizontal axis spatial frequency and vertical axis power spectral density values ​​of the smooth displacement power spectral density curve in a coordinate system after taking the logarithm of both. In the double logarithmic coordinate system, the road surface power spectral density curve approximately presents a linear relationship.

[0055] The fitted straight line refers to the straight line obtained by least squares linear fitting of the data points of the smooth displacement power spectral density curve in a double logarithmic coordinate system. The equation of the fitted straight line is that the logarithmic power spectral density value is equal to the slope multiplied by the logarithmic spatial frequency plus the intercept.

[0056] The frequency index is a parameter that describes the variation of the power spectral density of the road surface with spatial frequency. The frequency index is equal to the negative slope of the fitted straight line. In the standard road surface model, the frequency index is usually taken as 2, which means that the power spectral density decreases with the square of the spatial frequency.

[0057] The road surface roughness coefficient is the value of the road surface power spectral density at a reference spatial frequency, which is defined as 0.1. The road surface roughness coefficient reflects the overall intensity of road surface elevation undulations; a higher value indicates a more uneven road surface. The formula for calculating the road surface roughness coefficient is as follows: (Referencing a spatial frequency of 0.1...) By substituting the logarithm into the fitted linear equation, the logarithmic power spectral density value is calculated. Then, by taking the exponential operation on the logarithmic power spectral density value, the road surface roughness coefficient is obtained, with units of [unit missing]. .

[0058] The pose graph is a graph structure. The pose graph nodes represent the spatial position coordinates and attitude angles of the spherical mobile platform at critical moments. The pose graph edge constraints represent the relative motion transformation relationship between adjacent pose graph nodes or the absolute position constraints provided by the GNSS Beidou positioning module. By optimizing the position and attitude of all pose graph nodes to minimize the sum of squared residuals of the pose graph edge constraints, globally consistent trajectory estimation is achieved.

[0059] The critical moment refers to the moment during the movement of the spherical mobile platform that meets the time interval condition or the spatial interval condition. The time interval condition is that the time difference between adjacent critical moments is greater than or equal to a preset time interval threshold. The spatial interval condition is that the spatial coordinate distance between adjacent critical moments is greater than or equal to a preset spatial interval threshold. The preset time interval threshold ranges from 0.5s to 2s, and the preset spatial interval threshold ranges from 0.5m to 2m.

[0060] The loop closure constraint edge refers to the pose graph edge constraint added between the initial pose graph node and the return pose graph node when the spherical mobile platform's trajectory forms a closed loop, i.e., returns to a previously visited spatial position. The loop closure constraint edge requires that the relative pose transformation between the initial pose graph node and the return pose graph node should be close to a unit transformation. The accumulated trajectory drift error is distributed across the entire closed-loop path through the loop closure constraint edge. The criterion for determining the visited spatial position is that the Euclidean distance between the current spatial position coordinates and the spatial position coordinates of the historical pose graph node is less than a preset loop closure detection threshold. The preset loop closure detection threshold ranges from 1m to 5m.

[0061] The Levenberg-Marquardt optimization algorithm is a nonlinear least squares optimization method that combines the Gauss-Newton method and the gradient descent method. By introducing a damping factor, a diagonal term is added to the Hessian matrix of the objective function. When the solution is far from the optimal solution, the damping factor is large, which means that gradient descent guarantees convergence. When the solution is close to the optimal solution, the damping factor is small, which means that the Gauss-Newton method converges quickly.

[0062] The weighted residual sum of squares refers to the sum of squares of the residual vectors of all pose graph edge constraints. Each residual vector is weighted according to its corresponding information matrix. The information matrix is ​​the inverse of the covariance matrix and reflects the confidence of the pose graph edge constraints. The Levenberg-Marquardt optimization algorithm iteratively adjusts the position and orientation of the pose graph nodes to minimize the weighted residual sum of squares.

[0063] The globally consistent trajectory estimation refers to the motion trajectory of the spherical mobile platform obtained through pose graph optimization. Globally consistent trajectory estimation eliminates the cumulative drift over long periods of operation, ensuring that the entire trajectory remains consistent in the global coordinate system, and providing an accurate reference for road surface spatial positioning.

[0064] The road surface spectrum image refers to a two-dimensional image plotted with the smooth displacement power spectral density curve as the horizontal axis and the power spectral density value as the vertical axis. The road surface spectrum image intuitively reflects the vibration intensity distribution characteristics of the farmland road surface under different wavelengths.

[0065] The power spectral density curves for standard Grade E and standard Grade F pavements are standard curves corresponding to two common pavement roughness grades in farmland. The geometric mean of the pavement roughness coefficient for standard Grade E pavement is... The geometric mean of the road surface roughness coefficient for standard Class F road surface is: In the spatial frequency range of 0.011 Up to 2.83 The corresponding root mean square values ​​and geometric mean values ​​are respectively and .

[0066] The relationship between road surface power spectral density, road surface roughness coefficient, and frequency index is expressed by the following formula: the road surface power spectral density value equals the road surface roughness coefficient multiplied by the negative frequency exponent of the ratio of the spatial frequency to the reference spatial frequency. The conversion formula between velocity power spectral density and displacement power spectral density is expressed as follows: the velocity power spectral density value equals 2 multiplied by pi multiplied by the square of the spatial frequency, then multiplied by the displacement power spectral density value. The conversion formula between acceleration power spectral density and displacement power spectral density is expressed as follows: the acceleration power spectral density value equals 2 multiplied by pi multiplied by the fourth power of the spatial frequency, then multiplied by the displacement power spectral density value.

[0067] The specific implementation methods of the above steps are described in detail below.

[0068] The specific implementation of step S01 involves fixing a triaxial accelerometer, a dual-axis tilt sensor, and a triaxial gyroscope onto the inner platform of the spherical mobile platform, ensuring that the measurement axes of each sensor are aligned with the body coordinate system of the spherical mobile platform. A pressure sensor array or single-point pressure sensor is installed in the contact area between the bottom spherical shell of the spherical mobile platform and the ground to measure pressure changes along the longitudinal axis. An antenna for a GNSS BeiDou positioning module is installed on the top of the outer shell of the spherical mobile platform to obtain optimal satellite signal reception. A receiver for a remote control module is integrated on the control circuit board of the inner platform. Temperature and humidity sensors are installed at ventilation locations on the outer shell of the spherical mobile platform to monitor environmental parameters in real time. The purpose of this step is to construct a complete multi-sensor data acquisition system to provide a multi-dimensional physical quantity measurement basis for subsequent road surface spectrum detection.

[0069] The specific implementation of step S02 is as follows: the remote controller sends control commands to the remote control module, which drives the built-in motor of the spherical mobile platform to generate rolling motion of the spherical shell, thereby realizing the movement on the farmland road surface to be tested. The GNSS Beidou positioning module receives satellite signals in real time at a set sampling frequency and calculates the three-dimensional spatial coordinates of the spherical mobile platform, including longitude, latitude and altitude. The pressure sensor synchronously collects the pressure value of the bottom of the spherical mobile platform in the longitudinal direction, that is, the direction perpendicular to the ground, at the same or higher sampling frequency as the GNSS positioning module. The purpose of this step is to establish the spatiotemporal correspondence between the motion trajectory of the spherical mobile platform and the road surface excitation signal, so as to lay the foundation for the spatial distribution characteristics analysis of the road surface spectrum.

[0070] The specific implementation of step S03 involves the accelerometer acquiring the original acceleration values ​​of the inner platform along the three orthogonal axes and the quaternion parameters describing the spatial attitude of the inner platform at a high sampling frequency. The gyroscope simultaneously measures the instantaneous angular velocity values ​​of the spherical moving platform around the three axes to capture the dynamic characteristics of the rotational motion. The tilt sensor uses the principle of gravity sensing to measure the pitch and roll angles of the spherical moving platform relative to the horizontal plane. The temperature and humidity sensors record the ambient temperature and humidity values ​​during the test at a lower sampling frequency. The purpose of this step is to obtain complete descriptive information of the motion state of the spherical moving platform and environmental factor data to provide necessary input parameters for subsequent kinematic compensation and data quality assessment.

[0071] The specific implementation of step S04 is to perform numerical differentiation calculation on the pressure time series collected by the pressure sensor to obtain the rate of change of the pressure value with time, and compare the absolute value of the rate of change with the preset motion state threshold. When the absolute value is less than the preset motion state threshold, it indicates that the force on the bottom of the spherical moving platform is basically constant and is judged as a uniform motion state. At the beginning and end of the uniform motion state segment, a zero-speed correction operation is performed, that is, the speed deviation and displacement deviation accumulated since the last zero-speed correction moment are forcibly set to zero. The reference value of the preset motion state threshold is 5N / s. The purpose of this step is to identify the uniform motion range of the spherical moving platform and use the zero-speed correction technology to eliminate the drift error generated during the acceleration integration process and improve the accuracy of displacement estimation.

[0072] The specific implementation of step S05 is to calculate the influence of the rotational motion of the spherical moving platform on the acceleration sensor measurement value based on the rigid body kinematics principle. The first intermediate vector is obtained by cross-product operation of the three-axis angular velocity vector measured by the gyroscope and the position vector of the acceleration sensor relative to the center of mass of the spherical moving platform. Then, the centrifugal acceleration component vector is obtained by cross-product operation of the three-axis angular velocity vector and the first intermediate vector and taking the negative value. The angular acceleration vector is obtained by time differentiation of the three-axis angular velocity vector. The tangential acceleration component vector is obtained by cross-product operation of the angular acceleration vector and the position vector. The longitudinal linear acceleration is obtained by subtracting the longitudinal component of the centrifugal acceleration component, the longitudinal component of the tangential acceleration component, and the gravitational acceleration component from the original longitudinal acceleration collected by the acceleration sensor. The purpose of this step is to remove the interference of the rotational motion of the spherical moving platform itself on the acceleration measurement through kinematic compensation technology, thereby extracting the linear acceleration signal purely caused by the road surface unevenness excitation.

[0073] The specific implementation of step S06 involves constructing a state-space model of a Kalman filter. The state equation describes the evolution of the longitudinal acceleration over time, and the observation equation describes the relationship between the measured value and the system state. The kinematically compensated longitudinal linear acceleration and the system dynamics model are used to predict the state value at the next moment. The pressure value collected by the pressure sensor is divided by the mass of the spherical moving platform, and the gravitational acceleration component is added to obtain the estimated true value of the longitudinal acceleration. The Kalman gain is calculated based on the prediction error covariance matrix and the measurement error covariance matrix. The optimal longitudinal acceleration estimate is obtained by fusing the predicted value and the estimated true value of the longitudinal acceleration using the Kalman gain. The error covariance matrix is ​​updated for recursive calculation at the next moment. The purpose of this step is to achieve the optimal estimate of the longitudinal acceleration by fusing the measurement information of the acceleration sensor and the pressure sensor through the Kalman filter algorithm, effectively reducing the impact of measurement noise on subsequent integration calculations.

[0074] The specific implementation of step S07 is to perform the first integration of the optimal longitudinal axis acceleration estimate sequence output by the Kalman filter using the trapezoidal integral rule. In each small time interval, the velocity increment of the interval is obtained by multiplying the average of the optimal longitudinal axis acceleration estimates of two adjacent times by the time interval. All velocity increments are accumulated to obtain the velocity sequence. The velocity sequence is then integrated a second time using the trapezoidal integral rule to obtain the longitudinal axis displacement sequence. A high-pass filter is designed and the cutoff frequency is set to 0.05Hz. The high-pass filter is applied to the longitudinal axis displacement sequence to remove the low-frequency drift trend caused by the accumulation of integration error and retain the high-frequency vibration components caused by road surface excitation to obtain the vibration displacement sequence of the longitudinal axis of the spherical moving platform relative to the average position of the longitudinal axis. The purpose of this step is to convert the acceleration signal into a displacement signal and eliminate the integration drift through the high-pass filter to obtain pure vibration displacement data that can reflect the characteristics of road surface unevenness.

[0075] The specific implementation of step S08 involves dividing the vibration displacement sequence into multiple vibration displacement data segments according to a preset segment length. An overlap rate of 50% to 75% is set between adjacent vibration displacement data segments to increase the number of effective data segments. The reference value for the preset segment length is 1024 data points. A Hanning window function is applied to each vibration displacement data segment for weighting. The Hanning window function value is calculated by subtracting 0.5 multiplied by the cosine value from 0.5. The independent variable of the cosine value is 2 multiplied by pi multiplied by the current sample number divided by the total number of samples minus 1. The Welch method is used to perform a fast Fourier transform on the weighted vibration displacement data segments to convert the time-domain signal into a frequency-domain representation. The periodogram of each weighted vibration displacement data segment, i.e., the single-segment displacement power spectral density, is calculated. The arithmetic mean of the single-segment displacement power spectral density results of multiple vibration displacement data segments is used to obtain a smoothed displacement power spectral density curve. The purpose of this step is to reduce the variance and spectral leakage of the power spectral density estimation by combining the Welch method and the Hanning window function, thereby improving the stability and accuracy of the road surface power spectral density estimation.

[0076] The specific implementation of step S09 involves obtaining eight standard pavement power spectral density curves (grades A to H) from the standard pavement roughness classification table according to the international standard ISO 8608. A dynamic time warping algorithm is used to measure the similarity between the smooth displacement power spectral density curve and each standard pavement power spectral density curve. The dynamic time warping algorithm constructs a cumulative distance matrix and uses dynamic programming to search for the nonlinear alignment path that minimizes the sum of the distances between corresponding points on the two curves. The dynamic time warping distance between the smooth displacement power spectral density curve and each standard pavement power spectral density curve is calculated. The pavement roughness grade of the standard pavement power spectral density curve corresponding to the minimum dynamic time warping distance is selected as the pavement roughness grade of the farmland pavement to be tested. The purpose of this step is to achieve automated classification and determination of the roughness grade of the farmland pavement to be tested through the dynamic time warping matching algorithm, providing qualitative indicators for pavement quality assessment.

[0077] The specific implementation of step S10 involves logarithmically transforming both the horizontal axis spatial frequency and the vertical axis power spectral density values ​​of the smoothed displacement power spectral density curve into a double logarithmic coordinate system. In this double logarithmic coordinate system, the data points of the smoothed displacement power spectral density curve are fitted with a straight line using the least squares method to obtain a fitted straight line. The slope of the fitted straight line is the negative value of the frequency exponent. The reference spatial frequency is 0.1. The logarithmic power spectral density is obtained by substituting the logarithm into the fitted linear equation. The road surface roughness coefficient is then calculated by taking the exponent of the logarithmic power spectral density. The unit of the road surface roughness coefficient is... The purpose of this step is to extract the characteristic parameters of the road power spectral density curve, namely the frequency index and the road roughness coefficient, to provide standardized numerical indicators for the quantitative description of road roughness and the lateral comparison of different road sections.

[0078] The specific implementation of step S11 involves constructing a pose graph data structure using the spatial position coordinate sequence recorded by the GNSS BeiDou positioning module. At critical moments that meet the time interval or spatial interval conditions, the spatial position coordinates and attitude angles of the spherical mobile platform are added as pose graph nodes to the pose graph. The reference value for the preset time interval threshold is 1 second, and the reference value for the preset spatial interval threshold is 1 meter. The relative motion transformation between adjacent pose graph nodes is used as the pose graph edge constraint. When the Euclidean distance between the current spatial position coordinates and the spatial position coordinates of historical pose graph nodes is detected to be less than the preset loop closure detection threshold, it is determined that the spherical mobile platform has returned to normal operation. The spatial location is accessed, and the reference value for the preset loop closure detection threshold is 3 meters. Loop closure constraint edges are added between the starting pose graph node and the returning pose graph node. The Levenberg-Marquardt optimization algorithm is used to iteratively adjust the position and attitude of all pose graph nodes. The Levenberg-Marquardt optimization algorithm combines the advantages of the Gauss-Newton method and the gradient descent method by introducing a damping factor, and minimizes the weighted sum of squared residuals of all pose graph edge constraints to obtain a globally consistent trajectory estimate. The purpose of this step is to eliminate the cumulative drift error caused by long-term operation through pose graph optimization and loop closure detection technology, so as to provide a reliable globally consistent trajectory reference for accurate spatial positioning of the road surface spectrum.

[0079] The specific implementation of step S12 involves establishing a corresponding relationship between the road surface roughness grade, road surface roughness coefficient, frequency index, smooth displacement power spectral density curve, and spatial location coordinates in the globally consistent trajectory estimation. A road surface spectrum image is plotted with spatial frequency as the horizontal axis and power spectral density value as the vertical axis to visually display the vibration intensity distribution characteristics of the farmland road surface under test at different wavelengths. The ambient temperature and humidity values ​​recorded during the test are displayed next to the road surface spectrum image to reflect the influence of environmental factors on the test results. Simultaneously, the power spectral density curves of standard E-grade road surface and standard F-grade road surface are overlaid on the road surface spectrum image as reference benchmarks. The geometric mean of the road surface roughness coefficient of standard E-grade road surface is... The geometric mean of the road surface roughness coefficient of standard Class F road surface is: The purpose of this step is to generate a comprehensive pavement spectrum test report containing qualitative and quantitative parameters and spatial distribution information, providing comprehensive data support for farmland pavement quality assessment and agricultural machinery design optimization.

[0080] It should be noted that this invention also solves the following technical problem: In the process of farmland road surface spectrum detection, the long-term movement of the detection platform causes cumulative drift errors in the acceleration integral calculation, resulting in displacement calculation results deviating from the true value, thus affecting the accuracy of the power spectral density curve. This invention monitors the changes in bottom force in real time using pressure sensors, determines the motion state of the spherical moving platform based on the absolute value of the time derivative of the pressure value, and performs zero-velocity correction at the beginning and end of the state segment when uniform motion is detected, forcibly returning the cumulative velocity deviation and cumulative displacement deviation to zero, effectively eliminating the influence of integral drift on displacement calculation. In addition, this invention also solves the technical problem of inaccurate spatial positioning of the road surface spectrum caused by long-term drift of the detection trajectory. By constructing a pose graph and using the spatial position coordinates recorded by the GNSS Beidou positioning module as nodes, when the spherical moving platform is detected to return to a visited position, a closure constraint edge is added. The Levenberg-Marquardt optimization algorithm is used to minimize the weighted sum of squared residuals of all edge constraints, distributing the cumulative trajectory drift error to the entire closed-loop path, achieving globally consistent trajectory estimation, and providing accurate spatial positioning reference for the road surface spectrum.

[0081] Specifically, the principle of this invention is as follows: This invention solves the technical problem of kinematic interference introduced by the rotational motion of the sensor. Its core principle lies in establishing a complete separation mechanism from the original measurement value of the accelerometer to the pure vibration acceleration caused by the road surface excitation. When the spherical moving platform travels on the farmland road, the accelerometer not only senses the linear vibration acceleration caused by the road surface excitation, but also generates centrifugal acceleration and tangential acceleration at its installation position due to the platform's rotational motion. The magnitudes of these two additional accelerations depend on the angular velocity, angular acceleration, and the sensor's position offset relative to the center of mass. This invention accurately measures the three-axis angular velocity vector using a gyroscope, calculates the accurate values ​​of centrifugal and tangential acceleration using cross product operations and time derivatives, and then subtracts these calculable kinematic interference components from the original acceleration measurement value to obtain the pure linear acceleration generated solely by the road surface unevenness excitation. This kinematic compensation process eliminates the contamination of vibration measurement by rotational motion from a physical mechanism perspective, ensuring that subsequent power spectral density analysis can truly reflect the unevenness characteristics of the farmland road surface. Therefore, the technical solution of this invention conforms to the basic principles of rigid body kinematics and has sufficient logical rationality.

[0082] The following provides a specific embodiment 1 of the present invention. The specific implementation methods of steps S01, S02 and S03 in this embodiment 1 are the same as those described above, and will not be repeated in detail here. The specific implementation methods of other steps are described in detail below.

[0083] The specific implementation of step S04 involves calculating the absolute value sequence of the time derivative based on the pressure value sequence collected by the pressure sensor, determining the motion state of the spherical moving platform, and segmenting the sensor data stream. The formula for calculating the time derivative of the pressure value is as follows:

[0084] ;

[0085] In the formula, The time derivative of the pressure value, in units of ; For the first The pressure value collected by the pressure sensor at all times, in units of ; For the first The pressure value collected by the pressure sensor at all times, in units of ; The time interval between two consecutive samples, in units of The value is typically between 0.01 and 0.05. The motion state determination conditions are described as follows:

[0086] ;

[0087] In the formula, As a marker of motion state, Indicates a state of uniform motion. This indicates a state of non-uniform motion; The preset motion state threshold is in units of The value range is 1 to 10; For reference time derivative, the unit is . The empirical value is 5, used for dimensionless processing. When a uniform motion segment is detected, at the beginning of the uniform motion segment... and end time Perform zero-speed correction; the formula for zero-speed correction is as follows:

[0088] , ;

[0089] , ;

[0090] In the formula, Cumulative speed deviation, in units of ; Cumulative displacement deviation, in units of ; The starting time of the uniform motion segment, in units of ; The time at the end of the uniform motion segment, in units of .

[0091] The specific implementation of step S05 involves performing kinematic compensation on the raw longitudinal axis acceleration collected by the accelerometer to extract the longitudinal axis linear acceleration purely caused by road surface excitation. The calculation formula for the centrifugal acceleration component is as follows:

[0092] ;

[0093] In the formula, This is the centrifugal acceleration component vector, in units of... ; This is a three-axis angular velocity vector, in units of... ; The position vector of the accelerometer relative to the center of mass of the spherical moving platform is given by . ; This represents the cross product operation. The formula for calculating the tangential acceleration component is as follows:

[0094] ;

[0095] In the formula, This is the tangential acceleration component vector, in units of... ; This is the angular acceleration vector, in units of . The angular acceleration vector is obtained by taking the time derivative of the triaxial angular velocity vector, and the calculation formula is as follows: ,in This is the derivative of angular velocity with respect to time, in units of . The formula for calculating longitudinal linear acceleration is as follows:

[0096] ;

[0097] In the formula, The longitudinal linear acceleration is expressed in units of 1 / 2 Ω. ; The original acceleration along the vertical axis, in units of ; This represents the component of centrifugal acceleration along the vertical axis, with units of . ; This represents the component of the tangential acceleration along the vertical axis, with units of . ; This represents the component of gravitational acceleration, in units of... The experience value is 9.81.

[0098] The specific implementation of step S06 involves inputting the kinematically compensated longitudinal axis linear acceleration into a Kalman filter for noise reduction to obtain the optimal longitudinal axis acceleration estimate. The formula for calculating the true value of the longitudinal axis acceleration estimate is as follows:

[0099] ;

[0100] In the formula, Estimate the true value of the vertical axis acceleration, in units of ; The pressure value collected by the pressure sensor, in units of ; The mass of the spherical moving platform is expressed in units of . ; This represents the component of gravitational acceleration, in units of... The state equation and observation equation of the Kalman filter are expressed as follows:

[0101] ;

[0102] ;

[0103] In the formula, For the first The state vector at time t, which includes acceleration and velocity along the vertical axis; This is the state transition matrix; The process noise has a mean of zero and a covariance matrix of... Gaussian distribution; For the first The observed value at time t, i.e., the true value of the longitudinal acceleration estimate. The unit is ; The observation matrix; The observed noise follows a zero-mean covariance matrix. The Gaussian distribution is used. The formula for calculating the Kalman gain is as follows:

[0104] ;

[0105] In the formula, For the first Kalman gain at time step; The prediction error covariance matrix; Observation matrix Transpose of; The observation noise covariance matrix is ​​used. The formula for calculating the optimal vertical axis acceleration estimate is as follows:

[0106] ;

[0107] In the formula, For the first The optimal estimate of the vertical acceleration at time t, in units of ; To predict acceleration values, the units are... ; For the first The true value of the longitudinal acceleration estimate at time t, in units of .

[0108] The specific implementation of step S07 involves performing two trapezoidal integrals on the optimal longitudinal acceleration estimate sequence output by the Kalman filter to obtain the longitudinal displacement sequence, and then applying a high-pass filter to remove low-frequency drift. The trapezoidal integral formula is expressed as follows:

[0109] ;

[0110] ;

[0111] In the formula, For the first Speed ​​at time, in units of ; For the first Speed ​​at time, in units of ; For the first The optimal estimate of the vertical acceleration at time t, in units of ; For the first The optimal estimate of the vertical acceleration at time t, in units of ; The time interval is expressed in units of 1 / 2. ; For the first The vertical displacement at time t, in units of ; For the first The vertical displacement at time t, in units of The transfer function of a high-pass filter is expressed as follows:

[0112] ;

[0113] In the formula, For high-pass filters at frequency The transfer function at the location; The signal time frequency, in units of ; Cutoff frequency, unit: The value ranges from 0.01 to 0.1. The formula for calculating the vibration displacement sequence is as follows:

[0114] ;

[0115] In the formula, For the first Vibration displacement at time t, in units of ; The vertical displacement is the result of high-pass filtering, in units of... ; This represents the average value of the vertical axis displacement sequence, in units of... .

[0116] The specific implementation of step S08 involves dividing the vibration displacement sequence into segments according to a preset segment length and applying a Hanning window function, then using the Welch method to calculate the smoothed displacement power spectral density curve. The calculation formula for the Hanning window function is as follows:

[0117] ;

[0118] In the formula, For the first Hanning window function values ​​for each sample point; This is the current sample number, with a value ranging from 0 to... ; The preset segment length ranges from 256 to 2048 data points. Let π be the mathematical constant pi. The formula for calculating the weighted vibration displacement data segment is as follows:

[0119] ;

[0120] In the formula, For the first The weighted vibration displacement value of each sample point, in units of ; For the first Vibration displacement values ​​at each sample point, in units of The formula for converting between time frequency and spatial frequency is expressed as follows:

[0121] ;

[0122] In the formula, For the first Spatial frequencies, in units of ; For the first Time frequency, unit is ; The average operating speed of the spherical moving platform, in units of By analyzing the velocity sequence The arithmetic mean is used to obtain the power spectral density of a single displacement segment. The formula for calculating the power spectral density of a single displacement segment is as follows:

[0123] ;

[0124] In the formula, For the first The single-segment displacement power spectral density at a spatial frequency point, in units of ; The result is the Fast Fourier Transform of the weighted vibration displacement data segment; For the first 1 spatial frequency point, unit is ; Spatial frequency resolution, in units of The calculation formula is: ,in Spatial sampling interval, in units of Through average running speed Multiply by the time interval Calculated The formula for calculating the smooth displacement power spectral density curve is as follows:

[0125] ;

[0126] In the formula, For the first Smooth displacement power spectral density at spatial frequency points, in units of ; This represents the total number of vibration displacement data segments; For the first The vibration displacement data segment in the first The single-segment displacement power spectral density at a spatial frequency point, in units of .

[0127] The specific implementation of step S09 involves dynamically time-warping the smooth displacement power spectral density curve with the power spectral density curves of standard pavements from grades A to H in the standard pavement roughness classification table to determine the pavement roughness level. The formula for calculating the dynamic time-warped distance is as follows:

[0128] ;

[0129] In the formula, For dynamic time-normalized distance; This represents the total number of nodes in the optimal matching path, which is obtained through a dynamic programming algorithm. To smooth the displacement power spectral density curve at path 1 Spatial frequency index at each node, in units of ; The standard road surface power spectral density curve at the path number Spatial frequency index at each node, in units of ; The power spectral density curve of a certain grade of standard road surface at spatial frequency The value at the location, in units of ; This is the normalized power spectral density value, in units of Experience value This is used for dimensionless processing. The formula for determining the road surface roughness level is expressed as follows:

[0130] ;

[0131] In the formula, The roughness level of the farmland road surface to be tested; To smooth the displacement power spectral density curve and level Dynamic time-normalized distance of standard road surface power spectral density curve; Indicates taking Lowest level .

[0132] The specific implementation of step S10 is to calculate the road surface roughness coefficient and frequency index by fitting a double logarithmic coordinate system. The fitted linear equation is expressed as follows:

[0133] ;

[0134] In the formula, This is the logarithm of the power spectral density of the smooth displacement. For the first Smooth displacement power spectral density at spatial frequency points, in units of ; For the first 1 spatial frequency point, unit is ; The slope of the fitted line; This is the intercept of the fitted line. The formula for calculating the frequency exponent is as follows:

[0135] ;

[0136] In the formula, This is the frequency index. The formula for calculating the road surface roughness coefficient is as follows:

[0137] ;

[0138] In the formula, This is the road surface roughness coefficient, in units of... ; The reference spatial frequency is 0.1, and the unit is 0.1. ; It is an exponential function. The formula relating pavement power spectral density and pavement roughness coefficient is expressed as follows:

[0139] ;

[0140] In the formula, Spatial frequency The power spectral density value of the road surface at the location, in units of ; Spatial frequency, unit: ; The reference spatial frequency is 0.1, and the unit is 0.1. .

[0141] The specific implementation of step S11 is to obtain a globally consistent trajectory estimate using a pose graph optimization algorithm. The state vector representation of the pose graph nodes is as follows:

[0142] ;

[0143] In the formula, For the first The state vector of each pose graph node; , , The first Spatial coordinates of each pose graph node, in units of ; , , The first The pose angle of each pose graph node, in units of ; superscript This represents the vector transpose. The Euclidean distance for loop closure detection is expressed as follows:

[0144] ;

[0145] In the formula, The Euclidean distance between the current pose graph node and the historical pose graph nodes, in units of 1. ; , , The spatial coordinates of the current pose graph node, in units of ; , , The spatial coordinates of the nodes in the historical pose graph, in units of .when The time is determined to be a loop, where For reference distance, the unit is The empirical value is 1, used for dimensionless processing; The preset loop closure detection threshold is in units of The value ranges from 1 to 5. The formula for calculating the weighted sum of squared residuals is as follows:

[0146] ;

[0147] In the formula, This is the weighted sum of squared residuals; It is the set of edge constraints of the pose graph; For the first The first pose graph node and the first The residual vector between each pose graph node; The information matrix has elements whose dimensions are the reciprocals of the dimensions of the corresponding elements in the residual vector. Weighting is performed using this matrix to... Dimensionless; superscript This represents the transpose of a vector. The iterative formula for the Levenberg-Marquardt optimization algorithm is as follows:

[0148] ;

[0149] In the formula, This is the increment of the state vector of the pose graph node; The Hessian matrix is ​​calculated using the following formula: ,in Let be the Jacobian matrix, representing the partial derivative of the residual vector with respect to the state vector; It is the damping factor; It is an identity matrix, and its dimensions are the same as those of the identity matrix. same; Let be the gradient vector, and its calculation formula is: ,in The set of all residual vectors; superscript This indicates the matrix transpose.

[0150] The specific implementation of step S12 is the same as described above, and will not be repeated in detail here. The power spectral density curves of standard Class E and standard Class F pavements are in the spatial frequency range of 0.011... Up to 2.83 The corresponding root mean square values ​​and geometric mean values ​​are respectively and The conversion formula between velocity power spectral density and displacement power spectral density is expressed as follows:

[0151] ;

[0152] In the formula, Spatial frequency The velocity power spectral density value at that location, in units of ; Spatial frequency The displacement power spectral density value at the location, in units of ; Spatial frequency, unit: ; Let π be the mathematical constant pi. The conversion formula between acceleration power spectral density and displacement power spectral density is expressed as follows:

[0153] ;

[0154] In the formula, Spatial frequency The acceleration power spectral density value at the location, in units of ; Spatial frequency The displacement power spectral density value at the location, in units of ; Spatial frequency, unit: .

[0155] To better understand and implement this invention, the following is a specific application scenario example 2: A technical team conducted a comprehensive inspection and evaluation of the road surface conditions of a contiguous farmland area. This farmland area covers approximately 180 hectares and has a complex internal road network, including various road surface types. Traditional inspection methods are inefficient and struggle to obtain continuous road surface spectrum data. The technical team decided to use a farmland road surface spectrum detection method based on a spherical mobile platform for on-site testing.

[0156] The technical team first installed a three-axis accelerometer, a two-axis tilt sensor, and a three-axis gyroscope on the inner platform of the spherical mobile platform. The accelerometer's measurement range was set to ±16g with a sampling frequency of 200Hz, and the gyroscope's measurement range was ±2000° / s. A pressure sensor with a range of 0 to 500N was installed at the bottom of the spherical mobile platform, also with a sampling frequency of 200Hz. A GNSS / BeiDou dual-mode positioning module was installed on the outside of the spherical mobile platform, achieving a positioning accuracy of 0.8m and an update frequency of 10Hz. A remote control module operating at a frequency of 2.4GHz, as well as temperature and humidity sensors, were also installed. The temperature sensor's measurement range was -40℃ to 85℃ with an accuracy of 0.5℃, and the humidity sensor's measurement range was 0% to 100% relative humidity with an accuracy of 3%. The technical team measured the position vector of the accelerometer relative to the center of mass of the spherical moving platform using mechanical design drawings. The vector was 0.05m in the longitudinal direction, 0.02m in the transverse direction, and 0.08m in the vertical direction. The total mass of the spherical moving platform was measured to be 15kg.

[0157] The ambient temperature on the day of the test was 28℃, the relative humidity was 65%, and the weather conditions were favorable for road surface inspection. The technical team activated the spherical mobile platform and controlled it via a remote control module to move it across the farmland surface to be tested. The travel speed was set between 1.2 m / s and 1.5 m / s to ensure uniform movement. The GNSS Beidou positioning module recorded the spatial coordinates of the spherical mobile platform in real time, updating the position information 10 times per second. The pressure sensor synchronously collected the pressure value along the longitudinal axis of the bottom of the spherical mobile platform. In the initial static state, the pressure value stabilized at around 147 N, corresponding to the pressure generated by the weight of the spherical mobile platform itself.

[0158] Accelerometers continuously collect raw triaxial acceleration data from the inner platform and simultaneously output quaternion attitude data for attitude calculation. Gyroscopes synchronously measure the triaxial angular velocities of the spherical moving platform. On smooth sections, the angular velocity amplitude is generally less than 10° / s, while on bumpy sections, the angular velocity fluctuation increases significantly, reaching a maximum of 45° / s. Tilt sensors measure the pitch and roll angles of the spherical moving platform. On smooth sections, the pitch angle remains within the range of -2° to +2°, and the roll angle remains within the range of -1.5° to +1.5°. Temperature and humidity sensors record ambient temperature and humidity values ​​every 5 seconds. During the test, the temperature remained stable between 27°C and 29°C, and the humidity remained stable between 63% and 67%.

[0159] The technical team set a preset motion state threshold of 5 N / s and calculated the absolute value of its time derivative based on the pressure value collected by the pressure sensor. When the absolute value of the time derivative of the pressure value is less than 5 N / s, the spherical moving platform is determined to be in a uniform motion state. During the entire test, the system identified 23 uniform motion state segments, each lasting between 3 and 12 seconds. At the beginning and end of each uniform motion state segment, the system automatically performed a zero-speed correction operation, forcing the accumulated velocity deviation and accumulated displacement deviation to zero, effectively suppressing the accumulation of integral drift error.

[0160] The technical team performed kinematic compensation processing on the raw longitudinal acceleration collected by the accelerometer. Using the triaxial angular velocity vectors measured by the gyroscope and the position vector of the accelerometer relative to the center of mass of the spherical moving platform, they calculated the centrifugal and tangential acceleration components. In a typical bumpy road test, the amplitude of the raw longitudinal acceleration ranged from -3.2g to +2.8g. After kinematic compensation, the amplitude of the centrifugal acceleration component was approximately 0.3g, and the amplitude of the tangential acceleration component was approximately 0.2g. Subtracting the centrifugal, tangential, and gravitational acceleration components (1g) from the raw longitudinal acceleration yielded the longitudinal linear acceleration, whose amplitude range was reduced to -2.1g to +1.9g, more realistically reflecting the acceleration changes caused by road surface excitation.

[0161] The technical team inputs the kinematically compensated longitudinal axis linear acceleration into a Kalman filter for noise reduction. The state-space model of the Kalman filter includes a one-dimensional state equation and an observation equation, with the process noise covariance set to 0.05 and the measurement noise covariance set to 0.1. The support force acceleration is obtained by dividing the pressure value collected by the pressure sensor by the mass of the spherical moving platform (15 kg), and then adding the gravitational acceleration component (9.8). The Kalman filter is used to fuse the vertical axis linear acceleration with the true vertical axis acceleration estimate, outputting the optimal vertical axis acceleration estimate sequence. After filtering, the high-frequency noise of the acceleration signal is significantly reduced, and the signal-to-noise ratio is improved by approximately 15 dB.

[0162] The technical team performed a two-step integration transformation on the optimal longitudinal acceleration estimate sequence output by the Kalman filter using the trapezoidal integration rule. The first integration yielded the longitudinal velocity sequence, and the second integration yielded the longitudinal displacement sequence. Since the integration process introduces low-frequency drift, the team applied a high-pass filter to the longitudinal displacement sequence to remove this drift. The high-pass filter was a Butterworth design, fourth order, with a cutoff frequency set to 0.05Hz. After high-pass filtering, the vibration displacement sequence of the spherical moving platform relative to its average longitudinal position was obtained, as shown below. Figure 1 As shown, the amplitude range of the vibration displacement is from -18mm to +22mm, which clearly reflects the vibration characteristics caused by road surface unevenness.

[0163] The technical team segmented the vibration displacement sequence into 512 data points of a preset segment length, with an overlap rate of 50% between adjacent segments, resulting in 89 vibration displacement data segments. A Hanning window function was applied to each segment for weighting, effectively reducing spectral leakage. The Welch method was used to perform a Fast Fourier Transform (FFT) on each weighted segment to calculate the power spectral density of the single-segment displacement. The FFT achieved a frequency resolution of 0.39 Hz, corresponding to a spatial frequency resolution of 0.26 Hz. The power spectral density results of individual displacement data segments from 89 vibration displacement data points were averaged to obtain a smoothed displacement power spectral density curve. This smoothed displacement power spectral density curve operates within a spatial frequency range of 0.05. Up to 3.0 It exhibits a clear power-law decay characteristic.

[0164] The technical team performed dynamic time warping matching between the smooth displacement power spectral density curve and the power spectral density curves of standard pavements in grades A to H of the standard pavement roughness classification table. By constructing a cumulative distance matrix and using dynamic programming to search for the optimal path, the dynamic time warping distance between the smooth displacement power spectral density curve and each standard pavement power spectral density curve was calculated. The calculation results are shown in Table 1. The minimum dynamic time warping distance is 0.73, corresponding to the standard grade F pavement power spectral density curve. Therefore, the pavement roughness grade of the farmland pavement under test was determined to be grade F.

[0165] Table 1 Calculation results of dynamic time warping distance

[0166]

[0167] The technical team performed linear fitting based on the logarithmic coordinate relationship between the smooth displacement power spectral density curve and the spatial frequency. After taking the logarithm of both the spatial frequency (horizontal axis) and power spectral density values ​​(vertical axis), a scatter plot was plotted, and the least squares method was used for linear fitting. The slope of the fitted line was -1.98; taking the negative value yielded a frequency exponent of 1.98, close to the typical value of 2 for the frequency exponent in the standard pavement model. The reference spatial frequency was 0.1. Taking the logarithm and substituting it into the fitted linear equation, the calculated logarithmic power spectral density is -4.18. Taking the exponential operation on this value yields the road surface roughness coefficient. This value falls within the range of pavement roughness coefficient for standard Class F pavement. to This further verified the accuracy of the road surface grade determination results.

[0168] The technical team constructed a pose graph using spatial coordinate sequences recorded by a GNSS BeiDou positioning module. A preset time interval threshold of 1 second and a preset spatial interval threshold of 1.2 meters were set, and moments meeting either threshold were identified as critical moments. A total of 486 critical moments were identified during the entire test, and the spatial coordinates of the spherical mobile platform at these critical moments were used as pose graph nodes. The relative motion transformations between adjacent pose graph nodes were used as pose graph edge constraints, each containing relative displacement and relative attitude change information. During the test, the spherical mobile platform traveled along a network of farmland roads and repeatedly passed the same location. The system set a preset loop closure detection threshold of 2.5 meters. When the Euclidean distance between the current spatial coordinates and the spatial coordinates of historical pose graph nodes was less than 2.5 meters, it was determined to return to the visited spatial location, and a loop closure constraint edge was added between the initial pose graph node and the returned pose graph node. The system detected a total of 17 valid loop closures, such as... Figure 2 As shown, the addition of loop constraint edges significantly improves the global consistency of trajectory estimation.

[0169] The technical team minimized the weighted sum of squared residuals of all pose graph edge constraints using the Levenberg-Marquardt optimization algorithm. The initial damping factor of the optimization algorithm was set to 0.01, and it was adaptively adjusted in each iteration based on changes in the objective function. After 23 iterations, the weighted sum of squared residuals converged to a stable value, and the maximum positional deviation of the trajectory before and after optimization reached 3.7m, resulting in a globally consistent trajectory estimate. This globally consistent trajectory estimate eliminated the cumulative drift accumulated over long periods, ensuring consistency of the entire trajectory in the global coordinate system. The positional error at the trajectory loop closure point decreased from 3.2m before optimization to 0.6m after optimization, significantly improving the accuracy of road surface spectrum spatial positioning.

[0170] The internal structure of the small spherical mobile platform and the structure of various sensors in this embodiment are as follows: Figure 3 and Figure 4 As shown, the outer spherical shell consists of two perfectly mating hemispheres. The part in contact with the ground is covered with an anti-slip module with a higher coefficient of friction, designed with a grooved structure like a tire. The interior of the shell is designed with treads that fit the transmission device, ensuring smooth rolling. The motion parameters of the outer spherical shell traveling over the test surface are collected and processed by internal sensors and converted into shape parameters of the test surface. Combined with positioning and communication functions, a large-scale, intuitive representation of the road surface is generated in real time from a distance.

[0171] The technical team determined the road surface roughness level to be F and the road surface roughness coefficient to be... The frequency index (1.98), the smooth displacement power spectral density curve, and the spatial coordinates in the globally consistent trajectory estimation were correlated. A pavement spectral image was plotted, with spatial frequency on the horizontal axis and power spectral density on the vertical axis, clearly showing the vibration intensity distribution characteristics of the farmland pavement under test at different wavelengths. The ambient temperature (28℃) and humidity (65%) were displayed next to the pavement spectral image, and the power spectral density curves of standard Class E and Class F pavements were overlaid as reference benchmarks. Comparative analysis showed that the measured power spectral density curve of the farmland pavement under test was slightly higher than that of the standard Class F pavement in the low-frequency range and slightly lower in the high-frequency range, indicating that the pavement fluctuates significantly in the long wavelength range but is relatively smooth in the short wavelength range.

[0172] The technical team completed a comprehensive inspection of the road network across 180 hectares of farmland, accumulating a test mileage of 42 kilometers and collecting 84 million valid data points, identifying areas with different road surface grades. Test results showed that approximately 35% of the roads in the farmland area were Grade E, approximately 52% were Grade F, and approximately 13% were Grade G. This road surface spectrum data provides crucial information for optimizing the parameters of the adaptive suspension system for agricultural machinery. This allows the suspension system to dynamically adjust damping and stiffness coefficients based on real-time road conditions, significantly improving the ride comfort and operational stability of the agricultural machinery. Simultaneously, the precise correlation between the road surface spectrum data and spatial coordinates provides high-quality input for intelligent operation path planning algorithms, enabling the path planning system to comprehensively consider multiple objectives such as road surface smoothness, ride comfort, and operational efficiency, automatically generating the optimal operation path.

[0173] It should be noted that the variables involved in this invention are explained in detail in Tables 2, 3, and 4.

[0174] Table 2. Variable Explanation Table (Part 1)

[0175]

[0176] Table 3. Variable Explanation Table (Part Two)

[0177]

[0178] Table 4. Variable Explanation Table (Part 3)

[0179]

[0180] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for detecting the pavement spectrum of farmland based on a spherical moving platform, characterized in that, An accelerometer, tilt sensor, gyroscope, pressure sensor, GNSS / BeiDou positioning module, remote control module, temperature sensor, and humidity sensor are installed on a spherical mobile platform. The spherical mobile platform is controlled by the remote control module to move along the farmland surface under test and collect data synchronously. Kinematic compensation is performed on the raw longitudinal axis acceleration collected by the accelerometer to obtain the longitudinal axis linear acceleration. The longitudinal axis linear acceleration is input into a Kalman filter for noise reduction to obtain the optimal longitudinal axis acceleration estimate sequence. The optimal longitudinal axis acceleration estimate sequence is integrated twice and a high-pass filter is applied to obtain the vibration displacement sequence. The vibration displacement sequence is then... The smooth displacement power spectral density curve is obtained by segmented weighting and Welch method calculation; the smooth displacement power spectral density curve and the standard road surface power spectral density curve are dynamically time-warped to obtain the road surface roughness level; the frequency index and road surface roughness coefficient are obtained by linear fitting based on the double logarithmic coordinate relationship between the smooth displacement power spectral density curve and the spatial frequency; the pose map is constructed using the spatial position coordinate sequence recorded by the GNSS Beidou positioning module to obtain the globally consistent trajectory estimation; the road surface roughness level, road surface roughness coefficient, frequency index, smooth displacement power spectral density curve and spatial position coordinates are correlated to draw the road surface spectral image.

2. The method according to claim 1, characterized in that, The accelerometer, tilt sensor, and gyroscope are installed on the inner platform of the spherical mobile platform, the pressure sensor is installed at the bottom of the spherical mobile platform, and the GNSS Beidou positioning module, remote control module, temperature sensor, and humidity sensor are installed on the outside of the spherical mobile platform.

3. The method according to claim 2, characterized in that, The synchronously acquired data includes the GNSS Beidou positioning module recording the spatial coordinates of the spherical mobile platform in real time, the pressure sensor collecting the pressure value along the longitudinal axis of the bottom of the spherical mobile platform, the accelerometer collecting the three-axis raw acceleration and quaternion attitude data of the inner platform, the gyroscope measuring the three-axis angular velocity of the spherical mobile platform, the tilt sensor measuring the pitch and roll angles of the spherical mobile platform, and the temperature and humidity sensors recording the ambient temperature and humidity values.

4. The method according to claim 3, characterized in that, Before the kinematic compensation, the motion state of the spherical moving platform is determined based on the absolute value of the time derivative of the pressure value collected by the pressure sensor and the sensor data stream is segmented. When the absolute value of the time derivative of the pressure value is less than the preset motion state threshold, it is determined to be a uniform motion state. Zero speed correction is performed at the beginning and end of the uniform motion state segment to force the cumulative velocity deviation and cumulative displacement deviation to zero.

5. The method according to claim 4, characterized in that, The preset motion state threshold ∈ [1, 10] N / s, and the sensor data stream consists of the triaxial raw acceleration collected by the accelerometer, the triaxial angular velocity measured by the gyroscope, and the pitch and roll angles measured by the tilt sensor.

6. The method according to claim 5, characterized in that, The kinematic compensation specifically involves using the three-axis angular velocity vector measured by the gyroscope and the position vector of the accelerometer relative to the center of mass of the spherical moving platform to calculate the centrifugal acceleration component and the tangential acceleration component. The longitudinal axis linear acceleration is obtained by subtracting the centrifugal acceleration component, the tangential acceleration component, and the gravitational acceleration component from the original longitudinal axis acceleration.

7. The method according to claim 6, characterized in that, The calculation of the centrifugal acceleration component specifically involves performing a cross product operation between the three-axis angular velocity vector and the position vector to obtain a first intermediate vector, and then performing a cross product operation between the three-axis angular velocity vector and the first intermediate vector. The negative value of the result is the centrifugal acceleration component vector.

8. The method according to claim 7, characterized in that, The calculation of the tangential acceleration component specifically involves taking the time derivative of the three-axis angular velocity vector to obtain the angular acceleration vector, and then performing a cross product operation between the angular acceleration vector and the position vector to obtain the tangential acceleration component vector.

9. The method according to claim 8, characterized in that, The Kalman filter performs noise reduction processing. Specifically, it uses the pressure value collected by the pressure sensor divided by the mass of the spherical moving platform plus the gravitational acceleration component as the true value of the longitudinal axis acceleration estimation. The Kalman filter then fuses the true value of the longitudinal axis acceleration estimation with the longitudinal axis linear acceleration to output the optimal longitudinal axis acceleration estimation value.

10. The method according to claim 9, characterized in that, The two cumulative integrations are performed using the trapezoidal integration rule. The high-pass filter is used to remove the low-frequency drift trend of the vertical axis displacement sequence. The cutoff frequency of the high-pass filter is ∈ [0.01, 0.1] Hz.

Images