Non-decoupling rapid calculation method for seismic electromagnetic field caused by piezomagnetic effect
By solving the PSVSH-TMTE equations in a non-decoupled manner and utilizing surface harmonic coordinate transformation, the problem of low efficiency in seismic electromagnetic field calculation is solved, enabling rapid quantitative calculation of electric and magnetic fields in a three-dimensional model.
Patent Information
- Application Number
- CN202511436919.X
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-10-09
- Publication Date
- 2025-12-26
AI Technical Summary
Existing technologies are difficult to efficiently calculate the electromagnetic response generated by earthquakes, especially when considering dynamic problems, and the calculation efficiency is low. Furthermore, they only provide static magnetic field results while ignoring the quantitative calculation of the electric field.
A non-decoupled approach is adopted to jointly solve the PSVSH-TMTE equations in the frequency-wavenumber domain. Through surface harmonic coordinate transformation, three sets of FWD non-decoupled seismograph linear differential equations are obtained, and displacement, stress, electric field and magnetic field are solved respectively. Finally, the equations are transformed to the time-space domain to improve computational efficiency.
It enables rapid calculation of seismic electromagnetic fields, and can simultaneously provide quantitative results for both electric and magnetic fields, significantly improving calculation speed and making it suitable for three-dimensional models of complex crustal structures.
Smart Images

Figure CN121208922A_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application belongs to the technical field of geophysics, and particularly relates to a non-decoupling fast calculation method for seismic electromagnetic field caused by piezomagnetic effect. BACKGROUND
[0002] Sometimes, abnormal electromagnetic signals are generated along with earthquakes, which has been observed for a long time, but the specific mechanism is not clear. In the rock rich in ferromagnetic material, the change of stress will cause the change of magnetization intensity, thereby generating a magnetic field, i.e. piezomagnetic effect. The previous research is based on the elastic statics of rock, and the permanent static magnetic field caused by piezomagnetic effect is calculated, without considering the dynamic problem to quantitatively calculate the electromagnetic response generated by piezomagnetic effect in earthquake.
[0003] The unpaired extranuclear electrons of the micro-particles constituting the ferromagnetic material can form the intrinsic magnetic moment of the atom. Below the Curie temperature (the temperature at which the ferromagnetic and paramagnetic phases are converted to each other), the exchange of electron spins between adjacent atoms will cause the ordered arrangement of the intrinsic magnetic moment, i.e. spontaneous magnetization phenomenon. In addition, when subjected to stress, the deformation of the crystal lattice will increase the magnetoelastic property of the system. In order to reduce the total free energy, the magnetization intensity (magnetic moment per unit volume) of the material will change. This is the micro-mechanism of piezomagnetic effect, which will theoretically lead to the generation of seismic electromagnetic signals.
[0004] The piezomagnetic effect seismic-electric coupling control equation in the frequency-wave number domain can also be solved directly in the rectangular coordinate system, but when the result is transformed back to the frequency-space domain, two wave number integrals, i.e. k x and k y integrals, must be performed, and the calculation efficiency is very low.
[0005] To solve the above problems, the present application provides a non-decoupling fast calculation method for seismic electromagnetic field caused by piezomagnetic effect.
[0006] It should be noted that the above content belongs to the technical cognition range of the inventor, and does not necessarily constitute the prior art. SUMMARY
[0007] To solve the above problems, the purpose of the present application is to provide a non-decoupling fast calculation method for seismic electromagnetic field caused by piezomagnetic effect, which proposes that the piezomagnetic effect of ferromagnetic minerals in the crust will lead to the generation of seismic electromagnetic phenomenon, and gives a complete three-dimensional piezomagnetic effect seismic-electric coupling equation. The PSVSH-TMTE equation is solved in a non-decoupling manner, and the displacement, stress, electric field and magnetic field in the layered model can be obtained at the same time.
[0008] To achieve the above purpose, the present application provides a non-decoupling fast calculation method for seismic electromagnetic field caused by piezomagnetic effect, which comprises the following steps:
[0009] S1: Transform the three-dimensional seismic-electromagnetic coupling equation of the frequency-space domain to FWD, and solve the elastic dynamic PSVSH wave equation and the electromagnetic field TMTE Maxwell equation jointly by using a non-decoupling algorithm to obtain three groups of FWD non-decoupling seismic-electromagnetic linear differential equations.
[0010] S2: Fully solve the PSVSH-TMTE equation of the three groups of FWD non-decoupling seismic-electromagnetic linear differential equations in the FWD to obtain the seismic and electromagnetic wave fields in the corresponding 3DHLM.
[0011] S3: Transform the FWD to TSD by two independent domain transformation methods, i.e., ω→t and k→r.
[0012] Further, the S1 specifically includes the following:
[0013] The face harmonic coordinate base is introduced, the frequency-space domain seismic-electromagnetic coupling control equation is given, the SHC component of the magnetization intensity change J is derived, and then the CCS component is obtained from the SHC component. The displacement u(r, θ, z), the body force F(r, θ, z), the horizontal stress vector Q(r, θ, z), the electric field E(r, θ, z), the magnetic field H(r, θ, z) and the magnetization intensity change J(r, θ, z) under the CCS component are respectively applied to any vector ξ(r, θ, z) under the CCS component to obtain the corresponding SHC expansion formula. The SHC expansion formula of each physical quantity is substituted into the control equation, and then three groups of FWD non-decoupling seismic-electromagnetic linear differential equations are obtained through mathematical derivation:
[0014]
[0015] where V q (q = -1, 0, +1) is an unknown displacement-stress-electromagnetic vector,
[0016]
[0017] u S , u T , u R are face harmonic components of the displacement u, Q S , Q T , Q R are face harmonic components of the horizontal stress vector Q, and are horizontal electromagnetic fields generated by J q , is a source term,
[0018]
[0019] F S , F T , F Ris the surface harmonic component of the body force F(r, θ, z).
[0020] Further, the coefficient sub-matrix M, Γ and N in the differential equation q respectively
[0021]
[0022]
[0023] Further, the SHC component of the magnetization change amount J includes the following steps:
[0024] First, the CCS component of the initial magnetization J0 is given to obtain J 0x , J 0y and J 0z , J 0x , J 0y and J 0z are three components of J0 in the rectangular coordinate system, and is substituted into three components, and then the SHC expansion of the magnetization change amount J is finally obtained through mathematical derivation.
[0025] Further, the S2 specifically includes:
[0026] Using the boundary conditions in the 3D HLM and the contribution of the source, three sets of linear algebraic equations about the amplitude are constructed to obtain the displacement u, stress τ, electric field E q and magnetic field H q under FWD, the magnetization change amount J q is further obtained according to u and τ, and then the formula,
[0027] B q = μ (H q + J q ),
[0028] The magnetic induction B q is obtained.
[0029] Further, the S3 specifically includes:
[0030] First, the displacement, stress, electric field and magnetic induction of FWD are transformed into the frequency-space domain using the k→r transformation formula, and then the displacement, stress, electric field and magnetic induction in the frequency-space domain are transformed into the TSD using the ω→t transformation formula.
[0031] The non-decoupling fast calculation method for the magnetostrictive effect caused seismic electromagnetic field provided by the application can bring the following beneficial effects:
[0032] 1. The calculation method of the present application proposes that the piezomagnetic effect of ferromagnetic minerals in the crust will lead to the generation of seismoelectric phenomenon, and gives a complete three-dimensional piezomagnetic effect seismoelectric coupling equation, which is solved by using a non-decoupling method, and the PSVSH-TMTE equation is solved simultaneously to obtain displacement, stress, electric field and magnetic field in a layered model.
[0033] 2. The calculation method of the present application transforms the original equation to a surface harmonic coordinate system for solving, and three groups of first-order linear differential equations with different q values are obtained in the process, and the displacement, stress, electric field and magnetic induction intensity in the frequency-wave number domain are obtained by solving the three groups of equations respectively. At this time, only one wave number integral, i.e. k integral, is needed when transforming back to the time-space domain, compared with two integrals (k x and k y integral) in the rectangular coordinate system, which can greatly improve the calculation speed.
[0034] 3. The calculation method of the present application considers elastic dynamics for calculation, compared with only considering statics results (only magnetic field), the electric field and magnetic field induced by piezomagnetic effect to generate earthquake can be quantitatively given at the same time. DETAILED DESCRIPTION
[0035] The drawings described herein are used to provide further understanding of the present application, and form a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application, and do not constitute an improper limitation on the present application. In the drawings:
[0036] Figure 1 It is a flow chart of the non-decoupling fast calculation method of the piezomagnetic effect induced seismoelectric field of the present application;
[0037] Figure 2 It is a three-dimensional layered model calculation schematic diagram adopted by the present application;
[0038] Figure 3 It is a displacement, electric field and magnetic induction intensity schematic diagram of the piezomagnetic effect induced in the full space model in Example 1
[0039] Figure 4 It is a schematic diagram of the source electromagnetic wave radiation pattern induced by the piezomagnetic effect in the full space model in Example 1;
[0040] Figure 5 It is a displacement, electric field and magnetic induction intensity schematic diagram of the piezomagnetic effect induced in the half space model in Example 2. DETAILED DESCRIPTION
[0041] In order to more clearly explain the overall concept of the present application, the following will be described in detail with reference to the drawings of the specification in an exemplary manner.
[0042] The embodiment of the application provides a non-decoupling fast calculation method of a magnetostrictive effect caused seismic electromagnetic field. Figure 1 As shown in the figure, the calculation method comprises the following steps:
[0043] 3DHLM (three-dimensional layered model) establishment:
[0044] Whether it is a seismograph or an electromagnetic instrument, it is generally installed near the ground surface, and the real crust has a very complex structure, often containing different types of rocks, and containing structures such as pores, cracks and faults, the scales of these structures range from several microns (such as rock pores) to hundreds of kilometers (such as faults), and complex physical and chemical interactions occur between different structures, and it is difficult to establish an accurate model to simulate the actual underground structure. However, from a larger scale, the underground rock is roughly layered, and it is reasonable to approximate the crust as a layered structure.
[0045] Figure 2 A schematic diagram of a three-dimensional layered structure model is shown, which sets that each rock layer is composed of uniform and isotropic elastic material, and there is a difference between each rock layer in physical parameters, which is specifically manifested in the difference distribution of mechanical parameters (such as density, Young's modulus, etc.) and electromagnetic parameters (such as electrical conductivity, initial magnetization intensity, etc.). The model uses a Cartesian coordinate system for spatial representation, and the coordinate origin is set at the ground surface reference plane, wherein the x-axis points to the geographic north, the y-axis corresponds to the east direction, and the z-axis vertically points to the center of the earth direction, and the seismic source is set on the z-axis. The top air layer is treated as a vacuum medium.
[0046] S1: Transform the frequency space domain magnetostrictive effect three-dimensional seismic electromagnetic coupling equation to FWD (frequency wave number domain), and solve the elastic dynamics PSVSH wave equation and the electromagnetic field TMTE Maxwell equation to obtain three groups of FWD non-decoupling seismic electromagnetic linear differential equations by using a non-decoupling algorithm, and the specific calculation method is as follows:
[0047] The surface harmonic coordinate basis vector is introduced,
[0048]
[0049] Wherein, k is the wave number, i is the imaginary unit, J m (kr) is the m-order first kind Bessel function, e r , e θ , e z constitute a CCS (cylindrical coordinate system), and the parameter m depends on the seismic source. When the seismic source is a point force source: m = 0, ±1; when the seismic source is a seismic moment tensor source: m = 0, ±1, ±2. And constitute a surface harmonic coordinate system.
[0050] The transformation of the three-dimensional seismic-electric coupling equation with magnetostrictive effect in the frequency-space domain to the FWD is as follows:
[0051] The control equation of seismic-electric coupling with magnetostrictive effect in the frequency-space domain is given as follows:
[0052]
[0053] where, is the Hamiltonian operator, τ is the stress tensor, f is the body force density vector, I is the unit matrix, ρ is the medium density, ω is the angular frequency, λ and G are the Lame coefficients, E and H represent the electric field and magnetic field, respectively, B and J represent the magnetic induction intensity and the change of magnetization, respectively, represents the equivalent dielectric constant, i.e. where σ, ε, and μ represent the electrical conductivity, dielectric constant, and magnetic permeability of the rock medium, respectively, β represents the stress sensitivity coefficient, P represents the deviatoric stress tensor, and J0 represents the initial magnetization of the rock.
[0054] The SHC (surface harmonic coordinate) component of the change of magnetization J is derived as follows:
[0055] First, the CCS component of the initial magnetization J0 is given as follows:
[0056] J0= J 0r e r + J 0θ e θ + J 0z e z ,
[0057] J 0r = J 0x cosθ + J 0y sinθ,
[0058] J 0θ = -J 0x sinθ + J 0y cosθ.
[0059] where J 0x , J 0y , and J 0z are the three components of J0 in the rectangular coordinate system, J 0r , J 0θ are the corresponding cylindrical coordinate components, e r , e θ , and e z are the cylindrical coordinate base vectors, and sinθ and cosθ are the sine value and cosine value of the azimuth angle θ, respectively. For the SHC component of J, it can be obtained by the following integral,
[0060]
[0061] Substituting the above three components into the equation,
[0062] The mathematical derivation finally leads to the SHC expansion of the magnetization change J,
[0063]
[0064]
[0065] where,
[0066]
[0067] The next vector of the CCS, ξ(r,θ,z), can be expanded as follows:
[0068]
[0069]
[0070] where,
[0071] ξ S,m (z,k), ξ T,m (z,k) and ξ R,m (z,k) are functions of m and k, called SHC (surface harmonics coordinate) expansion coefficients, which can be calculated from the following integral expressions,
[0072]
[0073] [] * represents the conjugate.
[0074] Similarly, the CCS components can also be obtained from the SHC components,
[0075]
[0076] It is noted that the transformation of ξ from the wave-number domain to the spatial domain requires only one wave-number integration. For the displacement u(r,θ,z), the body force F(r,θ,z), the horizontal stress Q(r,θ,z), the electric field E(r,θ,z), the magnetic field H(r,θ,z) and the magnetization change J(r,θ,z) under the CCS, the following equations are applied respectively,
[0077] The corresponding SHC expansion can be obtained in turn by expanding,
[0078]
[0079] Substitute the SHC expansion of each physical quantity into the control equation, and obtain three groups of FWD non-decoupling piezomagnetic effect seismoelectric linear differential equations through mathematical derivation, which specifically include the following:
[0080] Substitute the SHC expansion of each physical quantity into the control equation
[0081]
[0082] Through mathematical derivation, three groups of FWD non-decoupling seismoelectric linear differential equations are obtained,
[0083]
[0084] where V q (q=-1,0,+1) is an unknown displacement-stress-electromagnetic vector,
[0085]
[0086] u S , u T , u R are the surface harmonic components of displacement u, Q S , Q T , Q R are the surface harmonic components of horizontal stress vector Q. and are the horizontal electromagnetic fields generated by J q , is a source term,
[0087] F=[0,0,0,-F S ,-F T ,-F R ,0,0,0,0] T .
[0088] F S , F T , F R are the surface harmonic components of body force F(r,θ,z).
[0089] The coefficient sub-matrices M, Γ and N q in the differential equation are respectively:
[0090]
[0091]
[0092] S2: completely solve the PSVSH-TMTE equations of the three groups of FWD non-decoupling seismoelectric linear differential equations in the FWD, and obtain the seismic and electromagnetic wave fields in the corresponding 3DHLM;
[0093] Use the boundary conditions in the 3DHLM,
[0094]
[0095] And the contribution of the epicenter,
[0096]
[0097]
[0098] in:
[0099]
[0100] M xx M yy M zz M xy M yx M xz M zx M yz M zy These are the nine components of the seismic moment tensor. δ m,0 δ m,-1 δ m,1 δ m,-2 δ m,2 It is the Dirac delta function, with a value of 1 when the subscripts are the same and a value of zero when they are different.
[0101] By constructing three sets of linear algebraic equations about the amplitude, the displacement u, stress τ, and electric field E under FWD are obtained. q and magnetic field H q Based on u and τ, the change in magnetization J is further obtained. q Then use the formula,
[0102] B q =μ(H q +J q ),
[0103] Magnetic induction intensity B q ;
[0104] S3: Transform the FWD to the TSD (time-space domain) through two independent domain transformation methods, namely ω→t and k→r, including;
[0105] First, use the k→r transformation formula:
[0106]
[0107] The displacement, stress, electric field, and magnetic flux density of the FWD are transformed into the frequency-space domain, and then the ω→t transformation formula is used.
[0108]
[0109] Transform the displacement, stress, electric field and magnetic induction in the frequency-space domain to TSD.
[0110] It is noted that each physical quantity in is a three-dimensional space field quantity, but only one k integration is needed when performing spatial transformation k→r, instead of two integrations (k x and k y ) in the rectangular coordinate system, so the efficiency of the algorithm can be greatly improved, and the calculation speed can be significantly improved.
[0111] Embodiment 1
[0112] Considering a homogeneous full-space model, as shown in Figure 2 (a), the electromagnetic signal generated by body waves is studied. The rock parameters used are: V P = 6055 m / s, V S = 3500 m / s, ρ = 2600 kg / m 3 , σ = 0.001 S / m, |J0| = 1 A / m, β = 1 × 10 -9 pa -1 , ε = 6.64 × 10 -11 F / m, μ = 4π × 10 -7 H / m, θ d = -6°, θ i = 49°. The source uses a double force couple point source, M xz = M zx = 1.2 × 10 18 Nm, which corresponds to a simulated M w 6.0 earthquake. The initial magnetization J0 is parallel to the geomagnetic field, and the declination and inclination angles are θ d and θ i , respectively, and the intensity is |J0|. β is the stress sensitivity coefficient of the rock.
[0113] Figure 3 The time-domain waveforms of the displacement field, electric field and magnetic induction at the measuring point (40 km, 75 km, -10 km) are given. From the figure, it can be seen that there are obvious P waves (arrival time 14.1 s) and S waves (arrival time 24.5 s), as well as the electric field and magnetic field appearing at the same time, called coseismic electromagnetic field, which is the black curve in (d)-(i). At this time, the displacement amplitude is about 0.01 m, the coseismic electric field reaches about 1, and the coseismic magnetic field reaches about 0.1 nT, and the amplitudes of the electric field and the magnetic field can be measured by existing electromagnetic instruments. It can be seen that the piezomagnetic effect can not only produce a magnetic field but also produce an electric field under the action of the seismic source, which is different from the result obtained based on elastic statics (only static magnetic field). Figure 3
[0114] Figure 3 It is shown that the earthquake can also generate early electromagnetic waves arriving earlier than P waves, with very small energy, the electric field is only 10 -5 μV / m order of magnitude, the magnetic field is only 10 -8 nT order of magnitude, which may be because the distance from the receiving point to the source is too far, the electromagnetic wave is greatly attenuated in the process of transmission. In order to further study the early electromagnetic wave generated by the earthquake, such as Figure 2 (a), 10201 receivers are arranged in the horizontal plane 10 km above the source, and the amplitude of the electromagnetic wave at each receiving point is extracted to obtain Figure 4 It can be seen that the energy of the electromagnetic wave is mainly concentrated near the epicenter, the electric field can reach about 0.01 μV / m, and the magnetic field can reach about 0.01 nT, which can be observed.
[0115] Example 2
[0116] The actual stratum always has an interface, such as the free surface. As Figure 2 (b), a two-layer model is used. The rock parameters are the same as in the full space. The double force couple source M xz = M zx = 1.2 x 10 18 N m at a depth of 10 km. The Curie isotherm depth H c is 15 km. The initial magnetization intensity of the rock above the isotherm depth is J0, and there is no piezomagnetic effect below the isotherm depth, and the initial magnetization intensity is 0. The electromagnetic response generated by the earthquake is calculated in this model.
[0117] Figure 5 The displacement and electromagnetic field of the observation point near the ground surface (40 km, 75 km, 0.1 m) are given. In the displacement field, the P wave, S wave and Rayleigh wave (R wave) can be seen Figure 5 a-c). In the electric field and magnetic field Figure 5 d-i), their co-seismic electromagnetic signals can be seen. Figure 5 (d) shows the electromagnetic signal E x before the arrival of the seismic wave, which arrives at about 2.8 s, which is the interface signal generated by the direct S wave vertically incident to the free surface excited by the source. It should also exist in the other components of the electric field and the magnetic field, but it is not obvious in Figure 5 (e)-(i) because it is much smaller than the co-seismic signal. The source can directly excite electromagnetic waves, but due to the conductivity of the stratum, the electromagnetic wave decays very quickly, and its amplitude is smaller than that of the interface electromagnetic wave, which is not shown in the figure.
[0118] The method provides a method for quickly calculating a magneto-piezo effect caused earthquake electromagnetic field. Considering that the crustal rock contains ferromagnetic material with a magneto-piezo effect, complete three-dimensional magneto-piezo effect earthquake-electricity coupling equations are given. A unique non-decoupling method is adopted to solve the three-dimensional earthquake-electricity coupling PSVSH-TMTE equations. In order to improve the calculation efficiency, the surface harmonic coordinate system is transformed to solve the first-order linear differential equations of three groups of different q values. By solving the three groups of two-dimensional equations, the displacement, stress, electric field and magnetic induction intensity in the three-dimensional space in the frequency-wave number domain can be obtained. At this time, the transformation back to the time-space domain only needs to perform a wave number integral, that is, k integral, compared with the two times of integral (k x and k y integral) in the rectangular coordinate, the calculation speed is significantly improved. Since the elastic dynamics is considered for calculation, compared with the statics results (only static magnetic field), the electric field and magnetic field caused by the magneto-piezo effect can be quantitatively given at the same time. The theoretical calculation results show that the earthquake early electromagnetic wave caused by the magneto-piezo effect can help to carry out earthquake early warning.
[0119] Each of the embodiments in the specification is described in a progressive manner, and the same and similar parts between the embodiments can be referred to each other, and each embodiment mainly describes the difference from other embodiments. Especially, for the system embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and the related parts can be referred to the part of the method embodiment.
[0120] The above only describes the embodiments of the present application and does not limit the present application. For those skilled in the art, the present application can have various modifications and changes. Any modification, equivalent replacement, improvement and the like within the spirit and principle of the present application shall be included in the scope of claims of the present application.
Claims
1. A method for fast non-decoupled calculation of magnetostrictive effect induced seismic electromagnetic fields, characterized in that, The calculation method comprises: S1: transforming the frequency space domain magnetostatic effect three-dimensional seismic-electric coupling equation to FWD, adopting a non-decoupling algorithm, jointly solving the elastic dynamics PSVSH wave equation and the electromagnetic field TMTE Maxwell equation to obtain three groups of FWD non-decoupling seismic-electric linear differential equations; S2: completely solving the PSVSH-TMTE equation of the three groups of FWD non-decoupling seismic-electric linear differential equations in FWD to obtain the seismic and electromagnetic wave field in the corresponding 3DHLM; S3: transforming FWD to TSD through two independent domain transformation methods, namely ω→t and k→r.
2. The method according to claim 1, wherein, The S1 specifically comprises the following: introducing a surface harmonic coordinate base vector, giving a frequency space domain magnetostatic effect seismic-electric coupling control equation, deducing the SHC component of the magnetization intensity change J, and then obtaining the CCS component from the SHC component; performing surface harmonic expansion on the displacement u, the body force F, the horizontal stress vector Q, the electric field E, the magnetic field H and the magnetization intensity change J vector under the CCS component to obtain corresponding SHC expansion formulas, substituting the SHC expansion formulas of each physical quantity into the control equation, and then obtaining three groups of FWD non-decoupling seismic-electric linear differential equations through mathematical derivation: where V q (q = -1, 0, +1) is an unknown displacement-stress-electromagnetic vector, u S , u T , u R is the surface harmonic component of the displacement u, Q S , Q T , Q R is the surface harmonic component of the horizontal stress vector Q, and is J q the horizontal direction electromagnetic field generated by, is the source term, F S , F T , F R are the surface harmonic components of the force F(r, θ, z).
3. The method of claim 2, wherein, The coefficient sub-matrices M, Γ and N in the differential equation q are respectively 4. The method of claim 3, wherein, The deducing of the SHC component of the magnetization intensity change J comprises the following steps: First, the CCS components of the initial magnetization J0 are given to obtain J 0x , J 0y , and J 0z , J 0x , J 0y , and J 0z are three components of J0 in the rectangular coordinate system, and are substituted into the three components, and then the SHC expansion of the magnetization change J is finally obtained through mathematical derivation.
5. The method of claim 4, wherein, The S2 specifically comprises: With the boundary conditions in 3D HLM and the contribution of the source, the displacement u, stress τ, electric field E and magnetic field H under FWD are solved by constructing three sets of linear algebraic equations about amplitude q q The change of magnetization J is further obtained from u and τ q Then the formula is used, B q = μ(H q + J q ), The magnetic induction B is obtained q .
6. The method of claim 5, wherein, The S3 specifically comprises: firstly, using the k→r transformation formula to transform the displacement, stress, electric field and magnetic induction intensity of FWD to the frequency-space domain, and then using the ω→t transformation formula to transform the displacement, stress, electric field and magnetic induction intensity of the frequency-space domain to TSD.