Traveling wave ranging method and device based on LM nonlinear least squares optimization algorithm

Abstract

The present invention discloses a traveling wave ranging method and device based on the LM (Levenberg‑Marquardt) nonlinear least squares optimization algorithm, and relates to the field of transmission line relay protection. The method includes obtaining the voltage and current values ​​at the protection installation locations on both sides, and using phase mode decomposition and wavelet decomposition to calculate the line mode high-frequency components of the current traveling wave; comprehensively utilizing multi-band information to calculate the product results of multiple levels of line mode high-frequency components; using the LM algorithm to fit the high-frequency component product results based on parabolas and sine functions, respectively, to determine fitting equations and parameters; calculating the peak point of the product result according to the fitting equation, determining the arrival time of the traveling wave of the protection on both sides, and completing accurate traveling wave ranging. The present invention uses the LM nonlinear least squares optimization algorithm, and can accurately capture the sampling points corresponding to the peak value of the wave head through theoretical calculations, thereby improving the accuracy of traveling wave ranging.

Description

Technical Field

[0001] The invention relates to the technical field of power transmission line relay protection, and more particularly to a traveling wave ranging method and device based on LM nonlinear least squares optimization algorithm. Background Art

[0002] With the continuous commissioning of ultra-high voltage long-distance transmission lines in my country, improving the accuracy and precision of fault distance measurement is of great significance to reducing the pressure of on-site maintenance and shortening power outage time. The traditional impedance distance measurement method is easily affected by factors such as system operation mode, short-circuit transition resistance and transformer transmission error, and it is difficult to meet today's precise positioning requirements. The modern traveling wave distance measurement method is affected by factors such as the uncertainty of the traveling wave source, the variable busbar wiring method and the difficulty in capturing the wave head peak, and its reliability and accuracy still need to be improved.

[0003] In response to the above problems, scholars at home and abroad have carried out a lot of research work, mainly including: related papers ("Research on combined traveling wave ranging method for transmission lines", Wang Kuixin, Zhu Cheng, Sun Jiajia, etc., "Power System Protection and Control", Issue 40, 2012, pp. 83-86), proposed a ranging scheme using the fusion of single-end traveling wave and double-end traveling wave ranging information, but this method does not consider the problem that the initial wave head and the reflected wave head peak of the fault are difficult to accurately locate; related papers ("Single-end fault location method based on multi-source transient information fusion", Deng Feng, Xu Fan, Zeng Zhe, etc., "Transactions of China Electrotechnical Society", Issue 37, 2022, pp. 3202-3212) A traveling wave distance measurement method based on the frequency-varying characteristics of transmission line parameters is proposed. This method uses the time domain characteristics of traveling waves or the difference between high and low frequency components to perform fault distance measurement. However, considering the complex nonlinear relationship between traveling wave steepness, traveling wave attenuation characteristics and fault distance, the distance measurement result is not accurate. A related paper ("Traveling wave fault distance measurement for transmission lines based on energy ratio method", Han Zhikun, Rao Shuyong, Jiang Yushan, etc., "Power System Technology", Issue 35, 2011, pp. 217-220) proposes a distance measurement method that uses traveling wave energy to resist interference. This method uses the square sum of traveling wave sampling points in a fixed time window to calculate the energy ratio to achieve a distance measurement function with noise resistance. However, it does not consider the peak point error problem caused by the interpolation deviation of protection sampling on both sides in the actual system.

[0004] In summary, for the traveling wave ranging method of AC transmission lines, the existing research is still imperfect in capturing the peak point of the wave head. Summary of the invention

[0005] In view of this, the present invention provides a traveling wave ranging method and device based on the LM nonlinear least squares optimization algorithm, which can accurately capture the sampling points corresponding to the wave head peak through theoretical calculation, thereby improving the accuracy of traveling wave ranging.

[0006] In order to achieve the above object, the present invention adopts the following technical solution:

[0007] S1: Get the current values ​​on the M and N sides of the protection installation;

[0008] S2: according to the current value, using phase mode decomposition and wavelet decomposition to calculate the line mode high frequency components of multiple levels of the current traveling wave;

[0009] S3: Comprehensively utilize multi-band information to calculate the product results of high-frequency components of multiple-level line modes;

[0010] S4: Using the LM algorithm, the high-frequency component product result is fitted based on the parabolic function;

[0011] S5: Using the LM algorithm, the high-frequency component product result is fitted based on the sine function;

[0012] S6: Determine the final fitting equations of the M side and the N side according to the variance of the fitting results in S4 and S5;

[0013] S7: Calculate the peak sampling point numbers of the protection on both sides according to the final fitting equation;

[0014] S8: Calculate the fault distance based on the peak sampling point numbers of the protection on both sides.

[0015] In one embodiment, step S2 specifically includes the following steps:

[0016] S21. Perform phase mode decomposition on the current value to obtain the line mode current. The specific formula is as follows:

[0017] i M1 (n)=[i Ma (n)-i Mb (n)] / 3;

[0018] Where: n is the current sampling point number; i Ma (n), i Mb (n) and i M1 (n) are respectively the A-phase current, B-phase current and line mode current of the M-side protection at the current sampling point number n;

[0019] S22, calculating the line mode current high frequency component and the line mode current low frequency component of the first layer of wavelet decomposition according to the line mode current, the specific calculation formula is as follows:

[0020]

[0021]

[0022] In the formula, i MH1 (n) and i ML1(n) respectively represent the line mode current high frequency component and line mode current low frequency component of the first layer of wavelet decomposition at the current sampling point number n; N H and N L are the lengths of the high-pass filter coefficient and the low-pass filter coefficient respectively; k represents the current filter coefficient point number; P H represents the high-pass filter coefficient of wavelet decomposition; P L Represents the low-pass filter coefficients of wavelet decomposition;

[0023] S23, using the wavelet decomposition of the line mode current low frequency component of the first layer to calculate the line mode current high frequency component and the line mode current low frequency component of the second layer of wavelet decomposition, the specific calculation formula is as follows:

[0024]

[0025]

[0026] In the formula, i MH2 (n) and i ML2 (n) respectively represent the line mode current high frequency component and line mode current low frequency component of the second layer of wavelet decomposition at the current sampling point number n;

[0027] S24, using the low-frequency component of the line mode current of the second layer of wavelet decomposition to calculate the high-frequency component of the line mode current of the third layer of wavelet decomposition, the specific calculation formula is as follows:

[0028]

[0029] Where: i MH2 (n) represents the high-frequency component of the line mode current in the third layer of wavelet decomposition at the current sampling point number n.

[0030] In one embodiment, in step S2, the filter coefficient P H and P L They are set to [-0.23040.7148-0.6309-0.02800.18700.0308-0.0329-0.0106] and [-0.01060.03290.0308-0.1870-0.02800.63090.71480.2304] respectively, with a coefficient length of N H and N L They are 8 respectively.

[0031] In one embodiment, in step S3, the multi-band information is comprehensively utilized to calculate the product result of the high-frequency components of the line modes of multiple levels, which specifically includes the following calculation formula:

[0032]

[0033] Where: i MH(n) and i MHj (n) respectively represent the comprehensive component of the line mode current at the current sampling point number n and the high-frequency component of the jth layer of wavelet decomposition; П is the standard symbol for polynomial multiplication; N lv is the number of frequency bands for calculating the integrated components.

[0034] In one embodiment, in step S3, the number of frequency bands N of the integrated component is calculated. lv Set to 3.

[0035] In one embodiment, in step S4, the LM algorithm is used to fit the high-frequency component product result based on the parabolic function, specifically including using the following fitting calculation formula:

[0036]

[0037] Where: Δa 1j , Δa 2j and Δa 3j are the regression step lengths of the three parameters to be fitted in the jth iteration; i MH (n) represents the comprehensive component of the line mode current at the current sampling point number n; y(n) is the parabola equation containing the coefficients to be fitted; N T is the number of data window points to be adjusted; a1, a2 and a3 are the three coefficients to be fitted respectively; J y is the Jacobian matrix of the parabolic function fitting; the superscripts “T” and “-1” are the standard symbols for matrix transposition and inversion, respectively; λ is the damping factor; I is the identity matrix; a 1set 、a 2set and a 3set are the initial setting values ​​of the three parameters to be fitted; N D Set an upper limit for the number of iterations.

[0038] Among them, the Jacobian matrix based on parabolic function fitting can be calculated by the following formula:

[0039]

[0040] Where: n is the current sampling point number; y(n) is the parabola equation containing the coefficients to be fitted; a1, a2 and a3 are the three coefficients to be fitted; y1(n), y2(n) and y3(n) are the partial derivative equations of the function with respect to the three coefficients to be fitted; Standard notation for partial derivatives; i MH is the comprehensive component of line mode current; N T is the number of data window points to be set; J y Jacobian matrix for the parabolic function fit.

[0041] In one embodiment, in step S5, the LM algorithm is used to fit the high-frequency component product result based on the sine function, specifically including using the following fitting calculation formula:

[0042]

[0043] Where: Δb 1j , Δb 2j , Δb 3j and Δb 4j are the regression step lengths of the four parameters to be fitted in the jth iteration; i MH (n) represents the comprehensive component of the line mode current at the current sampling point number n; z(n) is the sinusoidal function equation containing the coefficients to be fitted; N T is the number of data window points to be adjusted; b1, b2, b3 and b4 are the four coefficients to be fitted respectively; J z is the Jacobian matrix of the parabolic function fitting; the superscripts “T” and “-1” of the matrix are the standard symbols for matrix transposition and inversion, respectively; λ is the damping factor; I is the unit matrix; b 1set 、b 2set 、b 3set and b 4set are the initial setting values ​​of the four parameters to be fitted; N D Set an upper limit for the number of iterations.

[0044] Among them, the Jacobian matrix based on sine function fitting can be calculated by the following formula:

[0045]

[0046] Where: n is the current sampling point number; z(n) is the parabola equation containing the coefficients to be fitted; b1, b2, b3 and b4 are the four coefficients to be fitted; z1(n), z2(n), z3(n) and z4(n) are the partial derivative equations of the function with respect to the four coefficients to be fitted; Standard notation for partial derivatives; f samp is the sampling frequency; i MH is the comprehensive component of line mode current; N T is the number of data window points to be set; J y Jacobian matrix for the parabolic function fit.

[0047] In one embodiment, in steps S4 and S5, the number of data window points N T is set to 20, the damping factor λ is set to 3, and the upper limit of the number of iterations is N D Set to 25.

[0048] In one embodiment, in step S6, the final fitting equation of the M side is determined according to the variance of the fitting results in S4 and S5, which specifically includes the following steps:

[0049] If the variance of the fitting result satisfies If n is 0, the sine equation is selected as the final fitting equation, otherwise the parabola equation is selected as the final fitting equation, where y(n) is the parabola equation containing the coefficients to be fitted; z(n) is the sine function equation containing the coefficients to be fitted; i MH (n) represents the comprehensive component of the line mode current at the current sampling point number n; N T is the number of data window points to be adjusted.

[0050] In one embodiment, in step S7, the peak sampling point numbers of the protection on both sides are calculated according to the final fitting equation, which specifically includes the following steps:

[0051] If the parabola equation is determined to be the final fitting equation, the calculation formula for the peak sampling point number of the wave head is:

[0052]

[0053] Where: n max is the sampling point number of the wave head peak; a1 and a2 are the fitting coefficients of the two-term parabola equation respectively.

[0054] If the sine equation is determined to be the final fitting equation, the calculation formula for the peak sampling point number of the wave head is:

[0055]

[0056] Where: n max is the sampling point number of the wave head peak; b2 and b3 are the fitting coefficients of the two-term sine equation respectively; f samp is the sampling frequency.

[0057] In one embodiment, in step S8, the fault distance is calculated using the peak sampling point numbers of the wave heads on both sides, specifically including using the following calculation formula:

[0058]

[0059] Where: L F and L MN N are the traveling wave ranging distance between the fault point and the M-side protection and the total length of the line respectively; Mmax and N Nmax are the wave head peak sampling point numbers of the M-side and N-side protection calculated based on the LM fitting method; v is the wave velocity, and its value is Where L and C represent the positive sequence inductance and positive sequence capacitance of the line per unit length respectively; f samp is the sampling frequency.

[0060] It can be seen from the above technical solution that, compared with the prior art, the present invention discloses a traveling wave ranging method and device based on the LM nonlinear least squares optimization algorithm, which has the following beneficial effects:

[0061] The technical solution of the present invention is different from the prior art. In order to address the problem that the wave head peak point is difficult to accurately extract, the LM nonlinear least squares optimization algorithm is adopted to accurately capture the sampling points corresponding to the wave head peak through theoretical calculation, thereby improving the accuracy of traveling wave ranging. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0063] Figure 1 It is a schematic diagram of the method flow provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0064] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0065] Embodiment 1:

[0066] like Figure 1 As shown, the present invention discloses a traveling wave ranging method based on LM nonlinear least squares optimization algorithm, which includes the following steps:

[0067] S1: Get the current values ​​on the M and N sides of the protection installation;

[0068] S2: According to the above current value, the line mode high frequency components of multiple levels of the current traveling wave are calculated by using phase mode decomposition and wavelet decomposition; for a three-phase transmission line, the line mode component can be the AB phase current or the AC phase current. This embodiment takes the AB phase current as an example to introduce in detail the calculation process of the line mode high frequency components of multiple levels of the current traveling wave on the M side:

[0069] First, the current value is decomposed into phase mode to obtain the line mode current. The specific formula is as follows:

[0070] i M1 (n)=[i Ma (n)-iMb (n)] / 3, n is the current sampling point number; i M1 is the line mode current of the M side protection, which is 1 / 3 times of the AB phase current or AC phase current. Here, the AB phase current is taken as an example; i Ma and i Mb They are the A-phase current and B-phase current of the M-side protection respectively.

[0071] Then, wavelet decomposition is performed according to the line mode current to obtain multiple levels of line mode high frequency components. The specific formula is as follows:

[0072]

[0073] Where: i MH1 、i MH2 and i MH3 are the high-frequency components of the line-mode current in the first, second and third layers of wavelet decomposition, and their corresponding frequency bands are [f samp / 4,f samp / 2]、[f samp / 8,f samp / 4] and [f samp / 16,f samp / 8], where f samp is the sampling frequency; i ML1 and i ML2 are the low-frequency components of the line mode current in the first and second layers of wavelet decomposition, and their corresponding frequency bands are [0,f samp / 4] and [0,f samp / 8]; P H and P L are the high-pass and low-pass filter coefficients of wavelet decomposition; N H and N L are the lengths of high-pass and low-pass filter coefficients respectively; k is the current filter coefficient point number.

[0074] In this embodiment, the filter coefficient P H and P L Can be set to [-0.23040.7148-0.6309-0.02800.18700.0308-0.0329-0.0106] and [-0.01060.03290.0308-0.1870

[0075] -0.02800.63090.71480.2304], coefficient length N H and N L They are 8 respectively.

[0076] The above wavelet decomposition process can be expressed as:

[0077] According to the line mode current iM1 Calculate the high-frequency component i of the line mode current in the first layer of wavelet decomposition MH1 and the line mode current low frequency component i ML1 ;

[0078] Using the low-frequency component i of the line mode current in the first frequency band ML1 Calculate the high-frequency component i of the line mode current in the second layer of wavelet decomposition MH2 and the line mode current low frequency component i ML2 ;

[0079] The low-frequency component i of the line mode current in the second layer is decomposed by wavelet ML2 Calculate the high-frequency component i of the line mode current in the third layer of wavelet decomposition MH3 .

[0080] It should be noted that the calculation process of the N-side protection is the same as that of the M-side, and the specific calculation formula of the N-side protection is not listed here.

[0081] S3: Comprehensively utilize multi-band information to calculate the product results of high-frequency components of multiple-level line modes;

[0082] The calculation formula for comprehensive use of multi-band information is:

[0083]

[0084] Where: i MH and i MHj are the comprehensive component of line mode current and the high frequency component of the jth layer of wavelet decomposition respectively; П is the standard symbol for polynomial product; N lv is the number of frequency bands for calculating the integrated components.

[0085] Calculate the number of frequency bands N of the integrated component lv Set to 3.

[0086] S4: Using the LM algorithm, the high-frequency component product result is fitted based on the parabolic function;

[0087] The fitting calculation formula based on the parabolic function is:

[0088]

[0089] Where: Δa 1j , Δa 2j and Δa 3j are the regression step lengths of the three parameters to be fitted in the jth iteration; J y is the Jacobian matrix of the parabolic function fitting; the superscripts “T” and “-1” of the matrix are the standard symbols for matrix transposition and inversion, respectively; λ is the damping factor; I is the unit matrix; a 1set 、a 2set and a 3setare the initial setting values ​​of the three parameters to be fitted; N D Set an upper limit for the number of iterations.

[0090] Among them, the Jacobian matrix J based on parabolic function fitting y The calculation formula is:

[0091]

[0092] Where: n is the current sampling point number; y(n) is the parabola equation containing the coefficients to be fitted; a1, a2 and a3 are the three coefficients to be fitted; y1(n), y2(n) and y3(n) are the partial derivative equations of the function with respect to the three coefficients to be fitted; Standard notation for partial derivatives; i MH is the comprehensive component of line mode current; N T is the number of data window points to be set; J y Jacobian matrix for the parabolic function fit.

[0093] In the embodiment of the present invention, the number of data window points N T is set to 20, the damping factor λ is set to 3, and the upper limit of the number of iterations is N D Set to 25.

[0094] S5: Using the LM algorithm, the high-frequency component product result is fitted based on the sine function;

[0095] The fitting calculation formula based on the parabolic function is:

[0096]

[0097] Where: Δb 1j , Δb 2j , Δb 3j and Δb 4j are the regression step lengths of the four parameters to be fitted in the jth iteration; J y The Jacobian matrix of the sine function fit. The superscripts “T” and “-1” of the matrix are the standard symbols for matrix transposition and inversion, respectively; λ is the damping factor; I is the identity matrix; b 1set 、b 2sett 、b 3set and b 4set are the initial setting values ​​of the four parameters to be fitted; N D Set an upper limit for the number of iterations.

[0098] In the above calculation, the Jacobian matrix calculation formula based on sine function fitting is:

[0099]

[0100] Where: n is the current sampling point number; z(n) is the sine function equation containing the coefficients to be fitted; b1, b2, b3 and b4 are the four coefficients to be fitted; z1(n), z2(n), z3(n) and z4(n) are the partial derivative equations of the function with respect to the four coefficients to be fitted; Standard notation for partial derivatives; f samp is the sampling frequency; i MH is the comprehensive component of line mode current; N T is the number of data window points to be set; J y Jacobian matrix for the sine function fit.

[0101] In this embodiment, the number of data window points N T is set to 20, the damping factor λ is set to 3, and the upper limit of the number of iterations is N D Set to 25.

[0102] S6: Compare the final fitting variances in S4 and S5 to determine the fitting equation;

[0103] The criterion for determining the fitting equation is:

[0104] If the criteria are met The sine equation is judged as the final fitting equation, otherwise the parabola equation is judged as the final fitting equation.

[0105] It is worth noting that in the embodiments S3-S6 of the present invention, only the specific formula of the final fitting equation on the M side is described in detail, and the steps for obtaining the final fitting equation on the N side of the protection installation are the same as those on the M side, and the specific formula is no longer listed in detail in the present invention.

[0106] S7: Calculate the peak sampling point numbers of the protection on both sides according to the fitting equation;

[0107] If the parabola equation is determined to be the final fitting equation, the calculation formula for the peak sampling point number of the wave head is:

[0108]

[0109] Where: N max is the sampling point number of the wave head peak; a1 and a2 are the fitting coefficients of the two-term parabola equation respectively.

[0110] If the sine equation is determined to be the final fitting equation, the calculation formula for the peak sampling point number of the wave head is:

[0111]

[0112] Where: N max is the sampling point number of the wave head peak; b2 and b3 are the fitting coefficients of the two-term sine equation respectively; f samp is the sampling frequency.

[0113] S8: Calculate the fault distance based on the peak sampling point numbers of the protection on both sides.

[0114] The formula for calculating the fault distance using the peak sampling point numbers on both sides of the wave head is:

[0115]

[0116] Where: L F and L MN N are the traveling wave ranging distance between the fault point and the M-side protection and the total length of the line respectively; Mmax and N Nmax are the wave head peak sampling point numbers of the M-side and N-side protection calculated based on the LM fitting method; v is the wave velocity, and its value is Where L and C represent the positive sequence inductance and positive sequence capacitance of the line per unit length respectively; f samp is the sampling frequency.

[0117] Embodiment 2:

[0118] This embodiment discloses a device, which may be a computer device or other terminal. A computer program for implementing any method in Embodiment 1 is provided in the device.

[0119] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0120] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A traveling wave ranging method based on LM nonlinear least squares optimization algorithm, characterized in that: The following steps are included: S1: Get the current values ​​on the M and N sides of the protection installation; S2: according to the current value, using phase mode decomposition and wavelet decomposition to calculate the line mode high frequency components of multiple levels of the current traveling wave; S3: Comprehensively utilize multi-band information to calculate the product results of high-frequency components of multiple-level line modes; S4: Using the LM algorithm, the high-frequency component product result is fitted based on the parabolic function; S5: Using the LM algorithm, the high-frequency component product result is fitted based on the sine function; S6: Determine the final fitting equations of the M side and the N side according to the variance of the fitting results in S4 and S5; S7: Calculate the peak sampling point numbers of the protection on both sides according to the final fitting equation; S8: Calculate the fault distance based on the peak sampling point numbers of the protection on both sides.

2. The traveling wave ranging method based on the LM nonlinear least squares optimization algorithm according to claim 1 is characterized in that: The step S2 specifically includes the following steps: S21. Perform phase mode decomposition on the current value to obtain the line mode current. The specific formula is as follows: i M1 (n)=[i Ma (n)-i Mb (n)] / 3; Where: n is the current sampling point number; i Ma (n), i Mb (n) and i M1 (n) are respectively the A-phase current, B-phase current and line mode current of the M-side protection at the current sampling point number n; S22, calculating the line mode current high frequency component and the line mode current low frequency component of the first layer of wavelet decomposition according to the line mode current, the specific calculation formula is as follows: In the formula, i MH1 (n) and i ML1 (n) respectively represent the line mode current high frequency component and line mode current low frequency component of the first layer of wavelet decomposition at the current sampling point number n; N H and N L are the lengths of the high-pass filter coefficient and the low-pass filter coefficient respectively; k represents the current filter coefficient point number; P H represents the high-pass filter coefficient of wavelet decomposition; P L Represents the low-pass filter coefficients of wavelet decomposition; S23, using the wavelet decomposition of the line mode current low frequency component of the first layer to calculate the line mode current high frequency component and the line mode current low frequency component of the second layer of wavelet decomposition, the specific calculation formula is as follows: In the formula, i MH2 (n) and i ML2 (n) respectively represent the line mode current high frequency component and line mode current low frequency component of the second layer of wavelet decomposition at the current sampling point number n; S24, using the low-frequency component of the line mode current of the second layer of wavelet decomposition to calculate the high-frequency component of the line mode current of the third layer of wavelet decomposition, the specific calculation formula is as follows: Where: i MH3 (n) represents the high-frequency component of the line mode current in the third layer of wavelet decomposition at the current sampling point number n.

3. The traveling wave ranging method based on the LM nonlinear least squares optimization algorithm according to claim 2 is characterized in that: In step S3, the multi-band information is comprehensively utilized to calculate the product result of the high-frequency components of the line modes of multiple levels, which specifically includes the following calculation formula: Where: i MH (n) and i MHj (n) respectively represent the comprehensive component of the line mode current at the current sampling point number n and the high-frequency component of the jth layer of wavelet decomposition; П is the standard symbol for polynomial multiplication; N lv is the number of frequency bands for calculating the integrated components.

4. The traveling wave ranging method based on the LM nonlinear least squares optimization algorithm according to claim 3 is characterized in that: In step S3, the number of frequency bands N of the comprehensive component is calculated. lv Set to 3.

5. The traveling wave ranging method based on LM nonlinear least squares optimization algorithm according to claim 1 is characterized in that: In step S4, the LM algorithm is used to fit the high-frequency component product result based on the parabolic function, specifically including using the following fitting calculation formula: Where: Δa 1j , Δa 2j and Δa 3j are the regression step lengths of the three parameters to be fitted in the jth iteration; i MH (n) represents the comprehensive component of the line mode current at the current sampling point number n; y(n) is the parabola equation containing the coefficients to be fitted; N T is the number of data window points to be adjusted; a1, a2 and a3 are the three coefficients to be fitted respectively; J y is the Jacobian matrix of the parabolic function fitting; the superscripts "T" and "-1" are the standard symbols for matrix transposition and inversion, respectively; λ is the damping factor; I is the identity matrix; a 1set 、a 2set and a 3set are the initial setting values ​​of the three parameters to be fitted; N D Set an upper limit for the number of iterations.

6. The traveling wave ranging method based on LM nonlinear least squares optimization algorithm according to claim 1, characterized in that: In step S5, the LM algorithm is used to fit the high-frequency component product result based on the sine function, specifically including using the following fitting calculation formula: Where: Δb 1j , Δb 2j , Δb 3j and Δb 4j are the regression step lengths of the four parameters to be fitted in the jth iteration; i MH (n) represents the comprehensive component of the line mode current at the current sampling point number n; z(n) is the sinusoidal function equation containing the coefficients to be fitted; N T is the number of data window points to be adjusted; b1, b2, b3 and b4 are the four coefficients to be fitted respectively; J z is the Jacobian matrix of the parabolic function fitting; the superscripts "T" and "-1" of the matrix are the standard symbols for matrix transposition and inversion, respectively; λ is the damping factor; I is the identity matrix; b 1set 、b 2set 、b 3set and b 4set are the initial setting values ​​of the four parameters to be fitted; N D Set an upper limit for the number of iterations.

7. The traveling wave ranging method based on LM nonlinear least squares optimization algorithm according to claim 1 is characterized in that: In step S6, the final fitting equation of the M side is determined according to the variance of the fitting results in S4 and S5, which specifically includes the following steps: If the variance of the fitting result satisfies If n is 0, the sine equation is selected as the final fitting equation, otherwise the parabola equation is selected as the final fitting equation, where y(n) is the parabola equation containing the coefficients to be fitted; z(n) is the sine function equation containing the coefficients to be fitted; i MH (n) represents the comprehensive component of the line mode current at the current sampling point number n; N T is the number of data window points to be adjusted.

8. The traveling wave ranging method based on LM nonlinear least squares optimization algorithm according to claim 1, characterized in that: In step S7, according to the final fitting equation, the peak sampling point numbers of the protection on both sides are calculated, which specifically includes: If the parabola equation is selected as the final fitting equation, the calculation formula of the peak sampling point number of the wave head is: Where: N max is the sampling point number of the wave head peak; a1 and a2 are the fitting coefficients of the two-term parabola equation respectively; If the sine equation is selected as the final fitting equation, the calculation formula of the wave head peak sampling point number is: Where: N max is the sampling point number of the wave head peak; b2 and b3 are the fitting coefficients of the two-term sine equation respectively; f samp is the sampling frequency.

9. The traveling wave ranging method based on LM nonlinear least squares optimization algorithm according to claim 1, characterized in that: In step S8, the fault distance is calculated using the peak sampling point numbers of the wave heads on both sides, specifically including using the following calculation formula: Where: L F and L MN N are the traveling wave ranging distance between the fault point and the M-side protection and the total length of the line respectively; Mmax and N Nmax are the wave head peak sampling point numbers of the M-side and N-side protection calculated based on the LM fitting method; v is the wave velocity; f is the wave velocity; samp is the sampling frequency.

10. A traveling wave ranging device based on LM nonlinear least squares optimization algorithm, characterized in that: The device includes a computer system, and the computer system includes a computer program capable of implementing the traveling wave ranging method based on the LM nonlinear least squares optimization algorithm as described in any one of claims 1-9.

Images