[go: up one dir, main page]
More Web Proxy on the site http://driver.im/

CN115061083B - A method for estimating the direction of arrival of lightning signals based on energy distribution - Google Patents

A method for estimating the direction of arrival of lightning signals based on energy distribution Download PDF

Info

Publication number
CN115061083B
CN115061083B CN202210659040.1A CN202210659040A CN115061083B CN 115061083 B CN115061083 B CN 115061083B CN 202210659040 A CN202210659040 A CN 202210659040A CN 115061083 B CN115061083 B CN 115061083B
Authority
CN
China
Prior art keywords
frequency point
lightning signal
sub
band
critical frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210659040.1A
Other languages
Chinese (zh)
Other versions
CN115061083A (en
Inventor
陈建
李唯一
林琳
孙晓颖
燕学智
佴威至
王庆龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN202210659040.1A priority Critical patent/CN115061083B/en
Publication of CN115061083A publication Critical patent/CN115061083A/en
Application granted granted Critical
Publication of CN115061083B publication Critical patent/CN115061083B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/14—Systems for determining direction or deviation from predetermined direction

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种基于能量分布对雷电信号进行波达方向估计的方法，属于对雷电信号进行波达方向估计的方法。将频谱能量积累到90％处的频点确定为临界频率点，在临界频率点前，计算原偶数段子带的能量，根据各个原偶数子带的能量与其能量总和之比，确定新子带的长度，原偶数子带能量越高，新子带越长，能量越低，新子带越短；在临界频率点后，由于能量较低，直接选择原偶数段子带作为新子带进行处理，通过对各个子带进行波达方向估计，能够实现对雷电信号来波方向的精确估计。本发明优点是处理的数据较少，运算量相应减小，且由于有效带宽缩短，精度和成功率有提高，为雷电信号定位提供一种计算量较低、精度较高的波达方向估计方法。

The present invention relates to a method for estimating the direction of arrival of a lightning signal based on energy distribution, and belongs to a method for estimating the direction of arrival of a lightning signal. The frequency point where the spectrum energy accumulates to 90% is determined as a critical frequency point. Before the critical frequency point, the energy of the original even-numbered sub-band is calculated, and the length of the new sub-band is determined according to the ratio of the energy of each original even-numbered sub-band to its total energy. The higher the energy of the original even-numbered sub-band, the longer the new sub-band, and the lower the energy, the shorter the new sub-band; after the critical frequency point, due to the low energy, the original even-numbered sub-band is directly selected as the new sub-band for processing, and the direction of arrival of each sub-band is estimated, so that the accurate estimation of the direction of arrival of the lightning signal can be achieved. The advantages of the present invention are that the data to be processed is less, the amount of calculation is correspondingly reduced, and due to the shortening of the effective bandwidth, the accuracy and success rate are improved, and a method for estimating the direction of arrival with low calculation amount and high accuracy is provided for lightning signal positioning.

Description

Method for estimating direction of arrival of lightning signal based on energy distribution
Technical Field
The invention belongs to a method for estimating the direction of arrival of a lightning signal, and particularly relates to a method for improving a general broadband direction of arrival estimation method based on energy distribution in a frequency band to realize lightning signal positioning.
Background
The occurrence of lightning may generate many damages, and the physical effects of high temperature, violent shock wave, strong electromagnetic radiation and the like generated by the lightning cause the lightning to generate huge damage in the moment. The personal injury and death can be caused frequently, the building or the communication equipment is destroyed, forest fires are caused, the computer information system is interrupted, the combustion and even explosion of a warehouse, an oil refinery, an oil field and the like are caused, the property and the personal safety of people are endangered, and the threat to the vehicles such as aviation and aerospace is also great. The effective positioning method can provide valuable data for lightning monitoring and meteorological research, and also can provide important support for lightning early warning and disaster prevention and reduction of departments such as an aerospace system, an electric power system and the like. Lightning location is important for understanding lightning physics, and is helpful for revealing physical mechanism and understanding characteristic parameters of lightning discharge, thereby facilitating lightning protection.
The existing lightning positioning technology comprises a traditional magnetic direction finding technology, a very low frequency time arrival method, a very high frequency interferometer method and a positioning technology for processing a focus array signal. The magnetic direction measurement technology is easily affected by the topography of the measurement site or surrounding buildings, the requirements on the position of the antenna and the installation of the antenna are high, the requirements are difficult to meet in practical application, and the error is large. The principle of the time arrival method is that the lightning position is positioned by utilizing the time difference generated when the lightning arrives at different measuring stations, and single-station information has no use value and larger positioning error. The interference technology utilizes a plurality of periods of signals, the average time difference in the whole integration period is determined through the phase difference, the interferometer is divided into a narrow-band interferometer and a wide-band interferometer, the discharge process of thunder and lightning can be known while positioning, but the detection distance of the interference technology is shorter, the commercialization is difficult, more applications are realized in scientific research work, and less applications are realized in actual life. The lightning positioning method focusing on the array signal processing technology eliminates the limitation of the number of the antennas on the signal processing, and can improve the positioning precision by increasing the number of the antennas.
The method for estimating the direction of arrival is used as a main branch of array signal processing, and the principle is that signals are made to be incident on an antenna array which is arranged according to a certain rule, the spatial position of an incident signal source is determined by a certain method, and the method for estimating the direction of arrival can further improve the positioning precision and the noise suppression capability by increasing the number of antennas. The conventional wideband MUSIC method utilizes orthogonality of signal and noise subspace generated by covariance matrix decomposition to construct a signal space spectrum function, and obtains the space parameters of the incident signal by searching peak values. The method has been applied to lightning localization and proved to have high resolution and powerful sidelobe suppression capability. The traditional broadband MUSIC method divides the signal into equal-length subsections in the frequency domain, and estimates each subsection respectively to obtain the average value of the spatial spectrum function. When the sub-segments are segmented, equal-length segmentation is performed only according to the frequency range, and energy distribution in an actual signal frequency band is not considered. On one hand, the bandwidth of the lightning signal is wider, the traditional MUSIC method needs more segments and has larger calculated amount, and the accuracy is reduced along with the increase of the bandwidth, on the other hand, the spectrum of the lightning signal is unevenly distributed and is exponentially distributed, and the idea of evenly dividing cannot fully consider the change of the energy in the frequency band. The lightning signal direction of arrival estimation method based on energy distribution has the advantages of small calculated amount and high precision, and has practical application value.
The China patent 'a VHF frequency band cloud lightning detection and positioning system' (application number 201110204820.9) discloses a VHF frequency band cloud lightning detection and positioning system. The system comprises a central processing station and a plurality of detection processing sub-stations which are arranged at different positions. The central processing station comprises a router, a firewall, a communication and monitoring server, a cloud lightning database, a TDOA calculation server and a cloud lightning positioning and resolving server. Each detection processing sub-station comprises a time frequency unit, a VHF frequency band array antenna, a VHF frequency band array receiving unit, a high-speed synchronous acquisition unit, a VHF frequency band array calibration signal and local oscillator generating unit, a detection processing sub-station digital processing sub-system and an embedded controller. The detection processing sub-station and the central processing station are connected into the public wide area network in a wireless or wired mode; and the obtained cloud lightning detection data are uploaded to the central processing station by the detection processing sub-station through the wide area network, and the central processing station can remotely monitor the working state of each detection processing sub-station. The system can accurately estimate the incident directions of a plurality of simultaneous lightning radiation sources. The signal direction of arrival estimation unit in the digital processing system of the detection processing substation uses eight paths of lightning waveform data frames from the data caching and control unit, and performs spatial spectrum estimation by using a multiple signal classification method to obtain the lightning signal direction of arrival corresponding to the data frames. The invention provides a complete set of technology aiming at key links such as acquisition, detection, signal arrival direction estimation, signal arrival time difference estimation and the like of cloud lightning signals, but the arrival direction estimation link adopts a traditional multiple signal classification method, does not consider the frequency distribution of the lightning signal energy, and has lower precision.
Disclosure of Invention
The invention provides a method for estimating the direction of arrival of a lightning signal based on energy distribution, which aims to solve the problems of large calculated amount and low precision and success rate.
The technical scheme adopted by the invention is that the method comprises the following steps:
step 1: determining critical frequency points of lightning signals
According to the U.S. military standard MILs-STD-464, the general description of the lightning signal model is:
i(t)＝I0(e-αt-e-βt)
wherein the relevant parameter I 0＝218810A,α＝11354s-1,β＝647265s-1;
Carrying out Fourier transform on the lightning signals conforming to the model, analyzing the lightning signals, and accumulating energy to a frequency point of 90% to be used as a critical frequency point f h; the expression that the critical frequency point f h needs to satisfy is:
Wherein V (j 2 pi f) is Fourier transform of the lightning signal, j is an imaginary unit, f is a frequency point in the lightning signal spectrum, and h frequency point f h is a critical frequency point at which the energy of the lightning signal spectrum is accumulated to 90%;
Step 2: obtaining array receiving data generated after lightning signals pass through the uniform linear array;
Step 3: dividing lightning signal array received data X (f g) into a plurality of new sub-bands;
first, based on energy distribution, determining new sub-band length and number before and after critical frequency point,
Dividing the lightning signal array receiving data X (f g) into a plurality of new sub-bands, wherein the calculation method comprises the following steps:
W＝W1+W2
Step 4: performing broadband direction of arrival estimation on the divided W new subbands;
Step 5: and averaging the spatial spectrum functions of each new sub-band, and searching spectrum peaks to obtain the estimation result of the arrival direction of the lightning signal.
In the step (2), P lightning signals conforming to the description of the model in the step1 are incident into a uniform linear array formed by M omnidirectional array elements, wherein the array element distance of the uniform linear array is d;
Setting the 1 st array element as a reference array element, and receiving data of a time domain array of the m-th array element of the lightning signal as follows:
Wherein i p (t) is the incident p-th lightning signal, n m (t) is Gaussian noise added to the m-th array element at the moment of t, and the time difference t mp between the reception of the p-th lightning signal by the m-th array element and the reception of the p-th lightning signal by the first reference array element is (m-1) dsin theta p/c;
Fourier transforming the time domain array receiving data of the lightning signal, wherein the frequency domain array receiving data of the lightning signal is:
Wherein f g is the g-th frequency point in the lightning signal spectrum, θ p is the incident angle of the p-th lightning signal, a m(θp,fg) is the guiding vector of the p-th lightning signal, I p(fg) the fourier transform of the p-th lightning signal I p (t), and N m(fg) is the fourier transform of gaussian noise N m (t);
writing frequency domain array receiving data of lightning signals into a matrix form:
X(fg)＝A(θ,fg)I(fg)+N(fg)
expanding each part of the matrix respectively into the following parts:
A(θ,fg)＝[a(θ1,fg),a(θ2,fg),...,a(θP,fg)]
I(fg)＝[I1(fg),I2(fg),...,IP(fg)]T
N(fg)＝[N1(fg),N2(fg),...,NM(fg)]T
X(fg)＝[X1(fg),X2(fg),...,XM(fg)]T。
the determination method in the step 3 (one) is as follows:
Firstly, according to the frequency range of lightning signals, the frequency bands of the lightning signals before and after a critical frequency point are divided averagely, and the number of sub-bands such as the front and the back of the critical frequency point is determined;
Secondly, regarding each even number sub-band in front of the critical frequency point as an original even number sub-band in front of the critical frequency point, determining the length of each new sub-band in front of the critical frequency point according to the ratio of the energy of each original even number sub-band to the sum of the energy of each original even number sub-band, wherein the number of the new sub-bands in front of the critical frequency point is the same as the number of the original even number sub-bands in front of the critical frequency point, thus determining the number and the length of the new sub-bands in front of the critical frequency point, and directly selecting each even number sub-band after the critical frequency point as the new sub-band after the critical frequency point, thus determining the number and the length of the new sub-bands after the critical frequency point;
Equally dividing the lightning signal array received data before and after a critical frequency point to generate original even sub-bands, determining a new sub-band before the critical frequency point according to the duty ratio of the energy of each original even sub-band in the energy summation before the critical frequency point, and directly selecting each even sub-band after the critical frequency point as the new sub-band after the critical frequency point.
Before a critical frequency point of a lightning signal, the method for determining the length and the number of the new sub-band comprises the following steps:
(1) The frequency spectrum of the lightning signal is analyzed, the frequency band before the critical frequency point is divided into equal-length sub-bands, the number of the equal-length sub-bands is Q, each even sub-band is selected as an original even sub-band before the critical frequency point, and the number of the original even sub-bands before the critical frequency point is Q 1;
the calculation method of the number Q 1 of the original even number sub-bands before the critical frequency point of the lightning signal comprises the following steps:
n is a positive integer
(2) The calculation method of the total subband length L 1 of Q 1 original even subbands with the frequency band length L before the critical frequency point of the lightning signal is as follows:
L1＝L/Q×Q1
Wherein L is the length of a frequency band before a lightning signal critical frequency point, Q 1 is the number of original even sub-bands before the lightning signal critical frequency point, and Q is the number of equal-length sub-bands divided before the lightning signal critical frequency point;
(3) The calculation method of the energy E q of the q-th original even sub-band of the lightning signal before the critical frequency point comprises the following steps:
Wherein V (j 2 pi f r) is Fourier transform of a lightning signal, j is an imaginary unit, f r is an r-th frequency point in a lightning signal spectrum chart, a starting frequency point of an original even-numbered sub-band of energy to be calculated is an r 1 -th frequency point in the lightning signal spectrum, and a terminating frequency point of the original even-numbered sub-band of energy to be calculated is an r 2 -th frequency point in the lightning signal spectrum;
the energy of each original even number sub-band before the critical frequency point of the lightning signal is determined as
(4) The energy of Q 1 original even sub-bands in front of the calculated critical frequency point of the lightning signal is summed to obtain total energy E, and the calculation method comprises the following steps:
Wherein E q is the energy of the Q-th original even number sub-band before the critical frequency point of the lightning signal, and Q 1 is the number of the original even number sub-bands before the critical frequency point of the lightning signal;
(5) The length C q of the q new sub-band of the lightning signal before the critical frequency point is calculated by the following steps:
Cq＝Eq/E×L1
wherein E q is the energy of the Q-th original even sub-band before the critical frequency point of the lightning signal, E is the energy sum of the Q 1 original even sub-bands before the critical frequency point of the lightning signal, and L 1 is the sub-band length sum of the Q 1 original even sub-bands;
(6) The calculation method of the number W 1 of the new sub-bands before the critical frequency point of the lightning signal is as follows:
W1＝Q1
the number of new subbands is the same as the number of original even subbands.
After the lightning signal is at a critical frequency point, the length and the number of the new sub-bands are determined as follows:
(1) The frequency spectrum of the lightning signal is analyzed, the frequency band after the critical frequency point is divided into equal-length sub-bands, the number of the equal-length sub-bands is G, each even sub-band is selected as an original even sub-band after the critical frequency point, and the number of the original even sub-bands after the critical frequency point is G 1;
The calculation method of the number G 1 of the original even number sub-bands after the critical frequency point of the lightning signal comprises the following steps:
n is a positive integer
(2) The length of the frequency band after the critical frequency point of the lightning signal is L';
(3) The calculation method of the length B k of the kth new sub-band of the lightning signal after the critical frequency point is as follows:
Bk＝L'/G
Wherein L' is the frequency band length of the lightning signal after the critical frequency point, and G is the number of equal-length sub-bands of the lightning signal after the critical frequency point;
(4) The calculation method of the number W 2 of the new sub-bands after the critical frequency point of the lightning signal is as follows:
W2＝G1
the number of new sub-bands after the critical frequency point is the same as the number of original even sub-bands after the critical frequency point.
In the step 3 and the step (two), according to the principle of dividing the new sub-bands based on energy distribution, the critical frequency point of the lightning signal is f h, before the critical frequency point f h, the lightning signal array receiving data X (f g) is divided into W 1 new sub-bands, the length of the q new sub-band is C q, and the center point of the original even sub-band before the critical frequency point is taken as the center point of the new sub-band before the critical frequency point; after the critical frequency point f h, dividing the lightning signal array received data X (f g) into W 2 new sub-bands, wherein the length of the kth new sub-band is B k, and taking the central point of the original even sub-band after the critical frequency point as the central point of the new sub-band after the critical frequency point, so that the lightning signal array received data X (f g) can be accurately divided, and the lightning array received data is divided into W new sub-bands.
In step 4 of the present invention, a rotated signal subspace RSS method is used for the divided new subbands, and the specific steps are as follows:
(1) The incidence angle of the lightning signal is estimated in advance by adopting a conventional beam forming method CBF, and an estimated angle eta is obtained;
(2) Dividing array receiving data of each new sub-band into K sections in a time domain, and carrying out Fourier transform on each new sub-band to obtain array receiving data of the new sub-band in a frequency domain;
(3) Dividing array received data of a new sub-band into N shorter sub-bands in a frequency domain, and calculating to obtain covariance matrixes of the array received data of the shorter sub-bands, wherein the covariance matrixes are as follows:
R(fn,η)＝E[X(fn)XH(fn)]
Wherein η is a pre-estimated incident angle of the lightning signal, f n is an nth frequency point in a shorter subband, and X (f n) is array receiving data of the lightning signal at a frequency point f n;
(4) Determining the reference frequency as f c, and reconstructing array flow pattern matrixes A (f n, eta) and A (f c, eta) according to the estimated angle eta
(5) The focus matrix T (f n) is determined by the following expression:
T(fn)＝V(fn)UH(fn)
Wherein U (f n) and V (f n) are left and right singular vectors of A (f n,η)AH(fc, eta), respectively, (. Cndot.) H represents a conjugate transpose;
(6) The array received data on each frequency point is focused on the reference frequency f c through focusing operation, and the covariance matrix R (f c, eta) after focusing is as follows:
T(fn)R(fn,η)TH(fn)≈R(fc,η)
(7) Performing eigenvalue decomposition on the covariance matrix R (f c, eta) after focusing to obtain a spatial spectrum function of the covariance matrix, and obtaining a spatial spectrum function of a w new subband The calculation method comprises the following steps:
Wherein, The value between [0,180] is taken at intervals of 0.01 to realize the search of angles between [0,180], f c is the reference frequency, a (f c) is the guide vector of the signal subspace, and Z N is the noise subspace.
The step of estimating the incidence angle of the lightning signal in the step (1) is as follows:
a. the incident lightning signal passes through the uniform linear array, and the received data of neglecting noise is as follows:
II(t)＝a(θ)i(t)
wherein a (theta) is a guiding vector of a lightning signal, and i (t) is the lightning signal;
b. the received signals are subjected to phase compensation and summation, and an output signal y (t) obtained after the lightning signal received data are subjected to beam forming is:
y(t)＝γH·II(t)
Wherein γ is the phase compensation vector, (. Cndot.) H represents the conjugate transpose;
c. the output power P CBF of the conventional beamformer is:
PCBF＝γH·R·γ
wherein R is a covariance matrix of data II (t) received after the lightning signal passes through the uniform linear array;
d. And searching the peak value of the output power, wherein the angle eta corresponding to the peak value is the estimated value of the incident angle of the lightning signal.
In step5 of the invention, the spatial spectrum function of each new sub-band is averaged, and the spatial spectrum function of the lightning signal is calculatedThe calculation method is as follows:
Wherein W is the number of new sub-bands, A space spectrum function of a w new sub-band in a lightning signal spectrum;
Spatial spectral function of lightning signal Searching spectral peaks, and searchingAngle corresponding to maximum valueThe lightning signal arrival direction estimation result is obtained.
The invention has the advantages that: the invention fully considers the distribution of energy in the lightning signal spectrum, determines the frequency point where the energy of the lightning signal spectrum is accumulated to 90% as a critical frequency point, calculates the energy of the original even-numbered sub-bands before the critical frequency point, and determines the length of a new sub-band according to the ratio of the energy of each original even-numbered sub-band to the sum of the energy of each original even-numbered sub-band, wherein the higher the energy of the original even-numbered sub-band is, the longer the new sub-band is, the lower the energy is, and the shorter the new sub-band is; after the critical frequency point, the original even-numbered segment sub-band is directly selected as a new sub-band for processing due to lower energy. By carrying out the arrival direction estimation on each new sub-band, the accurate estimation on the arrival direction of the lightning signal is realized. Compared with the traditional broadband MUSIC method, the data processed by the method is reduced by half, the operand is also reduced by half, the accuracy and the success rate are also improved compared with the traditional broadband MUSIC method, and the method has better lightning signal arrival direction estimation performance.
Drawings
FIG. 1 is a far field one-dimensional linear array DOA estimation geometry model diagram;
FIG. 2 is a flow chart of the present invention;
FIG. 3 is a schematic diagram of an original even number of sub-bands generated after equal-length division of lightning signals before and after a critical frequency point;
FIG. 4 is a schematic diagram of a new sub-band determined from energy before and after a critical frequency point for a lightning signal;
FIG. 5 is a graph showing root mean square error of direction of arrival estimation for lightning signals using the method of the present invention and the conventional method, respectively;
FIG. 6 is a graph showing the success rate of estimating the direction of arrival of lightning signals by the method of the present invention and the conventional method, respectively.
Detailed Description
The invention is described in detail below in connection with the flow chart of fig. 2.
Step 1: determining critical frequency points of lightning signals
According to the U.S. military standard MILs-STD-464, the general description of the lightning signal model is:
i(t)＝I0(e-αt-e-βt)
wherein the relevant parameter I 0＝218810A,α＝11354s-1,β＝647265s-1;
Carrying out Fourier transform on the lightning signals conforming to the model to obtain a spectrogram of the lightning signals shown in the curve of fig. 3, analyzing the spectrogram, and taking a frequency point at which energy is accumulated to 90% as a critical frequency point f h;
In fig. 3, f h is a critical frequency point, and the expression to be satisfied by the critical frequency point f h is:
Wherein V (j 2 pi f) is Fourier transform of the lightning signal, j is an imaginary unit, f is a frequency point in the lightning signal spectrum, and h frequency point f h is a critical frequency point at which the energy of the lightning signal spectrum is accumulated to 90%;
Step 2: obtaining array reception data generated after passing lightning signals through uniform linear array
P lightning signals conforming to the model description in the step 1 are incident into a uniform linear array formed by M omnidirectional array elements, wherein the array element spacing of the uniform linear array is d;
Setting the 1 st array element as a reference array element, and receiving data of a time domain array of the m-th array element of the lightning signal as follows:
Wherein i p (t) is the incident p-th lightning signal, n m (t) is Gaussian noise added to the m-th array element at the moment of t, and the time difference t mp between the reception of the p-th lightning signal by the m-th array element and the reception of the p-th lightning signal by the first reference array element is (m-1) dsin theta p/c;
Fourier transforming the time domain array receiving data of the lightning signal, wherein the frequency domain array receiving data of the lightning signal is:
Wherein f g is the g-th frequency point in the lightning signal spectrum, θ p is the incident angle of the p-th lightning signal, a m(θp,fg) is the guiding vector of the p-th lightning signal, I p(fg) the fourier transform of the p-th lightning signal I p (t), and N m(fg) is the fourier transform of gaussian noise N m (t);
writing frequency domain array receiving data of lightning signals into a matrix form:
X(fg)＝A(θ,fg)I(fg)+N(fg)
expanding each part of the matrix respectively into the following parts:
A(θ,fg)＝[a(θ1,fg),a(θ2,fg),...,a(θP,fg)]
I(fg)＝[I1(fg),I2(fg),...,IP(fg)]T
N(fg)＝[N1(fg),N2(fg),...,NM(fg)]T
X(fg)＝[X1(fg),X2(fg),...,XM(fg)]T
Step 3: dividing the lightning signal array received data X (f g) into a plurality of new sub-bands
Based on energy distribution, determining the length and the number of new sub-bands before and after a critical frequency point, wherein the determining method comprises the following steps:
Firstly, according to the frequency range of lightning signals, the frequency bands of the lightning signals before and after a critical frequency point are divided averagely, and the number of sub-bands such as the front and the back of the critical frequency point is determined;
Secondly, regarding each even number sub-band in front of the critical frequency point as an original even number sub-band in front of the critical frequency point, determining the length of each new sub-band in front of the critical frequency point according to the ratio of the energy of each original even number sub-band to the sum of the energy of each original even number sub-band, wherein the number of the new sub-bands in front of the critical frequency point is the same as the number of the original even number sub-bands in front of the critical frequency point, thus determining the number and the length of the new sub-bands in front of the critical frequency point, and directly selecting each even number sub-band after the critical frequency point as the new sub-band after the critical frequency point, thus determining the number and the length of the new sub-bands after the critical frequency point;
Equally dividing the received data of the lightning signal array before and after a critical frequency point, wherein the generated original even sub-bands are shown in a shadow part of fig. 3, determining a new sub-band before the critical frequency point according to the duty ratio of the energy of each original even sub-band in the energy summation of the original even sub-bands before the critical frequency point, directly selecting each even sub-band after the critical frequency point as a new sub-band after the critical frequency point, and dividing the received data of the lightning signal array into new sub-bands shown in a shadow part of fig. 4;
before a critical frequency point of a lightning signal, the length and the number of the new sub-band are determined as follows:
(1) The frequency spectrum of the lightning signal is analyzed, the frequency band before the critical frequency point is divided into equal-length sub-bands, the number of the equal-length sub-bands is Q, each even sub-band is selected as an original even sub-band before the critical frequency point, and the number of the original even sub-bands before the critical frequency point is Q 1;
the calculation method of the number Q 1 of the original even number sub-bands before the critical frequency point of the lightning signal comprises the following steps:
n is a positive integer
(2) The calculation method of the total subband length L 1 of Q 1 original even subbands with the frequency band length L before the critical frequency point of the lightning signal is as follows:
L1＝L/Q×Q1
Wherein L is the length of a frequency band before a lightning signal critical frequency point, Q 1 is the number of original even sub-bands before the lightning signal critical frequency point, and Q is the number of equal-length sub-bands divided before the lightning signal critical frequency point;
(3) The calculation method of the energy E q of the q-th original even sub-band of the lightning signal before the critical frequency point comprises the following steps:
Wherein V (j 2 pi f r) is Fourier transform of a lightning signal, j is an imaginary unit, f r is an r-th frequency point in a lightning signal spectrum chart, a starting frequency point of an original even-numbered sub-band of energy to be calculated is an r 1 -th frequency point in the lightning signal spectrum, and a terminating frequency point of the original even-numbered sub-band of energy to be calculated is an r 2 -th frequency point in the lightning signal spectrum;
the energy of each original even number sub-band before the critical frequency point of the lightning signal is determined as
(4) The energy of Q 1 original even sub-bands in front of the calculated critical frequency point of the lightning signal is summed to obtain total energy E, and the calculation method comprises the following steps:
Wherein E q is the energy of the Q-th original even number sub-band before the critical frequency point of the lightning signal, and Q 1 is the number of the original even number sub-bands before the critical frequency point of the lightning signal;
(5) The length C q of the q new sub-band of the lightning signal before the critical frequency point is calculated by the following steps:
Cq＝Eq/E×L1
wherein E q is the energy of the Q-th original even sub-band before the critical frequency point of the lightning signal, E is the energy sum of the Q 1 original even sub-bands before the critical frequency point of the lightning signal, and L 1 is the sub-band length sum of the Q 1 original even sub-bands;
(6) The calculation method of the number W 1 of the new sub-bands before the critical frequency point of the lightning signal is as follows:
W1＝Q1
The number of the new sub-bands is the same as that of the original even sub-bands;
After the lightning signal is at the critical frequency point, the length and the number of the new sub-bands are determined as follows:
(1) The frequency spectrum of the lightning signal is analyzed, the frequency band after the critical frequency point is divided into equal-length sub-bands, the number of the equal-length sub-bands is G, each even sub-band is selected as an original even sub-band after the critical frequency point, and the number of the original even sub-bands after the critical frequency point is G 1;
The calculation method of the number G 1 of the original even number sub-bands after the critical frequency point of the lightning signal comprises the following steps:
n is a positive integer
(2) The length of the frequency band after the critical frequency point of the lightning signal is L';
(3) The calculation method of the length B k of the kth new sub-band of the lightning signal after the critical frequency point is as follows:
Bk＝L'/G
Wherein L' is the frequency band length of the lightning signal after the critical frequency point, and G is the number of equal-length sub-bands of the lightning signal after the critical frequency point;
(4) The calculation method of the number W 2 of the new sub-bands after the critical frequency point of the lightning signal is as follows:
W2＝G1
the number of the new sub-bands after the critical frequency point is the same as the number of the original even sub-bands after the critical frequency point;
(II) dividing the lightning signal array received data X (f g) into a plurality of new sub-bands
According to the principle of dividing the new sub-bands based on energy distribution, the critical frequency point of the lightning signal is f h, before the critical frequency point f h, the lightning signal array receiving data X (f g) is divided into W 1 new sub-bands, the length of the q new sub-bands is C q, and the center point of the original even sub-bands before the critical frequency point is taken as the center point of the new sub-bands before the critical frequency point; after a critical frequency point f h, dividing the lightning signal array received data X (f g) into W 2 new sub-bands, wherein the length of the kth new sub-band is B k, and taking the central point of an original even sub-band after the critical frequency point as the central point of a new sub-band after the critical frequency point, so that the lightning signal array received data X (f g) can be accurately divided, and the lightning array received data is divided into W new sub-bands, and the calculation method comprises the following steps:
W＝W1+W2
step 4: wideband direction of arrival estimation for divided W new subbands
The method for using the rotating signal subspace RSS for the divided new subbands comprises the following specific steps:
(1) Pre-estimating the incident angle of the lightning signal by adopting a conventional beam forming method (Conventional Beamforming, CBF) to obtain an estimated angle eta;
The lightning signal incidence angle estimation method comprises the following steps:
a. the incident lightning signal passes through the uniform linear array, and the received data of neglecting noise is as follows:
II(t)＝a(θ)i(t)
wherein a (theta) is a guiding vector of a lightning signal, and i (t) is the lightning signal;
b. the received signals are subjected to phase compensation and summation, and an output signal y (t) obtained after the lightning signal received data are subjected to beam forming is:
y(t)＝γH·II(t)
Wherein γ is the phase compensation vector, (. Cndot.) H represents the conjugate transpose;
c. the output power P CBF of the conventional beamformer is:
PCBF＝γH·R·γ
wherein R is a covariance matrix of data II (t) received after the lightning signal passes through the uniform linear array;
d. Searching a peak value of the output power, wherein an angle eta corresponding to the peak value is an estimated value of the incident angle of the lightning signal;
(2) Dividing array receiving data of each new sub-band into K sections in a time domain, and carrying out Fourier transform on each new sub-band to obtain array receiving data of the new sub-band in a frequency domain;
(3) Dividing array received data of a new sub-band into N shorter sub-bands in a frequency domain, and calculating to obtain covariance matrixes of the array received data of the shorter sub-bands, wherein the covariance matrixes are as follows:
R(fn,η)＝E[X(fn)XH(fn)]
Where η is a pre-estimated incident angle of the lightning signal, f n is the nth frequency point in the shorter subband, and X (f n) is the array received data of the lightning signal at frequency point f n.
(4) Determining the reference frequency as f c, and reconstructing array flow pattern matrixes A (f n, eta) and A (f c, eta) according to the estimated angle eta
(5) The focus matrix T (f n) is determined by the following expression:
T(fn)＝V(fn)UH(fn)
Wherein U (f n) and V (f n) are left and right singular vectors of A (f n,η)AH(fc, eta), respectively, (. Cndot.) H represents the conjugate transpose.
(6) The array received data on each frequency point is focused on the reference frequency f c through focusing operation, and the covariance matrix R (f c, eta) after focusing is as follows:
T(fn)R(fn,η)TH(fn)≈R(fc,η)
(7) Performing eigenvalue decomposition on the covariance matrix R (f c, eta) after focusing to obtain a spatial spectrum function of the covariance matrix, and obtaining a spatial spectrum function of a w new subband The calculation method comprises the following steps:
Wherein, Taking a value between [0,180] at intervals of 0.01 to realize the search of angles between [0,180], wherein f c is a reference frequency, a (f c) is a guide vector of a signal subspace, and Z N is a noise subspace;
step 5: averaging the spatial spectrum functions of each new sub-band, and searching through spectrum peaks to obtain the estimation result of the arrival direction of the lightning signal;
averaging the spatial spectrum functions of each new subband and the spatial spectrum functions of the lightning signals The calculation method is as follows:
Wherein W is the number of new sub-bands, A space spectrum function of a w new sub-band in a lightning signal spectrum;
Spatial spectral function of lightning signal Searching spectral peaks, and searchingAngle corresponding to maximum valueThe lightning signal arrival direction estimation result is obtained.
According to the invention, the energy change in the lightning signal spectrum is fully considered, the frequency point where 90% of the lightning signal spectrum energy is accumulated is taken as the critical frequency point according to the characteristics of the lightning signal, and the signals are processed differently before and after the critical frequency point. Before a critical frequency point, calculating the energy of an original even-numbered section sub-band, and determining the length of a new sub-band according to the ratio of the energy of each original even-numbered sub-band to the sum of the energy of each original even-numbered sub-band, wherein the higher the energy of the original even-numbered sub-band is, the longer the new sub-band is, the lower the energy is, and the shorter the new sub-band is; after the critical frequency point, the original even-numbered segment sub-band is directly selected as a new sub-band for processing due to lower energy. And using a general broadband direction-of-arrival estimation method for each new sub-band, and averaging the obtained results to obtain a direction-of-arrival estimation result of the lightning signal.
Compared with the existing method for estimating the direction of arrival of the lightning signal, the method has the advantages that the operation amount is greatly reduced, and the accuracy and the success rate are improved.
The experiment uses an i7-7500U processor, the main memory is 8.00GB, and the simulation platform is Matlab 2016a. To ensure accuracy of the experiment, 200 monte carlo experiments were performed.
FIG. 5 shows the root mean square error of the estimated angle compared with the true angle obtained by estimating the direction of arrival of the lightning signal by the method of the present invention and the conventional method within a certain signal-to-noise ratio after 200 Monte Carlo experiments. As shown in fig. 5, the accuracy of the conventional method is slightly higher than that of the method of the present invention at a low signal-to-noise ratio, but the accuracy of the method of the present invention is significantly improved as the signal-to-noise ratio is improved.
Fig. 6 shows the direction of arrival estimation of lightning signals in the range of a certain signal to noise ratio after 200 monte carlo experiments, and the success rate of the estimation is compared with that of the conventional method. As shown in fig. 6, the success rate of the method of the present invention is higher than that of the conventional method.
Through experimental measurement, the time required for estimating the direction of arrival of 2 lightning signals by the traditional method is 25.87s, the time required for estimating the direction of arrival of 2 lightning signals by the method is 12.48s, and the operation time is greatly reduced.
In summary, the method of the invention not only improves the accuracy of estimation, but also improves the success rate of estimation and reduces the calculated amount.
The foregoing description of the preferred embodiments of the invention is not intended to limit the invention to the particular embodiments disclosed, but on the contrary, the intention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention.

Claims (9)

1.一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于，包括下列步骤：1. A method for estimating the direction of arrival of a lightning signal based on energy distribution, comprising the following steps: 步骤1：确定雷电信号的临界频率点Step 1: Determine the critical frequency point of the lightning signal 根据美国军事标准MIL-STD-464，雷电信号模型的一般描述为:According to the US military standard MIL-STD-464, the general description of the lightning signal model is: i(t)＝I0(e-αt-e-βt)i(t)＝I 0 (e -αt -e -βt ) 其中，相关参数I0＝218810A，α＝11354s-1，β＝647265s-1；Among them, the relevant parameters are I 0 =218810A, α=11354s -1 , β=647265s -1 ; 对符合上述模型的雷电信号做傅里叶变换，对其进行分析，将能量积累到90％的频率点作为临界频率点fh；临界频率点fh需要满足的表达式为：Perform Fourier transform on the lightning signal that meets the above model, analyze it, and take the frequency point where the energy accumulates to 90% as the critical frequency point f h ; the expression that the critical frequency point f h needs to satisfy is: 其中，V(j2πf)是雷电信号的傅里叶变换，j为虚数单位，f为雷电信号频谱中的频率点，第h个频率点fh为雷电信号频谱能量累积到90％的临界频率点；Where V(j2πf) is the Fourier transform of the lightning signal, j is an imaginary unit, f is the frequency point in the lightning signal spectrum, and the hth frequency point fh is the critical frequency point where the lightning signal spectrum energy accumulates to 90%; 步骤2：获得雷电信号通过均匀线性阵列后产生的阵列接收数据；Step 2: Obtain array receiving data generated after the lightning signal passes through the uniform linear array; 步骤3：将雷电信号阵列接收数据X(fg)划分为若干个新子带Step 3: Divide the lightning signal array received data X(f g ) into several new sub-bands (一)、基于能量分布，确定临界频率点前后新子带长度和数量，(I) Based on the energy distribution, determine the length and number of new sub-bands before and after the critical frequency point. (二)、将雷电信号阵列接收数据X(fg)划分为若干个新子带，计算方法为：(ii) Divide the lightning signal array receiving data X(f g ) into several new sub-bands, and the calculation method is: W＝W1+W2 W＝W 1 +W 2 步骤4：对划分的W个新子带进行宽带波达方向估计；Step 4: perform broadband DOA estimation on the W new sub-bands; 步骤5：对各个新子带的空间谱函数取平均，通过谱峰搜索获取雷电信号的波达方向估计结果。Step 5: Average the spatial spectrum functions of each new sub-band and obtain the direction of arrival estimation result of the lightning signal through spectrum peak search. 2.根据权利要求1所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：所述步骤(2)中，将符合步骤1模型描述的P个雷电信号入射到由M个全向阵元组成的均匀线性阵列中，均匀线阵的阵元间距为d；2. The method for estimating the direction of arrival of lightning signals based on energy distribution according to claim 1, characterized in that: in the step (2), the P lightning signals conforming to the model description of step 1 are incident on a uniform linear array composed of M omnidirectional array elements, and the array element spacing of the uniform linear array is d; 设置第1个阵元为参考阵元，则雷电信号的第m个阵元的时域阵列接收数据为：Set the first array element as the reference array element, then the time domain array receiving data of the mth array element of the lightning signal is: 其中，ip(t)为入射的第p个雷电信号，nm(t)为t时刻加到第m个阵元上的高斯噪声，第m个阵元接收第p个雷电信号与第一个参考阵元接收第p个雷电信号的时间差tmp为(m-1)dsinθp/c；Wherein, i p (t) is the incident p-th lightning signal, nm (t) is the Gaussian noise added to the m-th array element at time t, and the time difference t mp between the m-th array element receiving the p-th lightning signal and the first reference array element receiving the p-th lightning signal is (m-1)dsinθ p /c; 将雷电信号的时域阵列接收数据做傅里叶变换，则雷电信号的频域阵列接收数据为：Perform Fourier transform on the lightning signal time domain array receiving data, and then the lightning signal frequency domain array receiving data is: 其中，fg为雷电信号频谱中第g个频率点，θp为第p个雷电信号的入射角度，am(θp,fg)为第p个雷电信号的导向矢量，Ip(fg)第p个雷电信号ip(t)的傅里叶变换，Nm(fg)为高斯噪声nm(t)的傅里叶变换；Wherein, fg is the gth frequency point in the lightning signal spectrum, θp is the incident angle of the pth lightning signal, am ( θp , fg ) is the steering vector of the pth lightning signal, Ip ( fg ) is the Fourier transform of the pth lightning signal ip (t), and Nm ( fg ) is the Fourier transform of the Gaussian noise nm (t); 将雷电信号的频域阵列接收数据写成矩阵形式：Write the frequency domain array receiving data of lightning signals into matrix form: X(fg)＝A(θ,fg)I(fg)+N(fg)X(f g )=A(θ,f g )I(f g )+N(f g ) 将上述矩阵各个部分展开，分别为：Expand the above matrix into: A(θ,fg)＝[a(θ1,fg),a(θ2,fg),...,a(θP,fg)]A(θ,f g )=[a(θ 1 ,f g ),a(θ 2 ,f g ),...,a(θ P ,f g )] I(fg)＝[I1(fg),I2(fg),...,IP(fg)]T I(f g )=[I 1 (f g ), I 2 (f g ),..., I P (f g )] T N(fg)＝[N1(fg),N2(fg),...,NM(fg)]T N(f g )＝[N 1 (f g ),N 2 (f g ),...,N M (f g )] T X(fg)＝[X1(fg),X2(fg),...,XM(fg)]T。X(f g )=[X 1 (f g ),X 2 (f g ),...,X M (f g )] T . 3.根据权利要求1所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：所述步骤3、(一)中确定方法如下：3. According to the method for estimating the direction of arrival of lightning signals based on energy distribution according to claim 1, it is characterized in that: the determination method in step 3 (i) is as follows: 首先，根据雷电信号频率范围，对临界频率点前后的雷电信号频段进行平均分割，确定临界频率点前后等长子带的数量；Firstly, according to the lightning signal frequency range, the lightning signal frequency band before and after the critical frequency point is evenly divided to determine the number of equal-length sub-bands before and after the critical frequency point; 其次，将临界频率点前的各个第偶数个子带看作临界频率点前原偶数子带，根据各个原偶数子带的能量与其能量总和之比，确定临界频率点前各个新子带的长度，临界频率点前新子带数量与临界频率点前原偶数子带的数量相同，这样就确定了临界频率点前新子带数量和长度，在临界频率点后直接选择各个第偶数个子带作为临界频率点后的新子带，这样就确定了临界频率点后的新子带的数量和长度；Secondly, each even-numbered subband before the critical frequency point is regarded as the original even-numbered subband before the critical frequency point. According to the ratio of the energy of each original even-numbered subband to the total energy thereof, the length of each new subband before the critical frequency point is determined. The number of new subbands before the critical frequency point is the same as the number of original even-numbered subbands before the critical frequency point. In this way, the number and length of new subbands before the critical frequency point are determined. After the critical frequency point, each even-numbered subband is directly selected as the new subband after the critical frequency point. In this way, the number and length of new subbands after the critical frequency point are determined. 对雷电信号阵列接收数据在临界频率点前后进行等长分割，产生的原偶数子带，在临界频率点前根据各个原偶数子带的能量在其能量总和中的占比，确定临界频率点前的新子带，在临界点后直接选择各个第偶数个子带作为临界频率点后的新子带。The lightning signal array received data is divided into equal lengths before and after the critical frequency point, and the original even sub-bands generated are used to determine the new sub-bands before the critical frequency point according to the proportion of the energy of each original even sub-band in its total energy. After the critical point, each even sub-band is directly selected as the new sub-band after the critical frequency point. 4.根据权利要求3所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：雷电信号在临界频率点前，新子带的长度和数量的确定方法如下：4. The method for estimating the direction of arrival of a lightning signal based on energy distribution according to claim 3 is characterized in that: before the critical frequency point of the lightning signal, the length and number of the new sub-bands are determined as follows: (1)通过对雷电信号的频谱进行分析，将临界频率点前的频段分割为等长子带，等长子带的数量为Q，在临界频率点前选择各个第偶数个子带作为临界频率点前的原偶数子带，临界频率点前原偶数子带的数量为Q1；(1) By analyzing the spectrum of the lightning signal, the frequency band before the critical frequency point is divided into equal-length sub-bands, the number of which is Q. Each even-numbered sub-band before the critical frequency point is selected as the original even-numbered sub-band before the critical frequency point, and the number of the original even-numbered sub-band before the critical frequency point is Q 1 ; 其中，雷电信号临界频率点前原偶数子带的数量Q1的计算方法为：Among them, the calculation method of the number of original even sub-bands Q1 before the critical frequency point of the lightning signal is: n为正整数 n is a positive integer (2)雷电信号临界频率点前的频段长度为L，Q1个原偶数子带的子带长度总和L1的计算方法为：(2) The length of the frequency band before the critical frequency point of the lightning signal is L. The calculation method of the sum of the sub-band lengths L1 of the original even sub-bands Q1 is: L1＝L/Q×Q1 L 1 = L/Q × Q 1 其中，L为雷电信号临界频率点前的频段长度，Q1为雷电信号临界频率点前原偶数子带的数量，Q为雷电信号临界频率点前分割的等长子带的数量；Among them, L is the frequency band length before the critical frequency point of the lightning signal, Q1 is the number of original even sub-bands before the critical frequency point of the lightning signal, and Q is the number of equal-length sub-bands divided before the critical frequency point of the lightning signal; (3)雷电信号在临界频率点前第q个原偶数子带的能量Eq的计算方法为：(3) The calculation method of the energy Eq of the qth original even sub-band of the lightning signal before the critical frequency point is: 其中，V(j2πfr)是雷电信号的傅里叶变换，j为虚数单位，fr为雷电信号频谱图中第r个频率点，需要计算能量的原偶数子带的起始频率点为雷电信号频谱中的第r1个频率点，需要计算能量的原偶数子带的终止频率点为雷电信号频谱中的第r2个频率点；Wherein, V( j2πfr ) is the Fourier transform of the lightning signal, j is an imaginary unit, fr is the rth frequency point in the lightning signal spectrum, the starting frequency point of the original even subband for which energy needs to be calculated is the r1th frequency point in the lightning signal spectrum, and the ending frequency point of the original even subband for which energy needs to be calculated is the r2th frequency point in the lightning signal spectrum; 经过计算，雷电信号临界频率点前各个原偶数子带的能量确定为 After calculation, the energy of each original even sub-band before the critical frequency point of the lightning signal is determined as (4)将上述计算出来的雷电信号临界频率点前的Q1个原偶数子带的能量求和，得到总能量E，计算方法为：(4) Sum the energies of the Q 1 original even sub-bands before the critical frequency point of the lightning signal calculated above to obtain the total energy E. The calculation method is: 其中，Eq为雷电信号在临界频率点前的第q个原偶数子带的能量，Q1为雷电信号在临界频率点前的原偶数子带的个数；Where, Eq is the energy of the qth original even sub-band of the lightning signal before the critical frequency point, and Q1 is the number of original even sub-bands of the lightning signal before the critical frequency point; (5)雷电信号在临界频率点前的第q个新子带的长度Cq计算方法为：(5) The length Cq of the qth new sub-band of the lightning signal before the critical frequency point is calculated as follows: Cq＝Eq/E×L1 C q = E q /E × L 1 其中，Eq为雷电信号临界频率点前第q个原偶数子带的能量，E为雷电信号在临界频率点前的Q1个原偶数子带的能量总和，L1为Q1个原偶数子带的子带长度总和；Where, Eq is the energy of the qth original even sub-band before the critical frequency point of the lightning signal, E is the sum of the energies of the Q1 original even sub-bands before the critical frequency point of the lightning signal, and L1 is the sum of the sub-band lengths of the Q1 original even sub-bands; (6)雷电信号临界频率点前新子带的数量W1的计算方法如下：(6) The calculation method of the number of new sub-bands W1 before the critical frequency point of the lightning signal is as follows: W1＝Q1 W 1 =Q 1 新子带个数与原偶数子带的个数相同。The number of new subbands is the same as the number of original even subbands. 5.根据权利要求3所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：雷电信号在临界频率点后，新子带的长度和数量确定方法如下：5. The method for estimating the direction of arrival of a lightning signal based on energy distribution according to claim 3 is characterized in that: after the lightning signal reaches a critical frequency point, the length and number of the new sub-bands are determined as follows: (1)通过对雷电信号的频谱进行分析，将临界频率点后的频段分割为等长子带，等长子带的数量为G，在临界频率点后选择各个第偶数个子带作为临界频率点后的原偶数子带，临界频率点后的原偶数子带的数量为G1；(1) By analyzing the spectrum of the lightning signal, the frequency band after the critical frequency point is divided into equal-length sub-bands, the number of which is G, and each even-numbered sub-band after the critical frequency point is selected as the original even-numbered sub-band after the critical frequency point, the number of which is G 1 ; 其中，雷电信号临界频率点后原偶数子带的数量G1的计算方法为：Among them, the calculation method of the number of original even sub-bands G1 after the critical frequency point of the lightning signal is: n为正整数 n is a positive integer (2)雷电信号临界频率点后的频段长度为L’；(2) The length of the frequency band after the critical frequency point of the lightning signal is L’; (3)雷电信号在临界频率点后的第k个新子带的长度Bk计算方法如下：(3) The length Bk of the kth new sub-band of the lightning signal after the critical frequency point is calculated as follows: Bk＝L'/GB k = L'/G 其中，L’为雷电信号在临界频率点后的频段长度，G为雷电信号在临界频率点后分割的等长子带的数量；Wherein, L’ is the frequency band length of the lightning signal after the critical frequency point, and G is the number of equal-length sub-bands into which the lightning signal is divided after the critical frequency point; (4)雷电信号临界频率点后新子带的数量W2的计算方法如下：(4) The calculation method of the number of new sub-bands W2 after the critical frequency point of the lightning signal is as follows: W2＝G1 W2 ＝ G1 临界频率点后新子带个数与临界频率点后原偶数子带的个数相同。The number of new subbands after the critical frequency point is the same as the number of original even subbands after the critical frequency point. 6.根据权利要求1所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：所述步骤3、(二)中，根据上述基于能量分布的划分新子带的原则，雷电信号临界频率点为fh，在临界频率点fh前，将雷电信号阵列接收数据X(fg)划分为W1个新子带，第q个新子带的长度为Cq，将临界频率点前原偶数子带的中心点作为临界频率点前新子带的中心点；在临界频率点fh后，将雷电信号阵列接收数据X(fg)划分为W2个新子带，第k个新子带的长度为Bk，将临界频率点后原偶数子带的中心点作为临界频率点后新子带的中心点，这样就可以对雷电信号阵列接收数据X(fg)进行准确划分，将雷电阵列接收数据划分为W个新子带。6. A method for estimating the direction of arrival of a lightning signal based on energy distribution according to claim 1, characterized in that: in the step 3, (ii), according to the principle of dividing new sub-bands based on energy distribution, the critical frequency point of the lightning signal is f h , before the critical frequency point f h , the lightning signal array received data X (f g ) is divided into W 1 new sub-bands, the length of the qth new sub-band is C q , and the center point of the original even sub-band before the critical frequency point is used as the center point of the new sub-band before the critical frequency point; after the critical frequency point f h , the lightning signal array received data X (f g ) is divided into W 2 new sub-bands, the length of the kth new sub-band is B k , and the center point of the original even sub-band after the critical frequency point is used as the center point of the new sub-band after the critical frequency point, so that the lightning signal array received data X (f g ) can be accurately divided, and the lightning array received data can be divided into W new sub-bands. 7.根据权利要求1所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：所述步骤4中，对划分的新子带使用旋转信号子空间RSS方法，具体步骤为：7. The method for estimating the direction of arrival of lightning signals based on energy distribution according to claim 1 is characterized in that: in step 4, the rotated signal subspace RSS method is used for the divided new subbands, and the specific steps are: (1)采用常规波束形成法CBF对雷电信号的入射角度进行预先估计，得到预估角度η；(1) The incident angle of the lightning signal is estimated in advance using the conventional beamforming method CBF to obtain the estimated angle η; (2)将各新子带的阵列接收数据在时域内分为K段，对每一新子带都做傅里叶变换，得到新子带在频域内的阵列接收数据；(2) Divide the array received data of each new sub-band into K segments in the time domain, perform Fourier transform on each new sub-band, and obtain the array received data of the new sub-band in the frequency domain; (3)将新子带的阵列接收数据在频域内划分为N个更短子带，计算得到各个更短子带阵列接收数据的协方差矩阵为：(3) The array received data of the new subband is divided into N shorter subbands in the frequency domain, and the covariance matrix of the array received data of each shorter subband is calculated as: R(fn,η)＝E[X(fn)XH(fn)]R(f n ,η)=E[X(f n )X H (f n )] 其中，η为预先估计出的雷电信号入射角，fn为更短子带中的第n个频点，X(fn)为雷电信号在频率点fn处的阵列接收数据；Wherein, η is the pre-estimated lightning signal incident angle, fn is the nth frequency point in the shorter sub-band, and X( fn ) is the array received data of the lightning signal at the frequency point fn ; (4)确定参考频率为fc，根据预估角度η重新构造阵列流型矩阵A(fn,η)和A(fc,η)(4) Determine the reference frequency as f c , and reconstruct the array flow matrix A(f n ,η) and A(f c ,η) according to the estimated angle η (5)聚焦矩阵T(fn)用以下表达式确定：(5) The focusing matrix T(f n ) is determined by the following expression: T(fn)＝V(fn)UH(fn)T(f n )＝V(f n )U H (f n ) 其中，U(fn)和V(fn)分别为A(fn,η)AH(fc,η)的左右奇异向量，(·)H代表共轭转置；where U(f n ) and V(f n ) are the left and right singular vectors of A(f n ,η)A H (f c ,η), respectively, and (·) H represents the conjugate transpose; (6)使各个频点上的阵列接收数据通过聚焦操作“聚焦”到参考频率fc上，聚焦后协方差矩阵R(fc,η)为：(6) The array receiving data at each frequency point is “focused” to the reference frequency f c through a focusing operation. After focusing, the covariance matrix R (f c , η) is: T(fn)R(fn,η)TH(fn)≈R(fc,η)T(f n )R(f n ,η)T H (f n )≈R(f c ,η) (7)对聚焦后协方差矩阵R(fc,η)做特征值分解，获得其空间谱函数，第w个新子带的空间谱函数计算方法如下：(7) Perform eigenvalue decomposition on the focused covariance matrix R(f c ,η) to obtain its spatial spectrum function. The spatial spectrum function of the w-th new subband is The calculation method is as follows: 其中，以0.01的间隔在[0,180]之间取值，以实现角度在[0,180]之间的搜索，fc为参考频率，a(fc)为信号子空间的导向矢量，ZN为噪声子空间。in, The value is taken between [0,180] with an interval of 0.01 to realize the search of angle between [0,180]. fc is the reference frequency, a( fc ) is the steering vector of the signal subspace, and ZN is the noise subspace. 8.根据权利要求7所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：所述(1)中雷电信号入射角预估的步骤如下：8. The method for estimating the direction of arrival of a lightning signal based on energy distribution according to claim 7, characterized in that: the step of estimating the incident angle of the lightning signal in (1) is as follows: a.入射雷电信号经过均匀线性阵列，忽略噪声的接收数据为：a. The incident lightning signal passes through a uniform linear array, ignoring the noise, and the received data is: II(t)＝a(θ)i(t)II(t)＝a(θ)i(t) 其中，a(θ)为雷电信号的导向矢量，i(t)为雷电信号；Among them, a(θ) is the guiding vector of the lightning signal, and i(t) is the lightning signal; b.对接收到的信号做相位补偿并求和，得到雷电信号接收数据经过波束形成后的输出信号y(t)为：b. Perform phase compensation on the received signal and sum it up to obtain the output signal y(t) of the lightning signal received data after beamforming: y(t)＝γH·II(t)y(t)＝ γH ·II(t) 其中，γ为相位补偿向量，(·)H代表共轭转置；Where γ is the phase compensation vector, (·) H represents the conjugate transpose; c.常规波束形成器的输出功率PCBF为：c. The output power P CBF of the conventional beamformer is: PCBF＝γH·R·γ PCBF = γH ·R·γ 其中，R为雷电信号通过均匀线性阵列后接收到的数据II(t)的协方差矩阵；Where R is the covariance matrix of the data II(t) received after the lightning signal passes through the uniform linear array; d.搜索输出功率的峰值，峰值对应的角度η即为雷电信号入射角的估计值。d. Search for the peak value of the output power. The angle η corresponding to the peak value is the estimated value of the incident angle of the lightning signal. 9.根据权利要求1所述的一种基于能量分布对雷电信号进行波达方向估计的方法，其特征在于：所述步骤5中，对各个新子带的空间谱函数取平均，雷电信号空间谱函数计算方式如下：9. The method for estimating the direction of arrival of a lightning signal based on energy distribution according to claim 1, characterized in that: in step 5, the spatial spectrum function of each new sub-band is averaged, and the spatial spectrum function of the lightning signal is averaged. The calculation is as follows: 其中，W为新子带个数，为雷电信号频谱中第w个新子带的空间谱函数；Where W is the number of new subbands, is the spatial spectrum function of the wth new subband in the lightning signal spectrum; 对雷电信号空间谱函数进行谱峰搜索，搜索到的极大值对应的角度即为雷电信号波达方向估计结果。Spatial spectrum function of lightning signal Perform peak search and find The angle corresponding to the maximum value This is the estimated result of the lightning signal arrival direction.
CN202210659040.1A 2022-06-10 2022-06-10 A method for estimating the direction of arrival of lightning signals based on energy distribution Active CN115061083B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210659040.1A CN115061083B (en) 2022-06-10 2022-06-10 A method for estimating the direction of arrival of lightning signals based on energy distribution

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210659040.1A CN115061083B (en) 2022-06-10 2022-06-10 A method for estimating the direction of arrival of lightning signals based on energy distribution

Publications (2)

Publication Number Publication Date
CN115061083A CN115061083A (en) 2022-09-16
CN115061083B true CN115061083B (en) 2024-11-22

Family

ID=83200859

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210659040.1A Active CN115061083B (en) 2022-06-10 2022-06-10 A method for estimating the direction of arrival of lightning signals based on energy distribution

Country Status (1)

Country Link
CN (1) CN115061083B (en)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106950531A (en) * 2016-11-29 2017-07-14 中国人民解放军理工大学 The thunder and lightning multiple source directional method inverted based on frequency domain time

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20150130845A (en) * 2014-05-14 2015-11-24 삼성전자주식회사 Apparatus and Device for Position Measuring of Electronic Apparatuses
CN114019449B (en) * 2022-01-10 2022-04-19 南京理工大学 Direction of arrival estimation method, device, electronic device and storage medium for signal source

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106950531A (en) * 2016-11-29 2017-07-14 中国人民解放军理工大学 The thunder and lightning multiple source directional method inverted based on frequency domain time

Also Published As

Publication number Publication date
CN115061083A (en) 2022-09-16

Similar Documents

Publication Publication Date Title
Dai et al. Spatial smoothing for direction of arrival estimation of coherent signals in the presence of unknown mutual coupling
CN103091661B (en) Broadband signal arriving direction estimation method based on iteration spectral reconfiguration
CN102013911A (en) Broadband signal direction of arrival (DOA) estimation method based on threshold detection
CN110412499B (en) Wideband DOA Estimation Algorithm Based on RSS Algorithm Based on Compressed Sensing Theory
Zhang et al. An IDFT approach for coprime array direction-of-arrival estimation
CN107092007A (en) A kind of Wave arrival direction estimating method of virtual second order array extension
CN110161489B (en) Strong and weak signal direction finding method based on pseudo frame
He et al. Simultaneous interference localization and array calibration for robust adaptive beamforming with partly calibrated arrays
Dong et al. DOA estimation for coprime array with mixed noise scenario via phased fractional low-order moment
Qi et al. Time-frequency DOA estimation of chirp signals based on multi-subarray
CN111580042B (en) A deep learning direction finding method based on phase optimization
CN115061083B (en) A method for estimating the direction of arrival of lightning signals based on energy distribution
CN113805139A (en) Broadband signal sparse representation direction-of-arrival estimation method based on focusing transformation
CN113625220B (en) A new method for fast estimation of direction of arrival and diffusion angle of multipath signals
CN108594165B (en) A Method for Estimating Direction of Arrival for Narrowband Signals Based on Expectation-Maximization Algorithm
Ganguly et al. An efficient DOA estimation and jammer mitigation method by means of a single snapshot compressive sensing based sparse coprime array
CN118501803A (en) Direction of Arrival Estimation Method for Wideband Signals Based on Multi-band Joint Sparse Bayesian Learning
Yu et al. A robust minimum variance beamformer with new constraint on uncertainty of steering vector
CN117932319A (en) Rapid gridless direction-of-arrival estimation method suitable for broadband linear frequency modulation signals
Tayem et al. Propagator rooting method direction of arrival estimation based on real data
CN113655443B (en) Low-frequency-band SAR radio frequency interference suppression method for broadband digital television signals
Qu et al. A new method for weak signals' DOA estimation in the presence of strong interferences
Sarac et al. Experimental analysis of detection and localization of multiple emitters in multipath environments
Chen et al. A PCA-BP fast estimation method for broadband two-dimensional DOA of high subsonic flight targets based on the acoustic vector sensor array
CN112666558B (en) Low-complexity MUSIC direction finding method and device suitable for automobile FMCW radar

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant