A method and system for fault protection of a multi-terminal hybrid DC power transmission system

Abstract

The application relates to a fault protection method and system of a multi-terminal hybrid DC power transmission system, which comprises the following steps: based on the Bergeron numerical calculation method, simulating the current anti-travel wave derivative of a fault protection device to obtain a reference waveform; when the multi-terminal hybrid DC power transmission system is faulty, calculating the measured current anti-travel wave derivative; adding disturbance noise to the current anti-travel wave derivative to obtain a measured waveform; respectively normalizing the reference waveform and the measured waveform to obtain a normalized reference waveform and a normalized measured waveform; calculating a sliding correlation coefficient sequence between the normalized reference waveform and the normalized measured waveform; judging a fault occurrence area according to the sliding correlation coefficient sequence; and judging whether to execute fault protection according to the fault occurrence area. The application uses the Bergeron numerical calculation method to analyze the characteristics of the first wave head of the travel wave felt by the protection device under each fault scene, so as to construct a fault protection scheme which can meet the requirements of the multi-terminal hybrid DC power transmission system.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of direct current transmission, in particular to a fault protection method and system for a multi-terminal hybrid direct current transmission system. BACKGROUND

[0002] In a multi-terminal LCC-MMC-HVDC system, i.e. a multi-terminal hybrid direct current transmission system, due to the lack of adaptable DCCB (direct current circuit breaker), when a fault occurs in the direct current line, regardless of which section of the line the fault occurs in, the fault handling mode is to switch the system control mode to realize fault ride-through. The specific method is: through emergency phase shifting of the rectifier side LCC, the inverter side MMC converter is blocked or actively controlled to realize the release and ionization of the line energy, and if it is a transient fault, the line can be restarted to restore. Therefore, the coverage range of the protection device of the system is different from that of the point-to-point two-terminal system relying on the circuit breaker to remove the fault, and the principle and setting method are different from those of the existing direct current transmission system. However, the existing protection principle suitable for the multi-terminal hybrid direct current transmission system has the following defects:

[0003] Firstly, the mechanism of the existing protection principle is simple and convenient to use, such as traveling wave protection and differential under-voltage protection, which mainly relies on simulation test in setting, lacks a solid theoretical basis, which may introduce uncertainty in actual application, limits the further improvement of its sensitivity, and the anti-transition resistance capability in engineering is only dozens of ohms, which cannot meet the requirement of 600 ohm anti-transition resistance capability of the 800 kV system; secondly, although the multi-terminal hybrid direct current transmission system does not require millisecond-level fault isolation, fast and communication-independent fault detection is still crucial, and if the rectifier side can detect the fault at a super-high speed within 1 ms, the overvoltage risk of the inverter side can be greatly reduced; thirdly, the hybrid direct current transmission system has a tree-like topology structure, and the traditional direction criterion has the risk of failure, and unlike the multi-terminal flexible direct current grid, the branch line of the multi-terminal hybrid direct current transmission system has only one side of the current limiting reactor, and the limited deployment of the current limiting reactor affects the adaptability of the traditional boundary protection and the direction criterion based on the boundary.

[0004] Therefore, in the multi-terminal hybrid direct current transmission system, the existing protection method cannot meet the demand of fault partition protection. SUMMARY

[0005] In order to solve the technical problems in the prior art that the fault partition protection of the multi-terminal hybrid direct current transmission system cannot meet the demand, the present application provides a fault protection method and system for a multi-terminal hybrid direct current transmission system.

[0006] The technical solution of the present application to solve the above technical problems is as follows:

[0007] A fault protection method for a multi-terminal hybrid direct current transmission system, comprising the following steps:

[0008] Based on the Berelon numerical calculation method, the operating parameters of the multi-terminal hybrid HVDC transmission system when no fault occurs and the metallic fault parameters of the multi-terminal hybrid HVDC transmission system outside the forward zone are used to simulate the protection device p within a preset time after the forward zone fault occurs. m The reverse wave derivative of the current at the position is used to obtain a reference waveform; wherein the preset time is specifically 0.2s;

[0009] When the multi-terminal hybrid HVDC transmission system fails, based on the Berelon numerical calculation method, the protection device p is activated within the preset time after the fault start criterion is activated. m The line mode voltage and line mode current at the position are used to calculate the measured current reverse wave derivative;

[0010] Add disturbance noise to the measured current reverse wave derivative to obtain the measured waveform;

[0011] Normalizing the reference waveform and the measured waveform respectively to obtain a normalized reference waveform and a normalized measured waveform respectively;

[0012] Calculating a sliding correlation coefficient sequence between a normalized reference waveform and a normalized measured waveform;

[0013] Determine the fault occurrence area based on the sliding correlation coefficient sequence;

[0014] Determine whether to execute fault protection based on the fault occurrence area.

[0015] The beneficial effects of the present invention are: first, the Berelon method is used to establish a numerical equivalent model for equipment (including inductors, DC filters, converter stations, etc.), so that the protection relay can mathematically calculate the theoretical electrical quantity changes (such as current and voltage) at different fault locations, thereby improving the accuracy and theoretical basis of fault detection; second, the protection principle is realized based on the calibration of the initial current reverse wave derivative curve, the core of which is to capture the rate of change of the initial reverse current wave at each relay position. By comparing with the actual measured waveform, the local sliding correlation coefficient index is used for quantification, achieving fast and accurate fault judgment; third, the proposed protection scheme is applicable to both LCC side protection and downstream line protection, solving the problem that the traditional protection principle applicable to multi-terminal hybrid DC transmission systems requires separate adjustment of LCC side protection and downstream line protection, and has high adaptability and fast response capabilities; fourth, a non-directional but extremely sensitive starting criterion is introduced to address the problem of possible failure or insufficient sensitivity of directional elements. At the same time, a weak noise identical to the template waveform is artificially added to the derivative curve of the current reverse wave, and after normalization, a virtual out-of-zone fault scenario can be simulated, which not only improves sensitivity but also effectively avoids reverse faults, thereby improving overall safety.

[0016] Further, the calculation formula of adding disturbance noise to the measured current anti-travel wave derivative is as follows:

[0017]

[0018] wherein, represents the measured waveform, I' mh represents the current anti-travel wave derivative, I' mH represents the reference waveform, and a represents a disturbance noise coefficient, and a takes a value in the range of 0-5.

[0019] When a reverse fault occurs, the protection device p m After the starting criterion is started, the current anti-travel wave is almost always 0 within 0.2 ms, that is, the current anti-travel wave derivative also tends to 0, and the normalization behavior in the protection flow will amplify the small noise detected by the protection during the reverse fault, causing the protection to malfunction.

[0020] Without using a directional element to avoid a reverse fault, in order to make the above main protection criterion reliably inactionable in the case of a reverse fault, before the starting criterion is started to enter the main protection flow, fixed noise needs to be added to the measured current anti-travel wave derivative I' mh at the protection device. In view of the feature that the protection detects that the current anti-travel wave tends to 0 during a reverse fault, a times the reference waveform I' mh is added to the measured current anti-travel wave derivative I' mH . The calculation formula of the above measured waveform I' is obtained by this method.

[0021] Further, the reference waveform and the measured waveform are normalized respectively, and normalized reference waveform i' mtem and normalized measured waveform are correspondingly obtained, including the following steps:

[0022] The reference waveform and the measured waveform are respectively data-extended;

[0023] The data-extended reference waveform and the data-extended measured waveform are respectively normalized, and the normalized reference waveform and the normalized measured waveform are correspondingly obtained;

[0024] The normalization formula is as follows:

[0025]

[0026] wherein, i' mtem represents the normalized reference waveform, i' mact represents the normalized measured waveform, I' mtem represents the data-extended reference waveform, and I' mact represents the data-extended measured waveform.

[0027] Further, data augmentation is performed on the reference waveform and the measured waveform respectively, including the following steps:

[0028] The linear interpolation method is used to interpolate the reference waveform and the measured waveform respectively; wherein the interpolation formula of the linear interpolation method is as follows:

[0029]

[0030] Wherein, x N represents the abscissa value of the Nth point inserted by linear interpolation, y N represents the ordinate value of the Nth point inserted by linear interpolation, x a represents the abscissa value of the starting point of the linear interpolation interval, y a represents the ordinate value of the starting point of the linear interpolation interval, x b represents the abscissa value of the terminal point of the linear interpolation interval, y b represents the ordinate value of the terminal point of the linear interpolation interval, g represents the number of points of linear interpolation between two adjacent data points; the abscissa value is a time sequence, and the ordinate value is a waveform sequence; the data in the first preset time window of the reference waveform after interpolation is set to zero to obtain the data augmented reference waveform; the data in the first preset time window of the reference waveform after interpolation is set to zero to obtain the data augmented reference waveform; wherein the first preset time window is a time window before the starting time, and the length of the first preset time window is less than the preset time. m The length of the first preset time window is less than the preset time.

[0031] Let the length of the time window for calculating the sliding correlation coefficient be t w , t w <T, T represents the preset time, and the sampling frequency of the protection device p m is f s In order to facilitate sliding correlation analysis, data in a time window before the starting time needs to be reserved, and since the fault traveling wave has not arrived in this period, the current anti-traveling wave data in this period is set to zero. Let the reference current anti-traveling wave derivative at the starting time after linear interpolation be I' m , and the sampling value of the measured current anti-traveling wave derivative be I' mH(0) . mh(0) The data in a time window before the starting time is set as follows:

[0032]

[0033]

[0034] Wherein, the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f s t w the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f mH(-1) the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f s t w the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f mh(-1) the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f

[0035] the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f mtem the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f mact as shown in the following formula:

[0036]

[0037]

[0038] wherein, I' mH(1) the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f mH((g+1)fsT) the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f s t w the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f mh the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f m the sampling value of the reference current anti-travel wave derivative at the data point (g+1)f

[0039] Further, the sliding correlation coefficient sequence between the normalized reference waveform and the normalized measured waveform is calculated, including the following steps:

[0040] the time when the protection device p m starts to protect the fault is t 0m ; the normalized reference waveform i' mref and the normalized measured waveform i' mact can be expressed as shown in the following formula:

[0041]

[0042]

[0043] in, for Normalized data, i' mH ( -1) For I' mH(-1) Normalized data, i' mH(0) For I' mH(0) Normalized data, i' mH(1) For I' mH(1) Normalized data, i' mH((g+1)fsT) for Normalized data; for Normalized data, i' mh(-1) For I' mh(-1) Normalized data, i' mh(0) For I' mh(0) Normalized data, i' mh(1) For I' mh(1) After normalization, for Normalized data.

[0044] S101, initialize the sliding coefficient i to i=0,

[0045] S102, with t w is the calculation time window length of the sliding correlation coefficient, respectively from the normalized reference waveform i' mref and normalized measured waveform i' mact The comparison fragment is intercepted in the , and the corresponding reference fragment H is obtained mi and the measured fragment h mi ; wherein the reference fragment H mi The number of data and the measured fragment h mi The number of data in is equal;

[0046] S103, respectively, the reference fragment H mi and the measured fragment h mi The data in is arranged in descending order, corresponding to the reference descending sequence A mi and the measured descending sequence B mi ;

[0047] S104, obtaining the reference descending sequence A mi With the measured descending sequence B mi The difference between the data of the same order in the , get the rank difference sequence d mi ;

[0048] S105, according to the rank difference sequence d mi Calculate the local correlation coefficient ρ from the data in m(i);

[0049] Among them, the local correlation coefficient ρ m(i) The calculation formula is as follows:

[0050]

[0051] Specifically, s Indicates the protection device p m The frequency of collecting line mode voltage and line mode current, d represents the rank difference sequence d mi The data in , q represents the rank difference sequence d mi The order value of the data in d q represents the rank difference sequence d mi The data with order q, Indicates d q the square of

[0052] S106, set the sliding coefficient i=i+1, and execute steps S102 to S105 cyclically until the sliding coefficient i satisfies i=(g+1)f s When T, execute step S107; wherein T represents the preset time;

[0053] S107, the local correlation coefficient ρ obtained each time m(i) The data are integrated into a data sequence to obtain the sliding correlation coefficient sequence.

[0054] Reference fragment H mi t 0m -t w Sequence value corresponding to the moment is the starting point, t 0m The sequence value i' corresponding to the moment mH(0+i) and i' mh(0+i) is the end point; the measured fragment h mi t 0m -t w Sequence value corresponding to the moment is the starting point, t 0m The sequence value i' corresponding to the moment mh(0+i) As the end point, cut out the reference segment H mi With the measured fragment h mi As shown below:

[0055]

[0056]

[0057] For reference fragment H mi With the measured fragment h mi The internal data are arranged in descending order to obtain the reference descending sequence A miand the measured descending sequence B mi As shown below:

[0058]

[0059]

[0060] Reference fragment H mi The zth element in the reference descending sequence A mi The position in is denoted as r z , which is called the rank of the z-th element point, so that the reference segment H can be obtained mi The rank sequence r corresponding to all elements in mi Similarly, record the measured segment h mi The zth element point in the measured descending sequence B mi The position in is denoted as s z , and get the measured fragment h mi The rank sequence s corresponding to all elements in mi . The sequence r mi With the sequence s mi Subtract each element in the rank difference sequence d mi As shown below:

[0061]

[0062] The length of the sliding time window t for maintaining the captured waveform w The sliding window slides back one data point each time, that is, i=i+1, and the local correlation coefficient ρ within each sliding window is calculated according to the above steps. m(i) , until the time window continuously slides through all the data within the preset time T=0.2ms, that is, i=(g+1)f s Stop at T, and finally we can get the reference fragment H in each sliding window mi With the measured fragment h mi The local correlation coefficient ρ m(i) The sliding correlation coefficient sequence ρ m , as shown below:

[0063] ρ m(0) , ρ m(1) as well as Both represent the sliding correlation coefficient sequence ρ m The data in .

[0064] Furthermore, judging the fault occurrence area based on the sliding correlation coefficient sequence includes the following steps:

[0065] Use the sliding correlation coefficient sequence as the horizontal coordinate and time as the horizontal coordinate to establish the correlation coefficient coordinate system;

[0066] obtaining an area surrounded by the sliding correlation coefficient sequence and the abscissa and ordinate of the correlation coefficient coordinate system in the preset time, to obtain a sliding correlation coefficient area;

[0067] The formula for calculating the sliding correlation coefficient area is as follows:

[0068] S m represents the sliding correlation coefficient area, T represents the preset time, t 0m represents the protection device p m the moment of starting fault protection, t represents the time point on the abscissa of the correlation coefficient coordinate system, p m (t) represents a function of calculating the local correlation coefficient at time point t; the value of T is 0.2;

[0069] The fault occurrence area is determined by comparing the sliding correlation coefficient area with a preset area threshold.

[0070] Further, the calculation formula of the preset area threshold is as follows:

[0071] S set = k rel.main S out ;

[0072] wherein, S set represents the preset area threshold, k rel.main represents the protection reliability coefficient, S out represents the sliding correlation coefficient area when the multi-terminal hybrid DC power transmission system occurs a forward external fault, the protection reliability coefficient k mrel The value range of k mrel is 0.8≤k mrel ≤0.9.

[0073] Since the reference segment is generated by the Berrethon method according to the simulation calculation of the forward external metallic fault parameters, therefore, when the forward external fault occurs, the measured segment and the reference segment should be highly consistent, so the sliding correlation coefficient sequence when the forward external fault occurs should be approximately a straight line parallel to the x-axis, and the intercept with the y-axis is 1, then the area surrounded by the sliding correlation coefficient sequence and the coordinate axis should be a rectangle with a length of 0.2 and a width of 1, so S out = 0.2 x 1 = 0.2.

[0074] Further, the specific steps of determining the fault occurrence area by comparing the sliding correlation coefficient area with the preset area threshold are as follows:

[0075] When the sliding correlation coefficient area is greater than or equal to the preset area threshold, the fault is determined as an out-zone fault; when the sliding correlation coefficient area is less than the preset area threshold, the fault is determined as an in-zone fault.

[0076] Further, the specific steps of judging whether to execute the fault protection according to the fault occurrence region are as follows:

[0077] If the fault is determined as an in-zone fault, the fault protection component is controlled to execute the fault protection operation; if the fault is determined as an out-zone fault, the fault protection component is controlled not to execute the fault protection operation.

[0078] In order to solve the above technical problems, the application further provides a multi-terminal hybrid DC power transmission system, and the specific technical content is as follows:

[0079] A multi-terminal hybrid DC power transmission system comprises,

[0080] A reference waveform calculation module is configured to simulate the occurrence of a forward out-zone fault in the protection device p m for a preset time after the occurrence of the forward out-zone fault based on the Berenger numerical calculation method according to the operating parameters of the multi-terminal hybrid DC power transmission system when no fault occurs and the forward out-zone fault parameters of the multi-terminal hybrid DC power transmission system, and obtain a reference waveform based on the current anti-travel wave derivative of the protection device p

[0081] A measured waveform calculation module is configured to calculate a measured current anti-travel wave derivative based on the Berenger numerical calculation method according to the line module voltage and the line module current of the protection device p m after the fault start criterion is started for the preset time when the multi-terminal hybrid DC power transmission system fails;

[0082] A noise adding module is configured to add disturbance noise to the measured current anti-travel wave derivative to obtain a measured waveform;

[0083] A normalization module is configured to normalize the reference waveform and the measured waveform respectively, and correspondingly obtain a normalized reference waveform and a normalized measured waveform;

[0084] A correlation coefficient calculation module is configured to calculate a sliding correlation coefficient sequence between the normalized reference waveform and the normalized measured waveform;

[0085] A fault judgment module is configured to judge the fault occurrence region according to the sliding correlation coefficient sequence;

[0086] A fault execution module is configured to judge whether to execute the fault protection according to the fault occurrence region. BRIEF DESCRIPTION OF DRAWINGS

[0087] Figure 1 The flowchart of the application is shown in the figure;

[0088] Figure 2 Topology diagram of Kunliu Long ± 800 kV three-terminal hybrid DC project in embodiment 1 of the present application;

[0089] Figure 3 Schematic diagram of inductance equivalent calculation circuit in embodiment 1 of the present application;

[0090] Figure 4 Schematic diagram of capacitance equivalent calculation circuit in embodiment 1 of the present application;

[0091] Figure 5 Schematic diagram of transmission line model in embodiment 1 of the present application;

[0092] Figure 6 Schematic diagram of transmission line equivalent calculation circuit in embodiment 1 of the present application;

[0093] Figure 7 Equivalent calculation circuit diagram of metal fault occurring at f8 in embodiment 1 of the present application;

[0094] Figure 8 Schematic diagram of reference waveform intercepted at p1 and p3 in embodiment 1 of the present application;

[0095] Figure 9 Line mode current detected at p1 and p3 after fault occurs in embodiment 1 of the present application;

[0096] Figure 10 Current reverse traveling wave and its derivative waveform detected at p1 after fault occurs in embodiment 1 of the present application;

[0097] Figure 11 Current reverse traveling wave and its derivative waveform detected at p3 after fault occurs in embodiment 1 of the present application;

[0098] Figure 12 Sliding correlation coefficient calculation process at p1 in embodiment 1 of the present application;

[0099] Figure 13 Sliding correlation coefficient calculation process at p3 in embodiment 1 of the present application;

[0100] Figure 14 Schematic diagram of area surrounded by sliding correlation coefficient sequence at p1 and coordinate axes in embodiment 1 of the present application;

[0101] Figure 15 Schematic diagram of area surrounded by sliding correlation coefficient sequence at p3 and coordinate axes in embodiment 1 of the present application;

[0102] Figure 16 Line mode current detected at p1 and p3 after fault occurs in embodiment 2 of the present application;

[0103] Figure 17The current anti-traveling wave and its derivative waveforms detected at p1 after the fault in Embodiment 2 of the present application;

[0104] Figure 18 The current anti-traveling wave and its derivative waveforms detected at p3 after the fault in Embodiment 2 of the present application;

[0105] Figure 19 The sliding correlation coefficient calculation process at p1 in Embodiment 2 of the present application;

[0106] Figure 20 The sliding correlation coefficient calculation process at p3 in Embodiment 2 of the present application;

[0107] Figure 21 The sliding correlation coefficient sequence and the coordinate axis surrounding area diagram at p1 in Embodiment 2 of the present application;

[0108] Figure 22 The sliding correlation coefficient sequence and the coordinate axis surrounding area diagram at p3 in Embodiment 2 of the present application; DETAILED DESCRIPTION

[0109] The principles and features of the present application are described below in conjunction with the accompanying drawings, and the examples are used only to explain the present application and are not intended to limit the scope of the present application.

[0110] Embodiment 1

[0111] As shown in Figure 1 , the present application provides a fault protection method for a multi-terminal hybrid DC power transmission system. The embodiment refers to the Qunliulong ± 800 kV three-terminal hybrid DC engineering topology as shown in Figure 2 , and builds an electromagnetic transient model in PSCAD / EMTDC simulation software to perform fault simulation detection to verify the effectiveness of the proposed method, wherein the MMC side adopts 70% full-bridge sub-modules and 30% half-bridge sub-modules. The overhead line in this embodiment adopts a frequency-varying parameter model to better analyze the transient state during the fault. The frequency of the protection sampling device is set to f s = 50 kHz, the simulation step is 0.02 ms, the linear interpolation point number g is 9, the sliding time window length t w = 10 μs, and the total window length T required after the start criterion is started is 0.2 ms. The protection device and the protection measuring points p1 and p3 are located at the first ends of the two DC lines.

[0112] When t = 0 ms, a positive 600 Ω ground fault occurs at the midpoint f6 of line 2, at this time, for the protection device p1 protecting the entire length of the transmission line and the protection device p3 protecting the single section of the line, it is a positive direction fault within the zone, and p1 and p3 are required to reliably act.

[0113] S1, based on the Bergeron numerical calculation method, according to the operating parameters of the multi-terminal hybrid DC power transmission system and the parameters of the external fault in the forward direction, the current anti-travel wave derivative collected at the protection device p1 and the protection device p3 within T=0.2ms after the external fault in the forward direction is simulated to obtain a reference waveform, and the specific steps are as follows:

[0114] The Bergeron numerical calculation method is a method for analyzing the multiple refraction or reflection process of waves by applying the concept of mixed waves, and its essence is still the solution of the wave equation. The core of the Bergeron numerical calculation method is to equivalent the distributed parameter components to lumped parameter components, so as to calculate the wave process on the line by using the commonly used solution method of lumped parameter components. Similarly, the inductance and capacitance of the lumped parameter components in the circuit also need to be converted into the corresponding numerical calculation circuit according to the requirements of numerical calculation. As shown in Figure 3 , L represents the lumped parameter inductance in the system parameters, and it is assumed that a voltage incident wave is injected from the m side at time t. At this time, the voltage drop of the inductance and the current flowing through the inductance have the following relationship:

[0115]

[0116] Where, u m (t) represents the voltage at the m side at time t, i m (t) represents the current at the m side at time t; u n (t) represents the voltage at the n side at time t, and i n (t) represents the current at the n side at time t.

[0117] Integrating formula (1) from t-Δt to t obtains the following formula:

[0118]

[0119] When Δt is small enough, formula (2) is approximately equal to the following formula:

[0120]

[0121] Where, I L (t-Δt) represents the equivalent current of the inductance L at (t-Δt) time.

[0122] Similarly, as shown in Figure 4 , the equivalent calculation formula of the capacitance is as follows:

[0123]

[0124] Where, I C (t-Δt) represents the equivalent current of the capacitance C at (t-Δt) time.

[0125] As shown in Figure 5 , let the wave impedance of the transmission line Line-x be Zx , length l x .

[0126] According to the law of mixed wave transmission, the mixed wave u x (t-τ m )+i x (t-τ mn ) emitted by the m side at t-τ x time will reach the n side after τ x milliseconds, where τ x represents the transmission time of the traveling wave in Line-x, and τ x = l / v x ms, where v x represents the wave speed of the traveling wave. Therefore, the voltage u n (t) and the current i nm (t) at t time on the n side can be represented by the voltage and the current at t-τ x time on the m side. The formula is as follows:

[0127]

[0128] where U nm (t-τ x ) is an equivalent voltage source determined by the voltage and the current at t-τ x time on the m side. Since the m side and the n side are mutually opposite, the voltage u n (t) and the current i mn (t) at t time on the m side can also be represented by the voltage and the current at t-τ x time on the n side. The equivalent result is shown in Figure 6 .

[0129] For a bipolar DC power transmission system, since there is a coupling relationship between the positive and negative lines, the following formula is usually used for decoupling processing:

[0130]

[0131] where u and i represent the positive / negative voltage / current of node j; u 0j and u 1j represent the zero-mode / line-mode voltage of node j; and i 0j and i 1j represent the zero-mode / line-mode current of node j.

[0132] Compared with the line-mode component, the zero-mode component is only generated during a ground fault, and the attenuation degree of the zero-mode component in the line is also greater than that of the line-mode, which will have a certain impact on the sensitivity of the protection. Therefore, under the premise of considering inter-pole faults and improving the sensitivity of the protection, the line-mode component is selected to construct the protection criterion.

[0133] In summary, distributed parameter components such as capacitors, inductors, and transmission lines are replaced with lumped parameter components, and the topology of the hybrid multi-terminal DC transmission system is reconstructed. The system operating parameters and the fault location parameters outside the forward zone of the system are substituted, and the node voltage method is used to simulate and calculate the protection device p within 0.2ms when the most serious fault outside the forward zone occurs. m The numerical solution of the current reverse wave derivative is loaded as a reference waveform at the line protection device. Taking the forward area protected by p1 and p3 as an example, the topology diagram when a metallic fault occurs outside the forward area is as follows: Figure 7 As shown, the reference waveform I' extracted from the simulation calculation results mH like Figure 8 shown.

[0134] S2. When a fault occurs in a DC line, the fault traveling wave will cause a change in the DC current when it is transmitted to the protection device. Taking advantage of the characteristics of line mode current, such as fast propagation speed and small fluctuation, the line mode current change rate is used as the protection triggering criterion. The formula for the line mode current change rate is as follows:

[0135]

[0136] Among them, I 1mset Indicates the position of each protection device when the system has the slightest fault in the forward zone m The detected line mode current mutation, k rel.start Represents the reliability coefficient of the starting component. When the starting criterion is met, the protection starts and the starting time is recorded as the arrival time of the fault traveling wave t 0m .

[0137] In order to ensure that the starting component based on the current change rate can start correctly under the slightest fault condition of the system, according to the PSCAD simulation model results of the three-terminal hybrid DC system, the threshold value is set to the line mode current change rate detected by the protection device under the scenario of a single-pole 600Ω grounding fault at the end of the line f7 in the area. 11set =123.82A, I 13set =89.19A. In order to ensure the sensitive start of the starting criterion, the reliability coefficient k of the starting component is set at the expense of selectivity in exchange for sensitivity. rel.start =0.8.

[0138] In this embodiment, the actual line mode current waveform I detected when p1 and p3 fail is 11 and I 13 like Figure 9 As shown, it can be seen that the startup components configured by p1 and p3 are respectively 01 =4.02ms, t 03 =Start at 0.91ms.

[0139] Calculate t 0m The measured current anti-travel wave I mh , and calculate the corresponding measured current anti-travel wave derivative I' mh , the calculation formula is as follows:

[0140]

[0141] I' mh = I mh (t)-I mh (t-1 / f s ), t ∈ [t 0m ,t 0m +0.2];

[0142] Where, U 1m represents the line mode voltage at the installation site p m of the protection device, I 1m represents the line mode current at the installation site p m of the protection device, and Z 1i is the line mode wave impedance of the line where the protection point is located. Figure 10 and Figure 11 are the current anti-travel wave and its derivative detected at p1 and p3, respectively.

[0143] S3, without using the directional element to avoid reverse faults, in order to make the main protection criterion reliably not act in the case of reverse faults, according to the characteristics that the current anti-travel wave and its derivative tend to 0 within T=0.2ms after the protection device starts the starting criterion in the case of reverse faults, a disturbance noise α times the same as the reference waveform is added to the measured current anti-travel wave derivative; α is defined as the noise addition multiple of the reference waveform, and in the embodiment, α is 1.

[0144] Through the above noise addition method, the system can make the fault protection criterion more reliable in the case of reverse faults without using the directional element to avoid reverse faults, that is, the fault protection device does not act in the case of reverse faults.

[0145] The current anti-travel wave derivative I' mh adds a disturbance noise to obtain the measured waveform The calculation formula is as follows:

[0146]

[0147] I′ mh(k) , I′ mH(k) , respectively represent the kth element in each sequence; represents the kth element data in the measured waveform, I' mh(k)represents the kth element data in the current anti-travel derivative, I' mH(k) represents the kth element data in the reference waveform.

[0148] S4, normalize the reference waveform and the measured waveform respectively, and correspondingly obtain a normalized reference waveform and a normalized measured waveform, including the following steps:

[0149] S401, before normalization, the waveform needs to be data augmented, first, the reference waveform I' mH and the measured waveform are linearly interpolated, and the formula is as follows:

[0150]

[0151] Where, x N represents the abscissa value of the Nth point inserted by linear interpolation, y N represents the ordinate value of the nth point inserted by linear interpolation, x a represents the abscissa value of the starting point of the linear interpolation interval, y a represents the ordinate value of the starting point of the linear interpolation interval, x b represents the abscissa value of the terminal point of the linear interpolation interval, y b represents the ordinate value of the terminal point of the linear interpolation interval.

[0152] In order to facilitate the sliding correlation analysis, data in a time window before the starting time needs to be reserved, and since the fault traveling wave has not arrived in this period, the current anti-travel wave data in this period is set to zero. Denote the linearly interpolated reference current anti-travel derivative of the protection device p m The reference current anti-travel derivative and the measured current anti-travel derivative at the starting time t 0m are I' mH(0) and I' mh(0) respectively, and the data from (t 0m -t w ) to t 0m has the following relationship.

[0153]

[0154]

[0155] Denote the augmented reference waveform I' mtem and the augmented measured waveform I' mact obtained by the above method, as shown below.

[0156]

[0157]

[0158] S402, respectively, the reference waveform I' after data expansion mtem And the measured waveform I' after data expansion mact Normalization is performed. The normalization process is as follows

[0159]

[0160] Get the normalized reference waveform i' mref and normalized measured waveform i' mact As shown below

[0161]

[0162]

[0163] S5. Calculating the sliding correlation coefficient between the normalized reference waveform and the normalized measured waveform to obtain a sliding correlation coefficient sequence, including the following steps:

[0164] S501, initialize the sliding coefficient to i=0, and use t w =10μs is the sliding correlation coefficient calculation window length, respectively from the normalized reference waveform i' mtem and normalized measured waveform i' mact The two fragments are respectively t 0m -t w Sequence value corresponding to the moment and is the starting point, t 0m The sequence value i' corresponding to the moment mH(0+i) and i' mh(0+i) As the end point, the window length is t w Reference fragment H mi With the measured fragment h mi As shown below

[0165]

[0166]

[0167] S502, reference segment H mi With the measured fragment h mi The internal data are arranged in descending order to obtain the reference descending sequence A mi and the measured descending sequence B mi As shown below

[0168]

[0169]

[0170] Reference fragment Hmi The zth element in the reference descending sequence A mi The position ranking in is recorded as r z , which is called reference segment H mi The rank of the z-th element point, so that the entire reference segment H can be obtained mi The rank sequence r of the internal data mi Similarly, record the measured segment h mi The zth element point in the measured descending sequence B mi The position ranking in is recorded as s z , we can get the measured fragment h mi The rank sequence s composed of all data in mi . The sequence r mi With the sequence s mi Subtract each element in the rank difference sequence d mi As shown below:

[0171]

[0172] Substituting it into the following formula, we can get the local correlation coefficient ρ when i=0 m(i) ;

[0173]

[0174] S503, maintaining the sliding time window length t of the intercepted waveform w The sliding window remains unchanged. Every time the sliding window slides back one data point, that is, i=i+1, steps S501 and S502 are executed once to calculate the local correlation coefficient ρ within the sliding window. m(i) , until the time window continues to slide through all data within T = 0.2ms, that is, i = (g + 1)f s Stop at T. The calculation process of sliding correlation coefficient at p1 and p3 is as follows Figure 12 and Figure 13 shown.

[0175] Finally, we can get the reference fragment H in each sliding window mi With the measured fragment h mi The sliding correlation coefficient sequence ρ composed of the local correlation coefficient m .

[0176]

[0177] S6. Determine the fault occurrence area based on the sliding correlation coefficient sequence;

[0178] Specifically, according to the sliding correlation coefficient sequence ρ m , determine the fault occurrence area, including the following steps:

[0179] S601, establish a coordinate system with the local correlation coefficient as the ordinate and time as the abscissa;

[0180] S602, introduce the sliding correlation coefficient sequence p mtem between the normalized reference waveform i' mact and the normalized measured waveform i' m into the coordinate system to obtain a sliding correlation coefficient curve;

[0181] S603, calculate the area enclosed by the sliding correlation coefficient curve and the coordinate axis in the coordinate system to obtain the sliding correlation coefficient area S m , and the calculation formula is as follows

[0182]

[0183] S604, judge the fault occurrence area by comparing the sliding correlation coefficient area S m with a preset area threshold S set , and the calculation formula of the preset area threshold is as follows:

[0184] S set =k rel.main S out ;

[0185] By introducing the protection reliability coefficient, the protection reliability coefficient value can be set as needed to improve the reliability of fault protection. The protection reliability coefficient k rel.main has a value range as follows: 0.8≤k rel.main ≤0.9. Specifically, the value of k rel.main is preferably 0.85.

[0186] Since the reference segment is generated by the Bergeron method according to the forward out-of-zone fault position parameter simulation calculation, the measured segment and the reference segment should be highly consistent when the actual forward out-of-zone fault occurs. Therefore, the sliding correlation coefficient sequence when the forward out-of-zone fault occurs should be a straight line, and the area enclosed by it and the coordinate axis should be a rectangle with a length of 0.2 and a width of 1, that is, S out =0.2×1=0.2, and S set =0.85×0.2=0.17.

[0187] When the sliding correlation coefficient area S m is greater than or equal to the preset area threshold S set , it is determined that the fault is an out-of-zone fault; when the sliding correlation coefficient area S m is less than the preset area threshold S set , it is determined that the fault is an in-zone fault.

[0188] In this embodiment, the sliding correlation coefficient sequence ρ of p1 and p3 is obtained by sliding window p1 and ρ p3 The area enclosed by the coordinate system is as follows Figure 14 and Figure 15 As shown, ρ is obtained by the area calculation formula of S603 sliding correlation coefficient p1 and ρ p3 The areas enclosed by the coordinate axes are S p1 =0.1273, S p1 =0.1200, both are less than the preset area threshold S set =0.17, so both p1 and p3 are determined to have an intra-zone fault.

[0189] S7. Determine whether to execute fault protection based on the fault occurrence area.

[0190] If the fault is determined to be an intra-zone fault, the fault protection component performs a fault protection operation; if the fault is determined to be an extra-zone fault, the fault protection component does not perform a fault protection operation.

[0191] In this embodiment, both p1 and p3 are determined to have an intra-zone fault, so the fault protection components at p1 and p3 both perform fault protection operations, which is consistent with theoretical analysis and meets the protection sensitivity requirement.

[0192] Example 2

[0193] Maintain the same system configuration as in Example 1. Assume that at t = 0ms, the current limiting reactor L A A positive metallic ground fault occurs at line f8, connected to the valve side of converter MMC-A. This fault is considered out of the forward zone for both protection device p1, which protects the entire length of the transmission line, and protection device p3, which protects a single section of the line. Both devices must reliably remain inoperative. Reference waveforms are obtained using the same method as in Example 1 and applied to protection devices p1 and p3.

[0194] At this time, the measured line mode current waveform I detected by p1 and p3 11 and I 13 like Figure 16 As shown. It can be seen that the startup components configured by p1 and p3 are respectively 01 =4.96ms, t 03 =Start at 1.84ms.

[0195] Intercept the line mode voltage and line mode current data within 0.2ms after the startup component is started, and calculate the measured current reverse wave I mh , and calculate the corresponding measured current reverse wave derivative I' mh , the calculation results at p1 and p3 are as follows Figure 17 and Figure 18 shown.

[0196] current anti-travel wave derivative I' mh adding disturbance noise to obtain the measured waveform to the measured waveform and the reference waveform I' mH performing the data expansion operation to obtain the data-expanded reference waveform I' mtem and the expanded measured waveform I' mact . Then the data-expanded waveforms are normalized respectively to obtain the normalized reference waveform i' mtem and the normalized measured waveform i' mact .

[0197] the normalized reference waveform i' mtem and the normalized measured waveform i' mact performing sliding correlation coefficient calculation, the calculation processes at p1 and p3 are respectively as shown in Figure 19 and Figure 20 . Then the sliding correlation curves at p1 and p3 are as shown in Figure 21 and Figure 22 , and ρ p1 and ρ p3 are calculated. p1 The area S p3 = 0.1981, both are greater than the setting value S set = 0.17, and the protections are not actuated, which is consistent with the theoretical analysis and meets the safety requirements of the protection.

[0198] The above only describes the preferred embodiments of the present application and is not used to limit the present application. Any modification, equivalent replacement, improvement, etc. within the concept and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A fault protection method for a multi-terminal hybrid direct current transmission system, characterized in that: The steps include: Based on the Berelon numerical calculation method, according to the operating parameters of the multi-terminal hybrid DC transmission system when no fault occurs and the metallic fault parameters of the multi-terminal hybrid DC transmission system outside the forward zone, the protection device is simulated within a preset time after the forward zone fault occurs. The reverse wave derivative of the current at is used to obtain the reference waveform; When the multi-terminal hybrid HVDC transmission system fails, based on the Berelon numerical calculation method, the protection device is activated within the preset time after the fault start criterion is activated. The line mode voltage and line mode current at the position are used to calculate the measured current reverse wave derivative; Add disturbance noise to the measured current reverse wave derivative to obtain the measured waveform; Normalizing the reference waveform and the measured waveform respectively to obtain a normalized reference waveform and a normalized measured waveform respectively; Calculating a sliding correlation coefficient sequence between a normalized reference waveform and a normalized measured waveform; Determine the fault occurrence area based on the sliding correlation coefficient sequence; Determine whether to execute fault protection based on the fault occurrence area.

2. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 1, characterized in that: Adding disturbance noise to the measured current reverse wave derivative specifically includes adding the current reverse wave derivative as disturbance noise to the measured current reverse wave derivative; the calculation formula of the measured waveform is as follows: ； in, represents the measured waveform, represents the current reverse wave derivative, represents the reference waveform, represents the disturbance noise coefficient, The value range is between 0 and 5.

3. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 1 or 2, characterized in that: Normalize the reference waveform and the measured waveform respectively to obtain the normalized reference waveform and normalizing the measured waveform, including the following steps: Perform data expansion on the reference waveform and the measured waveform respectively; Normalizing the reference waveform and the measured waveform after data expansion respectively, to obtain the normalized reference waveform and the normalized measured waveform respectively; The normalization formula is as follows: ； represents the normalized reference waveform, represents the normalized measured waveform, represents the reference waveform after data expansion, Represents the measured waveform after data expansion.

4. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 3, characterized in that: Data expansion is performed on the reference waveform and the measured waveform respectively, including the following steps: The reference waveform and the measured waveform are interpolated respectively using a linear interpolation method; wherein the interpolation formula of the linear interpolation method is as follows: ； Indicates the linear interpolation insertion The horizontal coordinate value of the point, Indicates the linear interpolation insertion The vertical coordinate value of the point, Indicates the horizontal coordinate value of the starting point of the linear interpolation interval, Indicates the vertical coordinate value of the starting point of the linear interpolation interval, Indicates the horizontal coordinate value of the end point of the linear interpolation interval, Indicates the ordinate value of the end point of the linear interpolation interval, Indicates the number of points of linear interpolation between two adjacent data points; the horizontal axis value is a time series, and the vertical axis value is a waveform sequence; The data in the first preset time window of the interpolated reference waveform is set to zero to obtain the reference waveform after data expansion; wherein the first preset time window is the protection device p m The time window before the start time, the length of the first preset time window is less than the preset time.

5. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 4, characterized in that: Calculating a sliding correlation coefficient sequence between a normalized reference waveform and a normalized measured waveform includes the following steps: S101, initialization sliding coefficient for , S102, with t w The calculation time window length of the sliding correlation coefficient is respectively from the normalized reference waveform and normalized measured waveform Extract the comparison fragment and get the corresponding reference fragment and measured fragments ; wherein the reference fragment The amount of data and the measured fragment The number of data in is equal; S103, respectively, the reference fragments and the measured fragment The data in is arranged in descending order, corresponding to the reference descending sequence and the measured descending sequence ; S104, finding the reference descending sequence With the measured descending sequence The difference between the data of the same order in the rank difference sequence ; S105, according to the rank difference sequence Calculate the local correlation coefficient using the data in ; Among them, the local correlation coefficient The calculation formula is as follows: ； Specifically, Indicates the protective device Collect the frequency of line mode voltage and line mode current, Represents the rank difference sequence The data in Represents the rank difference sequence The order value of the data in ; S106, set the sliding coefficient , loop through steps S102 to S105 until the sliding coefficient satisfy When , step S107 is executed; wherein, Indicates the preset time; S107, the local correlation coefficient obtained each time The data are integrated into a data sequence to obtain the sliding correlation coefficient sequence.

6. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 1, characterized in that: According to the sliding correlation coefficient sequence, the fault occurrence area is determined, including the following steps: Use the sliding correlation coefficient sequence as the horizontal coordinate and time as the horizontal coordinate to establish the correlation coefficient coordinate system; Calculating the area enclosed by the abscissa and the abscissa of the sliding correlation coefficient sequence and the correlation coefficient coordinate system within the preset time to obtain a sliding correlation coefficient area; The formula for calculating the sliding correlation coefficient area is as follows: ; represents the sliding correlation coefficient area, Indicates the preset time, Indicates the protective device When the fault protection is activated, represents the time point on the abscissa of the correlation coefficient coordinate system, Indicates the calculation time point A function of the local correlation coefficient at ; The fault occurrence area is determined by comparing the sliding correlation coefficient area with a preset area threshold.

7. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 6, characterized in that: The calculation formula of the preset area threshold is as follows: ； in, represents the preset area threshold, represents the protection reliability factor, It represents the sliding correlation coefficient area when a fault occurs outside the forward zone in the multi-terminal hybrid DC transmission system, and the protection reliability coefficient range is .

8. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 6, characterized in that: The specific steps for determining the fault area by comparing the sliding correlation coefficient area with the preset area threshold are as follows: When the sliding correlation coefficient area is greater than or equal to the preset area threshold, the fault is determined to be an out-of-zone fault; when the sliding correlation coefficient area is less than the preset area threshold, the fault is determined to be an in-zone fault.

9. The fault protection method for a multi-terminal hybrid direct current transmission system according to claim 8, characterized in that: The specific steps for determining whether to execute fault protection based on the fault occurrence area are as follows: If the fault is determined to be an internal fault, the fault protection component is controlled to perform a fault protection operation; if the fault is determined to be an external fault, the fault protection component is controlled not to perform a fault protection operation.

10. A multi-terminal hybrid direct current transmission system, characterized in that: include, The reference waveform calculation module is used to simulate the protection device within a preset time after the occurrence of a forward out-of-zone fault based on the operating parameters of the multi-terminal hybrid DC transmission system when no fault occurs and the forward out-of-zone metallic fault parameters of the multi-terminal hybrid DC transmission system based on the Berelon numerical calculation method. The reverse wave derivative of the current at is used to obtain the reference waveform; The measured waveform calculation module is used to, when the multi-terminal hybrid DC transmission system fails, calculate the protection device according to the protection device according to the Berelon numerical calculation method within the preset time after the fault start criterion is started. The line mode voltage and line mode current at the position are used to calculate the measured current reverse wave derivative; A noise adding module is used to add disturbance noise to the measured current reverse wave derivative to obtain the measured waveform; A normalization module is used to normalize the reference waveform and the measured waveform respectively, and obtain a normalized reference waveform and a normalized measured waveform respectively; A correlation coefficient calculation module, used to calculate a sliding correlation coefficient sequence between a normalized reference waveform and a normalized measured waveform; A fault judgment module is used to judge the fault occurrence area based on the sliding correlation coefficient sequence; The fault execution module is used to determine whether to execute fault protection according to the fault occurrence area.

Images