Electromagnetic-mechanical coupling analysis and optimization method of electromagnetic vibroseis

Electromagnetic vibroseis is an important seismic wave excitation device. At present, the research of electromagnetic vibroseis mainly focuses on application and prototype design, lacking research on electromagnetic-mechanical coupling. Analyzing the electromagnetic-mechanical coupling and optimizing the air gap magnetic field is of great significance to improve the excitation quality of low-frequency seismic waves. In this paper, the electromagnetic-mechanical coupling is analyzed by studying the parameter force factor, and two novel optimal designs of the air gap magnetic field are proposed. Firstly, the finite element model of the air gap magnetic field is established, and the influence of the force factor distribution curve on the vibration signal is analyzed by using the Gaussian function and parameter spline difference method. Further, an asymmetrical axial dual-magnet design is proposed to enhance the radial magnetic induction intensity Br in the air gap. The design of adding a radial magnet to the outer yoke is proposed to compensate for the nonlinearity of the Br in the air gap. The simulation results show that the asymmetric axial dual-magnet structure increases the Br by 22.2% compared with the axial single-magnet structure. Adding a radial magnet to the outer yoke reduces the amplitude ratio of the third harmonic to the fundamental wave from 23.24% to 4.66%. It is necessary to consider the influence of the height of the driving coils on the maximum displacement, force factor distribution curve, and harmonic distortion.


Introduction
Seismic exploration uses artificially excited seismic waves to infer the properties and morphology of underground medium, which plays an important role in the exploration * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. of oil, gas, and other mineral resources. Vibroseis radiates frequency-modulated sinusoidal sweep signals to the earth through continuous vibration, which has the characteristics of high efficiency, good repeatability, environmental protection, and safety. Land seismic exploration benefits from vibroseis technology [1][2][3]. Electromagnetic vibroseis is a kind of vibroseis driven by electromagnetic force. Land electromagnetic vibroseis is mainly used for high-resolution seismic exploration, including tunnel detection, mineral exploration, and detection of urban underground active structures [4][5][6]. Much research effort has been applied to the design of land electromagnetic vibroseis prototype. It is a common design method to build portable electromagnetic vibroseis prototypes based on commercial audio components, tactile sensors, or actuators as the vibration unit [7][8][9]. This design method can effectively reduce the cost. For example, Dean et al [9] constructed a simple portable electromagnetic vibroseis, which consists of four low-frequency actuators. However, the performance of electromagnetic vibroseis is limited by the vibration unit. The method of independently designing the vibration unit can flexibly adjust the performance of the electromagnetic vibroseis [10][11][12]. For example, electromagnetic vibroseis driven by linear synchronous motors consists of a U-shaped magnet track and a coil track and has been tested for performance and used in seismic exploration [13,14]. In recent years, the application of electromagnetic vibroseis is expanding to marine seismic exploration [15][16][17][18][19][20]. The research of marine electromagnetic vibroseis mainly focuses on prototype design and performance tests. The vibration core unit of the marine electromagnetic vibroseis is independently designed.
To conclude, the design of electromagnetic vibroseis is an important work. In electromagnetic vibroseis, the electromagnetic-mechanical coupling will affect the quality of low-frequency vibration signals. The low-frequency vibration signals could help to improve the depth of seismic exploration and the accuracy in the inversion of seismic data to acoustic impedance [21,22]. Therefore, it is important to analyze and suppress the nonlinearity of electromagneticmechanical coupling. However, there is a lack of research on electromagnetic-mechanical coupling in electromagnetic vibroseis. In other electromagnetic-mechanical fields. Klippel [23] studied the electromagnetic-mechanical coupling of loudspeakers. Saraswat and Tiwari [24] studied the electromagnetic-mechanical coupling in electrodynamic shakers. Kim et al [25] introduced an electrical, magnetic, and mechanical coupling analysis method in vibrating motors. Jiang et al [26] used the electromagnetic-mechanical-acoustic coupling method and resonance equations to study the sound pressure level and peak frequencies in micro-speakers.
Analysis and optimization of the air gap magnetic field is an important method to suppress the nonlinearity of electromagnetic-mechanical coupling. Merit et al [27] proposed a variety of iron-free magnetic circuit structures to improve the nonlinearity of air gap magnetic fields in loudspeakers. He et al [28] proposed that the uneven air gap structure can improve the magnetic induction intensity distribution of the air gap magnetic field in the horizontal electromagnetic vibration exciter. Sun et al [29] improved the performance of a permanent magnet synchronous motor by optimizing the air gap magnetic field. Qian et al [30] established the air gap magnetic field model of the permanent magnet synchronous linear motor and analyzed the relationship between the air gap magnetic field and the geometric size of the permanent magnet. Based on the magnetic equivalent circuit method, Ma et al [31] proposed a magnetic field analysis model for permanent magnet synchronous motors. Sun et al [32] improved the air gap magnetic field waveform in Hollow-cup motors through the eccentricity of the inner rotor.
In this paper, the electromagnetic-mechanical coupling of the electromagnetic vibroseis and its influence on the vibration signal is analyzed, and two novel methods for optimizing the air gap magnetic field are proposed. We study the electromagnetic-mechanical coupling of electromagnetic vibroseis by analyzing the parameter force factor. Firstly, the finite element model of the air gap magnetic field is established, and the B r and force factor distribution curve are calculated. Gaussian function and parameter spline difference (PSD) method are used to analyze the influence of force factor distribution curve on low-frequency vibration signal. To enhance the B r in the air gap, an asymmetric axial dualmagnet structure design is proposed. The design of adding a radial magnet to the outer yoke is proposed to compensate for the B r in the air gap region between the axial magnet and the outer yoke. The finite element simulation results show that the two methods proposed in this paper can effectively suppress the nonlinearity of the electromagnetic-mechanical coupling.
The rest of this paper is organized as follows. Section 2 introduces the electromagnetic vibroseis structure and electromagnetic-mechanical coupling. In section 3, a finite element model of the air gap magnetic field is established, and the influence of the nonlinear force factor on the vibration signal is analyzed. In section 4, two air gap magnetic field optimization models are proposed and established. Section 5 discusses the influence of the driving coils on the electromagnetic-mechanical coupling. Finally, section 6 outlines a conclusion.

Electromagnetic vibroseis structure
Electromagnetic vibroseis excites vibration signals through the electromagnetic force on the energized driving coils in the magnetic field, which is generally a magnet-driving coils structure. The structure of a land portable electromagnetic vibroseis is shown in figure 1, including the magnet, yoke, driving coils, guide rod, support spring, baseplate, etc. The axial magnet is used to excite the air gap magnetic field in land portable electromagnetic vibroseis. The yokes have high magnetic permeability, and the magnetic resistance can be approximately ignored. The driving coils made of enameled wire move in the air gap magnetic field. A guide rod is used to fix the driving coils to the baseplate. The magnet and yokes constitute the moving parts of the electromagnetic vibroseis, also known as the reaction mass. The support spring is used to balance the gravity of the reaction mass. The displacement of the reaction mass is limited by the height of the air gap and the height of the driving coils. The lowest operating frequency of the electromagnetic vibroseis is closely related to the maximum displacement.
The magnetic field structure of the electromagnetic vibroseis is generally cylindrical. In figure 1, a cylindrical coordinate system is established. The magnetic induction intensity B in the air gap can be expressed by B r , B z , B θ . When the alternating current i flows through the driving coils, only B r contributes to the vibration of the electromagnetic vibroseis. The electromagnetic force F z is expressed as follows: where r coils , h coils , and N are height, winding radius, and turns of the driving coils, respectively. B r,z is the radial magnetic induction intensity at coordinates are(r coils , z). L coils is the length of driving coils. B r,z is the average value of the B r,z integrated along the height of the driving coils.

Electromagnetic-mechanical coupling
The electromagnetic-mechanical coupling analysis mainly includes establishing the air gap magnetic field model, measuring the B r in the air gap, calculating the force factor distribution curve, fitting the force factor distribution curve, solving the electromagnetic-mechanical coupling motion equation, and analyzing the frequency domain acceleration signal. The flow chart of electromagnetic-mechanical coupling analysis of electromagnetic vibroseis is shown in figure 2. Firstly, the air gap magnetic field model of the electromagnetic vibroseis needs to be established. The air gap magnetic field satisfies the following formulas: where H and M are magnetic field intensity, and magnetization respectively. ∇× and ∇· are the curl and divergence operators of the vector. µ 0 is vacuum permeability. The method of solving the air gap magnetic field using Poisson equation or Laplace equation requires dividing the air domain, as shown in figure 3(a). Poisson equation based on magnetic vector potential A is:  where µ, J,and ∇ 2 are the magnetic permeability, current density and Laplace operator, respectively. The Poisson equation or Laplace equation requires boundary conditions that constrain regions I, II, and III to determine a unique solution. The boundaries of different regions are shown in the red line in figure 3(a). Figure 3(b) shows a lumped parameter magnetic circuit model of the air gap magnetic field, which can be solved by Kirchhoff law. However, factors such as magnetic flux leakage will bring additional difficulties to solving the model. With the rapid development of computer technology, the finite element method is widely used in engineering for magnetic field analysis. The finite element method provides numerical solutions of magnetic induction intensity. Therefore, we solve the air gap magnetic field based on the finite element method. Based on the established finite element model of the air gap magnetic field, the B r in the air gap is measured. The force factor BL is used to describe the electromagneticmechanical coupling and is defined as: The BL(x) is related to the displacement x of the driving coils in the air gap. When the electromagnetic vibroseis vibrates at high frequency, the displacement of the driving coils is very small, and the force factor is approximately constant. The displacement of the driving coils increases at low frequencies. The nonlinearity of the BL(x) is more obvious.
To establish the electromagnetic-mechanical coupling motion equation, it is necessary to function the BL(x). The BL(x) can be expressed in the form of a Gaussian function [24]: where b 1 and b 2 are the peak heights of the Gaussian function. c 1 and c 2 are the peak positions of the Gaussian function. d 1 and d 2 are the region widths of the Gaussian function. Based on the BL(x) after Gaussian function fitting, the motion equation of the electromagnetic vibroseis is established. The electromagnetic vibroseis system is simplified to a secondorder spring-mass-damper system. The motion equation of the reaction mass is as follows: where x,ẋ, andẍ represent the displacement, velocity, and acceleration of the reaction mass, respectively. M, D, and K are the mass, damping coefficient, and elastic coefficient of the electromagnetic vibroseis system, respectively. The motion equation of the reaction mass is transformed into a coupled nonlinear differential equation and solved by the PSD method [33]. The discrete-time iteration forms of reaction mass motion based on the PSD method are: where n is used to denote the time step, f represents the splinebased function, αand β are parameters of the spline-based function and ∆t is the discrete time interval. The sum of α and β is 0.5. If {α, β} is set to {0, 1/2}, the PSD method degenerates into a traditional finite-difference method. When {α, β} is set to {1/12, 5/12}, the error of the PSD method is much smaller than that of the traditional finite difference method. The spline-based function f is as follows: Through the iteration of the time step, the displacement, velocity, and acceleration of the reaction mass motion can be obtained. The influence of the electromagnetic-mechanical coupling nonlinearity on the vibration signal is obtained by analyzing the acceleration frequency domain signal.

Air gap magnetic field and vibration nonlinearity
The AC/DC module in COMSOL Multiphysics software is used to analyze the air gap magnetic field of the electromagnetic vibroseis. The steps of the finite element analysis method include establishing the geometric model, adding materials, meshing, calculating the magnetic field, and post-processing. The finite element model of the air gap magnetic field is a simplification of the electromagnetic vibroseis structure in figure 1, which only retains the axial magnet and yokes. Considering the axisymmetric nature of the air gap magnetic field, a 2D axisymmetric geometric model is selected in the AC/DC module, which can reduce the solution time. Both the axial magnet and the yokes are made of lossless soft iron material. The height and radius of the axial magnet are 20 mm and 60 mm, respectively. The magnetization direction of the axial magnet is the +z axis direction, and the residual magnetic flux density is 1.22 T. The height of the bottom yoke is 30 mm. The height of the inner yoke is 22 mm and the radius of the inner yoke is the same as that of the axial magnet. The thickness of the outer yoke is 28 mm. The height of the outer yoke is equal to the sum of the axial magnet height and the inner yoke height. The relative magnetic permeability of the yokes is set to 500. The air wraps the axial magnet and the yokes. The radius and height of the air domain are both 0.4 m. The air gap width is set to 4 mm. The finite element mesh is a free triangular mesh, and the maximum element mesh size is 1 mm. The simulation results of the air gap magnetic field model are shown in figure 4. In figure 4(a), the air gap region can be divided into three parts: region I between the axial magnet and the outer yoke, region II between the inner yoke and the outer yoke, and region III in the air domain. The red arrows indicate the B. A measurement path is set up to obtain the distribution curve of    In air gap regions I and III, the nonlinearity of the B r is obvious. In the air gap region II, the B r is maximum and stable. Due to the difference in magnetic permeability between region I and region III, the distribution of the B r in the air gap is not symmetrical. Figure 5(b) shows the distribution curve of the B z on the measurement path. The distribution curve of the B z is more complicated. Since the geometric model of the air gap magnetic field is cylindrical, B θ is 0.
Before calculating the force factor, it is necessary to determine the parameter information of the driving coils. The height and winding radius of the driving coils are set to 40 mm and 62 mm, respectively. The maximum displacement of the driving coils is set to ±10 mm. The turn number of the driving coils is set to 54. The force factor distribution curve is shown in figure 6. The force factor distribution curve has obvious nonlinearity. Set the +z direction as the positive direction of driving coils displacement. When the displacement of the driving coils is 0, the force factor is maximum. As the displacement of the driving coils increases, the force factor decreases gradually. When the driving coils displacement is −0.01 m and 0.01 m, the force factor is 90.84% and 87.01% of the maximum value, respectively. Therefore, the force factor distribution curve is asymmetric with the displacement of the driving coils. This is primarily caused by the asymmetry of the B r distribution curve. The green dotted line is the force factor distribution curve after Gaussian function fitting. The root mean square error of the fitting function is only 0.0126.
Before analyzing the influence of the force factor on the vibration signal, it is necessary to determine the parameters of the electromagnetic vibroseis. The spring coefficient is 50 000 N m −1 . The mass of the reaction mass is 57 kg. The damping coefficient is set as 1140 N m −1 s −1 . The response of the electromagnetic vibroseis system is shown in figure 7. The resonance frequency of the system is about 4.75 Hz. Figures 7(a) and (b) show the displacement and acceleration responses of the reaction mass, respectively. The movement of the reaction mass can be divided into two phases. When the vibration frequency is lower than the resonance frequency, the displacement response of the reaction mass is approximately constant. When the vibration frequency is greater than the resonance frequency, the acceleration response of the reaction mass is approximately constant.
The influence of force factor nonlinearity on the vibration signal is analyzed based on the PSD method. The force factor nonlinearity mainly affects low-frequency vibration signals. The lowest vibration frequency of the electromagnetic vibroseis is set to 2 Hz. The parameter information in the solution process of the PSD method is shown in figure 8. Figure 8(a) shows the alternating current flowing through the driving coils. The peak value of the alternating current is 35 A. Amplitude tapering is performed at the beginning and end of the alternating current, which can effectively suppress the Gibbs ringing effect and constrain the displacement of the reaction mass. The taper function selects the Hanning window, and the taper ratio is 30%. Figure 8(b) shows the electromagnetic force on the driving coils. The peak value of electromagnetic force is about 430 N. Figure 8(c) shows the fluctuation curve of the force factor with time. The force factor varies with the displacement of the driving coils. Figure 8  coincident. There is an obvious third harmonic in the acceleration signal. The amplitude ratio of the third harmonic to the fundamental wave is about 23.24%. Harmonics will reduce the fundamental wave energy of the electromagnetic vibroseis and cause interference to the signal recorded by the geophone.
The second and third harmonics of the acceleration at different vibration frequencies are analyzed. The vibration frequencies are set to 2 Hz, 3 Hz, and 4 Hz, respectively. When the vibration frequency is 2 Hz, 3 Hz, and 4 Hz, the maximum displacement of the reaction mass is set to 10 mm. The amplitude ratios of the multiple harmonics to the fundamental wave are shown in figure 9. The influence of the second harmonic on the fundamental wave is significantly less than the influence of the third harmonic on the fundamental wave. This is mainly because the asymmetry of the force factor leads to the second harmonic of the vibration signal, and the nonlinearity of the force factor leads to the third harmonic of the vibration signal. Compared with asymmetry, the nonlinearity of the force factor is more obvious. When the vibration frequency is 3 Hz and 4 Hz, the amplitude ratio of the third harmonic to the fundamental wave is almost the same. They are all significantly smaller than the amplitude ratio of the third harmonic to the fundamental when the vibration frequency is 2 Hz. This is due to the greater acceleration response near the resonant frequency. When the vibration frequency is 2 Hz, the second and third harmonic frequencies are 4 Hz and 6 Hz respectively, which are closer to the resonance frequency.

Optimization method of the air gap magnetic field
The force factor nonlinearity is mainly due to the uneven distribution of the B r in the air gap. Therefore, optimizing the air gap magnetic field can suppress the nonlinearity of the force factor. Analyzing the influencing factors is necessary to optimize the air gap magnetic field. The main parameter information of the air gap magnetic field includes the height and radius of the axial magnet, inner yoke height, air gap width, outer yoke thickness, and bottom yoke height. Based on the control variable method, a linear sweep is performed on every single variable. Figure 10(a) shows the effect of the axial magnet height on the B r in the air gap magnetic field. When the height of the axial magnet increases from 10 mm to 20 mm, the B r increases by 0.12 T. When the height of the axial magnet increases from 20 mm to 60 mm, the B r increases only by 0.05 T. As the height of the axial magnet increases, the enhancement effect of the B r decreases gradually. Increasing the height of the axial magnet not only increases the magnetomotive force but also increases the magnetic reluctance of the entire magnetic circuit. The increase of magnetic reluctance will weaken the enhancement effect of axial magnet height on the B r . Figure 10(b) shows the effect of the axial magnet radius on the B r . The magnetic reluctance of an axial magnet decreases as the radius of the axial magnet increases. Thus, the magnetic flux of the entire magnetic circuit increases. The increase of magnetic flux significantly increases the B r in the air gap. Figure 10(c) shows the effect of the inner yoke height on the B r . As the height of the inner yoke increases, the surface area of the air gap will increase. This reduces the magnetic flux density (B r ) in the air gap. Figure 10(d) shows the effect of the air gap width on the B r . As the air gap width increases, the air gap magnetic reluctance also increases. The increase of air gap magnetic reluctance reduces the B r . In figures 10(e) and (f), the height of the bottom yoke and the thickness of the outer yoke have little effect on the B r .

Asymmetric axial dual-magnet structure
In figure 5(a), the B r is evenly distributed in air gap region II (between the inner yoke and the outer yoke). Therefore, making the driving coils only move in region II can reduce force factor nonlinearity. This requires that the height of the air gap region II should be the sum of the height of the driving coils and twice the maximum displacement. Figure 11(a) shows the air gap magnetic field structure for an axial singlemagnet. As the height of the inner yoke increases, the B r in the air gap decreases gradually. To reduce the effect of increasing the height of the inner yoke on the B r , an air gap magnetic field structure for asymmetric axial dual-magnet is proposed, as shown in figure 11(b). The sum of the heights of the upper and lower axial magnets in figure 11(b) is equal to the height of a single axial magnet in figure 11(a). The magnetizing directions of the upper and lower axial magnets are indicated by red arrows. The heights of the upper and lower axial magnets are determined by linear sweep. Figure 12 shows the simulation results for two air gap magnetic field optimization schemes. The heights of the inner and outer yokes in figure 12 are both 62 mm. The width of the air gap, the thickness of the outer yoke, and the height of the bottom yoke are consistent with those in figure 4(a). Compared with the axial single-magnet structure, the asymmetric dual-magnet structure has two magnetic flux loops. Figure 13(a) shows the maximum B r in the air gap for different lower axial magnet heights when the sum of the upper and lower axial magnet heights is 20 mm. It can be found that when the height of the lower axial magnet is 9 mm, the B r in the air gap is the largest. Figure 13(b) shows the maximum B r in the air gap for different lower axial magnet heights when the sum of the upper and lower axial magnet heights is 30 mm. When the height of the lower axial magnet is 11 mm, the B r in the air gap is the largest. The asymmetry of the upper and lower axial magnet height is more obvious.
In the asymmetric axial dual-magnet structure model, the heights of the lower and upper axial magnets are set to 9 mm and 11 mm, respectively. The B r in the axial single-magnet and asymmetric axial dual-magnet models is shown in figure 14.
In the axial single-magnet model, the maximum B r in the air gap is about 0.472 T. In the asymmetric axial dual-magnet model, the maximum B r in the air gap is about 0.577 T. Compared with the axial single-magnet design, the asymmetric axial dual-magnet design increases the maximum B r by about 22.2%. Moreover, in the axial single-magnet model, as the height of the inner yoke increases, the nonlinearity of the B r in air gap region II also increases. However, this phenomenon in the asymmetric axial dual-magnet model is suppressed. Figure 15 shows the force factor distribution curves of the axial single-magnet and asymmetric axial dual-magnet models. The driving coils information is consistent with that in figure 6. In figure 15(a), the force factor is approximately constant within ±10 mm. The maximum force factor of the axial single-magnet model is 9.66 T m. The maximum force factor of the asymmetric axial dual-magnet model is 11.96 T m. Compared with the axial single-magnet design, the asymmetric axial dual-magnet design increases the maximum force factor by about 23.81%. Figure 15(b) shows the normalized force factor distribution curve. When the displacement range is from −20 mm to 12 mm, the normalized force factor distribution curves of the axial single-magnet and asymmetric axial dual-magnet models are approximately coincident. When the displacement is greater than 12 mm, the nonlinearity of the force factor distribution curve in the asymmetric axial dualmagnet model is more obvious. This is mainly because the upper axial magnet enhances the nonlinearity of the B r in the air gap region III. Therefore, the asymmetric axial dualmagnet model needs to ensure the movement of the driving coils in the air gap region II. Otherwise, it may cause more severe distortion of the low-frequency vibration signal than the axial single-magnet model.

Add a radial magnet to the outer yoke
Limiting the movement of the driving coils in the air gap region II can reduce the nonlinearity of the force factor distribution curve. Compensating the B r in air gap region I also reduces the nonlinearity of the force factor. We proposed an optimization scheme for the air gap magnetic field with a radial magnet added to the outer yoke, as shown in figure 16. Adding a radial magnet can compensate the B r in region I. Figure 17 shows the simulation results for the scheme of adding a radial magnet to the outer yoke. Different from the air gap magnetic field model in figure 11(a), the inner yoke height of the air gap magnetic field model in figure 16 is 42 mm. This is because compensating the nonlinearity of the B r in the region I can reduce the height of the inner yoke that needs to be increased.      The parameter information of the radial magnet includes residual magnetic flux density, height and thickness. The   residual magnetic flux density of the radial magnet is consistent with the residual magnetic flux density of the axial magnet.
To determine the optimal parameter information of the radial magnet, the height and thickness of the radial magnet are linearly swept. Figures 18(a) and (b) show the simulation results when the axial magnet height is 20 mm. Figures 18(c) and (d) show the simulation results when the axial magnet height is 30 mm. When the height of the radial magnet is higher than that of the axial magnet, the B r will have an obvious peak. When the height of the axial magnet is higher than that of the radial magnet, the B r will appear obvious depression. The greater the difference between the height of the radial magnet and the height of the axial magnet, the more obvious this phenomenon will be. Therefore, the height of the radial magnet should be consistent with the height of the axial magnet. When the height of the radial magnet is constant, the thickness of the radial magnet will also affect the distribution of the  B r in the air gap. As the thickness increases, the B r in region I increases, while that in region II decreases. Take the optimization model when the axial magnet height is 20 mm as an example. Figure 19 shows the force factor distribution curves for different radial magnet thickness models. It can be found that when the thickness of the radial magnet is 10 mm, the symmetry of the force factor distribution curve is better. Figure 20(a) shows the B r distribution in different air gap magnetic field models. The axial magnet height is 20 mm, and the inner yoke height is 22 mm, 42 mm, and 62 mm, respectively. As the height of the inner yoke increases, the B r decreases gradually. The height of the air gap region II increases with the height of the inner yoke. The purple curve represents the B r in the optimization model with a radial magnet added to the outer yoke when the height of the inner yoke is 42 mm. The height and thickness of the radial magnet are 20 mm and 10 mm, respectively. Although the addition of a radial magnet will reduce the B r in the region II, it can greatly increase the B r in the region I. Figure 20(b) shows the force factor distribution curves in different air gap magnetic field models. The height of the driving coils is 40 mm and the number of turns is 54. As the height of the inner yoke increases, the force factor decreases gradually. However, as the height of the inner yoke increases, the nonlinearity of the force factor first increases and then decreases. When the inner yoke height is 42 mm, the force factor nonlinearity in the air gap magnetic field model with a radial magnet added is smaller. The force factor when the displacement is 10 mm is about 97.59% of the force factor when the displacement is zero. Compared with the scheme of simply increasing the height of the inner yoke, the optimized scheme of adding a radial magnet increases the force factor and reduces the need to increase the height of the inner yoke. A Gaussian function is used to fit the force factor distribution curve for the optimization scheme of adding a radial magnet. The vibration frequency of the electromagnetic vibroseis is set at 2 Hz. The maximum displacement of the driving coils is ±10 mm. Based on the PSD method, the differential equation of the electromagnetic-mechanical coupling process is solved. The parameter information in the solution process is shown in figure 21. In figure 21(a), the force factor still fluctuates. Compared with figure 8(c), the fluctuation of the force factor is much smaller. Figures 21(c) and (d) show the time-domain and frequency-domain signals of the reaction mass acceleration. The amplitude ratio of the third harmonic to the fundamental wave decreased from 23.24% before optimization to 4.66% after optimization.

Discussion
The electromagnetic-mechanical coupling is not only affected by the B r in the air gap, but also by the height of the driving coils. According to formulas 1 and 6, the height of the driving coils will affect the maximum displacement, force factor, and harmonic distortion of the acceleration signal. Taking the air gap magnetic field model in figure 4 as an example, the influence of the driving coils height on the electromagneticmechanical coupling is discussed. Dashed lines of different colors represent driving coils of different heights, as shown in figure 22. The driving coils of different heights are symmetrical about the center of the inner yoke. The driving coils height is linear with the number of turns. For the convenience of calculation, when the height of the driving coils is 10 mm, the number of turns is 10. The driving coils heights are set to   Figure 23(a) shows the relationship between the height of the driving coils and the maximum displacement of the electromagnetic vibroseis. The height of the driving coils will limit the maximum displacement of the electromagnetic vibroseis. Figure 23(b) shows the force factor distribution curves in the air gap magnetic field model for different driving coils heights. As the height of the driving coils increases, the value of the force factor increases gradually. However, the nonlinearity of the force factor first increases and then decreases. The electromagneticmechanical coupling process of driving coils with different heights is solved based on the PSD method. The maximum displacement of the driving coils is ±10 mm. The vibration frequency is 2 Hz. Figure 23(c) shows the fluctuation of the force factor with time. Figure 23(d) shows the amplitude ratios of multiple harmonics and fundamental waves of acceleration in the air gap magnetic field model for different driving coils heights. The driving coils height has less influence on the second harmonic of acceleration. As the height of the driving coils increases, the amplitude ratio of the third harmonic to the fundamental wave first increases and then decreases.

Conclusion
In this paper, the electromagnetic-mechanical coupling is studied by analyzing the force factor in the finite element model of air gap magnetic field. The Gaussian fitting method and PSD method are used to analyze the influence of nonlinearity of force factor distribution curve on vibration signal. The nonlinearity of the radial magnetic induction intensity in the air gap is an important factor that causes low-frequency vibration distortion. Two novel optimization models for the air gap magnetic field are proposed.
An asymmetric axial dual-magnet design is proposed to enhance the radial magnetic induction intensity in the air gap. Compared with the axial single-magnet design, the radial magnetic induction intensity of the asymmetric axial dualmagnet design is increased by 22.2%. The design of adding a radial magnet to the outer yoke is proposed to compensate for the nonlinearity of the radial magnetic induction intensity in the air gap region between the axial magnet and the outer yoke. The proposed method can reduce the amplitude ratio of the third harmonic to the fundamental wave from 23.24% to 4.66%.
In the electromagnetic-mechanical coupling, the effect of the driving coils needs to be considered. The height of the driving coils will affect the maximum displacement of the electromagnetic vibroseis, the force factor distribution curve, and the vibration signal distortion.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).