CN102062868A

CN102062868A - Positioning and back-tracking method for earthquake electromagnetic wave source in ionized layer

Info

Publication number: CN102062868A
Application number: CN2009102380803A
Authority: CN
Inventors: 方广有; 曾中超; 王东峰; 张蓓
Original assignee: Institute of Electronics of CAS
Current assignee: Institute of Electronics of CAS
Priority date: 2009-11-18
Filing date: 2009-11-18
Publication date: 2011-05-18

Abstract

The invention discloses a positioning and back-tracking method for an earthquake electromagnetic wave source in an ionized layer, and relates to electromagnetic wave technology. By the method, a radiation source is positioned according to a received signal on a satellite. In the ionized layer, a route of an electromagnetic wave is determined by a Haselgrove differential equation at a similar high frequency; according to the equation, influences on the transmission route by the anisotropism and the dispersivity of the ionized layer are taken into consideration. The theory of the algorithm comprises that: an arithmetic solution of a ray equation is calculated by making use of the value technology; during calculation, the precision is controlled. In the actual operation, the actually measured electronic concentration of the ionized layer can be used. The effectiveness of the algorithm is verified by the actually measured data. The method can be applied to position of the radiation source in the ionized layer and also applied to position of a sound wave source in the ocean.

Description

Seismic Electromagnetic Wave Source Location and Backtracking Method in Ionosphere

技术领域technical field

本发明涉及电磁波技术领域，是一种电离层中的地震电磁波波源定位与后向追踪方法。The invention relates to the technical field of electromagnetic waves, and relates to a seismic electromagnetic wave source location and backward tracking method in the ionosphere.

背景技术Background technique

早在上世纪五六十年代，人们就开始了地震电磁研究，经过四十多年的不懈努力，人们已取得了大量与地震相关的电磁异常现象，这一课题的研究对地震的短期预测具有重要的意义。这些电磁前兆在地面上和高空中都有出现，但人们更感兴趣的主要还是电离层和磁层中的电磁异常，主要原因是高空中的异常现象所受的干扰和污染小，更容易借助于空间技术和遥感技术来探测和追踪。As early as the 1950s and 1960s, people began to study earthquake electromagnetics. After more than 40 years of unremitting efforts, people have obtained a large number of electromagnetic anomalies related to earthquakes. The research on this topic is of great significance to the short-term prediction of earthquakes. Significance. These electromagnetic precursors appear on the ground and in the high altitude, but people are more interested in the electromagnetic anomalies in the ionosphere and magnetosphere. Detection and tracking are based on space technology and remote sensing technology.

地球电离层是近地空间环境的重要组成部分，按照美国电器和电子工程师协会(IEEE)标准(1969)，电离层是“地球大气层的一部分，其中存在的粒子和电子数量多到足以影响无线电波的传播”。按照这个定义，电离层约是地面60km以上到磁层顶之间的整个空间。电磁波在电离层中的传播，主要受电离层电子浓度剖面、电离层的碰撞频率和电离层中地磁场的分布有关。The Earth's ionosphere is an important part of the near-Earth space environment. According to the Institute of Electrical and Electronics Engineers (IEEE) standard (1969), the ionosphere is "the part of the Earth's atmosphere in which there are enough particles and electrons to affect radio waves. Spread". According to this definition, the ionosphere is about the entire space between the ground 60km and the magnetopause. The propagation of electromagnetic waves in the ionosphere is mainly related to the electron concentration profile of the ionosphere, the collision frequency of the ionosphere and the distribution of the geomagnetic field in the ionosphere.

为了获得电磁波在空间中的传播方向，人们早已研发出能直接测量电磁场的各个分量的矢量天线，这些矢量天线实时记录电场和磁场的各个分量的幅度值和相位值。为了从这些原始测量数据中提取出波的传播信息(包括极化信息和波矢方向)，人们发展了好几种高效的波矢分析算法(如Means算法、波谱矩阵的SVD技术等)。In order to obtain the propagation direction of electromagnetic waves in space, people have developed vector antennas that can directly measure the various components of the electromagnetic field. These vector antennas record the amplitude and phase values of the various components of the electric field and magnetic field in real time. In order to extract wave propagation information (including polarization information and wave vector direction) from these raw measurement data, several efficient wave vector analysis algorithms (such as Means algorithm, SVD technology of wave spectrum matrix, etc.) have been developed.

发明内容Contents of the invention

本发明的目的是公开一种电离层中的地震电磁波波源定位与后向追踪方法，是一种高精度快速算法，以寻找空间无线电波的辐射源，该方法主要基于电离层的无线电波传播、射线微分方程的数值解法，并对追踪的误差进行了有效的控制，可精确定位地震区域。The purpose of the present invention is to disclose a seismic electromagnetic wave source location and backward tracking method in the ionosphere, which is a high-precision fast algorithm to find the radiation source of space radio waves. The method is mainly based on ionosphere radio wave propagation, Numerical solution of ray differential equation, and effectively control the tracking error, can accurately locate the earthquake area.

为达到上述目的，本发明的技术解决方案是：For achieving the above object, technical solution of the present invention is:

一种电离层中的地震电磁波波源定位与后向追踪方法，其利用接收电磁波的来波方向、频率、传播媒质的特性，进行电离层的射线追踪，对射线方程进行数值求解，来定位辐射源；其包括步骤：A seismic electromagnetic wave source location and backtracking method in the ionosphere, which uses the direction of arrival, frequency, and characteristics of the propagation medium of the received electromagnetic wave to perform ray tracing in the ionosphere, and numerically solve the ray equation to locate the radiation source ; which includes the steps of:

A)理论推导：A) Theoretical derivation:

(a)描述波的传播路径的方程是Haselgrove方程，将该方程改写为适合于数值求解的形式，求出折射率对各个变量的一阶偏导数；(a) The equation describing the propagation path of the wave is the Haselgrove equation, which is rewritten into a form suitable for numerical solution, and the first-order partial derivative of the refractive index to each variable is obtained;

(b)选择合适的数值方法来对射线方程求解，推导出具体的迭代格式；(b) Select an appropriate numerical method to solve the ray equation, and derive a specific iterative format;

B)计算机编程，即电离层射线追踪软件：B) Computer programming, i.e. ionospheric ray tracing software:

(a)根据A)步的理论推导结果，编制出计算电离层中波传播的计算机软件；(a) According to the theoretical derivation result of step A), the computer software for calculating the wave propagation in the ionosphere is compiled;

(b)由矢量天线的实测数据提取出波的初始传播方向；(b) Extract the initial propagation direction of the wave from the measured data of the vector antenna;

(c)由于波的传播路径具有可逆性，将接收点的初始方向反向，进行“后向”追踪，作为软件计算的初始条件之一；(c) Since the propagation path of the wave is reversible, the initial direction of the receiving point is reversed, and "backward" tracking is performed as one of the initial conditions for software calculation;

(d)设置软件计算的其他初始条件：接收点的位置、计算的容差、追踪的截止条件；(d) Set other initial conditions calculated by the software: the position of the receiving point, the tolerance of the calculation, and the cut-off condition of the tracking;

C)确定震中位置：C) Determine the location of the epicenter:

追踪的截止点就是辐射源在空间的波源，利用截止点确定震中位置。The cut-off point of the tracking is the wave source of the radiation source in space, and the epicenter position is determined by the cut-off point.

所述的方法，其所述B)步，直接对三维的Haselegrove方程进行数值求解，以避免地球磁场和电离层的色散性对波传播路径的影响，使路径更加精确。The method, the step B), directly solves the three-dimensional Haselegrove equation numerically, so as to avoid the influence of the dispersion of the earth's magnetic field and the ionosphere on the wave propagation path, and make the path more accurate.

所述的方法，其所述B)步，在由当前节点根据积分算法推进到下一个节点时，都对节点值进行校正，以避免误差的累积。In the method, in the step B), when the current node advances to the next node according to the integral algorithm, the node value is corrected to avoid the accumulation of errors.

所述的方法，其所述B)步，采用了变步长技术，该变步长技术充分考虑了媒质的折射率的空间梯度，在梯度大的地方步长小，梯度小的地方步长大，以保证追踪精度。Described method, its described B) step, has adopted variable step size technology, and this variable step size technology has fully considered the spatial gradient of the refractive index of medium, and the place step size is small in the big place of gradient, and the place step size of little gradient large to ensure tracking accuracy.

所述的方法，其所述B)步中的追踪截止条件，为在电离层的F2层峰值高度。In the described method, the tracking cut-off condition in the step B) is the peak height of the F2 layer in the ionosphere.

所述的方法，其所述C)步中，利用截止点确定震中位置，是在追踪到截止点后，以截止点沿着地球磁场的磁力线向地面投影，投影点就是震中位置。Described method, in its described C) step, utilize the cut-off point to determine the epicenter position, after tracking to the cut-off point, project to the ground with the cut-off point along the magnetic field line of the earth's magnetic field, and the projected point is exactly the epicenter position.

所述的方法，其前提是利用矢量天线获得多分量的电磁场数据，并根据波矢分析算法提取出波的频率、来波方向和极化特性信息。The premise of the method is to use the vector antenna to obtain multi-component electromagnetic field data, and extract the frequency, incoming wave direction and polarization characteristic information of the wave according to the wave vector analysis algorithm.

所述的方法，其所述积分算法，为Runge-Kutter算法、Adams-Moulton算法。Said method, said integration algorithm, is Runge-Kutter algorithm, Adams-Moulton algorithm.

所述的方法，其所述变步长技术，是由下式根据媒质的梯度自动调整步长：Described method, its described variable step size technique, is to be automatically adjusted step size according to the gradient of medium by following formula:

$step step = = \frac{kn k n}{| | \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; r r} + + \frac{11}{r r} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; θ θ} + + \frac{11}{r r sin sin θ θ} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; φ φ} | |}$

式中，k为一个比例常数，k的数值一般根据经验和追踪的精度来确定，在实际操作时，可以通过反复运行射线追踪程序来选择一个合适的值(一般介于10^-3和10^-5之间)。当折射率的沿着射线路径的梯度较小时，即折射率的变化较平缓，则步长较大；反之，当折射率沿着射线路径剧烈变化时，步长自动减小；In the formula, k is a proportional constant. The value of k is generally determined according to experience and tracking accuracy. In actual operation, a suitable value can be selected by repeatedly running the ray tracing program (generally between 10 ^-3 and 10 ^{- 5} ). When the gradient of the refractive index along the ray path is small, that is, the change of the refractive index is relatively gentle, the step size is large; on the contrary, when the refractive index changes sharply along the ray path, the step size is automatically reduced;

同时，对步长设定一个上限step_max和下限step_min，当按照上式计算出的步长超出这个范围时，就以上限步长或者下限步长作为实际的步长值。At the same time, set an upper limit step _max and a lower limit step _min for the step size. When the step size calculated according to the above formula exceeds this range, the upper limit step size or the lower limit step size is used as the actual step size value.

所述的方法，其所述波矢分析算法，为Means算法、波谱矩阵的SVD技术。In the method, the wave vector analysis algorithm is the Means algorithm and the SVD technology of the wave spectrum matrix.

本发明的方法，是一种高精度快速算法，对追踪的误差进行了有效的控制，可精确定位地震区域。The method of the invention is a high-precision fast algorithm, effectively controls the tracking error, and can accurately locate the earthquake area.

附图说明Description of drawings

图1为本发明一种电离层中的地震电磁波波源定位与后向追踪方法计算的坐标系示意图；其中：Fig. 1 is a coordinate system schematic diagram of seismic electromagnetic wave source location and backward tracking method calculation in a kind of ionosphere of the present invention; Wherein:

图1(a)为本发明计算所用的球坐标系；Fig. 1 (a) calculates the used spherical coordinate system for the present invention;

图1(b)为P点的局域坐标系的放大图；Figure 1(b) is an enlarged view of the local coordinate system of point P;

图2为本发明一种电离层中的地震电磁波波源定位与后向追踪方法的电离层射线追踪软件流程图。Fig. 2 is a flow chart of ionospheric ray tracing software for a seismic electromagnetic wave source location and backward tracing method in the ionosphere according to the present invention.

具体实施方式Detailed ways

本发明的一种电离层中的地震电磁波波源定位与后向追踪方法，是一种追踪高空电磁异常的方法，由此方法可以追踪得到产生这些异常的辐射源的空间位置。在利用本发明的方法之前，必须获得电离层的这些先验信息。借助于现有的波矢分析算法(如Means算法、波谱矩阵的SVD技术等)，从矢量天线的原始实测数据中提取出波的初始方向，是实施本发明方法的前提和基础，而这些已公开的波矢分析算法本身不是本发明的讨论范围，在此不再赘述。在应用本发明方法之前，已默认获得了下列信息：The invention relates to a seismic electromagnetic wave source positioning and backward tracking method in the ionosphere, which is a method for tracking high-altitude electromagnetic anomalies, whereby the method can track and obtain the spatial positions of radiation sources that produce these anomalies. Such a priori information on the ionosphere must be obtained prior to utilizing the method of the present invention. By means of existing wave vector analysis algorithms (such as Means algorithm, SVD technology of wave spectrum matrix, etc.), extracting the initial direction of wave from the original measured data of vector antenna is the premise and basis of implementing the method of the present invention, and these have been The disclosed wave vector analysis algorithm itself is not within the scope of the present invention, and will not be repeated here. Before applying the method of the present invention, the following information has been obtained by default:

(a)、已通过波矢分析算法获得了波的初始传播方向；(a), the initial propagation direction of the wave has been obtained through the wave vector analysis algorithm;

(b)、电离层的电子浓度的三维结构信息N(h，θ，

)是已知的，其中，N为电离层的电子浓度，单位为米^-3，h为高度(米)，θ和

分别为计算坐标系(见下文)中的极角和方位角，也就是地理上的余纬和经度。(b), three-dimensional structure information N(h, θ,

) is known, where N is the electron concentration of the ionosphere in m ^-3 , h is the height (m), θ and

are the polar and azimuth angles in the computational coordinate system (see below), that is, the geographic co-latitude and longitude, respectively.

见图1，为本发明一种电离层中的地震电磁波波源定位与后向追踪方法计算的坐标系示意图。图1(a)为本发明计算所用的球坐标系；图1(b)为P点的局域坐标系的放大图；See FIG. 1 , which is a schematic diagram of a coordinate system calculated by a seismic electromagnetic wave source location and backward tracking method in the ionosphere according to the present invention. Fig. 1 (a) is the spherical coordinate system that the present invention calculates used; Fig. 1 (b) is the enlargement figure of the local coordinate system of P point;

本发明方法计算的坐标系：The coordinate system that the inventive method calculates:

为了表示射线路径和波的传播方向，必须建立坐标系。以地心O为原点，以地球的自转轴为z轴建立一个计算坐标系(为球坐标系)，如图1(a)所示。在这个坐标系中，射线路径上的某点P用3个参数来表征：r，θ，φ。为了表征波在P点的传播方向(即波矢

的方向)，在P点建立一个局部直角坐标系

在这个坐标系中，方向表示OP连线的方向，

方向表示P点的正南方向，而

方向表示P点的正东方向，见图1(b)。在P点的波矢的方向既可以用局部直角坐标系中的两个参数(θ_k，φ_k)来表征，也可以用P点的仰角α和方位角β来表征(后文将用到这两个角度)。仰角α是从水平面旋转到波矢

的角度，范围为0°～90°；方位角β是从正北方向顺时针旋转到

在水平面的投影的角度，范围为0°～360°。因此，θ_k，φ_k和α，β之间的关系为θ_k＝π/2-α，φ_k＝π-β。In order to represent the ray path and the direction of wave propagation, a coordinate system must be established. A calculation coordinate system (spherical coordinate system) is established with the earth center O as the origin and the earth's rotation axis as the z-axis, as shown in Figure 1(a). In this coordinate system, a point P on the ray path is characterized by three parameters: r, θ, φ. In order to characterize the propagation direction of the wave at point P (namely, the wave vector

direction), establish a local Cartesian coordinate system at point P

In this coordinate system, The direction indicates the direction of the OP connection,

The direction represents the due south direction of point P, while

The direction represents the due east direction of point P, see Figure 1(b). Wave vector at point P The direction of can be characterized not only by two parameters (θ _k , φ _k ) in the local Cartesian coordinate system, but also by the elevation angle α and azimuth angle β of point P (these two angles will be used later). The elevation angle α is the rotation from the horizontal plane to the wave vector

The angle ranges from 0° to 90°; the azimuth β rotates clockwise from true north to

The projection angle on the horizontal plane ranges from 0° to 360°. Therefore, the relationship between θ _k , φ _k and α, β is θ _k =π/2-α, φ _k =π-β.

在高频近似下，波的传播路径用射线来描述。当考虑地球磁场对电磁波传播的影响时，电离层呈现出各向异性，此时射线的能量方向和波矢方向并不完全重合，而是存在一个夹角。1955年，J.Haselgrove导出了适合于计算机求解的射线方程：In the high-frequency approximation, the wave propagation path is described by rays. When considering the influence of the earth's magnetic field on the propagation of electromagnetic waves, the ionosphere presents anisotropy. At this time, the energy direction of the ray and the direction of the wave vector do not completely coincide, but there is an included angle. In 1955, J.Haselgrove derived the ray equation suitable for computer solution:

$\frac{dr dr}{dt dt} = = \frac{11}{{n no}^{22}} (({ρ ρ}_{r r} - - n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; {ρ ρ}_{r r}}))$

$\frac{dθ dθ}{dt dt} = = \frac{11}{r r {n no}^{22}} (({ρ ρ}_{θ θ} - - n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; {ρ ρ}_{θ θ}}))$

$\frac{dφ dφ}{dt dt} = = \frac{11}{{rn rn}^{22} sin sin θ θ} (({ρ ρ}_{φ φ} - - n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; {ρ ρ}_{φ φ}}))$ $((11 a a - - f f))$

$\frac{d d {ρ ρ}_{r r}}{dt dt} = = \frac{11}{n no} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; r r} + + {ρ ρ}_{θ θ} \frac{dθ dθ}{dt dt} + + {ρ ρ}_{φ φ} \frac{dφ dφ}{dt dt} sin sin θ θ$

$\frac{d d {ρ ρ}_{θ θ}}{dt dt} = = \frac{11}{r r} ((\frac{11}{n no} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; θ θ} - - {ρ ρ}_{θ θ} \frac{dr dr}{dt dt} + + r r {ρ ρ}_{φ φ} \frac{dφ dφ}{dt dt} cos cos θ θ))$

$\frac{d d {ρ ρ}_{φ φ}}{dt dt} = = \frac{11}{r r sin sin θ θ} ((\frac{11}{n no} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; φ φ} - - {ρ ρ}_{φ φ} \frac{dr dr}{dt dt} sin sin θ θ - - r r {ρ ρ}_{φ φ} \frac{dθ dθ}{dt dt} cos cos θ θ))$

其中，r，θ，φ是射线路径上的一点P在球坐标中的坐标，ρ_r，ρ_θ，ρ_φ是P点沿着波矢方向的折射率在球坐标系中的3个分量(满足)，自变量t可以为沿着射线路径单调增加的任意参数(因为等式右边不显含t)，n为媒质(主要是电离层)的折射率。由于r，θ，φ和ρ_r，ρ_θ，ρ_φ耦合在一起，无法单独求出r，θ，φ，而必须对上面的6个方程同时求解。Among them, r, θ, φ are the coordinates of a point P on the ray path in spherical coordinates, ρ _r , ρ _θ , ρ _φ are the three components of the refractive index of point P along the wave vector direction in spherical coordinates ( satisfy ), the argument t can be any parameter that monotonically increases along the ray path (because the right side of the equation does not contain t explicitly), and n is the refractive index of the medium (mainly the ionosphere). Since r, θ, φ and ρ _r , ρ _θ , ρ _φ are coupled together, r, θ, φ cannot be calculated separately, but the above six equations must be solved simultaneously.

方程组(1a～f)是由6个一阶微分方程构成的方程组，射线追踪的原理是对它们进行数值求解，最终获得完整的射线轨迹描述。在求解之前，必须将方程(1a～f)的右边写成待求量r，θ，φ、ρ_r，ρ_θ，ρ_φ和自变量t的显式函数。在数学上，也就是将方程组(1a～f)转化为下面的初值问题：The equations (1a~f) are composed of six first-order differential equations. The principle of ray tracing is to solve them numerically, and finally obtain a complete ray trajectory description. Before solving, the right side of the equation (1a～f) must be written as an explicit function of the quantity to be sought r, θ, φ, ρ _r , ρ _θ , ρ _φ and the independent variable t. Mathematically, it is to transform the equation system (1a～f) into the following initial value problem:

$\{\begin{matrix} y_{i}^{'} = f_{i} (x, y_{1}, y_{2}, . . ., y_{N}) \\ y_{i} (x_{0}) = y_{i}^{0} \end{matrix}$ (i＝1，2，…，N)(2a，b) $\{\begin{matrix} {the y}_{i}^{'} = f_{i} (x, {they}_{1}, {they}_{2}, . . ., {the y}_{N}) \\ {the y}_{i} (x_{0}) = {the y}_{i}^{0} \end{matrix}$ (i=1,2,...,N)(2a,b)

对于这类初值问题，在数学上已经发展了比较成熟的算法，如经典的四阶Runge-Kutter算法、Adams-Moulton预测校正系统等。这些算法都是公知的，本发明不对它们进行详细介绍。For this kind of initial value problem, relatively mature algorithms have been developed in mathematics, such as the classic fourth-order Runge-Kutter algorithm, Adams-Moulton predictive correction system, etc. These algorithms are well known, and the present invention does not introduce them in detail.

在电离层应用中，折射率n由Appleton-Hartree公式(后文简写为A-H公式)来表示：In ionospheric applications, the refractive index n is represented by the Appleton-Hartree formula (hereinafter abbreviated as A-H formula):

${n no}^{22} = = 11 - - \frac{X x ((11 - - X x))}{11 - - X x - - \frac{11}{22} {Y Y}_{}^{T T} &PlusMinus; &PlusMinus; \sqrt{\frac{11}{44} {Y Y}_{T T}^{44} + + {Y Y}_{L L}^{22} {((11 - - X x))}^{22}}}$

Y_T＝Y sin Θ (3)Y _T = Y sin Θ (3)

Y_L＝Y cos ΘY _L = Y cos Θ

其中in

$X x = = \frac{{f f}_{N N}^{22}}{{f f}^{22}} - - - - - - ((44))$

$Y Y = = \frac{{f f}_{H h}}{f f} - - - - - - ((55))$

f_N和f_H分别为电离层的等离子体频率和由地磁场因此的电子回旋频率，f为被追踪的波的频率，Θ为波矢方向和地磁场的夹角，Y_T和Y_L分别表示Y垂直于波矢的横向分量和平行于波矢的纵向分量。在实际应用时，f_N既可以由实测的电离层探测数据获得，也可以由电离层的经验模式获得；f_H由地磁场的经验模式(如国际地磁场参考模型IGRF、偶极子模型等)获得；f可以从矢量天线的实测数据中提取(根据已有的波矢分析算法)。而(3)式分母中的正负号分别对应电离层波传播中的寻常模和非常模。在空间中某个给定点的折射率只是方向的函数，这体现出媒质的各向异性，而电离层的各向异性是由地球磁场引起的。下面讲述如何将折射指数n及其偏导数

(η表示3个位置坐标r，θ，φ)表示成6个待求量r，θ，φ、ρ_r，ρ_θ，ρ_φ的显式函数。f _N and f _H are the plasma frequency of the ionosphere and the electron cyclotron frequency caused by the geomagnetic field, f is the frequency of the tracked wave, Θ is the angle between the wave vector direction and the geomagnetic field, Y _T and Y _L are respectively Denotes the transverse component of Y perpendicular to the wave vector and the longitudinal component parallel to the wave vector. In practical application, f _N can be obtained from the measured ionospheric detection data, or from the empirical model of the ionosphere; f _H can be obtained from the empirical model of the geomagnetic field (such as the international geomagnetic field reference model IGRF, dipole model, etc. ) to obtain; f can be extracted from the measured data of the vector antenna (according to the existing wave vector analysis algorithm). The positive and negative signs in the denominator of (3) respectively correspond to the ordinary mode and the extraordinary mode in the ionospheric wave propagation. The refractive index at a given point in space is only a function of direction, which reflects the anisotropy of the medium, and the anisotropy of the ionosphere is caused by the Earth's magnetic field. The following describes how the refractive index n and its partial derivatives

(η represents three position coordinates r, θ, φ) expressed as explicit functions of 6 quantities to be sought r, θ, φ, ρ _r , ρ _θ , ρ _φ .

1.折射率n1. Refractive index n

折射率n的方向代表了波矢方向，因此有下列表达式：The direction of the refractive index n represents the direction of the wave vector, so there are the following expressions:

${Y Y}_{}^{L L} = = \frac{{(({ρ ρ}_{r r} {Y Y}_{r r} + + {ρ ρ}_{θ θ} {Y Y}_{θ θ} + + {ρ ρ}_{φ φ} {Y Y}_{φ φ}))}^{22}}{{ρ ρ}_{r r}^{22} + + {ρ ρ}_{θ θ}^{22} + + {ρ ρ}_{φ φ}^{22}} - - - - - - ((66))$

${Y Y}_{}^{T T} = = {Y Y}^{22} - - {Y Y}_{L L}^{22} - - - - - - ((77))$

上式中Y_η表示Y在球坐标系中的3个分量，它们的表达式将在后文给出。将Y_T，Y_L的表达式带入A-H公式，则折射率可以表示为：In the above formula, Y _η represents the three components of Y in the spherical coordinate system, and their expressions will be given later. Putting the expressions of Y _T and Y _L into the AH formula, the refractive index can be expressed as:

${n no}^{22} = = 11 - - \frac{22 X x ((11 - - X x))}{D D.}$

D＝2(1-X)-Y_T ²+R (8a，b，c)D=2(1-X)-Y _T ² +R (8a,b,c)

$R R = = &PlusMinus; &PlusMinus; \sqrt{{Y Y}_{T T}^{44} + + 44 {Y Y}_{L L}^{22} {((11 - - X x))}^{22}}$

2、偏导数

2. Partial derivative

下面讲述如何将折射率的平方对位置变量的导数表示为待求量的显式函数：Here's how to express the derivative of the square of the index of refraction with respect to the position variable as an explicit function of the quantity to be sought:

$n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; η η} = = n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; X x} \frac{&PartialD; &PartialD; X x}{&PartialD; &PartialD; η η} + + n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; Y Y} \frac{&PartialD; &PartialD; Y Y}{&PartialD; &PartialD; η η} + + n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; Θ Θ} \frac{&PartialD; &PartialD; Θ Θ}{&PartialD; &PartialD; η η} - - - - - - ((99))$

又有there are

$n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; X x} = = - - \frac{((22 {((11 - - X x))}^{22} - - {Y Y}_{T T}^{22} ((11 - - 22 X x)) + + \frac{{Y Y}_{T T}^{44} ((11 - - 22 X x)) + + {44 Y Y}_{L L}^{22} {((11 - - X x))}^{33}}{R R}))}{{D D.}^{22}} - - - - - - ((1010))$

$n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; Y Y} = = \frac{22 X x ((11 - - X x))}{{D D.}^{22} Y Y} ((- - {Y Y}_{}^{T T} + + \frac{{Y Y}_{T T}^{44} + + 22 {Y Y}_{L L}^{22} {((11 - - X x))}^{22}}{R R})) - - - - - - ((1111))$

$\frac{n no}{{Y Y}_{T T} {Y Y}_{L L}} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; Θ Θ} = = \frac{22 X x ((11 - - X x)) ((- - 11 + + \frac{{Y Y}_{T T}^{22} - - 22 {((11 - - X x))}^{22}}{R R}))}{{D D.}^{22}} - - - - - - ((1212))$

$Y_{T} Y_{L} \frac{&PartialD; Θ}{&PartialD; η} = \frac{Y_{L}^{2}}{Y} \frac{&PartialD; Y}{&PartialD; η} - \underset{η}{Σ} ρ_{η} \frac{&PartialD; Y_{η}}{&PartialD; η} \frac{Y_{L}}{ρ}$ (η＝r，θ，φ) (13) $Y_{T} Y_{L} \frac{&PartialD; Θ}{&PartialD; η} = \frac{Y_{L}^{2}}{Y} \frac{&PartialD; Y}{&PartialD; η} - \underset{η}{Σ} ρ_{η} \frac{&PartialD; Y_{η}}{&PartialD; η} \frac{Y_{L}}{ρ}$ (η=r, θ, φ) (13)

而和

可以由电离层的电子浓度剖面和磁场模型求出，将在后文给出。and and

It can be obtained from the electron concentration profile of the ionosphere and the magnetic field model, which will be given later.

3、偏导数

3. Partial derivative

$n no \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; {ρ ρ}_{η η}} = = \frac{n no}{{Y Y}_{T T} {Y Y}_{L L}} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; Θ Θ} ((\frac{{ρ ρ}_{η η} {Y Y}_{}^{L L}}{{ρ ρ}^{22}} - - \frac{{Y Y}_{L L}}{ρ ρ} {Y Y}_{η η})) - - - - - - ((1414))$

4、电离层的电子浓度模型4. Electron concentration model of the ionosphere

在实际实施时，电子浓度模型采用Chapman函数：In actual implementation, the electron concentration model uses the Chapman function:

${f f}_{}^{N N} = = {f f}_{}^{c c} exp exp [[\frac{11}{22} ((11 - - z z - - exp exp ((- - z z))))]] - - - - - - ((1515))$

$z z = = \frac{h h - - {h h}_{max max}}{H h} - - - - - - ((1616))$

其中f_c为赤道上的临界频率，可取6.5MHz，可取300km，H可取62km，h表示离地面的高度。在计算坐标系中，高度可表示为：Among them, f _c is the critical frequency on the equator, which can be 6.5MHz, 300km, H can be 62km, and h indicates the height from the ground. In the computational coordinate system, height can be expressed as:

h＝r-R_e (17)h=rR _e (17)

r表示到地心的距离，R_e为地球半径，取为6370km。前文(9)式中出现的

可表示为：r represents the distance to the center of the earth, Re _is the radius of the earth, which is taken as 6370km. Appearing in formula (9) above

Can be expressed as:

$\frac{&PartialD; &PartialD; X x}{&PartialD; &PartialD; r r} = = \frac{X x}{22 H h} [[exp exp ((- - z z)) - - 11]] - - - - - - ((1818))$

5、地磁场模型5. Geomagnetic field model

地球磁场采用偶极子模型来近似，则电子回旋频率可表示为：The earth's magnetic field is approximated by a dipole model, and the electron cyclotron frequency can be expressed as:

${f f}_{H h} ((r r,, θ θ)) = = {f f}_{{H h}_{00}} {((\frac{{R R}_{e e}}{r r}))}^{33} {((11 + + 33 {cos cos}^{22} θ θ))}^{11 / / 22} - - - - - - ((2020))$

其中，

为赤道地面上的电子回旋频率，约为0.85MHz；θ为计算坐标系的极角，在偶极子模型下等效为地磁余纬。(9)式出现的与Y有关的偏导数为：in,

is the electron cyclotron frequency on the equatorial ground, about 0.85MHz; θ is the polar angle of the calculation coordinate system, which is equivalent to geomagnetic co-latitude under the dipole model. The partial derivative related to Y appearing in formula (9) is:

$\frac{&PartialD; &PartialD; Y Y}{&PartialD; &PartialD; r r} = = - - \frac{33 Y Y}{r r} - - - - - - ((21 twenty one))$

$\frac{&PartialD; &PartialD; Y Y}{&PartialD; &PartialD; θ θ} = = - - 33 Y Y \frac{sin sin θ θ cos cos θ θ}{11 + + 33 {cos cos}^{22} θ θ} - - - - - - ((22 twenty two))$

(13)式中出现了Y在计算坐标系中的各个分量及其导数：Each component of Y in the calculation coordinate system and its derivative appear in the formula (13):

${Y Y}_{r r} = = - - \frac{{f f}_{{H h}_{00}}}{f f} {((\frac{{R R}_{00}}{r r}))}^{33} 22 cos cos θ θ - - - - - - ((24 twenty four))$

${Y Y}_{θ θ} = = \frac{11}{22} {Y Y}_{r r} tan the tan θ θ - - - - - - ((2525))$

$\frac{&PartialD; &PartialD; {Y Y}_{r r}}{&PartialD; &PartialD; r r} = = - - \frac{33 {Y Y}_{r r}}{r r},, \frac{&PartialD; &PartialD; {Y Y}_{r r}}{&PartialD; &PartialD; θ θ} = = - - 22 {Y Y}_{θ θ} - - - - - - ((2626))$

$\frac{&PartialD; &PartialD; {Y Y}_{θ θ}}{&PartialD; &PartialD; r r} = = - - \frac{33 {Y Y}_{θ θ}}{r r},, \frac{&PartialD; &PartialD; {Y Y}_{θ θ}}{&PartialD; &PartialD; θ θ} = = \frac{{Y Y}_{r r}}{22} - - - - - - ((2727))$

6、追踪的初始条件的设定：6. Setting of initial conditions for tracking:

将方程组(1a～f)表示成(2a)的形式后，还必须给出形如(2b)式的初始条件，即r，θ，φ，ρ_r，ρ_θ，ρ_φ的初值。前3个量表示追踪的起始点，输入球坐标系中起始点P₀的坐标(r₀，θ₀，φ₀)即可；而后3个量表示在起始点沿着来波方向的反向的折射率n在局部坐标系(见图1)中的3个分量。假设在起始点的来波方向的仰角和方位角分别为α₀，β₀，为了逆向追踪波源(根据射线路径的可逆性)，将追踪的初始方向设为-α₀，-β₀，若将沿着该方向的折射率归一化到1，则3个分量可以表示为：After expressing the equations (1a～f) in the form of (2a), the initial conditions in the form of (2b) must also be given, namely the initial values of r, θ, φ, ρ _r , ρ _θ , ρ _φ . The first three quantities represent the starting point of the tracking, just input the coordinates (r ₀ , θ ₀ , φ ₀ ) of the starting point P ₀ in the spherical coordinate system; and the last three quantities represent the reverse direction of the incoming wave at the starting point The refractive index n has 3 components in the local coordinate system (see Figure 1). Assuming that the elevation angle and azimuth angle of the incoming wave direction at the starting point are α ₀ , β ₀ , in order to trace the wave source backwards (according to the reversibility of the ray path), the initial direction of tracing is set to -α ₀ , -β ₀ , if Normalize the refractive index along this direction to 1, then the three components can be expressed as:

ρ_r0＝sin(-α₀)，ρ_θ0＝-cos(-α₀)cos(-β₀)，ρ_φ0＝cos(-α₀)sin(-β₀) (28)ρ _r0 ＝sin(-α ₀ ), ρ _θ0 ＝-cos(-α ₀ )cos(-β ₀ ), ρ _φ0 ＝cos(-α ₀ )sin(-β ₀ ) (28)

由于这里只考虑了折射率的方向而没有考虑折射率的幅度，在实际计算时，应该对上述初值的幅度值作出修正：应用A-H公式计算出沿着(-α₀，-β₀)方向的折射率大小n，则修正后的初值为Since only the direction of the refractive index is considered here but not the magnitude of the refractive index, in the actual calculation, the magnitude of the above initial value should be corrected: apply the AH formula to calculate along the (-α ₀ ,-β ₀ ) direction The refractive index size n, then the corrected initial value is

$ρ_{η 0 correct} = ρ_{η 0} \frac{n}{ρ}$ (η＝r，θ，φ) (29) $ρ_{η 0 correct} = ρ_{η 0} \frac{no}{ρ}$ (η=r, θ, φ) (29)

实际上，在由当前节点按照Runge-Kutter积分算法(或者其他类似的积分算法)推进到下一个节点时，都应该对下一个节点的ρ_η值作这个校正，这样做可以避免误差的累积。In fact, when advancing from the current node to the next node according to the Runge-Kutter integral algorithm (or other similar integral algorithms), this correction should be made to the ρ _η value of the next node, which can avoid the accumulation of errors.

7、追踪中的变步长技术：7. Variable step size technology in tracking:

在很多已公开的射线追踪算法中，对(1)式求解时都是采用固定步长，当沿着射线路径电离层的电子浓度变化较大时(如在射线的反射点处)，固定步长容易导致射线路径和波矢方向的误差都较大，在严重时会出现折射率平方变负的情形(在A-H公式中，当波的频率小于电离层的等离子体频率时，波沿着某些方向是截止的，沿着这些不能传播的方向折射率的平方是负的。)为了解决这个问题，本发明方法引入了变步长技术，即在对微分方向组(1)的积分过程中，根据媒质的梯度自动调整步长：In many published ray tracing algorithms, a fixed step size is used when solving equation (1). It is easy to cause large errors in the ray path and wave vector direction, and in severe cases, the square of the refractive index will become negative (in the A-H formula, when the frequency of the wave is lower than the plasma frequency of the ionosphere, the wave along a certain These directions are cut-off, and the square of the refractive index along these directions that cannot propagate is negative.) In order to solve this problem, the method of the present invention introduces a variable step size technique, that is, in the integration process of the differential direction group (1) , automatically adjust the step size according to the gradient of the medium:

10、10.

$step step = = \frac{kn k n}{| | \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; r r} + + \frac{11}{r r} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; θ θ} + + \frac{11}{r r sin sin θ θ} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; φ φ} | |} - - - - - - ((3030))$

式中，k为一个比例常数，k的值一般根据经验和追踪的精度来确定，在实际操作时，可以通过反复运行射线追踪程序来选择一个合适的值(一般介于10^-3和10^-5之间)。当折射率的沿着射线路径的梯度较小时，即折射率的变化较平缓，则步长较大；反之，当折射率沿着射线路径剧烈变化时，步长自动减小；In the formula, k is a proportional constant. The value of k is generally determined according to experience and tracking accuracy. In actual operation, a suitable value can be selected by repeatedly running the ray tracing program (generally between 10 ^-3 and 10 ^{- 5} ). When the gradient of the refractive index along the ray path is small, that is, the change of the refractive index is relatively gentle, the step size is large; on the contrary, when the refractive index changes sharply along the ray path, the step size is automatically reduced;

8、追踪的截止条件：8. Cut-off conditions for tracking:

当追踪到电离层的F2层峰值高度(一般为300～350km)时则停止追踪。F2层的峰值高度的等离子体频率是一个极大值，这里的电动力学特性复杂，各种类型的等离子体不稳定性会激发出传播特性各异的电磁波。本发明方法中将F2层的峰值高度设为追踪的截止条件是合理的。Stop tracking when the peak height of the F2 layer of the ionosphere (generally 300-350 km) is tracked. The plasma frequency at the peak height of the F2 layer is a maximum value. The electrodynamic characteristics here are complex, and various types of plasma instabilities will excite electromagnetic waves with different propagation characteristics. In the method of the present invention, it is reasonable to set the peak height of the F2 layer as the cut-off condition for tracking.

9、地面波源或者震中的定位：9. Location of ground wave source or epicenter:

在F2层停止追踪后，沿着截止点的磁力线映射到地面，则该点即为震中或者地面的波源。大量的理论计算表明，地震震中区对电离层的影响是通过磁力线来映射的。反映磁力线的倾斜度的角度为磁倾角I，和地磁纬度Λ的关系为tanI＝2tanΛ，纬度越大，磁倾角也越大，意味着F2层的追踪截止点和震中区的偏移越小。而在低纬区或者赤道区，由于磁倾角较小，F2层的追踪截止点可能和实际的震中区有较大的偏移。关于这一课题的详细讨论可参考S.Pulinets和K.Boyarchuk合著的专著IonosphericPrecursors of Earthquakes，本发明方法的很多理论依据均取自其中，这里不再赘述。After the F2 layer stops tracking, the magnetic line along the cut-off point is mapped to the ground, and this point is the epicenter or the wave source of the ground. A large number of theoretical calculations have shown that the impact of the earthquake epicenter on the ionosphere is mapped through magnetic field lines. The angle reflecting the inclination of the magnetic field line is the magnetic inclination I, and the relationship with the geomagnetic latitude Λ is tanI=2tanΛ, the greater the latitude, the greater the magnetic inclination, which means that the offset between the tracking cut-off point of the F2 layer and the epicentral area is smaller. However, in low latitudes or equatorial regions, due to the small magnetic inclination, the tracking cut-off point of the F2 layer may have a large offset from the actual epicentral region. The detailed discussion about this topic can refer to the monograph Ionospheric Precursors of Earthquakes co-authored by S.Pulinets and K.Boyarchuk, and many theoretical basis of the inventive method are all taken from wherein, do not repeat them here.

综上所述，本发明一种电离层中的地震电磁波波源定位与后向追踪方法，其包括步骤：In summary, the present invention provides a seismic electromagnetic wave source location and backward tracking method in the ionosphere, which includes the steps of:

(a)描述波的传播路径的方程是Haselgrove方程，将该方程改写为适合于数值求解的形式，关键是求出折射率对各个变量的一阶偏导数。这一步属于理论推导工作。(a) The equation describing the propagation path of the wave is the Haselgrove equation. The equation is rewritten into a form suitable for numerical solution. The key is to obtain the first-order partial derivative of the refractive index with respect to each variable. This step belongs to the theoretical derivation work.

(b)选择合适的数值方法来对射线方程求解，推导出具体的迭代格式。这一步也属于理论分析工作。(b) Choose an appropriate numerical method to solve the ray equation and derive a specific iterative format. This step also belongs to theoretical analysis work.

(c)根据前几步的理论推导结果，编制出计算电离层中波传播的计算机软件。这一步属于计算机编程工作。(c) Compile computer software for calculating wave propagation in the ionosphere based on the theoretical derivation results of the previous steps. This step belongs to computer programming work.

(d)由矢量天线的实测数据提取出波的初始传播方向(利用已有的、已公开的波矢分析算法)。(d) Extract the initial propagation direction of the wave from the measured data of the vector antenna (using the existing and published wave vector analysis algorithm).

(e)由于波的传播路径(即射线路径)具有可逆性，将接收点的初始方向反向(体现出“后向”追踪的意义)，作为软件计算的初始条件之一。(e) Since the propagation path of the wave (that is, the ray path) is reversible, the initial direction of the receiving point is reversed (reflecting the meaning of "backward" tracking) as one of the initial conditions for software calculation.

(f)设置软件计算的其他初始条件：接收点的位置、计算的容差(反映了追踪的精度)、追踪的截止条件(在电离层的F2层峰值高度)等。(f) Set other initial conditions calculated by the software: the position of the receiving point, the calculation tolerance (reflecting the tracking accuracy), the tracking cut-off condition (the peak height of the F2 layer in the ionosphere), etc.

(g)追踪的截止点就是辐射源在空间的波源。对于实际的地震震中预报问题，则可以从追踪的截止点沿着地球磁场的磁力线向地面投影，则该投影点就是本发明方法计算出的可能的震中位置。(g) The cut-off point of the trace is the wave source of the radiation source in space. For the actual earthquake epicenter prediction problem, the cut-off point of tracking can be projected to the ground along the magnetic field lines of the earth's magnetic field, and then the projected point is the possible epicenter position calculated by the method of the present invention.

本发明的算法已用Fortran语言实现。用户输入接收点的位置(经度、纬度和高度)、频率、初始波矢方向、追踪容差，以及电子浓度剖面，程序就可以利用本发明方法提出的算法计算出射线路径和路径上的波矢方向。另外，该算法也具有很灵活的通用性，即用户可以自定义电离层的电子浓度剖面(利用以后的电离层模型或者实测数据)。The algorithm of the present invention has been implemented in Fortran language. The user inputs the position (longitude, latitude and height), frequency, initial wave vector direction, tracking tolerance, and electron concentration profile of the receiving point, and the program can use the algorithm proposed by the method of the present invention to calculate the ray path and the wave vector on the path direction. In addition, the algorithm also has very flexible versatility, that is, the user can customize the electron concentration profile of the ionosphere (by using the ionosphere model or measured data in the future).

由该算法计算出的电离层射线路径，并根据该射线路径定位波源的位置，已得到了实测数据的验证。本发明的方法不仅仅可以应用在电离层中的辐射源定位，也可以应用在海洋中的声波波源的定位。The ionospheric ray path calculated by this algorithm, and the location of the wave source according to the ray path, have been verified by the measured data. The method of the invention can not only be applied to the positioning of radiation sources in the ionosphere, but also can be applied to the positioning of acoustic wave sources in the ocean.

见图2，给出了电离层射线追踪软件的流程图，其中箭头表示调用关系。See Figure 2, which shows the flow chart of the ionospheric ray tracing software, where the arrows indicate the calling relationship.

Claims

1. Seismic electromagnetic wave source location and backward tracing method in ionosphere, it is characterized in that, utilize the incoming wave direction of receiving electromagnetic wave, the characteristic of frequency, propagation medium, carry out the ray tracing of ionosphere, carry out numerical solution to ray equation , to locate the radiation source; the steps include:

A) Theoretical derivation:

(a) The equation describing the propagation path of the wave is the Haselgrove equation, which is rewritten into a form suitable for numerical solution, and the first-order partial derivative of the refractive index to each variable is obtained;

(b) Select an appropriate numerical method to solve the ray equation, and derive a specific iterative format;

B) Computer programming, i.e. ionospheric ray tracing software:

(a) According to the theoretical derivation result of step A), the computer software for calculating the wave propagation in the ionosphere is compiled;

(b) Extract the initial propagation direction of the wave from the measured data of the vector antenna;

(c) Since the propagation path of the wave is reversible, the initial direction of the receiving point is reversed, and "backward" tracking is performed as one of the initial conditions for software calculation;

(d) Set other initial conditions calculated by the software: the position of the receiving point, the tolerance of the calculation, and the cut-off condition of the tracking;

C) Determine the location of the epicenter:

The cut-off point of the tracking is the wave source of the radiation source in space, and the epicenter position is determined by the cut-off point.

2. method as claimed in claim 1, it is characterized in that, described B) step, directly carry out numerical solution to the three-dimensional Haselegrove equation, to avoid the impact of the dispersive property of earth's magnetic field and ionosphere on wave propagation path, make path more precise.

3. The method according to claim 1, characterized in that, in the step B), when the current node advances to the next node according to the integration algorithm, the node value is corrected to avoid the accumulation of errors.

4. The method according to claim 1, characterized in that, said B) step adopts variable step size technology, which fully considers the spatial gradient of the refractive index of the medium, and steps in places where the gradient is large If the length is small, the step is large where the gradient is small to ensure tracking accuracy.

5. The method according to claim 1, characterized in that, said B) the tracking cut-off condition in the step is the F2 layer peak height in the ionosphere.

6. The method according to claim 1, characterized in that, in said C) step, using the cut-off point to determine the epicenter position is to project the cut-off point to the ground along the magnetic force lines of the earth's magnetic field after tracking the cut-off point, The projected point is the location of the epicenter.

7. The method as claimed in claim 1, wherein the premise is that the vector antenna is used to obtain multi-component electromagnetic field data, and the frequency, incoming wave direction and polarization characteristic information of the wave are extracted according to the wave vector analysis algorithm.

8. The method according to claim 3, wherein the integral algorithm is a Runge-Kutter algorithm or an Adams-Moulton algorithm.

9. method as claimed in claim 4, is characterized in that, described variable step size technique is to be automatically adjusted step size according to the gradient of medium by following formula:

step step = = \frac{kn k n}{| | \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; r r} + + \frac{11}{r r} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; θ θ} + + \frac{11}{r r sin sin θ θ} \frac{&PartialD; &PartialD; n no}{&PartialD; &PartialD; φ φ} | |}

In the formula, k is a proportional constant. When the gradient of the refractive index along the ray path is small, that is, the change of the refractive index is relatively gentle, the step size is larger; otherwise, when the refractive index changes sharply along the ray path, the step size Length automatically decreases;

At the same time, set an upper limit step _max and a lower limit step _min for the step size. When the step size calculated according to the above formula exceeds this range, the upper limit step size or the lower limit step size is used as the actual step size value.

10. The method according to claim 7, wherein the wave vector analysis algorithm is the SVD technique of Means algorithm and wave spectrum matrix.