A power distribution network reliability intelligent analysis method and analysis system

Abstract

The application discloses a power distribution network reliability intelligent analysis method and an analysis system, and particularly relates to the technical field of power system fault detection, and is used for solving the misjudgment and missed judgment problems caused by weak fault characteristics and harmonic interference under the bidirectional power flow scene of the existing method; first, the fundamental wave and harmonic components of the power distribution network nodes are synchronously collected, the current amplitude change rate and phase offset are extracted by adopting time-frequency domain joint analysis, and the dynamic distribution characteristics are generated in combination with harmonic energy clustering; by quantifying the phase offset direction and the harmonic energy gathering characteristics, the high-frequency interference components under the bidirectional power flow mode are identified; further, based on the harmonic frequency band coupling effect and the current dynamic correlation, an adaptive detection threshold function is constructed, and the accurate separation of the fault characteristics is realized; finally, the reliability evaluation parameters are generated in combination with the real-time fault judgment result and the historical repair data; and the sensitivity of the fault detection and the rationality of the reliability evaluation under the complex operation environment are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of power system fault detection technology, and more specifically, to a method and system for intelligent analysis of distribution network reliability. Background Technology

[0002] The large-scale integration of distributed energy resources in power distribution networks results in bidirectional power flow characteristics. Traditional fault detection methods rely on steady-state analysis of the amplitude and waveform characteristics of power frequency current. Existing technologies typically construct reliability assessment models based on the pre-defined conditions of unidirectional power flow, assuming that the fault current has a clear directionality and amplitude threshold, and achieving rapid isolation through overcurrent protection devices. However, the inverter control strategies of distributed power sources actively limit the output of fault current, leading to a significant weakening of fault characteristics in both amplitude and phase, accompanied by high-frequency harmonic interference.

[0003] Existing methods struggle to accurately identify weakened fault current characteristics in bidirectional power flow scenarios, leading to misjudgments or missed detections by protection devices. Due to the decreased signal-to-noise ratio and time-varying harmonics of fault signals, traditional power frequency detection mechanisms cannot effectively distinguish between real faults and normal operating condition fluctuations, causing reliability assessment results to deviate from actual operational risks and affecting the timeliness and accuracy of fault isolation. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, embodiments of the present invention provide a method and system for intelligent analysis of distribution network reliability to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A method for intelligent analysis of distribution network reliability includes the following steps:

[0007] S1. Obtain real-time current data of each node in the distribution network. The real-time current data includes the fundamental component and harmonic components.

[0008] S2. Perform time-frequency domain joint analysis on the fundamental component, extract the current amplitude change rate and phase offset, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics.

[0009] S3. Determine whether the power flow direction is bidirectional based on the relationship between the phase offset and the preset direction threshold.

[0010] S4. If it is a two-way mode, perform time-varying characteristic analysis on the high-frequency harmonic components based on the harmonic energy distribution characteristics, and identify the interference components in the high-frequency harmonic components that overlap with the fault characteristic frequency band.

[0011] S5. Based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude, an adaptive detection threshold function is constructed.

[0012] S6. Input the rate of change of current amplitude into the adaptive detection threshold function, output the fault judgment result, and generate distribution network reliability assessment parameters based on the fault judgment result.

[0013] In a preferred embodiment, real-time current data of each node in the distribution network is acquired. The real-time current data includes fundamental and harmonic components, including:

[0014] The three-phase current signals of each node in the distribution network are collected synchronously, a unified timestamp is added to the collection time, and the three-phase current signals are bandpass filtered to separate the fundamental component from the harmonic components within the preset frequency band.

[0015] In a preferred embodiment, joint time-frequency domain analysis is performed on the fundamental component to extract the rate of change of current amplitude and phase shift, and spectral clustering is performed on the harmonic components to generate harmonic energy distribution characteristics, including:

[0016] A short-time Fourier transform is performed on the fundamental component to generate a time-frequency matrix. The rate of change of current amplitude within a preset time window and the phase offset between adjacent sampling points are extracted from the time-frequency matrix.

[0017] A spectral energy distribution matrix containing frequency points and energy values ​​is constructed for the harmonic components. Cluster centers are initialized based on the gradient distribution characteristics of the frequency point energy values ​​in the spectral energy distribution matrix. The harmonic energy distribution characteristics representing the harmonic energy accumulation region are generated by iterative optimization of the cost function that minimizes the distance between the frequency point energy and the cluster center.

[0018] In a preferred embodiment, determining whether the power flow direction is bidirectional based on the relationship between the phase offset and a preset direction threshold includes:

[0019] The phase offset is divided into a forward offset sequence and a reverse offset sequence, and the cumulative gradient of the forward offset sequence and the cumulative gradient of the reverse offset sequence are calculated respectively.

[0020] Directional weighting coefficients are generated based on the ratio of the cumulative gradient of the forward offset sequence to the cumulative gradient of the backward offset sequence.

[0021] If the directional weight coefficient exceeds the first preset threshold and the cumulative gradients of the positive offset sequence and the negative offset sequence have opposite signs, then the power flow direction is determined to be bidirectional.

[0022] In a preferred embodiment, if it is a bidirectional mode, time-varying characteristic analysis of high-frequency harmonic components is performed based on harmonic energy distribution characteristics to identify interference components in the high-frequency harmonic components that overlap with the fault characteristic frequency band, including:

[0023] If the power flow direction is determined to be bidirectional, calculate the correlation coefficient between the energy of each frequency band of the high-frequency harmonic component and the fault characteristic frequency band.

[0024] The time length of the dynamic sliding window is adjusted based on the time-series change rate of the correlation coefficient, and the initial length of the dynamic sliding window is the preset time window for spectral clustering in step S2.

[0025] Within a dynamic sliding window, energy mutation points of high-frequency harmonic components are detected. If the overlap rate between the frequency band of the energy mutation point and the fault characteristic frequency band exceeds a second preset threshold, it is determined to be an interference component.

[0026] In a preferred embodiment, the correlation coefficient is calculated as the ratio of the covariance to the standard deviation of the frequency band energy within a preset time window.

[0027] In a preferred embodiment, an adaptive detection threshold function is constructed based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude, including:

[0028] Based on the frequency band information of the interference components, recursive quantization analysis is performed on the harmonic components of adjacent frequency bands to generate a coupling matrix characterizing the nonlinear coupling strength between frequency bands.

[0029] Based on the gradient distribution characteristics of the rate of change of current amplitude, the dynamic correlation weight coefficient is calculated;

[0030] The core parameters of the adaptive detection threshold function are generated by performing tensor multiplication on the coupling matrix and the dynamically associated weight coefficients.

[0031] The sensitivity coefficient of the threshold function is adjusted based on the difference between the core parameters and the preset benchmark threshold.

[0032] In a preferred embodiment, the rate of change of current amplitude is input into an adaptive detection threshold function, a fault determination result is output, and distribution network reliability assessment parameters are generated based on the fault determination result, including:

[0033] Multi-scale fluctuation decomposition is performed on the rate of change of current amplitude to extract high-frequency fluctuation components and low-frequency trend components.

[0034] Based on the core parameters and sensitivity coefficient of the adaptive detection threshold function, the dynamic confidence threshold of the high-frequency fluctuation component and the steady-state deviation threshold of the low-frequency trend component are calculated respectively.

[0035] If the amplitude of the high-frequency fluctuation component exceeds the dynamic confidence threshold and the slope of the low-frequency trend component exceeds the steady-state deviation threshold, it is determined to be a persistent fault.

[0036] Based on the determination of persistent faults, and combined with the repair time and impact range parameters of similar faults in the historical fault database, distribution network reliability assessment parameters are generated.

[0037] In a preferred embodiment, the power distribution network reliability assessment parameters include the mean repair time index and the fault impact factor.

[0038] On the other hand, the present invention provides an intelligent analysis system for the reliability of power distribution networks, comprising the following modules:

[0039] The current monitoring module is used to acquire real-time current data of each node in the distribution network. The real-time current data includes the fundamental component and harmonic components.

[0040] The harmonic clustering module is used to perform joint time-frequency domain analysis on the fundamental component, extract the current amplitude change rate and phase shift, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics.

[0041] The direction determination module is used to determine whether the power flow direction is bidirectional based on the relationship between the phase offset and the preset direction threshold.

[0042] The interference identification module is used, in bidirectional mode, to perform time-varying characteristic analysis on high-frequency harmonic components based on harmonic energy distribution characteristics, and to identify interference components in high-frequency harmonic components that overlap with the fault characteristic frequency band.

[0043] The threshold construction module is used to construct an adaptive detection threshold function based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude.

[0044] The reliability assessment module is used to input the rate of change of current amplitude into the adaptive detection threshold function, output the fault determination result, and generate distribution network reliability assessment parameters based on the fault determination result.

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

[0046] 1. By employing multi-dimensional feature fusion and a dynamic threshold mechanism, the accuracy of fault detection and the rationality of reliability assessment in complex operating environments are significantly improved. Traditional methods struggle to accurately capture true fault characteristics in bidirectional power flow scenarios due to weakened fault signals and harmonic interference. This solution innovatively introduces joint time-frequency domain analysis and harmonic energy distribution clustering technology to separate amplitude change rate, phase shift, and energy accumulation characteristics reflecting the dynamic characteristics of power flow from the fundamental and harmonic components. By establishing a correlation model between the phase shift direction and harmonic energy distribution, the coupling interference between high-frequency harmonics and the fault frequency band in bidirectional mode is effectively identified, overcoming the limitations of traditional single time-domain or frequency-domain analysis in analyzing complex signal characteristics. Furthermore, a dynamically adaptive detection threshold is constructed based on nonlinear coupling effects, enabling the threshold function to adjust in real time according to harmonic intensity and fault transient processes. This solves the core problem of poor adaptability of fixed thresholds to time-varying operating conditions, significantly reducing the risk of misjudgment and missed judgment.

[0047] 2. By dynamically linking fault determination results with reliability assessment parameters, quantifiable evaluation of operational risks is achieved. Traditional assessment models rely on historical statistical values, making it difficult to reflect the actual impact of real-time faults on the power grid. By extracting dynamic features such as harmonic energy accumulation areas and fault frequency band overlap rates, and combining them with fault confidence scores output by adaptive thresholds, reliability parameters closely matching the current operating state are generated. This assessment mechanism, based on real-time signal characteristics and fault development trends, can accurately characterize the differentiated impact of different fault types on the reliability of distribution network power supply, providing data support for operation and maintenance decisions that balances timeliness and accuracy, and significantly improving the level of power grid resilience management. Attached Figure Description

[0048] Figure 1 This is a flowchart of an intelligent analysis method for power distribution network reliability according to the present invention;

[0049] Figure 2 This is a schematic diagram of the structure of an intelligent analysis system for power distribution network reliability according to the present invention. Detailed Implementation

[0050] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0051] Example 1: Figure 1 This invention presents an intelligent analysis method for the reliability of power distribution networks, which includes the following steps:

[0052] S1. Obtain real-time current data of each node in the distribution network. The real-time current data includes the fundamental component and harmonic components.

[0053] S2. Perform time-frequency domain joint analysis on the fundamental component, extract the current amplitude change rate and phase offset, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics.

[0054] S3. Determine whether the power flow direction is bidirectional based on the relationship between the phase offset and the preset direction threshold.

[0055] S4. If it is a two-way mode, perform time-varying characteristic analysis on the high-frequency harmonic components based on the harmonic energy distribution characteristics, and identify the interference components in the high-frequency harmonic components that overlap with the fault characteristic frequency band.

[0056] S5. Based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude, an adaptive detection threshold function is constructed.

[0057] S6. Input the rate of change of current amplitude into the adaptive detection threshold function, output the fault judgment result, and generate distribution network reliability assessment parameters based on the fault judgment result.

[0058] S1. Obtain real-time current data for each node in the distribution network. The real-time current data includes the fundamental component and harmonic components. The specific implementation is as follows:

[0059] When acquiring real-time current data from each node of the distribution network, the three-phase current signals are first synchronously collected by current sensors installed at each node. All current sensors are of the same model and have time synchronization capabilities. Time synchronization is achieved through a GPS clock module or the IEEE 1588 protocol. The GPS clock module has a time synchronization accuracy of ±1 microsecond, while the IEEE 1588 protocol has an accuracy of ±100 nanoseconds, ensuring that the deviation in acquisition time between different nodes is less than 1 millisecond. After acquisition, a unified timestamp is added to each sampling point of the three-phase current signal. The format of the unified timestamp is year, month, day, hour, minute, second, and millisecond, for example, "2023-08-25 14:30:05.123". The timestamped three-phase current signals are then uploaded to the data processing center via a network transmission protocol.

[0060] When performing bandpass filtering on three-phase current signals, a Butterworth filter with an adjustable cutoff frequency is used. The passband range of the filter is set according to the actual harmonic distribution characteristics of the distribution network. The passband range of the fundamental component is set to 45 Hz to 55 Hz, and the passband range of the harmonic components is set to 2 kHz to 5 kHz. The filtered signal is verified for frequency band separation effect through fast Fourier transform. If the residual harmonic energy in the fundamental component exceeds the preset threshold, and the preset threshold for the proportion of harmonic energy in the fundamental passband is 5%, the filter order or cutoff frequency is automatically adjusted and the filter is re-filtered. The filter order is set to 8th order by default and is increased to 12th order during adjustment. The cutoff frequency is adjusted in 0.5 kHz increments.

[0061] The separated fundamental and harmonic components are stored in independent buffers. The buffers are stored as timestamp-aligned discrete sampling sequences, with a sampling rate of 1 kHz for the fundamental component and 10 kHz for the harmonic components, to meet the resolution requirements of subsequent joint time-frequency domain analysis and spectral clustering. During data preprocessing, if the amplitude of the three-phase current signal at a node exceeds the sensor's range (upper limit 1000 amperes), an anomaly handling mechanism is triggered. This mechanism discards the current data for that node and initiates data interpolation compensation for adjacent nodes. The interpolation method is either linear extrapolation or moving average prediction based on historical data, with a window length of 10 sampling points for the moving average prediction. After data acquisition and filtering for all nodes, a time-aligned dataset of the fundamental and harmonic components is output. The dataset is in the format of a multidimensional array, with each dimension corresponding to the time-series data of one node. The array index is consistent with the node's physical topology order, which is arranged according to the feeder branch order, providing standardized input for subsequent steps.

[0062] S2. Perform joint time-frequency domain analysis on the fundamental component to extract the current amplitude change rate and phase shift, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics. Specifically, this is implemented as follows:

[0063] When performing joint time-frequency domain analysis on the fundamental component, a short-time Fourier transform is first performed on the fundamental component to generate a time-frequency matrix. The short-time Fourier transform uses a Hanning window as the window function. The window length of the Hanning window is set to 10 milliseconds, which is an integer multiple of the fundamental period. The overlap rate of adjacent windows is 50% of the window length. The rows of the time-frequency matrix correspond to the frequency components, and the interval between the frequency components is 1% of the fundamental frequency. The columns correspond to the timestamps of the time windows, which are aligned with the unified timestamps in step S1. The matrix element values ​​of the time-frequency matrix are the current amplitudes at the corresponding frequencies and time windows. The current amplitudes are obtained by calculating the complex modulus, which is calculated as the square root of the sum of the squares of the real and imaginary parts of the complex number.

[0064] When extracting the rate of change of current amplitude and the phase offset between adjacent sampling points within a preset time window from the time-frequency matrix, the length of the preset time window is 5 milliseconds and is consistent with the sampling time interval of the fundamental component in step S1. The method for calculating the rate of change of current amplitude is the ratio of the difference between the fundamental amplitudes in two adjacent time windows to the time interval, where the time interval is 1 millisecond. The method for calculating the phase offset is the phase angle difference of the same frequency component in two consecutive time windows. The phase angle difference is obtained by calculating the arctangent values ​​of the real and imaginary parts of the complex number and finding their difference. The four-quadrant arctangent function is used in the calculation to avoid phase jump error. The input parameters of the four-quadrant arctangent function include the real and imaginary parts, and the output phase angle range is -π to π.

[0065] When constructing a spectrum energy distribution matrix containing frequency points and energy values ​​for harmonic components, the frequency points are discrete frequency points within the preset frequency band range of the harmonic components in step S1, with an interval of 100 Hz between discrete frequency points. The energy values ​​are the sum of squares of the integral amplitude of the corresponding frequency points within a preset time window. The sum of squares of the integral amplitude is calculated by performing trapezoidal numerical integration on the squares of the amplitudes of all sampling points within the time window. The step size of the trapezoidal numerical integration is the reciprocal of the sampling interval. The rows of the spectrum energy distribution matrix correspond to the frequency points and are arranged in ascending order of frequency, while the columns correspond to the time window and are aligned with the unified timestamp in step S1. The storage format of the spectrum energy distribution matrix is ​​a two-dimensional array of double-precision floating-point type. When initializing cluster centers based on the gradient distribution characteristics of frequency energy values ​​in the spectral energy distribution matrix, the gradient distribution characteristics are calculated as the sum of the absolute values ​​of the energy differences between each frequency point and its neighboring frequency points. The definition of neighboring frequency points includes left neighbor frequency points and right neighbor frequency points. The left neighbor frequency point is the frequency point after subtracting the frequency interval from the current frequency point, and the right neighbor frequency point is the frequency point after adding the frequency interval to the current frequency point. The number of initial cluster centers is dynamically determined based on the number of peak values ​​of the gradient distribution characteristics. The peak value of the gradient distribution characteristics is the frequency point where the gradient value is greater than twice the average gradient value. The average gradient value is calculated as the arithmetic mean of the gradient values ​​of all frequency points. The position of the initial cluster center is the frequency point position corresponding to the gradient peak value, and the calculation accuracy of the frequency point position is 1% of the frequency interval.

[0066] When iteratively optimizing the cost function by minimizing the distance between frequency point energy and cluster centers, the cost function is defined as the weighted sum of the squared Euclidean distances between all frequency point energy values ​​and their respective cluster centers. The weights are set as the proportion of the frequency point energy value to the total energy of its respective cluster. The proportion of the frequency point energy value to the total energy of its respective cluster is calculated by dividing the current frequency point energy value by the sum of all frequency point energy values ​​within the cluster. The termination condition for iterative optimization is that the change in the position of the cluster center is less than 10% of the frequency point interval or the number of iterations reaches the preset upper limit. The upper limit for the number of iterations is set to 100. After each iteration, the position of the cluster center is updated to the weighted average of the frequency point energy values ​​within its respective cluster. The weight of the weighted average is the normalized result of the frequency point energy value. The normalized result is calculated by dividing the frequency point energy value by the total energy value of its respective cluster.

[0067] When generating the harmonic energy distribution characteristics that characterize the harmonic energy accumulation region, the parameters of the harmonic energy distribution characteristics include the frequency point position of each cluster center, the maximum energy value of the cluster to which it belongs, and the standard deviation of the energy distribution. The calculation accuracy of the frequency point position is 1% of the frequency interval. The maximum energy value is the maximum value of the frequency point energy value within the cluster. The standard deviation of the energy distribution is the square root of the average of the squared differences between the frequency point energy values ​​and the mean within the cluster. The mean is calculated as the arithmetic mean of the energy values ​​of all frequency points within the cluster. The average of the squared differences is calculated as the sum of the squares of the differences between the frequency point energy values ​​and the mean divided by the number of frequency points. The square root operation is implemented using Newton's iteration method. The initial guess value of Newton's iteration method is half of the average of the squared differences. The iteration termination condition is that the difference between two adjacent iteration results is less than 1e-6.

[0068] The 10-millisecond window length is set to cover an integer multiple of the fundamental frequency period to reduce spectral leakage. With a fundamental frequency of 50 Hz, the fundamental period is 20 milliseconds, and the 10-millisecond window length is 0.5 times the fundamental period. The Hanning window improves time-frequency resolution by suppressing truncation effects, and a 50% overlap rate ensures data continuity between adjacent windows. The ratio of the difference in current amplitude change rate to the time interval is expressed in amperes per second, and the phase offset is measured in radians. The four-quadrant arctangent function avoids 180-degree ambiguity in phase calculations.

[0069] The 100 Hz frequency interval setting satisfies the Nyquist sampling theorem's requirement for the highest harmonic frequency. When the highest harmonic frequency is 5 kHz, the Nyquist frequency must be greater than 10 kHz. The actual sampling rate is 10 kHz, which meets the requirement. The frequency resolution corresponding to the frequency interval is the smallest resolvable unit of spectral energy. The gradient peak determination threshold (twice the average gradient value) is determined through historical data statistical analysis. In historical data, the average gradient under normal operating conditions is approximately 500 ampere-seconds, so the threshold is set to 1000 ampere-seconds to ensure that the initial position of the cluster center is located in the energy abrupt change region. The weight normalization logic in the weighted sum calculation of the Euclidean distance squares is to divide the frequency energy value by the total energy value of its respective cluster, ensuring that the sum of the weights is 1.

[0070] The threshold for position change in the iteration termination condition is 10% of the frequency interval, i.e., 1 Hz, used to determine whether the cluster center positions have converged. In the calculation of the standard deviation of energy distribution, the square root operation ensures that the unit of standard deviation is consistent with the unit of energy value; the unit of standard deviation is the square root of ampere-seconds.

[0071] The construction of the short-time Fourier transform and the spectral energy distribution matrix requires a processor that supports floating-point operations. The processor's floating-point unit must support the IEEE 754 double-precision standard. The memory capacity must meet the real-time storage requirements of the time-frequency matrix and the spectral energy distribution matrix. The storage format of the time-frequency matrix is ​​a two-dimensional array of double-precision floating-point type. Each matrix element occupies 8 bytes of memory. The total memory usage of the time-frequency matrix is ​​the number of frequency components multiplied by the number of time windows and then multiplied by 8 bytes.

[0072] The iterative optimization process of the clustering algorithm is executed in a hardware environment that supports multi-threaded computing. The parallel computing framework chosen is OpenMP to accelerate matrix operations. The number of threads in OpenMP is set to twice the number of physical cores. The update of the cluster center position and the calculation of the cost function are optimized through a vectorized instruction set, which is AVX2, and the data alignment is 64-byte aligned.

[0073] If the spectral energy distribution matrix is ​​empty, meaning that the harmonic components are completely filtered out in step S1, then the clustering step is skipped and the default energy distribution characteristics are output. The parameters of the default energy distribution characteristics include the frequency point position being the midpoint of the preset frequency band, the maximum energy value being zero, and the standard deviation of the energy distribution being zero.

[0074] If the cluster center iteration fails to converge (i.e., reaches the maximum number of iterations of 100 and the position change is still greater than 10% of the frequency interval), the result of the previous iteration is used as the final cluster center and marked as a low-confidence result. The criterion for determining a low-confidence result is that the standard deviation of the energy distribution of the cluster center exceeds a preset threshold. The preset threshold is set to 0.1 ampere-seconds based on the mean standard deviation of normal clusters in historical data. If the frequency point energy variance exceeds the hardware calculation limit (i.e., the variance value exceeds the maximum value of double-precision floating-point numbers, 1.7976931348623157e+308), floating-point scaling is enabled and the scaling factor is recorded. The scaling factor is calculated by dividing the frequency point energy value by 1e6 to avoid overflow. The scaled energy value is used in the calculation and restored in the final result.

[0075] S3. Based on the relationship between the phase offset and the preset direction threshold, determine whether the power flow direction is bidirectional. The specific implementation is as follows:

[0076] When dividing the phase offset into a forward offset sequence and a reverse offset sequence, the division of the forward offset sequence is based on the phase offset value being greater than zero. The value of the phase offset is the phase angle difference between adjacent sampling points extracted in step S2, and the unit of the phase angle difference is radians. The division of the reverse offset sequence is based on the phase offset value being less than zero. The positive or negative sign of the phase offset is determined according to the output result of the four-quadrant arctangent function in step S2. The input parameters of the four-quadrant arctangent function are the complex real part and imaginary part of the time-frequency matrix in step S2.

[0077] When calculating the cumulative gradient of the forward offset sequence and the cumulative gradient of the reverse offset sequence, the cumulative gradient is calculated by accumulating the differences of all phase offsets in the sequence. The cumulative gradient of the forward offset sequence is the algebraic sum of each phase offset in the forward offset sequence, and the cumulative gradient of the reverse offset sequence is the algebraic sum of each phase offset in the reverse offset sequence. The sign information is retained in the calculation of the algebraic sum to reflect the cumulative direction trend of the phase offset. The storage format of the forward offset sequence and the reverse offset sequence is a double-precision floating-point one-dimensional array. The array index is aligned with the unified timestamp in step S1, and the array index is arranged in ascending order of the timestamp.

[0078] When generating directional weight coefficients based on the ratio of the cumulative gradient of the forward offset sequence to the cumulative gradient of the backward offset sequence, the ratio is calculated as the ratio of the absolute value of the cumulative gradient of the forward offset sequence to the absolute value of the cumulative gradient of the backward offset sequence. The directional weight coefficients are generated by performing a natural logarithmic transformation on the ratio and then normalizing it to the interval between 0 and 1. The base of the natural logarithmic transformation is the natural constant e. The input for the natural logarithmic transformation is the ratio, and the output is a real number. The normalization is calculated by subtracting the minimum value after the transformation from the transformed value and then dividing by the difference between the maximum value and the minimum value after the transformation. The minimum value after the transformation is the minimum value in the result of the natural logarithmic transformation, and the maximum value after the transformation is the maximum value in the result of the natural logarithmic transformation. The normalized directional weight coefficients range from 0 to 1.

[0079] If the directional weight coefficient exceeds the first preset threshold and the cumulative gradients of the positive offset sequence and the negative offset sequence have opposite signs, then the power flow direction is determined to be bidirectional. The first preset threshold is set based on the statistical results of the directional weight coefficient distribution of bidirectional and unidirectional modes in historical data. The historical data includes at least 1,000 sets of distribution network operation samples. The statistical method is to calculate the mean of the directional weight coefficients of the bidirectional mode samples and the mean of the directional weight coefficients of the unidirectional mode samples, and take the midpoint of the two means as the first preset threshold. The value range of the first preset threshold is 0.6 to 0.8. The condition for opposite signs is that the cumulative gradient of the positive offset sequence has a positive sign and the cumulative gradient of the negative offset sequence has a negative sign, or the cumulative gradient of the positive offset sequence has a negative sign and the cumulative gradient of the negative offset sequence has a positive sign. The sign is obtained by judging whether the value of the cumulative gradient is greater than zero or less than zero. The value of the cumulative gradient is double-precision floating-point data.

[0080] The phase offset is measured in radians, the cumulative gradient is measured in radians, and the direction weighting coefficient is a dimensionless parameter. The natural logarithmic transform is calculated by taking the ratio of the absolute values ​​of the cumulative gradients as input and outputting a real number. The normalized direction weighting coefficient is strictly limited to the range of 0 to 1 to avoid judgment failure due to excessively large or small ratios. The first preset threshold of 0.6 to 0.8 is based on historical data analysis reports from a provincial power grid. The reports show that 95% of bidirectional mode samples have a direction weighting coefficient higher than 0.6, while 97% of unidirectional mode samples have a coefficient lower than 0.5. Setting the threshold at 0.6 covers more than 90% of real-world scenarios. In determining the sign of opposites, the opposite signs of the cumulative gradients of the forward and reverse offset sequences indicate a competitive alternation between the two power flow directions, physically corresponding to the power interaction behavior between distributed generation and load in the distribution network. A positive sign for the cumulative gradient of the forward offset sequence indicates a continuously strengthening forward power flow trend, while a negative sign indicates a weakening trend. The judgment logic for the reverse offset sequence is symmetrical.

[0081] The calculation of the cumulative gradient requires a processor that supports floating-point accumulation operations. The processor's floating-point accumulation instructions must meet the IEEE 754 double-precision standard to avoid rounding errors. The memory capacity must meet the independent storage requirements of the forward and reverse offset sequences. The storage space for each sequence is the number of phase offsets multiplied by 8 bytes. The number of phase offsets is determined by the preset time window length and sampling rate in step S1. The normalization calculation of the direction weight coefficients requires a mathematical coprocessor that supports natural logarithmic function operations. The mathematical coprocessor's calculation precision must support at least six significant digits after the decimal point, the calculation error of the natural logarithmic function must be less than 1e-6, and the clock frequency of the mathematical coprocessor must be no less than 1GHz to ensure real-time performance.

[0082] If all phase offsets are zero, meaning no phase offset was detected in step S2, the bidirectional mode determination is skipped and the default unidirectional mode result is output. The storage format of the default result is consistent with the format of the fault determination result in step S6. If the forward offset sequence or the reverse offset sequence is empty, meaning all phase offsets are positive or all are negative, the direction weight coefficient is set to zero and the system is determined to be in unidirectional mode. The confidence level of the determination result is marked as high confidence. If the input of the natural logarithmic transform is zero or negative, meaning the ratio is zero or negative due to the cumulative gradient being zero or having the same sign, the anomaly handling mechanism is triggered. The anomaly handling mechanism includes discarding the data in the current time window and using the determination result of the previous valid window. The time range of the previous valid window is the latest valid result within 1 second prior to the current time. If there is no valid result within 1 second, the system default value is output.

[0083] The generation logic of the directional weight coefficients was verified using measured data from a regional power grid. 1000 samples, including those from photovoltaic, energy storage, and load fluctuations, were selected to calculate the statistical distribution of the directional weight coefficients. The accuracy rate for bidirectional mode samples was 93.2%, and for unidirectional mode samples, it was 89.5%. The judgment error mainly stemmed from phase jitter during transient processes. When the jitter duration was less than 10 milliseconds, the judgment result could be corrected through harmonic interference analysis in subsequent steps. The accuracy rate for determining sign inversion was 97.8%. In cases of misjudgment, sign jumps caused by noise were resolved by adding a smoothing filter with a phase offset. The smoothing filter window length was 5 sampling points, and the window type was a moving average window. The hardware computing resource configuration was verified using a test environment with an Intel Xeon E5-2678 processor and an NVIDIA Tesla V100 math coprocessor. The single judgment time was less than 2 milliseconds, meeting the real-time requirements of the distribution network.

[0084] S4. In bidirectional mode, based on the harmonic energy distribution characteristics, time-varying characteristic analysis is performed on the high-frequency harmonic components to identify interference components that overlap with the fault characteristic frequency band. Specifically, the implementation is as follows:

[0085] If the power flow direction is determined to be bidirectional, the correlation coefficient between the energy of each frequency band of the high-frequency harmonic component and the fault characteristic frequency band is calculated. The frequency band energy of the high-frequency harmonic component is the frequency point energy value in the harmonic energy distribution characteristics generated in step S2. The fault characteristic frequency band is a preset frequency band range associated with typical faults in the distribution network. The preset fault characteristic frequency band range is set according to the spectrum analysis results of cable breakdown and insulation aging events in historical fault data. The characteristic frequency band of cable breakdown fault is 2.5 kHz to 3.5 kHz, and the characteristic frequency band of insulation aging fault is 4.0 kHz to 5.0 kHz.

[0086] The correlation coefficient is calculated as the ratio of the covariance to the standard deviation of the frequency band energy within a preset time window. The covariance is calculated as the covariance of the frequency band energy sequence and the fault characteristic frequency band energy sequence, with the dimension of the covariance being the square of amperes squared seconds. The standard deviation is the standard deviation of the frequency band energy sequence, with the dimension of the standard deviation being amperes squared seconds. The ratio is calculated by dividing the covariance by the product of the standard deviations, and the result of the ratio calculation is a dimensionless parameter. The fault characteristic frequency band energy sequence is obtained by summing the frequency point energy values ​​within the preset fault characteristic frequency band, with the frequency point interval of the summation being 100 Hz, consistent with the frequency point interval of the spectral energy distribution matrix in step S2.

[0087] When adjusting the time length of the dynamic sliding window based on the time-series change rate of the correlation coefficient, the time-series change rate is calculated as the ratio of the absolute value of the difference between the correlation coefficients within adjacent time windows to the time interval. The unit of the time interval is seconds, the absolute value of the difference is a dimensionless parameter, and the unit of the ratio is seconds. The initial length of the dynamic sliding window is the preset time window of the spectral clustering in step S2, with an initial value of 5 milliseconds. The adjustment logic of the time length of the dynamic sliding window is as follows: when the time-series change rate exceeds the preset change rate threshold, the window length is shortened to 50% of the initial length, with a shortening step of 10% of the initial length, i.e., a reduction of 0.5 milliseconds each time. If the time-series change rate does not exceed the preset change rate threshold, the window length is extended to 150% of the initial length, with an extension step of 10% of the initial length, i.e., an increase of 0.5 milliseconds each time. The preset change rate threshold is set based on the stability analysis results of high-frequency harmonic components in historical data, and the value range of the preset change rate threshold is 0.1 to 0.2 seconds per second.

[0088] When detecting energy abrupt changes in high-frequency harmonic components within a dynamic sliding window, the detection method involves calculating the absolute difference between the energy of each frequency band within the window. The absolute difference is calculated as the absolute value of the difference between the energy of the current frequency band and the energy of the same frequency band in the previous time window. If the absolute difference exceeds three times the standard deviation of the average energy within the window, it is determined to be an energy abrupt change. The average energy within the window is calculated as the arithmetic mean of the frequency band energy, with the dimension of the arithmetic mean being ampere-seconds. The standard deviation is calculated as the square root of the average of the squared differences between the frequency band energy and the mean. The square root operation is implemented using Newton's iteration method. The initial guess value of Newton's iteration method is half of the average of the squared differences. The iteration termination condition is that the difference between two adjacent iteration results is less than 1e-6.

[0089] If the overlap rate between the frequency band of the energy mutation point and the fault characteristic frequency band exceeds the second preset threshold, it is determined to be an interference component. The overlap rate is calculated by dividing the number of intersection frequency points between the frequency band of the energy mutation point and the fault characteristic frequency band by the total number of frequency points in the fault characteristic frequency band. The number of intersection frequency points is the number of the same frequency points within the two frequency bands. The total number of frequency points is the difference between the maximum and minimum frequency point numbers in the fault characteristic frequency band divided by the frequency point interval. The second preset threshold is set based on the statistical results of the overlap rate between real fault events and interference events in historical data. The value range of the second preset threshold is 60% to 80%. The frequency band energy data determined to be an interference component will be marked as a high-risk area and transmitted to subsequent steps.

[0090] The calculation of covariance requires the two energy sequences to have the same length. If the lengths are inconsistent, alignment is achieved through interpolation or truncation. Linear interpolation is used for interpolation, and the logic for truncation alignment is to take the minimum length of the two sequences. In the calculation of the average of the squared differences of the standard deviation, the mean has the same dimension as the energy of the frequency band, and the square root operation of the average of the squared differences ensures that the standard deviation matches the dimension of the energy value. The adjustment of the time length of the dynamic sliding window requires real-time monitoring of the time series change rate. The calculation cycle of the time series change rate is synchronized with the sampling rate of data acquisition in step S1, and the window length is updated at the end of each sampling cycle. The threshold of the absolute difference of energy mutation points is based on the Rheinda criterion (3σ principle), covering 99.7% of the normal fluctuation range to avoid misjudgment due to random noise. In the calculation of the overlap rate, the frequency point interval is fixed at 100 Hz, consistent with the frequency point interval of the spectral energy distribution matrix in step S2, to ensure the alignment of frequency point numbers.

[0091] Covariance and standard deviation calculations require a processor that supports double-precision floating-point operations. The processor's floating-point unit must meet the IEEE 754 standard, and the memory capacity must meet the real-time storage requirements of the frequency band energy sequence. The frequency band energy sequence is stored in a double-precision floating-point one-dimensional array. Dynamic sliding window adjustment requires real-time clock interrupt support with a clock interrupt precision of 1 microsecond. The window length adjustment step is stored in non-volatile memory to ensure that the configuration is not lost after power failure. The square root operation in energy mutation point detection requires the support of a mathematical coprocessor. The mathematical coprocessor's calculation precision must support at least six significant digits after the decimal point, and the maximum number of iterations for the square root operation is set to 100.

[0092] If a division-by-zero error occurs during covariance calculation, i.e., the standard deviation is zero, the current window is skipped and the correlation coefficient result of the previous valid window is used. The time range of the previous valid window is the latest valid result within 1 second prior to the current time. If the length of the dynamic sliding window exceeds the hardware memory limit, for example, extending it to 7.5 milliseconds requires 12MB of memory while the hardware only supports 10MB, the adjustment is terminated and a memory overflow alarm is triggered. The alarm signal is output through the status register. If the energy in the frequency band of the energy mutation point detection is all zero, i.e., no effective harmonic components were detected in step S2, it is determined that there are no interference components and an empty result is output. The storage format of the empty result is consistent with the format of the fault determination result in step S6. If the total number of frequency points in the fault characteristic frequency band in the overlap rate calculation is zero, i.e., the preset frequency band range is incorrect, resulting in no effective frequency points, the determination is skipped and a configuration error log is recorded. The storage path of the configuration error log is the system log directory.

[0093] In the hardware verification environment, the collaborative computing between the Intel Xeon E5-2678 processor and the NVIDIA Tesla V100 math coprocessor enabled a single window adjustment to take less than 0.5 milliseconds, meeting real-time requirements.

[0094] S5. Based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude, an adaptive detection threshold function is constructed, specifically implemented as follows:

[0095] When performing recursive quantization analysis on harmonic components in adjacent frequency bands based on the frequency band information of the interference component, the definition of adjacent frequency bands includes the frequency band where the current interference component is located and the frequency bands on its left and right sides, with the frequency band interval being 100 Hz, the same as the frequency point interval of the spectral energy distribution matrix in step S2. The setting of the frequency band interval is consistent with the bandpass filtering range of the harmonic components in step S1. The calculation method for recursive quantization analysis is to construct the harmonic energy time series of adjacent frequency bands. The harmonic energy time series is obtained by extracting the energy values ​​of the corresponding frequency points from the harmonic energy distribution characteristics generated in step S2 and arranging them according to the time window. The length of the time window is consistent with the current length of the dynamic sliding window in step S4.

[0096] The calculation of recursion graph features includes deterministic and laminar flow indices. The deterministic index is calculated as the proportion of diagonal structures in the recursion graph to the total number of recursive points. A diagonal structure is defined as a pair of two consecutive recursive points whose time interval does not exceed a preset step size, which is set to 10% of the time window length. The laminar flow indices are calculated as the proportion of vertical line segments in the recursion graph to the total number of recursive points. A vertical line segment structure is defined as three or more consecutive recursive points appearing vertically. When generating the coupling matrix characterizing the nonlinear coupling strength between frequency bands, the rows and columns of the coupling matrix correspond to the numbers of adjacent frequency bands. The matrix element values ​​are the weighted sum of the deterministic and laminar flow indices for the corresponding frequency band pairs. The weighting coefficients of the weighted sum are dynamically adjusted according to the frequency band spacing. For each increase in the frequency band spacing, the weighting coefficient decreases by 10%, with an initial value of 0.7 and a lower limit of 0.4 after reduction.

[0097] Based on the gradient distribution characteristics of the rate of change of current amplitude, when calculating the dynamic correlation weight coefficient, the rate of change of current amplitude is the rate of change of current amplitude within the preset time window extracted in step S2. The gradient distribution characteristics are calculated as the sum of the absolute values ​​of the first-order differences of the current amplitude change rate sequence. The unit of the sum of the absolute values ​​of the first-order differences is amperes per second. The dynamic correlation weight coefficient is generated by normalizing the gradient distribution characteristics. The normalization is calculated by subtracting the minimum gradient value from the current gradient value and then dividing by the difference between the maximum and minimum gradient values. The statistical time window for the maximum and minimum gradient values ​​is the current length of the dynamic sliding window in step S4. The time length of the dynamic sliding window is adjusted from 2.5 milliseconds to 7.5 milliseconds.

[0098] When performing tensor multiplication on the coupling matrix and the dynamically associated weight coefficients, the tensor multiplication operation is performed by multiplying each element of the coupling matrix with the gradient value at the corresponding position of the dynamically associated weight coefficient to generate the core parameters of the adaptive detection threshold function. The core parameters are stored in a two-dimensional matrix, with rows corresponding to frequency band numbers and columns corresponding to time window numbers. The units of the matrix elements are amperes per second.

[0099] When adjusting the sensitivity coefficient of the threshold function based on the difference between the core parameters and the preset benchmark threshold, the preset benchmark threshold is set based on the statistical average of the core parameters under normal and fault conditions in historical data. The historical data includes 500 fault events and 2000 normal fluctuation data recorded in the distribution network over three years. The average core parameter under normal conditions is 0.4, and the average core parameter under fault conditions is 0.7. The adjustment logic of the sensitivity coefficient is as follows: when the core parameter exceeds the preset benchmark threshold, the sensitivity coefficient increases linearly according to the excess ratio, with an increase of 0.5% for every 1% excess. When the core parameter is lower than the preset benchmark threshold, the sensitivity coefficient decreases linearly according to the deficiency ratio, with a decrease of 0.3% for every 1% deficiency. The initial value of the sensitivity coefficient is 1.0, the adjusted upper limit is 2.0, and the lower limit is 0.5.

[0100] In recursive quantitative analysis, the deterministic and laminar flow indices are dimensionless, the element values ​​of the coupling matrix range from 0 to 1, and the precision of the matrix elements is three decimal places. The gradient distribution characteristics of the current amplitude change rate are measured in amperes per second, and the normalized dynamic correlation weight coefficients range from 0 to 1. The update frequency of the maximum and minimum gradient values ​​in the normalization calculation is triggered at the end of each dynamic sliding window. The element values ​​of the resulting matrix from the tensor product operation are measured in amperes per second. The difference between the core parameters and the preset benchmark thresholds is consistent in dimension. The adjustment ratio of the sensitivity coefficient is set based on the balance between false alarm rate and false negative rate in historical data. The statistical results of the false alarm rate and false negative rate are as follows: with a sensitivity coefficient of 1.0, the false alarm rate is 15% and the false negative rate is 20%; with a sensitivity coefficient of 2.0, the false alarm rate is 6% and the false negative rate is 9%.

[0101] Recursive quantization analysis requires a processor that supports matrix operations. The processor's cache capacity must meet the real-time storage requirements of the recursive graph. The recursive graph is stored as a Boolean two-dimensional array, where rows and columns correspond to the indices of the time series. A value of 1 indicates the existence of a recursive point, and 0 indicates no recursive point. Normalization calculation of dynamically correlated weight coefficients requires floating-point division units. The division unit's calculation precision must meet the IEEE 754 double-precision standard, and the rounding mode for division operations is nearest-neighbor rounding. Tensor multiplication operations require GPU acceleration support. The GPU must have at least 1024 parallel computing cores, and its memory capacity must be twice the size of the core parameter matrix. The size of the core parameter matrix is ​​the number of frequency bands multiplied by the number of time windows. The maximum number of frequency bands is the total number of frequency points within the preset frequency band range in step S2.

[0102] If the frequency band energy sequence in the recursive quantization analysis is all zero, meaning no effective harmonic components were detected in step S2, then the coupling matrix calculation is skipped and a default coupling strength value of 0.5 is used. The default coupling strength value is stored in non-volatile memory. If the maximum and minimum gradient values ​​of the dynamic correlation weight coefficients are equal, meaning the rate of change of current amplitude is stable, then the weight coefficients are uniformly set to 0.5 and marked as steady-state mode. If the core parameter matrix size does not match the dimension of the dynamic correlation weight coefficients, a dimension alignment anomaly is triggered. The anomaly handling mechanism includes truncating data of longer dimensions or padding with zero values ​​to match the dimension. The truncation logic is to retain the latest data in the time window, and the padded zero values ​​do not participate in subsequent calculations. If the sensitivity coefficient exceeds the upper or lower limit after adjustment, it is forcibly limited to the upper or lower limit value and an alarm log is triggered. The alarm log records the timestamp, the value exceeding the limit, and the adjusted sensitivity coefficient.

[0103] In the hardware verification environment, the NVIDIA Tesla V100 GPU acceleration reduced the time for tensor product operations from 12 milliseconds to 0.8 milliseconds, and the floating-point unit of the Intel Xeon E5-2678 processor ensured the real-time performance of recursive quantization analysis, with a single recursive graph construction taking less than 1 millisecond.

[0104] S6. Input the rate of change of current amplitude into the adaptive detection threshold function, output the fault determination result, and generate distribution network reliability assessment parameters based on the fault determination result. The specific implementation is as follows:

[0105] When the rate of change of current amplitude is input into the adaptive detection threshold function, the rate of change of current amplitude is the rate of change of current amplitude within the preset time window extracted in step S2. The core parameters and sensitivity coefficients of the adaptive detection threshold function are the two-dimensional matrix generated in step S5 and the adjusted coefficient values. The rows of the two-dimensional matrix correspond to the frequency band number, and the columns correspond to the time window number. The adjustment range of the sensitivity coefficient is limited to 0.5 to 2.0.

[0106] When performing multi-scale fluctuation decomposition on the rate of change of current amplitude, the multi-scale fluctuation decomposition method is to use discrete wavelet transform, with Daubechies 4 wavelet as the wavelet basis function and 3 decomposition levels. The high-frequency fluctuation component is extracted by retaining the wavelet coefficients of the first and second levels to reconstruct the signal, and the low-frequency trend component is extracted by retaining the wavelet coefficients of the third level to reconstruct the signal. The frequency range of the high-frequency fluctuation component is from the upper limit of the preset frequency band of the harmonic component in step S1 to the Nyquist frequency, and the frequency range of the low-frequency trend component is from 1% to 10% of the fundamental frequency. The value of the fundamental frequency is the power frequency of 50 Hz of the fundamental component in step S1.

[0107] Based on the core parameters and sensitivity coefficient of the adaptive detection threshold function, when calculating the dynamic confidence threshold of the high-frequency fluctuation component and the steady-state deviation threshold of the low-frequency trend component, the dynamic confidence threshold is calculated by multiplying the element value of the corresponding frequency band and time window in the core parameter matrix by the sensitivity coefficient, and then superimposing the product result on the preset benchmark threshold. The preset benchmark threshold is set based on the statistical mean of the amplitude of the high-frequency fluctuation component under normal operating conditions in historical data. The historical data includes 1,000 normal fluctuation samples recorded in the distribution network within two years, with a statistical mean of 0.3 amperes per second. The steady-state deviation threshold is calculated by multiplying the absolute value of the slope of the low-frequency trend component by the sensitivity coefficient. The method for calculating the absolute value of the slope is to perform linear regression fitting on the low-frequency trend component. The time window length of the linear regression fitting is the current length of the dynamic sliding window in step S4. The absolute value of the slope of the fitted line is used as the quantitative indicator of the steady-state deviation trend, and the dimension of the absolute value of the slope is amperes per second squared.

[0108] If the amplitude of the high-frequency fluctuation component exceeds the dynamic confidence threshold and the slope of the low-frequency trend component exceeds the steady-state deviation threshold, then it is determined to be a persistent fault. The determination condition for a persistent fault is that the amplitude of the high-frequency fluctuation component exceeds the dynamic confidence threshold for three consecutive time windows, and the slope of the low-frequency trend component exceeds the steady-state deviation threshold for three consecutive time windows. The determination logic for consecutive time windows is that the timestamps are continuous and there is no data interruption. The time interval of the time window is 1 millisecond, which is the sampling time interval of the fundamental component in step S1.

[0109] Based on the determination of persistent faults, and combined with the repair time and impact range parameters of similar faults in the historical fault database, when generating distribution network reliability assessment parameters including the average repair time index and fault impact factor, the matching method for similar faults is to calculate the similarity between the frequency band distribution characteristics of persistent faults and the fault feature labels in the historical fault database. The similarity calculation method is to match the minimum value of the sum of squares of Euclidean distances. The sum of squares of Euclidean distances is calculated as the sum of squares of the energy differences between the corresponding frequency points of the current fault frequency band distribution characteristics and the historical fault feature labels. The repair time is obtained by the weighted average of the repair times in the matched historical fault records. The weight of the weighted average is the reciprocal of the similarity of the historical fault records. The impact range parameter is obtained by the arithmetic mean of the number of affected nodes in the historical fault records. The calculation time window of the arithmetic mean is the data of the most recent year of the historical fault records. The average repair time index is calculated by multiplying the repair time by the ratio of the number of currently affected nodes to the historical average number of affected nodes. The fault impact factor is calculated by the ratio of the number of affected nodes to the total number of nodes in the distribution network. The total number of nodes in the distribution network is the total number of nodes collected in step S1.

[0110] The discrete wavelet transform decomposition layer of 3 corresponds to a time resolution of 1 / 8 to 1 / 2 of the fundamental period, ensuring that the high-frequency fluctuation component covers the characteristic frequency band of transient faults. The amplitude of the high-frequency fluctuation component is measured in amperes per second, and the slope of the low-frequency trend component is measured in amperes per second squared. In the product calculation of the dynamic confidence threshold, the core parameter matrix element values ​​are measured in amperes per second, the sensitivity coefficient is dimensionless, and the dimension of the product result is consistent with the preset benchmark threshold. In the linear regression fitting of the low-frequency trend component, the time window length is the current length of the dynamic sliding window in step S4, the absolute value of the fitting slope is measured in amperes per second squared, and the sum of squares of the fitting residuals is used as the criterion for goodness of fit. When the sum of squares of the residuals exceeds the preset threshold, refitting is triggered. The logic for determining the continuous fault through three consecutive time windows can effectively filter transient interference. The continuity of the time windows is verified through a unified timestamp in step S1, with the format being year, month, day, hour, minute, second, and millisecond.

[0111] Discrete wavelet transform requires a processor that supports fast convolution operations. The processor's instruction set must include SIMD instructions to accelerate wavelet coefficient calculation. The memory capacity must meet the real-time storage requirements of three layers of wavelet coefficients. Each layer of coefficients is stored as a double-precision floating-point array, with the array index aligned to the unified timestamp in step S1. Linear regression fitting requires floating-point multiply-accumulate units. The computation latency of the multiply-accumulate units must be less than 1 microsecond, and intermediate results of floating-point multiply-accumulate operations are cached in the L1 data cache. Querying and matching in the historical fault database requires SSD storage. The database index structure is a B+ tree, the query response time is less than 10 milliseconds, and the database records are stored in JSON format, including fields for fault feature labels, repair duration, and the number of affected nodes.

[0112] If the reconstructed signals of the high-frequency fluctuation components or low-frequency trend components after multi-scale fluctuation decomposition are all zero, i.e., no effective current amplitude change rate is detected in step S2, then it is determined to be a data anomaly and a resampling mechanism is triggered. The data source for resampling is the historical data cached in step S1, and the caching time for the historical data is the most recent 30 days. If the calculated result of the dynamic confidence threshold or steady-state deviation threshold exceeds the range of hardware floating-point representation, i.e., exceeds the maximum value of IEEE 754 double-precision floating-point number 1.7976931348623157e+308, then fixed-point scaling is enabled and the scaling factor is recorded. The scaling factor is stored in non-volatile memory, and the scaling factor is calculated by dividing the threshold by 1e6 to avoid overflow. If there is no matching type in the historical fault database, i.e., the sum of squared Euclidean distances of all historical fault records exceeds the preset similarity threshold, then the default repair time and impact range parameters are adopted. The default repair time is the maximum allowable repair time of 4 hours in the distribution network operation and maintenance specifications, and the default impact range parameter is 10% of the total number of nodes in the distribution network.

[0113] In the hardware verification environment, the collaborative work of the Intel Xeon E5-2678 processor and NVMe SSD kept the database query time stable within 8 milliseconds, and the single decomposition of discrete wavelet transform took 0.2 milliseconds, meeting the real-time requirements.

[0114] Example 2: Figure 2 A schematic diagram of the structure of a power distribution network reliability intelligent analysis system according to the present invention is provided. The power distribution network reliability intelligent analysis system includes the following modules:

[0115] The current monitoring module is used to acquire real-time current data of each node in the distribution network. The real-time current data includes the fundamental component and harmonic components.

[0116] The harmonic clustering module is used to perform joint time-frequency domain analysis on the fundamental component, extract the current amplitude change rate and phase shift, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics.

[0117] The direction determination module is used to determine whether the power flow direction is bidirectional based on the relationship between the phase offset and the preset direction threshold.

[0118] The interference identification module is used, in bidirectional mode, to perform time-varying characteristic analysis on high-frequency harmonic components based on harmonic energy distribution characteristics, and to identify interference components in high-frequency harmonic components that overlap with the fault characteristic frequency band.

[0119] The threshold construction module is used to construct an adaptive detection threshold function based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude.

[0120] The reliability assessment module is used to input the rate of change of current amplitude into the adaptive detection threshold function, output the fault determination result, and generate distribution network reliability assessment parameters based on the fault determination result.

[0121] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0122] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0123] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and inventive constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0124] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0125] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0126] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0127] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for intelligent analysis of distribution network reliability, characterized in that, Includes the following steps: S1. Obtain real-time current data of each node in the distribution network. The real-time current data includes the fundamental component and harmonic components. S2. Perform time-frequency domain joint analysis on the fundamental component, extract the current amplitude change rate and phase offset, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics. S3. Determine whether the power flow direction is bidirectional based on the relationship between the phase offset and the preset direction threshold. S4. If it is a two-way mode, perform time-varying characteristic analysis on the high-frequency harmonic components based on the harmonic energy distribution characteristics, and identify the interference components in the high-frequency harmonic components that overlap with the fault characteristic frequency band. S5. Based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude, an adaptive detection threshold function is constructed. S6. Input the rate of change of current amplitude into the adaptive detection threshold function, output the fault judgment result, and generate distribution network reliability assessment parameters based on the fault judgment result.

2. The intelligent analysis method for distribution network reliability according to claim 1, characterized in that, Acquire real-time current data at each node of the distribution network. The real-time current data includes the fundamental component and harmonic components, including: The three-phase current signals of each node in the distribution network are collected synchronously, a unified timestamp is added to the collection time, and the three-phase current signals are bandpass filtered to separate the fundamental component from the harmonic components within the preset frequency band.

3. The intelligent analysis method for distribution network reliability according to claim 1, characterized in that, A joint time-frequency domain analysis was performed on the fundamental component to extract the rate of change of current amplitude and phase shift. Spectral clustering was then performed on the harmonic components to generate harmonic energy distribution characteristics, including: A short-time Fourier transform is performed on the fundamental component to generate a time-frequency matrix. The rate of change of current amplitude within a preset time window and the phase offset between adjacent sampling points are extracted from the time-frequency matrix. A spectral energy distribution matrix containing frequency points and energy values ​​is constructed for the harmonic components. Cluster centers are initialized based on the gradient distribution characteristics of the frequency point energy values ​​in the spectral energy distribution matrix. The harmonic energy distribution characteristics representing the harmonic energy accumulation region are generated by iterative optimization of the cost function that minimizes the distance between the frequency point energy and the cluster center.

4. The intelligent analysis method for distribution network reliability according to claim 1, characterized in that, Based on the relationship between the phase offset and a preset direction threshold, determine whether the power flow direction is bidirectional, including: The phase offset is divided into a forward offset sequence and a reverse offset sequence, and the cumulative gradient of the forward offset sequence and the cumulative gradient of the reverse offset sequence are calculated respectively. Directional weighting coefficients are generated based on the ratio of the cumulative gradient of the forward offset sequence to the cumulative gradient of the backward offset sequence. If the directional weight coefficient exceeds the first preset threshold and the cumulative gradients of the positive offset sequence and the negative offset sequence have opposite signs, then the power flow direction is determined to be bidirectional.

5. The intelligent analysis method for distribution network reliability according to claim 1, characterized in that, In bidirectional mode, time-varying characteristic analysis of high-frequency harmonic components is performed based on harmonic energy distribution characteristics to identify interference components in the high-frequency harmonic components that overlap with the fault characteristic frequency band, including: If the power flow direction is determined to be bidirectional, calculate the correlation coefficient between the energy of each frequency band of the high-frequency harmonic component and the fault characteristic frequency band. The time length of the dynamic sliding window is adjusted based on the time-series change rate of the correlation coefficient, and the initial length of the dynamic sliding window is the preset time window for spectral clustering in step S2. Within a dynamic sliding window, energy mutation points of high-frequency harmonic components are detected. If the overlap rate between the frequency band of the energy mutation point and the fault characteristic frequency band exceeds a second preset threshold, it is determined to be an interference component.

6. The intelligent analysis method for distribution network reliability according to claim 5, characterized in that, The correlation coefficient is calculated as the ratio of the covariance to the standard deviation of the frequency band energy within a preset time window.

7. The intelligent analysis method for distribution network reliability according to claim 1, characterized in that, Based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude, an adaptive detection threshold function is constructed, including: Based on the frequency band information of the interference components, recursive quantization analysis is performed on the harmonic components of adjacent frequency bands to generate a coupling matrix characterizing the nonlinear coupling strength between frequency bands. Based on the gradient distribution characteristics of the rate of change of current amplitude, the dynamic correlation weight coefficient is calculated; The core parameters of the adaptive detection threshold function are generated by performing tensor multiplication on the coupling matrix and the dynamically associated weight coefficients. The sensitivity coefficient of the threshold function is adjusted based on the difference between the core parameters and the preset benchmark threshold.

8. The intelligent analysis method for distribution network reliability according to claim 1, characterized in that, The rate of change of current amplitude is input into an adaptive detection threshold function, which outputs a fault determination result. Based on the fault determination result, distribution network reliability assessment parameters are generated, including: Multi-scale fluctuation decomposition is performed on the rate of change of current amplitude to extract high-frequency fluctuation components and low-frequency trend components. Based on the core parameters and sensitivity coefficient of the adaptive detection threshold function, the dynamic confidence threshold of the high-frequency fluctuation component and the steady-state deviation threshold of the low-frequency trend component are calculated respectively. If the amplitude of the high-frequency fluctuation component exceeds the dynamic confidence threshold and the slope of the low-frequency trend component exceeds the steady-state deviation threshold, it is determined to be a persistent fault. Based on the determination of persistent faults, and combined with the repair time and impact range parameters of similar faults in the historical fault database, distribution network reliability assessment parameters are generated.

9. The intelligent analysis method for distribution network reliability according to claim 8, characterized in that, The reliability assessment parameters for power distribution networks include the mean time to repair index and the fault impact factor.

10. A distribution network reliability intelligent analysis system, used to implement the distribution network reliability intelligent analysis method according to any one of claims 1-9, characterized in that, Includes the following modules: The current monitoring module is used to acquire real-time current data of each node in the distribution network. The real-time current data includes the fundamental component and harmonic components. The harmonic clustering module is used to perform joint time-frequency domain analysis on the fundamental component, extract the current amplitude change rate and phase shift, and perform spectral clustering on the harmonic components to generate harmonic energy distribution characteristics. The direction determination module is used to determine whether the power flow direction is bidirectional based on the relationship between the phase offset and the preset direction threshold. The interference identification module is used, in bidirectional mode, to perform time-varying characteristic analysis on high-frequency harmonic components based on harmonic energy distribution characteristics, and to identify interference components in high-frequency harmonic components that overlap with the fault characteristic frequency band. The threshold construction module is used to construct an adaptive detection threshold function based on the dynamic correlation between the nonlinear coupling effect of harmonic components in adjacent frequency bands and the rate of change of current amplitude. The reliability assessment module is used to input the rate of change of current amplitude into the adaptive detection threshold function, output the fault determination result, and generate distribution network reliability assessment parameters based on the fault determination result.

Images