Broadband lightning very high frequency radiation source positioning method based on ESPRIT algorithm

Abstract

The invention discloses a broadband lightning very high frequency radiation source positioning method based on an ESPRIT algorithm, and the method comprises the steps: 1, constructing an orthogonal uniform linear very high frequency antenna array, and collecting broadband lightning very high frequency signals; 2, discrete Fourier transform is carried out on the broadband thunder and lightning very-high-frequency signals, the broadband thunder and lightning very-high-frequency signals are divided into J frequency points, and the frequency domain value of each frequency point signal is worked out; step 3, dividing the orthogonal uniform linear very high frequency antenna array into two sub-arrays in an x-axis direction and a y-axis direction, positioning each frequency point value of the signals in the two directions by using an ESPRIT algorithm, and respectively obtaining a calculation incident angle of each frequency point in the two directions; step 4, carrying out phase ambiguity resolution on the calculated incident angles of each frequency point in the two directions so as to obtain a real incident angle; and 5, calculating the azimuth angle and the elevation angle of the radiation source according to the real incident angle. According to the invention, rapid and refined positioning imaging of the lightning radiation source can be realized.

Description

Technical Field

[0001] The invention belongs to the technical field of signal processing, and in particular relates to a broadband lightning very high frequency radiation source positioning method based on an ESPRIT algorithm. Background Art

[0002] The precise positioning of lightning very high frequency (VHF) radiation sources can effectively help analyze the development process of lightning channels, especially the capture and positioning of some weak radiation processes. It is very beneficial to the analysis of the mechanism of lightning occurrence and development, and also promotes lightning protection and early warning.

[0003] The electromagnetic time reversal method (EMTR), multiple signal classification algorithm (MUSIC) and orthonormal propagator method (OPM) and other radiation source direction of arrival (DOA) estimation algorithms have been verified by experiments and still have high positioning accuracy under low signal-to-noise ratio, and have received extensive attention and research. Under the same conditions, they have higher positioning performance than traditional methods based on time of arrival (TOA) and interferometer methods, and can effectively capture more weak radiation sources.

[0004] The limitation of the above-mentioned high-precision positioning algorithms is that, in order to obtain a refined lightning development positioning channel, it is necessary to continuously use a sliding window to extract data, and then for the positioning of each sliding window data, a global maximum search is required, which leads to a slow calculation speed and a long time for the overall positioning process of the algorithm. Experiments have shown that for a 500ms lightning VHF data positioning, the above-mentioned high-precision positioning algorithms need to consume hundreds of hours to obtain a more refined lightning development channel map, which is less timely and has limited application value for practical research. At present, in the field of lightning VHF radiation source positioning, DOA estimation methods based on array signal processing technology are rarely used, especially methods based on direct solution for rapid and refined DOA estimation, such as the estimation signal parameter via rotational invariance techniques (ESPRTI) algorithm, which are rarely reported. Summary of the invention

[0005] The purpose of the present invention is to provide a broadband lightning very high frequency radiation source positioning method based on the ESPRIT algorithm, and to apply the ESPRIT method based on array signal processing technology that can directly solve DOA parameters to the field of broadband lightning very high frequency radiation source positioning to achieve rapid and refined positioning.

[0006] In order to achieve the above object, the technical solution adopted by the present invention is as follows:

[0007] A broadband lightning very high frequency radiation source location method based on ESPRIT algorithm comprises the following steps:

[0008] Step 1: construct an orthogonal uniform linear VHF antenna array to collect broadband lightning VHF signals;

[0009] Step 2: Perform discrete Fourier transform on the broadband lightning VHF signal and divide the broadband lightning VHF signal into J frequency points. The frequency domain expression of each frequency point signal is as follows:

[0010] X(f j )=A(f j )S(f j )+N(f j ),j∈[1,J]

[0011] Among them, f j represents the jth frequency point, X(f j ),S(f j ),N(f j ) represent the frequency domain expressions of the received signal, radiation source and noise after Fourier transformation, respectively, and A represents the manifold vector of the array antenna;

[0012] Step 3: Divide the orthogonal uniform linear VHF antenna array into two sub-arrays in the x-axis direction and the y-axis direction, and use the ESPRIT algorithm to locate the frequency domain signals of each frequency point in the x-axis direction and the y-axis direction, and obtain the calculated incident angle v of the jth frequency point in the two directions respectively. j and u j ,

[0013] Assume that the azimuth and elevation angles of the far-field broadband incident signal are θ, then the calculated incident angle v for each frequency point in the x-axis direction j and The relationship of θ is expressed as:

[0014]

[0015] Similarly, the calculated incident angle u for each frequency point in the y-axis direction is j and The relationship of θ is expressed as:

[0016]

[0017] Among them, J represents the number of frequency points;

[0018] Step 4: Perform phase deambiguation on the calculated incident angles of each frequency point in both the x-axis and y-axis directions to obtain the true incident angle, specifically:

[0019] In the phase difference calculation process of the ESPRIT algorithm, phase ambiguity will occur because the phase difference is calculated by the phase angle. It can be seen that:

[0020]

[0021] Among them, arg() represents the phase calculation of the complex number, z represents the complex number result calculated by the ESPRIT algorithm, represents the calculated phase difference, while φ represents the actual phase difference, and there is a phase difference of several 2π between the two;

[0022] For each different frequency point f j , the calculated phase difference arg(z j ) satisfies the following relationship:

[0023]

[0024] c represents the speed of light, d represents the spacing of the array antenna, and α j represents the unresolved ambiguous incident angle of the jth frequency point. Due to the phase ambiguity problem, the α calculated by the above formula j Not the same, only when the phase difference of each frequency point compensates the correct phase ambiguity term, all α j It should be consistent, that is, as follows:

[0025]

[0026] Define a slope value sl as follows:

[0027]

[0028] Define a slope difference, the expression is as follows:

[0029]

[0030] When the phase difference of each frequency point correctly compensates the fuzzy phase, the overall slope difference will reach the minimum value, that is:

[0031]

[0032] At this time, the deblurred incident angle α' after compensationj It is expressed as:

[0033]

[0034] For broadband signals, the incident angle is determined by the following multi-frequency deambiguation incident angle weighted average method:

[0035]

[0036] Among them, ρ j represents the maximum eigenvalue of the signal subspace of the jth frequency point, which represents the strength of the signal; in the two-dimensional positioning of the present invention, v j and u j represent the calculated incident angle of the jth frequency point in the x-axis and y-axis directions, respectively. However, this angle has not been defuzzified, so for v j and u j , respectively use the method in step 4 to perform the defuzzification operation, that is, let α j Equal to v j and u j , and then perform phase deambiguation to finally obtain the calculated incident angles v and u of the broadband signal in the x-axis and y-axis directions;

[0037] Step 5: Calculate the azimuth and elevation angles of the radiation source from the actual incident angles in two directions, specifically:

[0038]

[0039] Furthermore, the frequency domain signals of each frequency point on the x-axis and y-axis in step 3 are located using the ESPRIT algorithm to obtain the calculated incident angle v of the jth frequency point in the two directions of the two axes. j and u j , specifically:

[0040] For a given uniform linear array of M antennas, assuming that there are N narrowband radiation sources in the far field transmitting to the array antenna, the signal received by the array antenna is expressed as follows:

[0041] X(t)=AS(t)+N(t)

[0042] Where X(t), S(t), N(t) represent the arrival signal, radiation source signal and noise received by the array antenna respectively; A represents the array manifold vector;

[0043] A=[a(β1),a(β2),...,a(β N )]

[0044]

[0045] βn represents the incident angle of the nth radiation source, λ n represents the wavelength of the signal, and d represents the spacing of the array antenna;

[0046] make Then A in the above formula can be rewritten as:

[0047]

[0048] Divide this linear array into two sub-arrays X1 and X2, which are the 1st to M-1st and 2nd to Mth items of array X respectively. Then the signals received by the two sub-arrays are expressed as follows:

[0049] X1(t)=A1S(t)+N1(t)

[0050] X2(t)=A2S(t)+N2(t)

[0051] Therefore, the following relationship exists:

[0052] A2=A1Φ

[0053]

[0054] The above formula shows that the manifold vectors of the two linear subarrays differ by only a fixed phase difference, and this phase difference represents the direction of the incoming wave of the radiation source. Therefore, the DOA parameter of the radiation source can be obtained by solving Φ.

[0055] The covariance matrix of the signal received by the array antenna is expressed as follows:

[0056] R X =E[XX H ]=AR S A H +R N

[0057] Among them, R X and R S Represent the covariance matrix of the received signal and the radiation source respectively, H represents the conjugate transpose, and R X Perform eigenvalue decomposition to obtain the following formula:

[0058]

[0059] Among them U S and U N Represent the signal subspace of the radiation source and the noise respectively. Since the signal subspace and the space spanned by the array manifold vector are the same, there exists a non-singular matrix T that satisfies the following formula:

[0060] U S =AT

[0061] The corresponding expressions for the two subarrays are written as follows:

[0062] U1=A1T

[0063] U2=A2T=A1ΦT

[0064] Therefore, we can get:

[0065] U1T -1 ΦT=A1ΦT=U2

[0066] Let Ψ = T -1 ΦT, we get:

[0067] U1Ψ=U2

[0068] By calculating Ψ and then performing eigenvalue decomposition on it, we can obtain its diagonal matrix value Φ, thereby solving the DOA parameter of the radiation source. The incident angle of the radiation source is obtained by the following formula:

[0069]

[0070] For a single radiation source, let N = 1, that is, only β1 is taken as the incident angle of the radiation source;

[0071] Calculate the incident angle v in the x-axis direction j When the M antennas of the uniform linear array are replaced by the antennas on the x-axis, X(t) is replaced by X(f j ) is used as the replacement, and the remaining calculation steps are the same as described above.

[0072]

[0073] Among them, d x Indicates the spacing between antennas on the x-axis.

[0074] Similarly, when calculating the incident angle u in the y-axis direction j When the M antennas of the uniform linear array are replaced by the antennas on the y-axis, X(t) is replaced by X(f j ) is used as the replacement, and the remaining calculation steps are the same as described above.

[0075]

[0076] Among them, d y Indicates the spacing of the antennas on the y-axis.

[0077] Furthermore, the orthogonal uniform linear VHF antenna array in step 1 is an L-shaped uniform linear VHF antenna array, which has 7 antennas in total, one located at the origin of the coordinate system, and the other 3 antennas are evenly distributed on the orthogonal x-axis and y-axis.

[0078] Furthermore, the interval between adjacent antennas on the x-axis and the y-axis is 9 m.

[0079] Compared with the prior art, the present invention has the following effects:

[0080] (1) The present invention applies the ESPRIT algorithm to the field of broadband lightning very high frequency radiation source positioning, which can realize rapid and refined positioning and imaging of lightning radiation sources.

[0081] (2) The present invention provides a fast and efficient method for resolving multi-frequency phase difference ambiguity, which can achieve accurate estimation of the actual phase difference. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 is a schematic diagram of a uniform linear array antenna;

[0083] Figure 2 is a schematic diagram of an L-shaped uniform linear lightning VHF antenna array used in the present invention;

[0084] Figure 3 This is a positioning imaging result diagram of the method of the present invention on the very high frequency data of an artificial lightning induction;

[0085] Figure 4 This is the positioning imaging result of the broadband interferometer positioning algorithm on the very high frequency data of the same artificial lightning. DETAILED DESCRIPTION

[0086] The implementation of the present invention is described in detail below in conjunction with specific embodiments and drawings.

[0087] A broadband lightning very high frequency radiation source location method based on ESPRIT algorithm comprises the following steps:

[0088] Step 1: construct an orthogonal uniform linear VHF antenna array to collect broadband lightning VHF signals;

[0089] Step 2: Perform discrete Fourier transform on the broadband lightning VHF signal and divide the broadband lightning VHF signal into J frequency points. The frequency domain expression of each frequency point signal is as follows:

[0090] X(f j )=A(f j )S(f j )+N(f j ),j∈[1,J]

[0091] Among them, X(f j ),S(f j ),N(f j ) represents the frequency domain expression of the Fourier transformed received signal, radiation source and noise. A represents the manifold vector of the array antenna.

[0092] Step 3: Divide the orthogonal uniform linear VHF antenna array into two sub-arrays in the x-axis direction and the y-axis direction, and use the ESPRIT algorithm to locate the frequency domain signals of each frequency point in the x-axis direction and the y-axis direction, and obtain the calculated incident angle v of the jth frequency point in the two directions respectively. j and u j , specifically:

[0093] Assume that the azimuth and elevation angles of the far-field broadband incident signal are θ, then the calculated incident angle v for each frequency point in the x-axis direction j and The relationship of θ is expressed as:

[0094]

[0095] Similarly, the calculated incident angle u for each frequency point in the y-axis direction is j and The relationship of θ is expressed as:

[0096]

[0097] Wherein, J represents the number of frequency points.

[0098] Step 4: Perform phase deambiguation on the calculated incident angles of each frequency point in both directions to obtain the true incident angle, specifically:

[0099] In the process of phase difference calculation, the ESPRIT algorithm will have phase ambiguity problem because the phase difference is calculated by phase angle. It can be seen that:

[0100]

[0101] Among them, arg() represents the phase calculation of the complex number, and z represents the complex number result calculated by the ESPRIT algorithm. represents the calculated phase difference, while φ represents the actual phase difference, and there is a phase difference of several 2π between the two;

[0102] In order to solve the phase ambiguity problem, the present invention proposes a minimum slope difference criterion. The specific principle is as follows: for each different frequency point f j , the calculated phase difference arg(z j ) satisfies the following relationship:

[0103]

[0104] c represents the speed of light, d represents the spacing of the array antenna, and α jrepresents the unresolved ambiguous incident angle of the jth frequency point. Due to the phase ambiguity problem, the α calculated by the above formula j Not the same, only when the phase difference of each frequency point compensates the correct phase ambiguity term, all α j It should be consistent, that is, as follows:

[0105]

[0106] Define a slope value sl as follows:

[0107]

[0108] Define a slope difference, the expression is as follows:

[0109]

[0110] When the phase difference of each frequency point correctly compensates the fuzzy phase, the overall slope difference will reach the minimum value, that is:

[0111]

[0112] At this time, the deblurred incident angle α' after compensation j It is expressed as:

[0113]

[0114] For broadband signals, the incident angle is determined by the following multi-frequency deambiguation incident angle weighted average method:

[0115]

[0116] Among them, ρ j represents the maximum eigenvalue of the signal subspace of the jth frequency point, which represents the strength of the signal; in the two-dimensional positioning of the present invention, v j and u j represent the calculated incident angle of the jth frequency point in the x-axis and y-axis directions, respectively. However, this angle has not been defuzzified, so for v j and u j , respectively use the method in step 4 to perform the defuzzification operation, that is, let α j Equal to v j and u j , and then perform phase deambiguation to finally obtain the calculated incident angles v and u of the broadband signal in the x-axis and y-axis directions;

[0117] Step 5: Calculate the azimuth and elevation angles of the radiation source from the actual incident angles in two directions, specifically:

[0118]

[0119] Preferably, the frequency domain signals of each frequency point on the x-axis and the y-axis in step 3 are positioned using the ESPRIT algorithm to obtain the calculated incident angle v of the jth frequency point in the two directions of the two axes respectively. j and u j , specifically:

[0120] like Figure 1 As shown in , for a given uniform linear array of M antennas, assuming that there are N narrowband radiation sources in the far field transmitting to the array antenna, the signal received by the array antenna is expressed as follows:

[0121] X(t)=AS(t)+N(t)

[0122] Where X(t), S(t), N(t) represent the arrival signal, radiation source signal and noise received by the array antenna respectively; A represents the array manifold vector;

[0123] A=[a(β1),a(β2),...,a(β N )]

[0124]

[0125] β n represents the incident angle of the nth radiation source, λ n represents the wavelength of the signal, and d represents the spacing of the array antenna;

[0126] make Then A in the above formula can be rewritten as:

[0127]

[0128] Divide this linear array into two sub-arrays X1 and X2, which are the 1st to M-1st and 2nd to Mth items of array X respectively. Then the signals received by the two sub-arrays are expressed as follows:

[0129] X1(t)=A1S(t)+N1(t)

[0130] X2(t)=A2S(t)+N2(t)

[0131] Therefore, the following relationship exists:

[0132] A2=A1Φ

[0133]

[0134] The above formula shows that the manifold vectors of the two linear subarrays differ by only a fixed phase difference, and this phase difference represents the direction of the incoming wave of the radiation source. Therefore, the DOA parameter of the radiation source can be obtained by solving Φ.

[0135] The covariance matrix of the signal received by the array antenna can be expressed as follows:

[0136] R X =E[XX H ]=AR S A H +R N

[0137] Among them, R X and R S Represent the covariance matrix of the received signal and the radiation source respectively, H represents the conjugate transpose, and R X Perform eigenvalue decomposition to obtain the following formula:

[0138]

[0139] Among them U S and U N Represent the signal subspace of the radiation source and the noise respectively. Since the signal subspace and the space spanned by the array manifold vector are the same, there exists a non-singular matrix T that satisfies the following formula:

[0140] U S =AT

[0141] The corresponding expressions of the two subarrays can be written as follows:

[0142] U1=A1T

[0143] U2=A2T=A1ΦT

[0144] Therefore, we can get:

[0145] U1T -1 ΦT=A1ΦT=U2

[0146] Let Ψ = T -1 ΦT, we can get:

[0147] U1Ψ=U2

[0148] By calculating Ψ and then performing eigenvalue decomposition on it, we can obtain its diagonal matrix value Φ, thereby solving the DOA parameter of the radiation source. The incident angle of the radiation source is obtained by the following formula:

[0149]

[0150] Calculate the incident angle v in the x-axis direction j When the M antennas of the uniform linear array are replaced by the antennas on the x-axis, X(t) is replaced by X(f j ) is used as the replacement, and the remaining calculation steps are the same as described above.

[0151]

[0152] Among them, d x Indicates the spacing between antennas on the x-axis.

[0153] Similarly, when calculating the incident angle u in the y-axis direction j When the M antennas of the uniform linear array are replaced by the antennas on the y-axis, X(t) is replaced by X(f j ) is used as the replacement, and the remaining calculation steps are the same as described above.

[0154]

[0155] Among them, d y Indicates the spacing of the antennas on the y-axis.

[0156] Preferably, combined Figure 2 The orthogonal uniform linear VHF antenna array in step 1 is an L-shaped uniform linear VHF antenna array, which has 7 antennas in total, one located at the origin of the coordinate system, and the other 3 evenly distributed on the orthogonal x-axis and y-axis, and the interval between adjacent antennas is 9m.

[0157] Example

[0158] The first step is to use a sliding window to extract lightning VHF data in sequence. The VHF data used in the present invention is Figure 2 The array antenna shown in the figure has a system sampling rate of 500 MHz and a sliding window length of 512 sampling points, i.e. 1.024 μs. The sliding window data is subjected to discrete Fourier transform to obtain the frequency point value.

[0159] In the second step, the 20 frequency points with the highest energy values ​​are selected and the above-mentioned broadband ESPRIT algorithm is used to calculate and obtain the final positioning imaging result.

[0160] Figure 3 and Figure 4 The positioning imaging results of the method of the present invention and the existing broadband interferometer method on the very high frequency data of an artificial lightning strike are respectively demonstrated. Figure 3For example, the figure is divided into two parts, the lower part is the elevation angle-azimuth coordinate diagram, the horizontal axis represents the azimuth, the angle range is [0,360°], the vertical axis represents the elevation angle, the angle range is [0,90°]. This sub-figure depicts the position and change of the radiation source in the two-dimensional coordinate system, and the color change indicates the time sequence, which corresponds to the color bar on the right; the upper part is the elevation angle-time coordinate diagram, the horizontal axis represents the time, the range is [0,500ms], and the vertical axis has the same meaning as the vertical axis in the lower part. This sub-figure depicts a trend of the elevation angle of the radiation source, that is, the height, over time. The number in the lower left box indicates the number of radiation source positioning points depicted in the figure. Figure 4 The subgraph structure and meaning in Figure 3 Be consistent.

[0161] Figure 3 The method of the present invention located 49,181 radiation sources, and Figure 4 The medium-wideband interferometer method located 30,256 radiation sources. In contrast, the lightning channel development map obtained by the method of the present invention is more detailed and contains more lightning development branches. This shows the superiority of the positioning performance of the present method. This section of VHF data is 500ms in total, and the method of the present invention takes a total of 30 minutes, which is comparable to the broadband interferometer method in terms of time.

[0162] The present invention combines the array signal processing technology ESPRIT algorithm with the lightning very high frequency radiation source positioning, which can directly solve the two-dimensional DOA parameters without the need for full-space maximum value search, and has a faster calculation speed while finely positioning the radiation source.

[0163] The above shows and describes the basic principles, main features and advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments. The above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention to be protected. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A broadband lightning VHF radiation source location method based on ESPRIT algorithm, characterized in that: The following steps are involved: Step 1: construct an orthogonal uniform linear VHF antenna array to collect broadband lightning VHF signals; Step 2: Perform discrete Fourier transform on the broadband lightning VHF signal and divide the broadband lightning VHF signal into J frequency points. The frequency domain expression of each frequency point signal is as follows: X(f j )=A(f j )S(f j )+N(f j ),j∈[1,J] Among them, f j represents the jth frequency point, X(f j ),S(f j ),N(f j ) represent the frequency domain expressions of the received signal, radiation source and noise after Fourier transformation, respectively, and A represents the manifold vector of the array antenna; Step 3: Divide the orthogonal uniform linear VHF antenna array into two sub-arrays in the x-axis direction and the y-axis direction, and use the ESPRIT algorithm to locate the frequency domain signals of each frequency point in the x-axis direction and the y-axis direction, and obtain the calculated incident angle v of the jth frequency point in the two directions respectively. j and u j , Assume that the azimuth and elevation angles of the far-field broadband incident signal are θ, then the calculated incident angle v for each frequency point in the x-axis direction j and The relationship of θ is expressed as: Similarly, the calculated incident angle u for each frequency point in the y-axis direction is j and The relationship of θ is expressed as: Among them, J represents the number of frequency points; Step 4: Perform phase deambiguation on the calculated incident angles of each frequency point in both the x-axis and y-axis directions to obtain the true incident angle, specifically: In the phase difference calculation process of the ESPRIT algorithm, phase ambiguity will occur because the phase difference is calculated by the phase angle. It can be seen that: Among them, arg() represents the phase calculation of the complex number, and z represents the complex number result calculated by the ESPRIT algorithm. represents the calculated phase difference, while φ represents the actual phase difference, and there is a phase difference of several 2π between the two; For each different frequency point f j , the calculated phase difference arg(z j ) satisfies the following relationship: c represents the speed of light, d represents the spacing of the array antenna, and α j represents the unresolved ambiguous incident angle of the jth frequency point. Due to the phase ambiguity problem, the α calculated by the above formula j Not the same, only when the phase difference of each frequency point compensates the correct phase ambiguity term, all α j It should be consistent, that is, as follows: Define a slope value sl as follows: Define a slope difference, the expression is as follows: When the phase difference of each frequency point correctly compensates the fuzzy phase, the overall slope difference will reach the minimum value, that is: At this time, the deblurred incident angle α' after compensation j It is expressed as: For broadband signals, the incident angle is determined by the following multi-frequency deambiguation incident angle weighted average method: Among them, ρ j represents the maximum eigenvalue of the signal subspace at the jth frequency point, representing the strength of the signal; v j and u j represent the calculated incident angle of the jth frequency point in the x-axis and y-axis directions, respectively. However, this angle has not been defuzzified, so for v j and u j , respectively use the method in step 4 to perform the defuzzification operation, that is, let α j Equal to v j and u j , and then perform phase deambiguation to finally obtain the calculated incident angles v and u of the broadband signal in the x-axis and y-axis directions; Step 5: Calculate the azimuth and elevation angles of the radiation source from the actual incident angles in two directions, specifically:

2. The broadband lightning VHF radiation source location method based on the ESPRIT algorithm according to claim 1 is characterized in that: In step 3, the frequency domain signals of each frequency point on the x-axis and y-axis are located using the ESPRIT algorithm to obtain the calculated incident angle v of the jth frequency point in both directions of the two axes. j and u j , specifically: For a given uniform linear array of M antennas, assuming that there are N narrowband radiation sources in the far field transmitting to the array antenna, the signal received by the array antenna is expressed as follows: X(t)=AS(t)+N(t) Where X(t), S(t), N(t) represent the arrival signal, radiation source signal and noise received by the array antenna respectively; A represents the array manifold vector; A=[a(β1),a(β2),...,a(β N )] β n represents the incident angle of the nth radiation source, λ n represents the wavelength of the signal, and d represents the spacing of the array antenna; make Then A in the above formula can be rewritten as: Divide this linear array into two sub-arrays X1 and X2, which are the 1st to M-1st and 2nd to Mth items of array X respectively. Then the signals received by the two sub-arrays are expressed as follows: X1(t)=A1S(t)+N1(t) X2(t)=A2S(t)+N2(t) Therefore, the following relationship exists: A2=A1Φ The above formula shows that the manifold vectors of the two linear subarrays differ by only a fixed phase difference, and this phase difference represents the direction of the incoming wave of the radiation source. Therefore, the DOA parameter of the radiation source can be obtained by solving Φ. The covariance matrix of the signal received by the array antenna is expressed as follows: R X =E[XX H ]=AR S A H +R N Among them, R X and R S Represent the covariance matrix of the received signal and the radiation source respectively, H represents the conjugate transpose, and R X Perform eigenvalue decomposition to obtain the following formula: Among them U S and U N Represent the signal subspace of the radiation source and the noise respectively. Since the signal subspace and the space spanned by the array manifold vector are the same, there exists a non-singular matrix T that satisfies the following formula: U S =AT The corresponding expressions for the two subarrays are written as follows: U1=A1T U2=A2T=A1ΦT Therefore, we can get: U1T -1 ΦT=A1ΦT=U2 Let Ψ = T -1 ΦT, we get: U1Ψ=U2 By calculating Ψ and then performing eigenvalue decomposition on it, we can obtain its diagonal matrix value Φ, thereby solving the DOA parameter of the radiation source. The incident angle of the radiation source is obtained by the following formula: Calculate the incident angle v in the x-axis direction j When the M antennas of the uniform linear array are replaced by the antennas on the x-axis, X(t) is replaced by X(f j ) to obtain: Among them, d x Indicates the spacing of antennas on the x-axis; Similarly, when calculating the incident angle u in the y-axis direction j When the M antennas of the uniform linear array are replaced by the antennas on the y-axis, X(t) is replaced by X(f j ) to obtain: Among them, d y Indicates the spacing of the antennas on the y-axis.

3. The broadband lightning VHF radiation source location method based on the ESPRIT algorithm according to claim 2 is characterized in that: The orthogonal uniform linear VHF antenna array in step 1 is an L-shaped uniform linear VHF antenna array, which has 7 antennas in total.

4. The broadband lightning VHF radiation source location method based on ESPRIT algorithm according to claim 3 is characterized in that: Among the 7 antennas, one is located at the origin of the coordinate system, and the other three are evenly distributed on the orthogonal x-axis and y-axis.

5. The broadband lightning VHF radiation source location method based on ESPRIT algorithm according to claim 3 is characterized in that: The distance between adjacent antennas on the x-axis and y-axis is 5-10m.

6. The broadband lightning VHF radiation source location method based on ESPRIT algorithm according to claim 5 is characterized in that: The distance between adjacent antennas on the x-axis and y-axis is 9m.

Images