A method for identifying current characteristics in a transformer substation based on correlation peak judgment method

Abstract

The present invention discloses a method for identifying the current characteristics of a transformer substation area based on the correlation peak decision method, including: collecting the synchronous phase information of the power grid, generating and injecting a unique phase-coded characteristic signal; collecting the current signal at the key nodes of the transformer substation area and extracting the phase information; constructing a dynamic adaptive phase relationship network to extract the feature vectors with topological invariance; performing correlation analysis on the extracted feature vectors and the standard template, and determining the identification result through adaptive correlation peak decision; dynamically adjusting the system parameters according to the identification result. The method of the present invention combines phase coding with topological invariant features, and has the characteristics of strong anti-interference ability, good environmental adaptability and high recognition accuracy.

Description

Technical Field

[0001] The present invention belongs to the field of power information acquisition, and in particular, to a method for identifying the current characteristics of a distribution transformer area based on the correlation peak decision method. Background Art

[0002] Identifying the current characteristics of a distribution transformer area is a key technology for the safe operation and energy management of smart grids, and has important research significance. With the large-scale access of distributed energy and new loads to the power grid, traditional power grid monitoring and management methods are facing huge challenges. Accurately identifying the current characteristics of a distribution transformer area can achieve power grid fault warning, load identification, unauthorized power consumption detection, and power quality assessment, which has important value for improving the safety, reliability, and economy of the power grid, and is a basic support technology for realizing smart grids and the energy Internet.

[0003] Currently, the identification of the current characteristics of a distribution transformer area mainly relies on frequency analysis and time-domain feature extraction methods. Frequency analysis methods usually use technologies such as Fourier transform and wavelet analysis, and achieve identification by injecting characteristic signals in a specific frequency band and performing spectrum analysis at the receiving end. Time-domain feature extraction methods are based on statistical characteristics, and construct feature vectors by extracting statistical quantities such as the peak value, mean value, and variance of the current waveform. Some studies also use machine learning-based methods to train models through a large amount of historical data to achieve automatic feature recognition. These methods perform well in ideal experimental environments, but still have limitations in actual complex power grid environments.

[0004] However, the existing identification methods have obvious deficiencies in complex power grid environments. First, the frequency feature-based methods are vulnerable to harmonic interference and frequency overlap, and the identification accuracy decreases significantly in areas with dense nonlinear loads; second, the existing methods are highly sensitive to environmental changes, and small changes in line parameters caused by seasonal temperature and humidity changes will cause the preset model to fail; third, the existing methods generally use absolute measurement values as criteria and lack the excavation of topological relationships. When the data of some monitoring points are lost or distorted, the overall identification ability drops sharply. These problems seriously restrict the popularization and application of the technology for identifying the current characteristics of a distribution transformer area in actual power grids, and there is an urgent need to develop more robust and adaptive identification methods. Summary of the Invention

[0005] The object of the invention is to provide a method for identifying the current characteristics of a distribution transformer area based on the correlation peak decision method, in order to solve at least one technical problem existing in the prior art.

[0006] Technical solution: A method for identifying the current characteristics of a distribution transformer area based on the correlation peak decision method includes:

[0007] Collect the synchronous phase information of the power grid, generate a characteristic signal with phase encoding and inject it into the power line to form an injected current signal;

[0008] Collect line current signals containing injected current signals at a predetermined key node in the substation area, preprocess them and extract phase information to obtain a phase time series; accordingly, extract the current characteristic phase difference pattern and generate a topological invariant feature vector.

[0009] Perform correlation analysis on the topological invariant feature vector and a preset standard feature template, calculate the correlation peak value and establish a decision criterion to determine the identification result and the identification credibility index; accordingly, dynamically adjust the system parameters, update the phase network structure and the decision threshold.

[0010] Advantageous effects: The present invention can effectively avoid harmonic interference and frequency overlap problems, and has strong resistance to local data loss or distortion. It establishes a complete adaptive optimization and parameter update mechanism, enabling the system to continuously learn and adapt to environmental changes, and has the characteristics of strong anti-interference ability, good environmental adaptability, and high recognition accuracy. Description of the Drawings

[0011] Figure 1 It is a step flow chart of a method for identifying the current characteristics of a substation area based on the correlation peak decision method provided by an embodiment of the present application.

[0012] Figure 2 It is a step flow chart of the steps for forming an injected current signal provided by an embodiment of the present application.

[0013] Figure 3 It is a step flow chart of the steps for obtaining a phase time series provided by an embodiment of the present application.

[0014] Figure 4 It is a step flow chart of the steps for generating a topological invariant feature vector provided by an embodiment of the present application.

[0015] Figure 5 It is a step flow chart of the steps for determining the identification result and the identification credibility index provided by an embodiment of the present application. Detailed Embodiments

[0016] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0017] It should be specifically noted that for clearly showing the step flow of this application, serial numbers are marked for each step in the specification. These serial numbers are only for the convenience of explanation and do not limit the execution order of the steps. In actual operation, according to the technical requirements of specific implementation scenarios, the steps can be executed in an order different from that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

[0018] As Figure 1 shown, the method for identifying the current characteristics of a power distribution area based on the correlation peak judgment method includes the following steps:

[0019] S1. Collect the grid synchronous phase information to generate a characteristic signal with phase encoding; inject the characteristic signal into the power line to form an injected current signal;

[0020] Specifically, the grid synchronous phase information can be the time characteristics of the current and voltage of the grid, generating a unique characteristic signal, which has a specific phase encoding and can be regarded as a digital "tag" or "mark"; inject the generated characteristic signal into the power line (that is, let this "tag signal" enter the grid). After injection, an injected current signal will be formed, which can be used as a basis for analyzing and identifying specific current characteristics.

[0021] S2. Collect the line current signals containing the injected current signal at a predetermined key node in the power distribution area, preprocess the line current signals and extract the phase information to obtain a phase time series;

[0022] Specifically, the line current signal can be a "mixed signal" on the power line, and processes such as filtering and denoising are performed to extract clearer and more meaningful content.

[0023] S3. Based on the phase time series, construct a dynamically adaptive phase relationship network, extract the current characteristic phase difference pattern, and generate a topologically invariant feature vector;

[0024] Specifically, the dynamically adaptive phase relationship network will automatically optimize according to the changes in the data, similar to building a relationship graph that can adapt to changes in real time. The phase difference characteristic pattern can be the "behavior rule" between currents, reflecting the core information of the current characteristics. The topologically invariant feature vector is stable and can remain consistent even when the network structure changes, so as to be used for subsequent analysis or applications.

[0025] S4. Perform correlation analysis on the topologically invariant feature vector and a preset standard feature template, calculate the correlation peak and establish a judgment criterion to determine the identification result and the identification credibility index;

[0026] Specifically, the relevant peak value, i.e., the level of similarity, represents the matching degree between the eigenvector and the standard template. The higher the peak value, the more similar the two are. The decision criterion can ensure the scientificity and consistency of the identification process. For example, reaching a certain threshold can be recognized as a match.

[0027] S5. Dynamically adjust the system parameters according to the identification result and the identification credibility index, update the phase network structure and the decision threshold, and maintain a high recognition rate of the system under different conditions.

[0028] Specifically, according to the analysis result, automatically adjust some key parameters in the system. For example, adjust some technical details that affect the current feature identification to enable the system to better adapt to the current conditions. Update the phase network structure to enable it to more accurately capture the features of the current. At the same time, also update the "threshold" of the decision rule (such as increasing or decreasing the standard for matching features) to ensure more accurate classification or identification.

[0029] In this embodiment, by dynamically adjusting the system parameters and updating the network structure, the adaptability of the solution under different operating conditions is ensured; by using the topological invariant eigenvector to resist the interference caused by the change of the network topological structure, the stability of the analysis result is ensured; it has the characteristics of strong anti-interference ability, good environmental adaptability, and high recognition accuracy.

[0030] As Figure 2 shown, according to one aspect of the present application, the steps of forming the injected current signal include:

[0031] S11. Read the grid voltage signal, and obtain the grid synchronization phase information through zero-crossing detection and phase-locked loop technology;

[0032] S12. Based on the grid synchronization phase reference, design a phase offset sequence and generate a rotating phase encoding sequence;

[0033] S13. Based on a preset fixed carrier frequency, map the rotating phase encoding sequence onto the carrier signal to generate a phase-modulated carrier signal;

[0034] S14. Segment the phase-modulated carrier signal at pseudo-random time intervals to form an intermittent pulse sequence; inject the intermittent pulse sequence into the power line through a dedicated current to generate an injected current signal.

[0035] Specifically, read the grid voltage signal, determine the grid reference phase through zero-crossing detection technology, and use phase-locked loop technology to track the grid frequency and phase changes in real time to obtain the grid synchronization phase reference φ0. Based on the grid synchronization phase reference φ0, design a phase offset sequence [Δφ1, Δφ2,..., Δφ n, each offset value ranges from 0° to 360°, and the difference between adjacent offsets is not less than 45° to ensure anti-interference ability, forming a rotation phase encoding sequence Φ e . Select a fixed carrier frequency f k (usually 700 Hz), map the rotation phase encoding sequence Φ e onto the carrier signal, so that the signal has a specific phase offset Δφ i within each time interval t i , generating a phase-modulated carrier signal S k (t). The specific modulation method is: S k (t) = A·sin(2πf k t + φ0 + Δφ i ), t ∈ [t i , t i+1 , where A is the amplitude and is adaptively adjusted according to the line impedance and noise level. Segment the phase-modulated carrier signal S k (t) at pseudo-random time intervals T r , each segment lasting for a duration T p (typical value 100 ms), forming an intermittent pulse sequence P k (t). Inject the intermittent pulse sequence P k (t) into the power line through a dedicated current injection device to generate a line injection current signal I k (t).

[0036] According to one aspect of the present application, the steps of generating a rotation phase encoding sequence include:

[0037] Read the grid synchronous phase information, construct a phase offset basic unit based on the Gray code principle, and form a coding base sequence;

[0038] Apply the Hamming code principle to the coding base sequence to increase coding redundancy and generate a redundant coding sequence;

[0039] Analyze the phase jump characteristics in the redundant coding sequence, adjust the phase value through a phase mapping function to generate an optimized redundant sequence; convert the optimized redundant sequence into a differential phase coding form to form a differential phase sequence;

[0040] Combine the differential phase sequence with a preset absolute phase reference point to construct a rotation phase encoding sequence, ensuring that the rotation phase encoding sequence has the characteristics of forward error correction ability and a single sharp peak autocorrelation function.

[0041] Specifically, read the grid synchronous phase reference φ0, design a phase offset basic unit based on the Gray code principle, including four basic phase states of 0°, 90°, 180°, and 270°, to form a coding base sequence Φ β. The Gray code is adopted to ensure that only one bit changes between adjacent states, improving the anti-interference ability. For the coding base sequence Φ β Apply the Hamming code principle to increase coding redundancy to achieve forward error correction ability and generate a redundant coding sequence Φ r . The specific operation is to insert 3 parity bits after every 4-bit base sequence, and the values of the parity bits are calculated according to the following rules: P1 = B1 Θ B2 Θ B4; P2 = B1 Θ B3 Θ B4; P3 = B2 Θ B3 Θ B4; where B1 to B4 are the 4-bit values of the base sequence, P1 to P3 are the parity bits, and Θ represents the exclusive OR operation.

[0042] Analyze the phase jump characteristics in the redundant coding sequence Φ r and calculate the jump size between adjacent phases: ΔΦ i = |Φ r (i + 1) - Φ r (i)|, i ∈ [1, n - 1]; if there is a case where ΔΦ i < 45°, then adjust the phase value through the phase mapping function: Φ r '(i) = Φ r (i) + Δ i , where Δ i is the optimization increment to ensure that all ΔΦ i ≥ 45° after adjustment, and generate an optimized redundant sequence Φ r '. To enhance the anti-environmental interference ability, convert the optimized redundant sequence Φ r ' into a differential phase coding form: ΔΦ e (i)=Φ r '(i + 1) - Φ r '(i), i ∈ [1, n - 1]; where ΔΦ e (i) represents the phase change amount between adjacent time points, forming a differential phase sequence ΔΦ e . Combine the differential phase sequence ΔΦ e with the absolute phase reference point to construct the final rotational phase coding sequence Φ e : Φ e (1) = Φ r '(1); Φ e (i + 1) = Φ e (i) + ΔΦ e (i), i ∈ [1, n - 1]; ensure that the coding sequence satisfies: the minimum phase difference is not less than 45°; has forward error correction ability; the periodic autocorrelation function has a single sharp peak; the non-periodic autocorrelation sidelobe peak does not exceed 0.3 times of the main peak.

[0043] According to one aspect of the present application, the steps of forming an intermittent pulse sequence include:

[0044] Generate a pseudo-random sequence based on the linear congruence method, map the pseudo-random sequence to the time interval domain to generate a pseudo-random time interval sequence for determining the interval time between pulses.

[0045] Analyze the current characteristics of the substation area and the interference spectrum characteristics, and determine the optimal pulse duration based on the principle of maximizing the signal-to-noise ratio to control the duration of each pulse.

[0046] Optimize the pulse edge based on the optimal pulse duration to generate an optimized window function for smoothing the pulse edge.

[0047] Obtain the line impedance characteristics and background noise level, dynamically adjust the pulse energy to generate an adaptive amplitude sequence for modulating the pulse amplitude.

[0048] Combine the phase-modulated carrier signal with the pseudo-random time interval sequence, the adaptive amplitude sequence, and the optimized window function to construct an intermittent pulse sequence.

[0049] Specifically, generate a pseudo-random sequence based on the linear congruence method to ensure that it cannot be predicted by external interference. The calculation formula is: X_{i + 1} = (a·X i + c) mod m, where a, c, and m are carefully selected parameters to ensure that the sequence has good statistical characteristics. Map the random sequence to the time interval domain to generate a pseudo-random time interval sequence T r , whose value range is [0.5s, 2.5s]. Analyze the current characteristics of the substation area and the interference spectrum characteristics, and determine the optimal pulse duration T p : T p = arg max{SNR(T)}, where SNR(T) is the estimated signal-to-noise ratio when the pulse duration is T, usually taking values in the range of 80ms - 120ms, and output the optimal pulse duration T p . According to the line impedance characteristics and background noise level, dynamically adjust the pulse energy to balance the detection performance and the interference to the power grid: A(t) = A0·[1 + α·sin(2πf0t)]·[1 +β·N(t)], where A0 is the reference amplitude, α is the power frequency modulation coefficient (typical value 0.1), f0 is the power frequency, β is the noise response coefficient, and N(t) is the normalized noise index, to generate an adaptive amplitude sequence A(t). To reduce spectral leakage and adjacent frequency point interference, optimize the pulse edge design and apply an improved Hamming window function: w(t) = 0.54 - 0.46·cos(2πt / T p ) +γ·t / T p ·(1 - t / T p), where γ is an adaptive parameter that is dynamically adjusted according to the characteristics of the grid noise to generate an optimized window function w(t). Combining the phase-modulated carrier signal S k (t), the pseudo-random time interval sequence T r , the adaptive amplitude sequence A(t), and the optimized window function w(t), a complete discontinuous pulse sequence P k (t) is constructed: P k (t) = A(t)·w(t - t j )·sin(2πf k t + φ0 + Φ e (j)), t ∈ [t j , t j +T p ; P k (t) = 0, t ∈ (t j +T p , t j +T r (j)); where t j is the starting time of the j-th pulse, satisfying t_{j + 1} = t j + T p + T r (j).

[0050] As Figure 3 shown, according to one aspect of the present application, the steps of obtaining the phase time series include:

[0051] S21. Select a predetermined number of key nodes in the substation area, and use a high-precision current sensor to collect the line current signal, and set the sampling frequency to a value higher than the carrier frequency;

[0052] S22. Segment the line current signal into synchronous signal segments with the same length to ensure the time alignment of signals at different nodes;

[0053] S23. Filter the synchronous signal segments through a band-pass filter with a center frequency of the carrier frequency to remove power frequency and its harmonic interference, and obtain a filtered signal;

[0054] S24. Apply the Hilbert transform to the filtered signal to construct an analytic signal, extract the instantaneous phase information, and obtain the phase time series of each node.

[0055] Specifically, select N key nodes in the substation area (such as the low-voltage side of the transformer, important branch lines, terminal user access points, etc.), and use a high-precision current sensor to collect the line current signal I i (t) (i = 1, 2,... N), and set the sampling frequency to 10 kHz to ensure the accurate capture of the 700 Hz carrier signal. Based on the zero-phase mark of the power grid, the collected line current signal Ii (t) is segmented into synchronous signal segments I with consistent lengths i，j (t), where each segment corresponds to a complete injection pulse period, ensuring the time alignment of signals at different nodes. For the synchronous signal segment I i，j (t), a band - pass filter with a center frequency of f k and a bandwidth of 20 Hz is designed for preliminary filtering to remove power frequency and its harmonic interferences, obtaining the filtered signal I i，j '(t). For the filtered signal I i，j '(t), the Hilbert transform is applied to construct an analytic signal and extract the instantaneous phase information, obtaining the phase time series θ i，j (t) of each node. The extraction process is expressed as: I i，j '(t) + j·H{I i，j '(t)} = A i，j (t)·exp[j·θ i，j (t)], where H{} represents the Hilbert transform, A i，j (t) is the instantaneous amplitude, and θ i，j (t) is the instantaneous phase.

[0056] According to one aspect of the present application, the steps of extracting the instantaneous phase information to obtain the phase time series of each node include:

[0057] Perform zero - mean processing on the filtered signal and remove abnormal spikes to generate a pre - processed signal;

[0058] Use the Hilbert transform to calculate the orthogonal component of the pre - processed signal to generate the Hilbert transform result;

[0059] Based on the pre - processed signal and its Hilbert transform result, construct an analytic signal, extract the instantaneous phase, and obtain the original phase sequence;

[0060] Perform phase unwrapping on the original phase sequence to eliminate the discontinuity caused by phase jumps, obtaining the unwrapped phase sequence;

[0061] Perform trend analysis on the unwrapped phase sequence, extract and remove the linear phase growth component, and perform phase normalization to obtain the final phase time series.

[0062] Specifically, perform zero - mean processing on the filtered signal I i，j '(t) to remove the DC component and low - frequency drift, and at the same time apply the local polynomial fitting method to remove abnormal spikes, generating the pre - processed signal I i，j ''(t). Use an improved Hilbert transform to calculate the orthogonal component of the pre - processed signal I i，j ''(t), generating the imaginary part of the analytic signal: H{I i，j''(t)} = 1 / π∫ -∞ ∞ (I i，j ''(τ) / (t - τ))dτ; where t is the current time point, τ is the integration variable; to reduce the computational complexity, it is indirectly implemented using FFT: Calculate the FFT of the preprocessed signal I i，j ''(t) to obtain F(ω); Construct the frequency response H(ω) of the Hilbert transformer: H(ω) = -j·sgn(ω) = {-j, ω > 0; 0, ω = 0; j, ω < 0}; where ω is the angular frequency and sgn() is the sign function; Multiply in the frequency domain: G(ω) = F(ω)·H(ω), calculate the IFFT of G(ω) to obtain the Hilbert transform result H{I i，j ''(t)}.

[0063] Based on the preprocessed signal I i，j ''(t) and its Hilbert transform result H{I i，j ''(t)}, construct the analytic signal: z(t) = I i，j ''(t) + j·H{I i，j ''(t)}; Extract the instantaneous phase of the analytic signal: θ i，j _raw(t) = atan2(H{I i，j ''(t)}, I i，j ''(t)), where atan2 is the four - quadrant arctangent function, to obtain the original phase sequence θ i，j _raw(t). Due to the range limitation of the atan2 function in [-π, π], there are phase jumps in the original phase sequence θ i，j _raw(t), and phase unwrapping processing is required. Define the cumulative phase difference: Δθ(k) = θ i，j _raw(k + 1) - θ i，j _raw(k); When |Δθ(k)| > π, it is determined that phase wrapping occurs, and an appropriate 2π offset is added: θ i，j _unwrap(k + 1) = θ i，j _unwrap(k) + Δθ(k) + 2π·sgn(-Δθ(k)), where sgn is the sign function, and the unwrapped phase sequence θ i，j _unwrap(t) is output. Perform trend analysis on the unwrapped phase sequence θ i，j _unwrap(t) to extract and remove the linear phase growth component (the basic phase growth caused by the carrier frequency): θ_trend(t) = 2π·f k ·t + θ0, where f kis the carrier frequency, and θ0 is the initial phase. Calculate the detrended phase: θ i，j _detrend(t) = θ i，j _unwrap(t) - θ_trend(t); Finally, perform phase normalization and map it to the interval [0, 2π) to obtain the final phase time series θ i，j (t): θ i，j (t) = mod(θ i，j _detrend(t), 2π).

[0064] As Figure 4 shown, according to one aspect of the present application, the steps of generating the topological invariant feature vector include:

[0065] S31. For the phase time series, calculate the phase difference between nodes to form a phase difference matrix;

[0066] S32. Based on the phase difference matrix, apply the self-organizing network construction algorithm to establish the phase network topology of the substation area;

[0067] S33. Perform time series analysis on the phase network topologies obtained in a continuous predetermined number of cycles, extract the dynamic change features, and obtain the network evolution feature matrix;

[0068] S34. Extract the topological invariant feature vector with environmental invariance from the network evolution feature matrix as the core basis for feature identification.

[0069] Specifically, for the extracted phase time series θ i，j (t), calculate the phase difference between nodes to form a phase difference matrix D j . For N nodes, the phase difference matrix is N×N-dimensional, and the elements are expressed as: D j (i, k) = θ i，j (t) - θ k，j (t); i, k ∈ [1, N]; the phase difference values are uniformly mapped to the interval [-180°, 180°]. Based on the phase difference matrix D j , apply the new self-organizing network construction algorithm to establish the phase network topology G j of the substation area. The key steps of this algorithm include: setting the phase difference threshold λ p (the typical value is 15°); when |D j (i, k)| < λ p , a strong connection is established between nodes i and k; when λ p ≤ |D j (i, k)| < 2λ p , a weak connection is established between nodes i and k; when |D j (i, k)| ≥ 2λ pWhen there is no connection established between nodes i and k; construct a phase network topological structure G represented by a weighted adjacency matrix based on the connection strength j 。

[0070] For the phase network topological structure G obtained for consecutive multiple cycles j (j = 1, 2,..., M), perform time series analysis, extract its dynamic change characteristics, and obtain a network evolution feature matrix E t 。It mainly includes: node connection stability characteristics, community structure change characteristics, connection strength fluctuation characteristics, and network clustering coefficient change characteristics. Specifically: for the phase network topological structure G at each time point j Extract basic topological parameters, including: node degree distribution D_deg(j), clustering coefficient C_clus(j), path length distribution L_path(j), centrality index B_cent(j), and spectral eigenvalue E_spec(j), and combine them to form a network parameter time series P(j). Analyze the time stability of the network parameter time series P(j), and calculate the coefficient of variation of each parameter: CV_P(i) = std(P i (j)) / mean(P i (j)); j ∈ [1, M]; where P i represents the i-th network parameter, and generate a parameter stability vector S_P. Apply an improved Louvain algorithm to the phase network topological structure G j for community partitioning: in the initial state, each node is an independent community; iteratively optimize the community structure to maximize the modularity function Q: Q = 1 / 2m∑ i，j [W'(i, j) - k i k j / 2m]δ(c i , c j ); where m is the total edge weight, k i is the degree of node i, c i is the community to which node i belongs, δ is the Kronecker function, and W'(i, j) is the edge weight between node i and node j; track the evolution relationship of communities at different time points through a matching algorithm; generate a community structure evolution sequence C_comm(j). Based on a time sliding window, calculate the dynamic characteristic indicators of the network structure: structure change rate: a normalized measure of the network topology difference between adjacent time points; community stability: the proportion of community members remaining unchanged; critical node migration: the proportion of nodes with a significant change in centrality ranking; link volatility: the proportion of newly added / removed links, and combine them to form a network dynamic characteristic vector D n et. Integrate the network parameter time series P(j), the parameter stability vector S_P, the community structure evolution sequence C_comm(j), and the network dynamic characteristic vector D net, construct the complete network evolution feature matrix E t : E t = [P(j); S_P; C_comm(j); D n et], where each column of the matrix corresponds to a time point and each row corresponds to a feature dimension.

[0071] Apply topological data analysis methods to extract topological invariant feature vectors V with environmental invariance from the network evolution feature matrix E t . The key extraction steps are: constructing a multi-scale simplicial complex (from the phase network to topological features); calculating persistent homology features (reflecting the network structure stability); extracting the Betti number sequence (the algebraic representation of the network cycle structure); generating topological spectral features (the spectral representation of the network topological structure). Combine the above features to form the topological invariant feature vector V t , which is used as the core basis for feature identification. t

[0072] According to one aspect of the present application, the steps of applying the self-organizing network construction algorithm to establish the phase network topological structure of the substation area include:

[0073] Evaluate the time dimension stability of each element of the phase difference matrix, calculate the standard deviation, and generate the phase difference stability matrix;

[0074] Generate an adaptive threshold matrix based on the phase difference stability matrix;

[0075] Combine the adaptive threshold matrix with the phase difference stability matrix to establish the connection strength of continuous weights and generate the connection weight matrix;

[0076] Evaluate the time stability of each connection in the connection weight matrix based on historical data, calculate the connection reliability index, and generate the connection reliability matrix;

[0077] Combine the connection reliability matrix with the connection weight matrix to construct a comprehensive weight matrix and form the final phase network topological structure.

[0078] Specifically, evaluate the time dimension stability of each element of the phase difference matrix D j , calculate the standard deviation: σ_D(i ， k) = std(D j (i, k, t)), t ∈ [t1, t2], where [t1, t2] is the analysis time window, and generate the phase difference stability matrix σ_D. Based on the phase difference stability matrix σ_D, adopt an adaptive threshold design method. First, calculate the overall distribution characteristics of the stability matrix to obtain the global statistical parameters: μ_σ = mean(σ_D); σ_σ = std(σ_D). Then design the adaptive threshold function: λ​p (i ， k) = λ0 + α·σ_D(i ， k) - β·exp(-σ_D(i ， k) / γ), where λ0 is the base threshold (typical value 15°), and α, β, and γ are parameters automatically adjusted according to the current environmental noise, generating an adaptive threshold matrix λ p .

[0079] Introduce a multi-layer weighted connection mechanism. Instead of using the traditional binary connection (connected / unconnected), establish the connection strength of continuous weights according to the relationship between the phase difference and the threshold: W(i, k) = {1, |D j (i, k)| < λ p (i, k); exp(-(|D j (i, k)| - λ p (i, k)) 2 / μ), λ p (i, k) ≤ |D j (i, k)| < 2λ p (i, k); 0, |D j (i, k)| ≥ 2λ p (i, k)}, where μ is the attenuation coefficient, generating a connection weight matrix W. Based on historical data, evaluate the time stability of each connection, calculate the connection reliability index: R(i, k) = 1 - var(W_hist(i, k, t)) / mean(W_hist(i, k, t)), where W_hist is the historical connection weight record, and var and mean represent variance and mean respectively, generating a connection reliability matrix R. Combine the connection weight matrix W and the connection reliability matrix R to construct a comprehensive weight matrix: W'(i, k) = W(i, k)·R(i, k); Apply the minimum spanning tree algorithm in graph theory to ensure network connectivity and remove possible noisy connections at the same time. Specific steps: Construct a complete connection graph based on the comprehensive weight matrix W'; Apply the improved Kruskal algorithm to extract the minimum spanning tree; On the basis of the minimum spanning tree, supplement other edges with weights greater than the threshold η; Form the final phase network topology structure G j , represented by an adjacency matrix.

[0080] According to one aspect of the present application, the steps of extracting topologically invariant feature vectors with environmental invariance from the network evolution feature matrix include:

[0081] Construct a multi-scale simplicial complex based on the network evolution feature matrix, forming a sequence of simplicial complexes at different distance thresholds;

[0082] Apply the persistent homology algorithm to the sequence of simplicial complexes, calculate the homology groups, track the birth and death of homology features, and generate a persistent homology feature set;

[0083] Construct a persistence diagram based on the persistent homology feature set, extract the statistical features of the persistence diagram, and form a persistence diagram feature vector;

[0084] Extract the sequence of Betti numbers from the persistent homology feature set and generate a Betti number sequence vector;

[0085] Calculate the topological spectral features based on the network evolution feature matrix, including spectral energy, spectral gap, and spectral moments, and generate a topological spectral feature vector;

[0086] Integrate the persistence diagram feature vector, the Betti number sequence vector, and the topological spectral feature vector to construct the final topological invariant feature vector.

[0087] Specifically, based on the network evolution feature matrix E t Construct multi-scale simplicial complexes. First, set a series of distance thresholds {ε1, ε2,..., ε n}, and for each threshold ε i : Calculate the pairwise distance matrix D_E between feature vectors; when d(x, y) ≤ ε i , establish a connection between x and y; when all edges (x, y, z) exist, form a 2-simplex (triangle); and so on to construct higher-order simplices; obtain the sequence of simplicial complexes K(ε i ) at different thresholds. Apply the persistent homology algorithm to the sequence of simplicial complexes K(ε i ): For each threshold ε i calculate the homology groups H0(K(ε i ), H1(K(ε i ), H2(K(ε i ), where H0 corresponds to the number of connected components; H1 corresponds to the number of one-dimensional holes; H2 corresponds to the number of two-dimensional cavities. Track the "birth" and "death" of homology features to generate the persistence interval of homology features: [ε βi rth , ε death ; Calculate the persistence duration: pers = ε death - ε βi rth ; Generate the persistent homology feature set PH.

[0088] Construct a persistence diagram based on the persistent homology feature set PH: Using (ε βi rth , ε death) Represent each persistence interval with coordinates in a two-dimensional plane; calculate the statistical features of the persistence diagram: the number of persistence intervals, the maximum duration, the duration distribution, and the persistence interval density; extract the core features of the persistence diagram to form the persistence diagram feature vector PD. Extract the sequence of Betti numbers from the persistence homology feature set PH: β*(ε i ) = rank(H*(K(ε i ))), k = 0, 1, 2; where β* represents the k-dimensional Betti number, rank represents the rank of the homology group, H* represents the homology group, to generate the Betti number sequence vector β. Based on the network evolution feature matrix E t Calculate the topological spectral features: construct the Laplacian matrix L = D - A, where D is the degree matrix and A is the adjacency matrix; calculate the eigenvalue spectrum {λ1, λ2,..., λ n} of the Laplacian matrix; extract the spectral features: spectral feature energy: ∑λ i 2 ; spectral gap: λ2 - λ1; spectral moment: ∑λ i k / n; eigenvalue distribution entropy; generate the topological spectral feature vector S_topo. Integrate the persistence diagram feature vector PD, the Betti number sequence vector β, and the topological spectral feature vector S_topo to construct the final topological invariant feature vector V t : V t = [w1·PD; w2·β; w3·S_topo], where w1, w2, and w3 are feature weights, which are adaptively adjusted according to the current environment and historical recognition results. The feature vector design ensures: being robust to topological structure deformation; being insensitive to noise and outliers; retaining the essential topological properties of the network; being able to effectively distinguish the topological patterns of different feature signals.

[0089] As Figure 5 shown, according to one aspect of the present application, the steps for determining the identification result and the identification credibility index include:

[0090] S41. Through the injection experiment in an ideal environment, obtain the topological invariant feature vector under standard conditions as the standard feature template;

[0091] S42. Conduct a multi-dimensional correlation analysis on the standard feature template and the topological invariant feature vector, introduce the topological weighted correlation function, and obtain the correlation peak sequence;

[0092] S43. Based on the pre-stored environmental noise level and historical recognition results, construct an adaptive decision threshold function, compare the correlation peak sequence with the threshold, and output the decision result and the index of the best matching template;

[0093] S44. Calculate the decision confidence index based on the adaptive decision threshold function; combine the decision confidence index with the consistency of the continuous predetermined round of decision results to generate the final identification result and the identification credibility index.

[0094] Specifically, through the injection experiment in the ideal environment, obtain the topological invariant feature vector V under standard conditions s , which is used as the standard feature template for subsequent comparison. To adapt to different environmental conditions, multiple templates can be established to form a template library V s k (k = 1, 2,..., L). Compare the extracted topological invariant feature vector V t with the standard feature template V s k through multi-dimensional correlation analysis. The key is to introduce the topological weighted correlation function: R(k, τ) = Σ i w(i)·V t (i + τ)·V s k (i), where w(i) is the topological importance weight and τ is the relative time shift. Calculate the maximum correlation value for each template k: R max (k) = max{R(k, τ)}, and obtain the correlation peak sequence R max . Conduct adaptive correlation peak decision. Based on the environmental noise level and historical recognition results, design the adaptive decision threshold function η(σ), where σ represents the current environmental noise index. When the condition: max{R max (k)} > η(σ)·R0 is satisfied, it is determined that the characteristic signal exists, where R0 is the reference correlation value in the pure noise environment. Output the decision result D (0 indicates no characteristic signal, 1 indicates the existence of a characteristic signal), and the corresponding best matching template index k_best. Conduct confidence evaluation and multi-round decision fusion. To improve the decision reliability, conduct multi-round decisions and fuse the results. Calculate the confidence index: C = [max{R max (k)} - η(σ)·R0] / [η(σ)·R0]; generate the final identification result Y and the identification credibility index Q according to the confidence C and the consistency of continuous multi-round decisions.

[0095] Among them, the specific process of conducting adaptive correlation peak decision is as follows: Based on the characteristics of the currently collected signal, calculate the multi-dimensional environmental noise index: signal average power: P_avg = mean(|I i，j '(t)| 2 ); signal volatility: F_sig = std(I i，j '(t)) / mean(|I i，j '(t)|); harmonic interference degree: H_dist = ∑ iP(f_harm, i) / P_total, where P(f_harm, i) is the harmonic frequency point power and P_total is the total power; Phase stability: S_phase = std(θ i，j (t)); Generate the ambient noise vector σ, which includes all the above indicators. Read the historical decision data, extract the optimal thresholds and corresponding performance indicators under different environmental conditions: Establish the environment-threshold-performance dataset {(σ_hist, η_hist, Perf_hist)}; Use the K-nearest neighbor algorithm to find the K historical environmental conditions that are most similar to the current environment σ; Extract the average optimal threshold under these K environments: η_KNN = ∑ i=1 K w i ·η_hist, i / ∑ i=1 K w i , where w i is the weight based on environmental similarity; Generate the historical reference threshold η_KNN. Based on the ambient noise vector σ and the historical reference threshold η_KNN, construct the adaptive threshold function: η(σ) = η_0 + α·exp(β·||σ||) + γ·η_KNN, where η_0 is the base threshold, and α, β, γ are system parameters optimized offline through historical data, and generate the adaptive threshold η_cur for the current environment. Based on the adaptive threshold η_cur, the correlation peak sequence R max and the peak significance vector SIG, construct the multi-dimensional decision criterion: D_Peak(k) = {1, if R max (k) > η_cur·R0 && SNR(k) > SNR_min && Sharp(k) > Sharp_min; 0, otherwise}; where R0 is the reference correlation value in a pure noise environment, and SNR_min and Sharp_min are the minimum thresholds of signal-to-noise ratio and sharpness respectively, and generate the peak decision vector D_Peak. Analyze the peak decision vector D_Peak to determine the final decision result: If all D_Peak(k) = 0, it is determined that there is no characteristic signal, and the decision result D = 0; If there exists D_Peak(k) = 1, it is determined that there is a characteristic signal, and the decision result D = 1; At the same time, determine the best matching template: k_best = arg max{R max (k) | D_Peak(k) = 1}; Output the decision result D and the best matching template index k_best.

[0096] Among them, the confidence evaluation and multi-round decision fusion are specifically as follows: Based on the correlation peak sequence R max, the adaptive threshold η_cur and the pure noise benchmark related value R0 are used to calculate the decision confidence: C_single = [max{R max (k)} - η_cur·R0] / [η_cur·R0]; when the decision result D = 0, the confidence is negative, indicating the degree of certainty of the tendency to have no feature signal, and the single-round confidence C_single is generated. According to the single-round confidence C_single and the peak significance vector SIG, the weighted confidence is calculated: C_weight = C_single·[α·SNR(k_best) + β·Sharp(k_best)] / (α + β), where α and β are weight coefficients, and the weighted confidence C_weight is generated. The decision result D, the best matching template index k_best, and the weighted confidence C_weight of the current round are stored in the decision history cache, retaining the decision results of the most recent M rounds to form the decision history sequence H_D: H_D ={(D1, k_best1, C_weight1), (D2, k_best2, C_weight2),..., (D m , k_best m , C_weight m )}. Analyze the consistency in the decision history sequence H_D: Calculate the decision result consistency rate: Con_D = count(D i = D1) / M; calculate the best template consistency rate: Con_k = count(k_best i = k_best1) / count(D i = 1); calculate the mean and variance of the confidence: Mean_C = mean(C_weight i ); Var_C = var(C_weight i ); generate the decision consistency vector Con. Based on the decision consistency vector Con and the current decision result D, multi-round fusion decision-making is performed: Calculate the fusion confidence: C_fusion = Mean_C·(1 + γ·Con_D - Δ·sqrt(Var_C)), where γ and Δ are weighted coefficients; set the fusion threshold λ_fusion; when C_fusion > λ_fusion, the identification result Y = 1 is output, otherwise the identification result Y = 0 is output; calculate the final identification credibility index Q: Q = |C_fusion|·(1 + ε·Con_k), where ε is the template consistency weight coefficient.

[0097] According to one aspect of the present application, the steps of obtaining the relevant peak sequence include:

[0098] Calculate the preliminary similarity between the topological invariant feature vector and each standard feature template to generate a preliminary similarity vector;

[0099] Based on the preliminary similarity vector, evaluate the discrimination ability of each feature dimension and generate a feature weight vector;

[0100] Based on the principle of topological data analysis, construct a topological importance weighting function, assign weights to persistent homology features, Betti number features, and spectral features, and generate a topological importance function;

[0101] Combine the feature weight vector with the topological importance function to construct the final weight, calculate the topologically weighted correlation function, determine the maximum correlation value, and generate a correlation peak sequence;

[0102] Perform a peak significance evaluation on the correlation peak sequence, calculate the peak signal-to-noise ratio and sharpness index, and generate a peak significance vector.

[0103] Specifically, read the topological invariant feature vector V t and the template library V s k , and calculate the preliminary similarity of each template: S i nit(k) = cos(V t， V s k ) = (V t ·V s k ) / (||V t ||·||V s k ||), where cos represents the cosine similarity, and generate the preliminary similarity vector S i nit. Apply the feature sensitivity analysis method to evaluate the importance of each feature dimension: construct a historical recognition data set containing correctly recognized and misrecognized samples; for each feature dimension i, calculate its discrimination ability: d i = |μ i + - μ i - | / sqrt((σ i + ) 2 + (σ i - ) 2 ), where μ i + and σ i + are the mean and standard deviation of the correctly recognized samples in dimension i, respectively, and μ i - and σ i - are the statistics of the misrecognized samples; normalize the discrimination ability to weights: w(i) = di / ∑ j d j , generate the feature weight vector w. Based on the principle of topological data analysis, design a topological importance weighting function: for persistent homology features, the weight is proportional to the duration; for Betti number features, the weight is proportional to the topological stability; for spectral features, the weight is proportional to the eigenvalue stability, and comprehensively generate the topological importance function w_topo(i). Combine the feature weight vector w with the topological importance function w_topo(i) to construct the final weight: w_final(i) = w(i)·w_topo(i); based on the final weight, calculate the topological weighted correlation function: R(k, τ) = ∑ i w_final(i)·V t (i + τ)·V s k (i), where τ is the relative time shift, and calculate the maximum correlation value for each template k: R max (k) = max{R(k, τ)}, τ ∈ [-τ max , τ max , generate the correlation peak sequence R max . For the correlation peak sequence R max , conduct peak significance evaluation: calculate the background noise level of the correlation function: R n oise = median(R(k, τ)); calculate the peak signal-to-noise ratio: SNR(k) = R max (k) / R n oise; calculate the peak sharpness (steepness): Sharp(k) =R max (k) / (0.5·[R(k, τ max -Δ) + R(k, τ max +Δ)]), where Δ is a small time shift offset, usually taking 1 - 3 sampling points; generate the peak significance vector SIG, which includes the signal-to-noise ratio and sharpness indicators.

[0104] According to one aspect of the present application, the steps of updating the phase network structure and the decision threshold include:

[0105] S51. Real-time monitor the power grid parameters and environmental data of the substation area, calculate the environmental change index vector, and use the unsupervised clustering method to classify the current environmental state to obtain the environmental category identifier;

[0106] S52. According to the environmental category identifier and the identification credibility index, adjust the key parameters for constructing the phase network, generate the updated parameter set, and trigger the network reconstruction mechanism under extreme environmental conditions;

[0107] S53. Optimize the decision threshold function using the Bayesian optimization method based on the statistical results of decision outcomes under different environmental conditions in historical data, and generate an optimized threshold parameter set;

[0108] S54. After obtaining a high-confidence identification result under specific environmental conditions, extract the corresponding topological invariant feature vector, evaluate its difference from the existing template, and achieve incremental update of the template library;

[0109] S55. Regularly evaluate the overall performance of the system, generate a performance evaluation report, determine the next optimization focus, form a closed-loop feedback mechanism, and continuously improve the system performance.

[0110] Specifically, real-time monitor parameters such as the voltage, current, and power of the distribution network in the substation area, combine with environmental data such as temperature and humidity, and calculate the environmental change index vector E. According to the index vector, use the unsupervised clustering method to classify the current environmental state into predefined environmental categories, and obtain the environmental category identifier C_env. According to the environmental category identifier C_env and the identification confidence index Q, adjust the key parameters for constructing the phase network, including the phase difference threshold λ p and the connection strength weight, etc., to generate an updated parameter set P n et. Under extreme environmental conditions, trigger the network reconstruction mechanism to ensure the stability of the topological structure. Based on the statistical results of decision outcomes under different environmental conditions in historical data, optimize the decision threshold function η(σ) to make it more accurately adapt to various noise environments. The specific optimization process uses the Bayesian optimization method to balance the false alarm rate and the miss detection rate, and obtain the optimized threshold parameter set P_th. After obtaining a high-confidence identification result under specific environmental conditions, extract the corresponding topological invariant feature vector V t , evaluate its difference from the existing template. If the difference is significant, add it as a new template to the template library V s k , achieve incremental update of the template library, and improve the adaptability of the system to environmental changes. Regularly evaluate the overall performance of the system, including indicators such as recognition accuracy, response time, and environmental adaptability, and generate a performance evaluation report R_perf. Based on the evaluation results, determine the next optimization focus, form a closed-loop feedback mechanism, and continuously improve the system performance.

[0111] In a specific embodiment of the present application, a method for identifying the current characteristics of a distribution network area based on the correlation peak decision method realizes high-precision identification of the current characteristics of the distribution network area in a complex power grid environment through technical means combining phase coding and topological invariants. The following will specifically describe the implementation process of the method with numerical examples.

[0112] Step 1. Generation and injection of the rotated phase-coded feature signal.

[0113] 1.1. Acquisition of grid synchronization phase reference. In this embodiment, the grid voltage signal is collected by a voltage acquisition device installed on the low-voltage side of the substation transformer, and the sampling frequency is set to 20 kHz. The grid reference phase is determined by the zero-crossing detection technique, and the phase-locked loop technique is used to track the changes in the grid frequency and phase in real time. The actually measured grid voltage frequency is 49.92 Hz, and the phase reference φ0 = 32.7°. Through the phase-locked loop technique, the tracking frequency drift rate is about ±0.05 Hz / min, and the phase tracking accuracy is better than ±0.5°.

[0114] 1.2. Design of rotating phase coding sequence.

[0115] 1.2.1. Generation of coding base sequence. Based on the Gray code principle, the basic phase shift unit is designed to form the coding base sequence Φ β , where Φ β represents the basic phase state sequence. In this embodiment, four basic phase states are selected: 0°, 90°, 180°, 270°, to form the basic Gray code unit, and the coding base sequence Φ β = [0°, 90°, 270°, 180°].

[0116] 1.2.2. Phase coding redundancy design. Applying the Hamming code principle, coding redundancy is added to achieve the forward error correction ability. Three parity bits are inserted after every 4-bit base sequence, and the values of the parity bits are calculated according to the following rules: P1 = B1 Θ B2 Θ B4; where P1 is the first parity bit; B1, B2, B4 are the values of the 1st, 2nd, and 4th bits of the base sequence respectively; Θ represents the exclusive OR operation. P2 = B1 Θ B3 Θ B4; where P2 is the second parity bit; B1, B3, B4 are the values of the 1st, 3rd, and 4th bits of the base sequence respectively; Θ represents the exclusive OR operation. P3 = B2 Θ B3 Θ B4; where P3 is the third parity bit; B2, B3, B4 are the values of the 2nd, 3rd, and 4th bits of the base sequence respectively; Θ represents the exclusive OR operation.

[0117] Taking the base sequence Φ β = [0°, 90°, 270°, 180°] as an example, mapping the phase values to 2-bit coding: 0° = 00, 90° = 01, 180° = 11, 270° = 10, then: B1 = 00, B2 = 01, B3 = 10, B4 = 11; calculating the parity bits: P1 = 00 Θ 01 Θ 11 = 10 = 270°; P2 = 00 Θ 10 Θ 11 = 01 = 90°; P3 = 01 Θ 10 Θ 11 = 00 = 0°; generating the redundant coding sequence Φ r = [0°, 90°, 270°, 180°, 270°, 90°, 0°], where Φ rRepresents a phase-encoded sequence containing redundant check bits.

[0118] 1.2.3. Phase transition optimization. Analyze the phase transition characteristics in the redundant encoded sequence and calculate the transition magnitude between adjacent phases: ΔΦ i = |Φ r (i + 1) - Φ r (i)|; where ΔΦ i represents the absolute difference between the i-th phase value and the (i + 1)-th phase value; Φ r (i) represents the i-th phase value of the redundant encoded sequence; i represents the sequence index, with a value range of [1 ， n - 1], and n is the sequence length.

[0119] For Φ r = [0°, 90°, 270°, 180°, 270°, 90°, 0°], calculate the adjacent phase transitions: ΔΦ1 = |90° - 0°| = 90°; ΔΦ2 = |270° - 90°| = 180°; ΔΦ3 = |180° - 270°| = 90°; ΔΦ4 = |270° - 180°| = 90°; ΔΦ5 = |90° - 270°| = 180°; ΔΦ6 = |0° - 90°| = 90°. All ΔΦ i are not less than 45°, so no phase transition optimization is required. Optimize the redundant sequence Φ r ' = Φ r = [0°, 90°, 270°, 180°, 270°, 90°, 0°], where Φ r ' represents the phase-encoded sequence after transition optimization.

[0120] 1.2.4. Differential phase encoding conversion. Convert the optimized redundant sequence into differential phase encoding form: ΔΦ e (i) = Φ r '(i + 1) - Φ r '(i); where ΔΦ e (i) represents the phase change amount between adjacent time points; Φ r '(i) represents the i-th phase value of the optimized redundant sequence; i represents the sequence index, with a value range of [1, n - 1], and n is the sequence length.

[0121] Calculate the differential phase sequence: ΔΦ e (1) = 90° - 0° = 90°; ΔΦ e (2) = 270° - 90° = 180°; ΔΦ e(3) = 180° - 270° = -90°; ΔΦ e (4) = 270° - 180° = 90°; ΔΦ e (5) = 90° - 270° = -180°; ΔΦ e (6) = 0° - 90° = -90°; Obtain the differential phase sequence ΔΦ e = [90°, 180°, -90°, 90°, -180°, -90°], where ΔΦ e represents the differential phase coding sequence.

[0122] 1.2.5. Generation of the final coding sequence. Combine the differential phase sequence with the absolute phase reference point to construct the final rotational phase coding sequence: Φ e (1) = Φ r '(1) = 0°; where Φ e (1) represents the first phase value of the final coding sequence; Φ r '(1) represents the first phase value of the optimized redundant sequence. Φ e (i + 1) = Φ e (i) + ΔΦ e (i); where Φ e (i + 1) represents the (i + 1)-th phase value of the final coding sequence; Φ e (i) represents the i-th phase value of the final coding sequence; ΔΦ e (i) represents the i-th value of the differential phase sequence; i represents the sequence index, with a value range of [1, n - 1], and n is the sequence length.

[0123] Calculate the final coding sequence: Φ e (1) = 0°; Φ e (2) = 0° + 90° = 90°; Φ e (3) = 90° + 180° = 270°; Φ e (4) = 270° + (-90°) = 180°; Φ e (5) = 180° + 90° = 270°; Φ e (6) = 270° + (-180°) = 90°; Φ e (7) = 90° + (-90°) = 0°; Obtain the final rotational phase coding sequence Φ e = [0°, 90°, 270°, 180°, 270°, 90°, 0°], where Φ eRepresents the final rotation phase encoding sequence. Examine the characteristics of the final encoding sequence: The minimum phase difference is not less than 45°: The minimum phase difference is 90°, which meets the condition; It has forward error correction ability: Achieved through Hamming code check bits, and can correct 1-bit errors; Autocorrelation function characteristics: Verified by calculation, the periodic autocorrelation function has a single sharp peak, and the non-periodic autocorrelation sidelobe peak does not exceed 0.3 times the main peak.

[0124] 1.3. Carrier signal phase modulation. Select a fixed carrier frequency f k = 700Hz, map the rotation phase encoding sequence Φ e onto the carrier signal, so that the signal has a specific phase offset within each time interval t i . The modulation formula is: S k (t) =A·sin(2πf k t + φ0 + Φ e (i)); where S k (t) represents the phase-modulated carrier signal; A represents the amplitude, which is set to 1.5A in this embodiment; f k represents the carrier frequency; t represents time; φ0 represents the power grid synchronization phase reference; Φ e (i) represents the i-th phase value of the final encoding sequence; t∈[t i , t i+1 indicates that t is within the i-th time interval. In this embodiment, each phase duration is 20ms. For t = 0 - 20ms, the phase is 0°, and the signal expression is: S k (t) = 1.5·sin(2π·700·t + 32.7° + 0°) = 1.5·sin(2π·700·t + 32.7°); For t = 20 - 40ms, the phase is 90°, and the signal expression is: S k (t) = 1.5·sin(2π·700·t + 32.7° + 90°) = 1.5·sin(2π·700·t +122.7°); And so on, to generate a complete phase-modulated carrier signal.

[0125] 1.4. Intermittent pulse injection implementation.

[0126] 1.4.1. Pseudo-random time interval generation. Generate a pseudo-random sequence based on the linear congruence method: X i+1 = (a·X i +c) mod m; where X i+1 represents the (i + 1)-th value of the random sequence; X irepresents the i-th value of the random sequence; a represents the multiplier parameter, set to 1664525; c represents the increment parameter, set to 1013904223; m represents the modulus parameter, set to 2 32 ; mod represents the modulo operation. Set the initial value X0 = 12345, and calculate the pseudo-random sequence: X1 = (1664525·12345 + 1013904223) mod 2 32 = 1666278168; X2 = (1664525·1666278168 + 1013904223) mod 2 32 = 1253251530; X3 = (1664525·1253251530 + 1013904223) mod 2 32 = 3777254875... Map the random sequence to the time interval domain: T r (i) = 0.5 + 2·X i / 2 32 ; where T r (i) represents the i-th pseudo-random time interval, in seconds; X i represents the i-th value of the random sequence. Calculate the pseudo-random time interval sequence: T r (1) = 0.5 + 2·1666278168 / 2 32 = 0.5 + 2·0.388 = 1.276s; T r (2) = 0.5 + 2·1253251530 / 2 32 = 0.5 + 2·0.292 = 1.084s; T r (3) = 0.5 + 2·3777254875 / 2 32 = 0.5 + 2·0.879 = 2.258s... Obtain the pseudo-random time interval sequence T r = [1.276s, 1.084s, 2.258s,...], where T r represents the pseudo-random time interval sequence.

[0127] 1.4.2. Pulse duration optimization. By experimentally analyzing the current characteristics and interference spectrum characteristics of the substation area, the signal-to-noise ratio is measured at different pulse durations, and the optimal pulse duration is selected: when the pulse duration (ms) is 60, the signal-to-noise ratio (dB) is 12.3; when the pulse duration (ms) is 80, the signal-to-noise ratio (dB) is 14.8; when the pulse duration (ms) is 100, the signal-to-noise ratio (dB) is 16.5; when the pulse duration (ms) is 120, the signal-to-noise ratio (dB) is 15.7; when the pulse duration (ms) is 140, the signal-to-noise ratio (dB) is 14.2. Based on the principle of maximizing the signal-to-noise ratio, the optimal pulse duration T p = 100 ms, where T p represents the pulse duration.

[0128] 1.4.3. Adaptive adjustment of pulse energy. According to the line impedance characteristics and background noise level, the pulse energy is dynamically adjusted. The adjustment formula is: A(t) = A0·[1 + α·sin(2πf0t)]·[1 + β·N(t)]; where A(t) represents the adaptive amplitude; A0 represents the reference amplitude, set to 1.5 A; α represents the power frequency modulation coefficient, set to 0.1; f0 represents the power frequency, which is 50 Hz; β represents the noise response coefficient, set to 0.2; N(t) represents the normalized noise index, which is obtained through real-time measurement. Suppose at a certain moment t = 0.05 s, the measured normalized noise index N(t) = 0.15, then: A(0.05) = 1.5·[1 + 0.1·sin(2π·50·0.05)]·[1 + 0.2·0.15] = 1.5·[1 + 0.1·sin(π / 2)]·[1 + 0.03] = 1.5·[1 + 0.1·1]·1.03 = 1.5·1.1·1.03 = 1.698 A. At different moments, the amplitude is dynamically adjusted according to the noise level to generate the adaptive amplitude sequence A(t), where A(t) represents the adaptive amplitude sequence.

[0129] 1.4.4. Optimal design of pulse waveform. To reduce spectral leakage and adjacent frequency point interference, the pulse edge is optimized, and the improved Hamming window function is applied: w(t) = 0.54 - 0.46·cos(2πt / T p ) + γ·t / T p ·(1 - t / T p ); where w(t) represents the optimized window function; t represents time, and the value range is [0, T p ; T prepresents the pulse duration; γ represents the adaptive parameter, set to 0.05. Calculate the window function values. For example: w(0) = 0.54 - 0.46·cos(0) + 0.05·0 / 0.1·(1 - 0 / 0.1) = 0.54 - 0.46·1 + 0 = 0.08; w(0.05) = 0.54 - 0.46·cos(π) + 0.05·0.05 / 0.1·(1 - 0.05 / 0.1) = 0.54 - 0.46·(-1) + 0.05·0.5·0.5 = 0.54 + 0.46 + 0.0125 = 1.0125; w(0.1) = 0.54 - 0.46·cos(2π) + 0.05·0.1 / 0.1·(1 - 0.1 / 0.1) = 0.54 - 0.46·1 + 0 = 0.08. In this way, the complete optimized window function w(t) is generated, where w(t) represents the optimized window function.

[0130] 1.4.5, Generation of discontinuous pulse sequence. Combine the phase - modulated carrier signal S k (t), the pseudo - random time - interval sequence T r , the adaptive amplitude sequence A(t) and the optimized window function w(t) to construct the complete discontinuous pulse sequence: P k (t) = A(t)·w(t - t j )·sin(2πf k t + φ0 + Φ e (j)); where P k (t) represents the discontinuous pulse sequence; t ∈ [t j , t j + T p represents that t is within the j - th pulse time; A(t) represents the adaptive amplitude sequence; w(t - t j ) represents the translated optimized window function; f k represents the carrier frequency; φ0 represents the power - grid synchronization phase reference; Φ e (j) represents the j - th phase value of the final coding sequence; t j represents the starting time of the j - th pulse. P k (t) = 0; where t ∈ (t j + T p , t j + T p + T r (j)) represents that t is within the interval from the end of the j - th pulse to the start of the j + 1 - th pulse.

[0131] Calculation of pulse starting time: t j+1 = t j+ T p + T r (j); where t j+1 represents the start time of the (j + 1)-th pulse; t j represents the start time of the j-th pulse; T p represents the pulse duration; T r (j) represents the j-th pseudo-random time interval. Assuming t1 = 0, then: t2 = 0 + 0.1 + 1.276 = 1.376 s; t3 = 1.376 + 0.1 + 1.084 = 2.56 s; t4 = 2.56 + 0.1 + 2.258 = 4.918 s.... A discontinuous pulse sequence is injected into the power line through a dedicated current injection device to generate a line injection current signal.

[0132] Step 2: Multi-point phase signal acquisition and preprocessing.

[0133] 2.1. Key node signal acquisition. Select N = 5 key nodes in the substation area, including the low-voltage side of the transformer, the access points of the main branch lines, and the access points of the end users. Use high-precision current sensors to collect the line current signals I i (t) (i = 1, 2,... N), and set the sampling frequency to 10 kHz to ensure the accurate capture of the 700 Hz carrier signal. The acquisition time for each node is 5 seconds, and the number of sampling points obtained is 5 × 10000 = 50000 points.

[0134] 2.2. Current signal segmentation and synchronization. Based on the zero-phase mark of the power grid, segment the collected line current signals I i (t) into synchronous signal segments I i，j (t) with the same length. Each segment corresponds to a complete injection pulse period to ensure the time alignment of signals at different nodes. Using the pulse start time calculated above as the mark, the segmentation is as follows: Segment 1: t ∈ [0, 0.1]; Segment 2: t ∈ [1.376, 1.476]; Segment 3: t ∈ [2.56, 2.66]...; The number of sampling points included in each segment is 0.1 × 10000 = 1000 points.

[0135] 2.3. Band-pass filtering and enhancement processing. For the synchronous signal segments I i，j (t), design a band-pass filter with a center frequency of f k = 700 Hz and a bandwidth of 20 Hz for preliminary filtering to remove power frequency and its harmonic interference, and obtain the filtered signal I i，j '(t). Design a second-order Butterworth band-pass filter, and its transfer function is: H(s) = (s / Q) / (s 2+ s / (Q + 1); where s represents the Laplace variable; Q represents the quality factor, set to f k / BW = 700 / 20 = 35. Apply this filter to each synchronization signal segment to obtain a filtered signal.

[0136] 2.4. Phase information extraction.

[0137] 2.4.1. Sampling signal preprocessing. For the filtered signal I i，j '(t), perform zero-mean processing to remove the DC component and low-frequency drift: I i，j ''(t) = I i，j '(t) - mean(I i，j '(t)); where I i，j ''(t) represents the preprocessed signal; I i，j '(t) represents the filtered signal; mean(I i，j '(t)) represents the mean of the filtered signal. At the same time, apply the local polynomial fitting method to remove abnormal spikes and generate the preprocessed signal I i，j ''(t).

[0138] 2.4.2. Analytic signal construction. Indirectly implement the Hilbert transform using FFT, and calculate the quadrature component of the preprocessed signal I i，j ''(t): Calculate the FFT of the preprocessed signal I i，j ''(t) to obtain F(ω). Construct the frequency response H(ω) of the Hilbert transformer: H(ω) = -j·sgn(ω); where j represents the imaginary unit; sgn(ω) represents the sign function, which takes the value 1 when ω > 0, 0 when ω = 0, and -1 when ω < 0. Multiply in the frequency domain: G(ω) = F(ω)·H(ω); where G(ω) represents the frequency-domain product; F(ω) represents the Fourier transform of the preprocessed signal; H(ω) represents the frequency response of the Hilbert transformer. Calculate the IFFT of G(ω) to obtain the Hilbert transform result H{I i，j ''(t)}.

[0139] 2.4.3. Instantaneous phase calculation. Based on the preprocessed signal I i，j ''(t) and its Hilbert transform result H{I i，j ''(t)}, construct the analytic signal: z(t) = I i，j ''(t) + j·H{I i，j ''(t)}; where z(t) represents the analytic signal; j represents the imaginary unit; I i，j ''(t) represents the preprocessed signal; H{I i，j''(t)} represents the result of the Hilbert transform of the preprocessed signal. Extract the instantaneous phase of the analytic signal: θ i，j _raw(t) = atan2(H{I i，j ''(t)}, I i，j ''(t)); where θ i，j _raw(t) represents the original phase sequence; atan2 represents the four-quadrant arctangent function; H{I i，j ''(t)} represents the result of the Hilbert transform; I i，j ''(t) represents the preprocessed signal.

[0140] 2.4.4. Phase unwrapping. Due to the range limitation of the atan2 function within [-π, π], there are phase jumps in the original phase sequence θ i，j _raw(t), and phase unwrapping processing is required. Define the cumulative phase difference: Δθ(k) = θ i，j _raw(k + 1) - θ i，j _raw(k); where Δθ(k) represents the phase difference between the k-th sampling point and the (k + 1)-th sampling point; θ i，j _raw(k) represents the original phase value of the k-th sampling point; θ i，j _raw(k + 1) represents the original phase value of the (k + 1)-th sampling point. When |Δθ(k)| > π, it is determined that phase wrapping occurs, and an appropriate 2π offset is added: θ i，j _unwrap(k + 1) = θ i，j _unwrap(k) + Δθ(k) + 2π·sgn(-Δθ(k)); where θ i，j _unwrap(k + 1) represents the unwrapped phase value of the (k + 1)-th sampling point; θ i，j _unwrap(k) represents the unwrapped phase value of the k-th sampling point; Δθ(k) represents the phase difference; sgn represents the sign function.

[0141] 2.4.5. Phase detrending and normalization. Perform trend analysis on the unwrapped phase sequence θ i，j _unwrap(t), extract and remove the linear phase growth component (the basic phase growth caused by the carrier frequency): θ_trend(t) = 2π·f k ·t + θ0; where θ_trend(t) represents the trend phase; f k represents the carrier frequency; t represents time; θ0 represents the initial phase. The initial phase θ0 = 0.58 rad is obtained by least squares fitting. Calculate the detrended phase: θ i，j _detrend(t) = θ i，j_unwrap(t) - θ_trend(t); where θ i，j _detrend(t) represents the detrended phase; θ i，j _unwrap(t) represents the unwrapped phase sequence; θ_trend(t) represents the trend phase. Finally, phase normalization is performed and mapped to the interval [0, 2π): θ i，j (t) = mod(θ i，j _detrend(t), 2π); where θ i，j (t) represents the final phase time series; mod represents the modulo operation; θ i，j _detrend(t) represents the detrended phase.

[0142] Step 3: Self-organizing phase network construction and feature extraction.

[0143] 3.1. Calculation of the phase difference matrix. For the extracted phase time series θ i，j (t), calculate the phase difference between nodes to form the phase difference matrix D j . For N = 5 nodes, the phase difference matrix is 5×5 dimensional, and the elements are expressed as: D j (i, k) = θ i，j (t) - θ k，j (t); where D j (i, k) represents the phase difference between node i and node k; θ i，j (t) represents the phase time series of node i at the j-th time segment; θ k，j (t) represents the phase time series of node k at the j-th time segment; i, k ∈ [1 ， N] represents the node index. The phase difference values are uniformly mapped to the interval [-180°, 180°], for example: D j (1, 2) = 78.3°; D j (1, 3) = -45.7°; D j (1, 4) = 125.4°; D j (1, 5) = -102.1°; D j (2, 3) = -124.0°...

[0144] 3.2. Construction of the phase network topology structure.

[0145] 3.2.1. Evaluation of the phase difference stability. Perform a time-dimensional stability evaluation on each element of the phase difference matrix D j , and calculate the standard deviation: σ_D(i, k) = std(D j(i, k, t)); where σ_D(i, k) represents the standard deviation of the phase difference between node i and node k; D j (i, k, t) represents the phase difference between node i and node k at different times t; std represents the standard deviation calculation function; t ∈ [t1, t2] represents the analysis time window. For example, for nodes 1 and 2, within the time window t ∈ [0, 100 ms], the standard deviation is calculated from the phase differences at 10 time points: D j (1, 2, t1) = 78.3°; D j (1, 2, t2) = 77.9°;...; D j (1, 2, t 10 ) = 79.1°; σ_D(1, 2) = 0.5°. Similarly, the phase difference stability of all node pairs is calculated to generate the phase difference stability matrix σ_D.

[0146] 3.2.2. Adaptive connection threshold design. Based on the phase difference stability matrix σ_D, an adaptive threshold design method is adopted. First, calculate the overall distribution characteristics of the stability matrix: μ_σ = mean(σ_D) = (0.5 + 3.2 + 1.8 + 2.1 + 1.5 + 0.8 + 4.3 + 2.7 + 1.2 + 0.9) / 10 = 1.9°; where μ_σ represents the average value of the stability matrix; mean represents the mean calculation function. σ_σ = std(σ_D) = 1.2°; where σ_σ represents the standard deviation of the stability matrix; std represents the standard deviation calculation function. Design the adaptive threshold function: λ p (i, k) = λ0 + α·σ_D(i, k) - β·exp(-σ_D(i, k) / γ); where λ p (i, k) represents the adaptive threshold between node i and node k; λ0 represents the base threshold, set to 15°; α represents the stability weighting coefficient, set to 1.5; β represents the exponential decay coefficient, set to 5; γ represents the exponential proportionality coefficient, set to 1; σ_D(i, k) represents the standard deviation of the phase difference between node i and node k. For example, for nodes 1 and 2: λ p (1, 2)= 15 + 1.5·0.5 - 5·exp(-0.5 / 1) = 15 + 0.75 - 5·0.607 = 15 + 0.75 - 3.035 =12.715°. Calculate the adaptive thresholds for all node pairs to generate the adaptive threshold matrix λ p .

[0147] 3.2.3. Establishment of multi-layer weighted connections. Introduce a multi-layer weighted connection mechanism. According to the relationship between the phase difference and the threshold, establish the connection strength of continuous weights: W(i, k) = 1; where |D j (i, k)| < λ p (i, k) represents the case where the absolute value of the phase difference is less than the threshold; W(i, k) represents the connection weight between node i and node k. W(i, k) = exp(-(|D j (i, k)| - λ p (i, k)) 2 / μ); where λ p (i, k) ≤ |D j (i, k)| < 2λ p (i, k) represents the case where the absolute value of the phase difference is greater than or equal to the threshold and less than twice the threshold; μ represents the attenuation coefficient, set to 100. W(i, k) = 0; where |D j (i, k)| ≥ 2λ p (i, k) represents the case where the absolute value of the phase difference is greater than or equal to twice the threshold.

[0148] For example, for node 1 and node 2, |D j (1, 2)| = 78.3° > 2·λ p (1, 2) = 2·12.715° = 25.43°, so W(1, 2) = 0. For node 1 and node 3, |D j (1, 3)| = 45.7° > 2·λ p (1, 3) = 2·21.8° = 43.6°, so W(1, 3) = 0. For node 2 and node 4, |D j (2, 4)| = 47.1° > λ p (2, 4) = 15.7° and |D j (2, 4)| < 2·λ p (2, 4) = 31.4°, so: W(2, 4) = exp(-(|47.1 - 15.7|) 2 / 100) = exp(-(31.4) 2 / 100) = exp(-9.86) = 0.0005. Calculate the connection weights of all node pairs to generate the connection weight matrix W.

[0149] 3.2.4. Connection reliability assessment. Based on historical data, evaluate the time stability of each connection and calculate the connection reliability index: R(i, k) = 1 - var(W_hist(i, k, t)) / mean(W_hist(i, k, t)); where R(i, k) represents the connection reliability index between node i and node k; var represents the variance calculation function; mean represents the mean calculation function; W_hist(i, k, t) represents the historical connection weight record. For example, for nodes 2 and 4, the historical connection weight record is: [0.0005, 0.0008, 0.0006, 0.0004, 0.0007], then: var(W_hist(2, 4, t)) = 2.5e-8 mean(W_hist(2, 4, t)) = 0.0006 R(2, 4) = 1 - 2.5e-8 / 0.0006 = 1 - 4.17e-5 = 0.99996. Calculate the connection reliability of all node pairs and generate the connection reliability matrix R.

[0150] 3.2.5. Topology optimization. Combine the connection weight matrix W and the connection reliability matrix R to construct the comprehensive weight matrix: W'(i, k) = W(i, k)·R(i, k); where W'(i, k) represents the comprehensive weight between node i and node k; W(i, k) represents the connection weight; R(i, k) represents the connection reliability. Apply the improved Kruskal algorithm to extract the minimum spanning tree to ensure network connectivity and remove possible noisy connections at the same time: construct a complete connection graph based on W'; sort all edges in descending order of W' values; select the edge with the largest weight and add it to the spanning tree if no loop is formed; repeat the above steps until all nodes are connected; supplement other edges with weights greater than the threshold η = 0.5; form the final phase network topology structure G j , represented by the adjacency matrix.

[0151] 3.3. Phase network time evolution analysis.

[0152] 3.3.1. Network structure parameter extraction. For the phase network topology structure G at each time point jExtract basic topological parameters: Node degree distribution D_deg(j): The distribution of the number of connections of each node; for example, at a certain moment j, the degrees of each node are: [2, 3, 2, 3, 2]. Clustering coefficient C_clus(j): Represents the degree of connection tightness between node neighbors; for example, at a certain moment j, the clustering coefficients of each node are: [0.5, 0.67, 0.5, 0.67, 0.5]. Path length distribution L_path(j): The distribution of the shortest paths between nodes; for example, the average shortest path length of the network is 2.1. Centrality index B_cent(j): A measure of the importance of a node in the network; for example, at a certain moment j, the betweenness centralities of each node are: [0.2, 0.5, 0.1, 0.4, 0.1]. Spectral feature value E_spec(j): The eigenvalue of the network adjacency matrix; for example, at a certain moment j, the eigenvalues are: [3.5, 1.2, 0.5, -0.8, -1.5]. Combine to form the network parameter time series P(j), which represents the topological characteristics of the network at different time points j.

[0153] 3.3.2. Structural stability analysis. Analyze the temporal stability of the network parameter time series P(j), and calculate the coefficient of variation of each parameter: CV_P(i) = std(P i (j)) / mean(P i (j)); where CV_P(i) represents the coefficient of variation of the i-th network parameter; std represents the standard deviation calculation function; mean represents the mean calculation function; P i (j) represents the change sequence of the i-th network parameter over time j. For example, for the node degree distribution: mean(D_deg(j)) = 2.4, std(D_deg(j)) = 0.5, CV_P(1) = 0.5 / 2.4 = 0.21; Calculate the coefficient of variation of all network parameters to generate the parameter stability vector S_P.

[0154] 3.3.3. Community structure identification and tracking. Apply the improved Louvain algorithm to the phase network topology G j for community partitioning: Initially, each node is an independent community; Iteratively optimize the community structure to maximize the modularity function Q: Q = (1 / 2m)·∑(i ， j)[W'(i, j) - (k i ·k j ) / (2m)]·Δ(c i , c j ); where m represents the total edge weight; k i represents the degree of node i; c i represents the community to which node i belongs; Δ represents the Kronecker function, when c i =c jIt takes the value of 1 at a certain time, otherwise 0. For example, in a network of 5 nodes, it may be divided into 2 communities: Community 1: Nodes 1, 3, 5; Community 2: Nodes 2, 4. By using a matching algorithm to track the evolution relationship of communities at different time points, a community structure evolution sequence C_comm(j) is generated.

[0155] 3.3.4. Extraction of network dynamic characteristics. Based on a time-sliding window, calculate the dynamic characteristic indicators of the network structure: Structure change rate: A normalized measure of the difference in network topology between adjacent time points. For example, the structure change rate at times t and t+1 is 0.15; Community stability: The proportion of community members remaining unchanged. For example, the community stability between two adjacent time points is 0.8; Critical node migration: The proportion of nodes with a significant change in centrality ranking. For example, the critical node migration proportion is 0.2; Link volatility: The proportion of newly added / removed links. For example, the link volatility is 0.1. Combine them to form the network dynamic characteristic vector D n et.

[0156] 3.3.5. Construction of the network evolution feature matrix. Integrate the time series of network parameters P(j), the parameter stability vector S_P, the community structure evolution sequence C_comm(j), and the network dynamic characteristic vector D n et to construct the complete network evolution feature matrix E t : E t = [P(j); S_P; C_comm(j); D n et]; where E t represents the network evolution feature matrix; each column corresponds to a time point, and each row corresponds to a feature dimension.

[0157] 3.4. Extraction of topological invariant features.

[0158] 3.4.1. Construction of multi-scale simplicial complexes. Based on the network evolution feature matrix E t Construct multi-scale simplicial complexes. Set a series of distance thresholds {ε1, ε2,..., ε n} = {0.1, 0.2, 0.3, 0.4, 0.5}, for each threshold ε i : Calculate the pairwise distance matrix D_E of feature vectors. For example, for feature vectors v1 and v2, the Euclidean distance is 0.25. When d(x, y) ≤ ε iWhen a connection is established between x and y. For example, at the threshold ε2 = 0.2, no connection is established between v1 and v3 because d(v1, v3) = 0.32 > 0.2. When all edges (x, y, z) exist, a 2-simplex (triangle) is formed. For example, at the threshold ε3 = 0.3, a triangle is formed among v1, v2, and v4 because d(v1, v2) = 0.25; d(v1, v4) = 0.18; d(v2, v4) = 0.22, all of which are less than 0.3. Higher-order simplices are constructed in this way to obtain a sequence of simplicial complexes K(ε i ).

[0159] 3.4.2 Persistent homology calculation. Apply the persistent homology algorithm to the sequence of simplicial complexes K(ε i ): For each threshold ε i calculate the homology groups H0(K(ε i )), H1(K(ε i )), H2(K(ε i )); H0 represents the number of connected components; H1 represents the number of one-dimensional holes; H2 represents the number of two-dimensional cavities. For example, at the threshold ε1 = 0.1, H0 = 5 (5 independent connected components), H1 = 0 (no holes), H2 = 0 (no cavities); at the threshold ε3 = 0.3, H0 = 1 (all points are connected), H1 = 2 (there are 2 holes), H2 = 0 (no cavities). Track the "birth" and "death" of homology features. For example, a certain hole appears (is born) at ε2 = 0.2 and disappears (dies) at ε4 = 0.4. Generate the persistence interval of homology features: [ε βi rth , ε death ; for example, [0.2, 0.4] represents the persistence interval of a hole. Calculate the persistence duration: pers = ε death - ε βi rth ; for example, pers = 0.4 - 0.2 = 0.2. Generate the persistent homology feature set PH, which contains information about all persistent homology features.

[0160] 3.4.3 Persistent diagram construction and feature extraction. Construct a persistent diagram based on the persistent homology feature set PH: Using (ε βi rth , ε deathThe coordinates represent each persistence interval on a two-dimensional plane; for example, the point (0.2, 0.4) represents the aforementioned annulus persistence interval. Calculate the statistical features of the persistence diagram: the number of persistence intervals: for example, there are 10 persistence intervals in total; the maximum persistence time: for example, max(pers) = 0.3; the persistence time distribution: for example, [0.1, 0.1, 0.2, 0.2, 0.3, 0.1, 0.2, 0.1, 0.1, 0.2]; the density of persistence intervals: for example, the density of persistence intervals in the region [0.1, 0.2]×[0.3, 0.4] is 0.4. Extract the core features of the persistence diagram to form the persistence diagram feature vector PD.

[0161] 3.4.4. Extraction of the Betti number sequence. Extract the Betti number sequence from the persistent homology feature set PH: β*(ε i ) = rank(H*(K(ε i ))); where β* represents the k-dimensional Betti number; rank represents the rank of the homology group; k takes values of 0, 1, or 2. For example: β0(0.1)= 5, β0(0.2) = 3, β0(0.3) = 1, β0(0.4) = 1, β0(0.5) = 1; β1(0.1) = 0, β1(0.2) = 1, β1(0.3) = 2, β1(0.4) = 1, β1(0.5) = 0; β2(0.1) = 0, β2(0.2) = 0, β2(0.3) = 0, β2(0.4) =1, β2(0.5) = 1. Generate the Betti number sequence vector β = [β0, β1, β2].

[0162] 3.4.5. Calculation of topological spectrum features. Based on the network evolution feature matrix E t Calculate the topological spectrum features: Construct the Laplacian matrix L = D - A; where L represents the Laplacian matrix; D represents the degree matrix, the diagonal elements are the degrees of each node, and the other elements are 0; A represents the adjacency matrix. Calculate the eigenvalue spectrum {λ1, λ2,..., λ n} of the Laplacian matrix; for example, the eigenvalue spectrum is {0, 0.8, 1.5, 2.3, 3.4}. Extract the spectrum features: the energy of the eigenvalue spectrum: ∑λ i 2 = 0 2 + 0.8 2 + 1.5 2 + 2.3 2 +3.4 2 = 19.74; the spectral gap: λ2 - λ1 = 0.8 - 0 = 0.8; the spectral moment: ∑λ i k / n, for example, the first moment = (0 + 0.8 + 1.5 + 2.3 + 3.4) / 5 = 1.6; eigenvalue distribution entropy: -∑(λ i / ∑λ i )·log(λ i / ∑λ i ) = 1.37; generate the topological spectral feature vector S_topo.

[0163] 3.4.6. Topological invariant feature integration. Integrate the persistent diagram feature vector PD, the Betti number sequence vector β, and the topological spectral feature vector S_topo to construct the final topological invariant feature vector V t : V t = [w1·PD; w2·β; w3·S_topo]; where w1, w2, and w3 represent feature weights and are adaptively adjusted according to the current environment and historical recognition results. In this embodiment, according to the historical recognition performance, w1 = 0.5, w2 = 0.3, and w3 = 0.2 are set. The feature vector design ensures: robustness to topological structure deformation; insensitivity to noise and outliers; retention of the essential topological characteristics of the network; and ability to effectively distinguish the topological patterns of different feature signals.

[0164] Step Four: Phase difference correlation peak decision and identification decision.

[0165] 4.1. Standard feature template generation. Through the injection experiment under ideal conditions, obtain the topological invariant feature vector V s under standard conditions as the standard feature template for subsequent comparison. To adapt to different environmental conditions, multiple templates are established to form a template library V s k (k = 1, 2,..., L), where L = 5, representing the standard templates under 5 typical environmental conditions. For example, the template library contains the standard feature vectors under the following environmental conditions: V s 1 : Standard load condition (residential area, weekday daytime); V s 2 : Heavy load condition (commercial area, weekday daytime); V s 3 : Light load condition (residential area, weekday night); V s 4 : Fluctuating load condition (industrial area, weekday daytime); V s 5 : Mixed load condition (comprehensive area, weekend).

[0166] 4.2. Phase difference multi-dimensional correlation calculation.

[0167] 4.2.1. Preliminary evaluation of multi-template similarity. Read the topological invariant feature vector V t and the template library V s k , and calculate the preliminary similarity of each template: S i nit(k) = cos(V t , V s k ) = (V t ·V s k ) / (‖V t ‖·‖V s k ‖); where S i nit(k) represents the preliminary similarity with the k-th template; cos represents the cosine similarity; V t represents the currently extracted topological invariant feature vector; V s k represents the k-th standard feature template; · represents the vector dot product; ‖·‖ represents the vector norm.

[0168] For example, calculate the preliminary similarity with 5 templates: S i nit(1) = cos(V t , V s 1 ) = 0.85; S i nit(2)= cos(V t , V s 2 ) = 0.62; S i nit(3) = cos(V t , V s 3 ) = 0.73; S i nit(4) = cos(V t , V s 4 ) = 0.58; S i nit(5) = cos(V t , V s 5 ) = 0.67. Generate the preliminary similarity vector S i nit = [0.85, 0.62, 0.73, 0.58, 0.67], where S i nit represents the preliminary similarity vector.

[0169] 4.2.2. Calculation of Feature Dimension Importance Weights. Apply the feature sensitivity analysis method to evaluate the importance of each feature dimension: construct a historical recognition dataset that includes correctly recognized and misrecognized samples; for each feature dimension i, calculate its discrimination ability: d i = |μ i + - μ i - | / sqrt((σ i + ) 2 + (σ i - ) 2 );where d i represents the discrimination ability of the i-th feature dimension; μ i + represents the mean of the correctly recognized samples in dimension i; μ i - represents the mean of the misrecognized samples in dimension i; σ i + represents the standard deviation of the correctly recognized samples in dimension i; σ i - represents the standard deviation of the misrecognized samples in dimension i. For example, for the first feature dimension: μ_1 + = 0.82, μ_1 - = 0.45, σ_1 + = 0.12, σ_1 - = 0.15; d_1 = |0.82 - 0.45| / sqrt(0.12 2 + 0.15 2 ) = 0.37 / 0.19 = 1.95. Normalize the discrimination ability to weights: w(i) = d i / ∑ j d j ; where w(i) represents the weight of the i-th feature dimension; d i represents the discrimination ability of the i-th feature dimension; ∑ j d j represents the sum of the discrimination abilities of all feature dimensions. Assume the feature vector has 20 dimensions, and the calculated sum of discrimination abilities ∑ j d j = 30, then: w(1) = 1.95 / 30 = 0.065. Generate the feature weight vector w.

[0170] 4.2.3. Design of Topological Importance Weighting Function. Based on the principle of topological data analysis, design a topological importance weighting function: for persistent homology features, the weight is proportional to the duration. For example, for a feature with a duration of 0.3, the weight is 0.3 / 0.5 = 0.6, where 0.5 is the maximum duration. For Betti number features, the weight is proportional to the topological stability. For example, if β1 remains unchanged within the threshold interval [0.3, 0.4], its stability weight is 0.7. For spectral features, the weight is proportional to the eigenvalue stability. For example, if the eigenvalue stability is 0.85, the weight is 0.85. Combine them to generate the topological importance function w_topo(i).

[0171] 4.2.4. Calculation of Multidimensional Weighted Correlation Function. Combine the feature weight vector w with the topological importance function w_topo(i) to construct the final weight: w_final(i) = w(i)·w_topo(i); where w_final(i) represents the final weight of the i-th feature dimension; w(i) represents the i-th element of the feature weight vector; w_topo(i) represents the i-th element of the topological importance function. For example, for the first feature dimension: w_final(1) = 0.065·0.6 = 0.039. Based on the final weight, calculate the topological weighted correlation function: R(k, τ) = ∑ i w_final(i)·V t (i + τ)·V s k (i); where R(k, τ) represents the correlation value with the k-th template at time shift τ; w_final(i) represents the final weight of the i-th feature dimension; V t (i + τ) represents the i-th element of the current feature vector at time shift τ; V s k (i) represents the i-th element of the k-th standard feature template; ∑ i represents the sum over all feature dimensions. Calculate the maximum correlation value for each template k: R max (k) = max{R(k, τ)}; where R max (k) represents the maximum correlation value with the k-th template; max represents the maximum value function; R(k, τ) represents the correlation value with the k-th template at time shift τ; τ ∈ [-τ max , τ max represents the time shift range, and τ max is set to 10. For example, calculate the maximum correlation values with 5 templates: R max (1) = 0.92; R max (2) = 0.65; R max (3) = 0.78; R max(4) = 0.61; R max (5) = 0.70. Generate the relevant peak sequence R max = [0.92, 0.65, 0.78, 0.61, 0.70], where R max represents the relevant peak sequence.

[0172] 4.2.5. Peak significance evaluation. For the relevant peak sequence R max conduct peak significance evaluation: Calculate the background noise level of the correlation function: R n oise = median(R(k, τ)); where R n oise represents the background noise level; median represents the median function; R(k, τ) represents the correlation values of all templates at all time shifts. For example, R n oise = 0.35. Calculate the peak signal-to-noise ratio: SNR(k) = R max (k) / R n oise; where SNR(k) represents the peak signal-to-noise ratio of the k-th template; R max (k) represents the maximum correlation value with the k-th template; R n oise represents the background noise level. For example, calculate the peak signal-to-noise ratios of 5 templates: SNR(1) = 0.92 / 0.35 = 2.63; SNR(2) = 0.65 / 0.35 = 1.86; SNR(3) = 0.78 / 0.35 = 2.23; SNR(4) = 0.61 / 0.35 = 1.74; SNR(5) = 0.70 / 0.35 = 2.00. Calculate the peak sharpness (steepness): Sharp(k) = R max (k) / (0.5·[R(k, τ max - Δ) + R(k, τ max + Δ)]); where Sharp(k) represents the peak sharpness of the k-th template; R max (k) represents the maximum correlation value with the k-th template; R(k, τ max - Δ) represents the correlation value at Δ points to the left of the maximum correlation value; R(k, τ max + Δ) represents the correlation value at Δ points to the right of the maximum correlation value; Δ represents the offset, set to 2. For example, for template 1, assume R(1, τ max - Δ) = 0.75, R(1, τ max + Δ) = 0.70: Sharp(1) = 0.92 / (0.5·[0.75 + 0.70]) = 0.92 / 0.725 = 1.27. Generate the peak significance vector SIG, including the signal-to-noise ratio and sharpness metrics.

[0173] 4.3. Adaptive correlation peak decision.

[0174] 4.3.1. Environmental noise index calculation. Based on the characteristics of the currently collected signal, calculate the multi-dimensional environmental noise index: Signal average power: P_avg = mean(|I i，j '(t)| 2 ); where P_avg represents the signal average power; mean represents the mean function; I i，j '(t) represents the filtered signal; |·| represents the absolute value. For example, P_avg = 1.25. Signal volatility: F_sig = std(I i，j '(t)) / mean(|I i，j '(t)|); where F_sig represents the signal volatility; std represents the standard deviation function; mean represents the mean function; I i，j '(t) represents the filtered signal; |·| represents the absolute value. For example, F_sig = 0.35. Degree of harmonic interference: H_dist = ∑ i P(f_harm, i) / P_total; where H_dist represents the degree of harmonic interference; P(f_harm, i) represents the power of the i-th harmonic frequency point; P_total represents the total power; ∑ i represents the summation over all harmonic frequency points. For example, H_dist = 0.15. Phase stability: S_phase = std(θ i，j (t)); where S_phase represents the phase stability; std represents the standard deviation function; θ i，j (t) represents the phase time series. For example, S_phase = 0.08. Generate the environmental noise vector σ = [P_avg, F_sig, H_dist, S_phase] = [1.25, 0.35, 0.15, 0.08], where σ represents the environmental noise vector.

[0175] 4.3.2. Historical Threshold Data Analysis. Read historical decision data, extract the optimal thresholds and corresponding performance metrics under different environmental conditions: establish an environment-threshold-performance dataset {(σ_hist, η_hist, Perf_hist)}. Use the K-Nearest Neighbor algorithm to find the K = 3 historical environmental conditions that are most similar to the current environment σ. For example, the three most similar historical environments and their corresponding optimal thresholds are: (σ_hist1, η_hist1) = ([1.20, 0.32, 0.14, 0.07], 0.75); (σ_hist2, η_hist2) = ([1.28, 0.38, 0.17, 0.09], 0.80); (σ_hist3, η_hist3) = ([1.18, 0.33, 0.13, 0.06], 0.73).

[0176] Calculate the environmental similarity weights: w i = 1 / d(σ, σ_hist, i); where w i represents the weight of the i-th historical environment; d(σ, σ_hist, i) represents the distance between the current environment σ and the i-th historical environment σ_hist, i. For example, calculate the weights of the three historical environments: d(σ, σ_hist1) = 0.08, w1 = 12.5; d(σ, σ_hist2) = 0.05, w2 = 20.0; d(σ, σ_hist3) = 0.10, w3 = 10.0.

[0177] Extract the average optimal threshold under these K environments: η_KNN = ∑ i=1 K w i ·η_hist, i / ∑ i=1 K w i ; where η_KNN represents the historical reference threshold based on K-Nearest Neighbor; w i represents the weight of the i-th historical environment; η_hist, i represents the optimal threshold of the i-th historical environment; ∑ i=1 K represents the sum over the K nearest neighbor environments. Calculate the historical reference threshold: η_KNN =(12.5·0.75 + 20.0·0.80 + 10.0·0.73) / (12.5 + 20.0 + 10.0) = 34.075 / 42.5 = 0.802.

[0178] 4.3.3. Construction of Adaptive Threshold Function. Based on the environmental noise vector σ and the historical reference threshold η_KNN, construct the adaptive threshold function: η(σ) = η0 + α·exp(β·‖σ‖) + γ·η_KNN; where η(σ) represents the adaptive threshold function; η0 represents the basic threshold, set to 0.5; α represents the exponential coefficient, set to 0.1; β represents the environmental sensitivity coefficient, set to -0.5; ‖σ‖ represents the norm of the environmental noise vector; γ represents the historical reference weight, set to 0.3; η_KNN represents the historical reference threshold. Calculate the adaptive threshold of the current environment: ‖σ‖ = sqrt(1.25 2 + 0.35 2 + 0.15 2 + 0.08 2 ) = sqrt(1.5625 + 0.1225 + 0.0225 + 0.0064) = sqrt(1.7139) = 1.31 η(σ) = 0.5 + 0.1·exp(-0.5·1.31) + 0.3·0.802 = 0.5 + 0.1·0.52 + 0.24 = 0.5 + 0.052 + 0.24 =0.792. Generate the adaptive threshold η_cur = 0.792 of the current environment, where η_cur represents the adaptive threshold of the current environment.

[0179] 4.3.4. Design of Peak Decision Criterion. Based on the adaptive threshold η_cur, the correlation peak sequence R max and the peak significance vector SIG, construct the multi-dimensional decision criterion: D_Peak(k) = 1; where R max (k) > η_cur·R0 && SNR(k) >SNR_min && Sharp(k) > Sharp_min means that the k-th template meets the decision condition; R max (k) represents the maximum correlation value with the k-th template; η_cur represents the adaptive threshold of the current environment; R0 represents the reference correlation value in a pure noise environment, set to 0.1; SNR(k) represents the peak signal-to-noise ratio of the k-th template; SNR_min represents the minimum signal-to-noise ratio threshold, set to 1.5; Sharp(k) represents the peak sharpness of the k-th template; Sharp_min represents the minimum sharpness threshold, set to 1.2. D_Peak(k) =0; where it means that the k-th template does not meet the decision condition.

[0180] For example, determine whether 5 templates meet the judgment condition: η_cur·R0 = 0.792·0.1 = 0.0792; D_Peak(1): 0.92 > 0.0792 && 2.63 > 1.5 && 1.27 > 1.2, which meets the condition, so D_Peak(1) = 1; D_Peak(2): 0.65 > 0.0792 && 1.86 > 1.5 && (assuming Sharp(2) = 1.15) < 1.2, which does not meet the condition, so D_Peak(2) = 0; D_Peak(3): 0.78 > 0.0792 && 2.23 > 1.5 && (assuming Sharp(3) = 1.25) > 1.2, which meets the condition, so D_Peak(3) = 1; D_Peak(4): 0.61 > 0.0792 && 1.74 > 1.5 && (assuming Sharp(4) = 1.10) < 1.2, which does not meet the condition, so D_Peak(4) = 0; D_Peak(5): 0.70 > 0.0792 && 2.00 > 1.5 && (assuming Sharp(5) = 1.22) > 1.2, which meets the condition, so D_Peak(5) = 1. Generate the peak judgment vector D_Peak = [1, 0, 1, 0, 1], where D_Peak represents the peak judgment vector.

[0181] 4.3.5. Comprehensive judgment of multiple templates. Analyze the peak judgment vector D_Peak to determine the final judgment result: If all D_Peak(k) = 0, it is determined that there is no characteristic signal, and the judgment result D = 0; If there exists D_Peak(k) = 1, it is determined that there is a characteristic signal, and the judgment result D = 1; At the same time, determine the best matching template: k_best = arg max{R max (k) | D_Peak(k) = 1}; where k_best represents the index of the best matching template; arg max represents taking the independent variable that makes the function reach the maximum value; R max (k) represents the maximum correlation value with the k-th template; D_Peak(k)=1 means that the k-th template meets the judgment condition. In this example, the template indices that meet the conditions are 1, 3, and 5, and the corresponding R max values are 0.92, 0.78, and 0.70 respectively. Therefore, k_best = 1. Output the judgment result D = 1 and the best matching template index k_best = 1.

[0182] 4.4. Confidence evaluation and multi-round judgment fusion.

[0183] 4.4.1. Single-round decision confidence calculation. Based on the correlation peak sequence R max , the adaptive threshold η_cur, and the pure noise reference correlation value R0, calculate the decision confidence: C_single = [max{R max (k)} - η_cur·R0] / [η_cur·R0]; where C_single represents the single-round confidence; max{R max (k)} represents the maximum correlation peak; η_cur represents the adaptive threshold of the current environment; R0 represents the reference correlation value in the pure noise environment. Calculate the single-round confidence: C_single = [0.92 - 0.792·0.1] / [0.792·0.1] = [0.92 - 0.0792] / 0.0792 = 0.8408 / 0.0792 = 10.6. When the decision result D = 0, the confidence is negative, indicating the degree of certainty of the tendency to have no feature signal. Generate the single-round confidence C_single = 10.6.

[0184] 4.4.2. Confidence weighted allocation. According to the single-round confidence C_single and the peak significance vector SIG, calculate the weighted confidence: C_weight = C_single·[α·SNR(k_best) + β·Sharp(k_best)] / (α + β); where C_weight represents the weighted confidence; C_single represents the single-round confidence; α represents the signal-to-noise ratio weight, set to 0.7; β represents the sharpness weight, set to 0.3; SNR(k_best) represents the peak signal-to-noise ratio of the best-matched template; Sharp(k_best) represents the peak sharpness of the best-matched template. Calculate the weighted confidence: C_weight = 10.6·[0.7·2.63 + 0.3·1.27] / (0.7 + 0.3) = 10.6·[1.841 + 0.381] / 1 = 10.6·2.222 = 23.55. Generate the weighted confidence C_weight = 23.55.

[0185] 4.4.3. Multi-round decision result caching. Store the decision result D, the best-matched template index k_best, and the weighted confidence C_weight of the current round in the decision history cache, and retain the decision results of the most recent M = 5 rounds to form the decision history sequence H_D: H_D = {(D1, k_best1, C_weight1), (D2, k_best2, C_weight2),..., (D m , k_best m , C_weight m )}; where H_D represents the decision history sequence; Di Represents the judgment result of the i-th round; k_best i Represents the index of the best matching template in the i-th round; C_weight i Represents the weighted confidence in the i-th round. For example, the judgment history sequence after adding the current round result is: H_D = {(1, 1, 23.55), (1, 1, 18.72), (1, 3, 15.30), (0, -1, -8.45), (1, 1, 20.18)}.

[0186] 4.4.4. Judgment Consistency Analysis. Analyze the consistency in the judgment history sequence H_D: Calculate the judgment result consistency rate: Con_D = count(D i = D1) / M; where Con_D represents the judgment result consistency rate; count(D i = D1) represents the number of historical judgments that are consistent with the latest judgment result; M represents the length of the historical sequence. For example, Con_D = 4 / 5 = 0.8. Calculate the best template consistency rate: Con_k = count(k_best i = k_best1) / count(D i = 1); where Con_k represents the best template consistency rate; count(k_best i = k_best1) represents the number of historical best templates that are consistent with the latest best template; count(D i = 1) represents the number of historical judgments with a judgment result of 1. For example, Con_k = 3 / 4 = 0.75. Calculate the mean and variance of the confidence: Mean_C = mean(C_weight i ); where Mean_C represents the mean of the weighted confidence; mean represents the mean function; C_weight i represents the weighted confidence in the i-th round. For example, Mean_C = (23.55 + 18.72 + 15.30 - 8.45 + 20.18) / 5 = 69.3 / 5 = 13.86 Var_C = var(C_weight i ); where Var_C represents the variance of the weighted confidence; var represents the variance function; C_weight i represents the weighted confidence in the i-th round. For example, Var_C = 171.25. Generate the judgment consistency vector Con = [Con_D, Con_k, Mean_C, Var_C] = [0.8, 0.75, 13.86, 171.25], where Con represents the judgment consistency vector.

[0187] 4.4.5. Multi - round decision fusion. Based on the decision consistency vector Con and the current decision result D, perform multi - round fusion decision: Calculate the fusion confidence: \(C_{fusion}=Mean\_C\cdot(1 + \gamma\cdot Con\_D-\Delta\cdot\sqrt{Var\_C})\); where \(C_{fusion}\) represents the fusion confidence; \(Mean\_C\) represents the mean of the weighted confidence; \(\gamma\) represents the decision consistency weight, set to 0.5; \(Con\_D\) represents the decision result consistency rate; \(\Delta\) represents the variance penalty coefficient, set to 0.01; \(\sqrt{Var\_C}\) represents the standard deviation of the weighted confidence. Calculate the fusion confidence: \(C_{fusion}=13.86\cdot(1 + 0.5\cdot0.8 - 0.01\cdot\sqrt{171.25}) = 13.86\cdot(1 + 0.4 - 0.01\cdot13.09)=13.86\cdot(1.4 - 0.1309)=13.86\cdot1.2691 = 17.59\). Set the fusion threshold \(\lambda_{fusion}=5\); when \(C_{fusion}>\lambda_{fusion}\), output the identification result \(Y = 1\), otherwise output the identification result \(Y = 0\). In this example, \(C_{fusion}=17.59>\lambda_{fusion}=5\), so the identification result \(Y = 1\). Calculate the final identification credibility index Q: \(Q = |C_{fusion}|\cdot(1+\varepsilon\cdot Con\_k)\); where Q represents the identification credibility index; \(|C_{fusion}|\) represents the absolute value of the fusion confidence; \(\varepsilon\) represents the template consistency weight coefficient, set to 0.3; \(Con\_k\) represents the best template consistency rate. Calculate the identification credibility index: \(Q = 17.59\cdot(1 + 0.3\cdot0.75)=17.59\cdot(1 + 0.225)=17.59\cdot1.225 = 21.55\). Output the final identification result \(Y = 1\) and the identification credibility index \(Q = 21.55\).

[0188] Step Five. Adaptive Optimization and Parameter Update.

[0189] 5.1. Environmental Change Monitoring and Classification. Real - time monitor parameters such as the voltage, current, and power of the distribution network in the substation area, combine with environmental data such as temperature and humidity, and calculate the environmental change index vector E: Grid voltage volatility: \(\Delta V=\frac{max(V)-min(V)}{mean(V)}\), for example, \(\Delta V = 0.05\); Load change rate: \(\Delta P=\frac{std(P)}{mean(P)}\), for example, \(\Delta P = 0.12\); Total harmonic distortion rate: \(THD=\sqrt{\sum\frac{V_h}{V_1}}\), for example, \(THD = 0.03\); Temperature index: T 2 / V_1 2 ) norm = (T - T_min) / (T_max - T_min), for example T n orm = 0.65; Humidity index: H n orm = (H - H_min) / (H_max - H_min), for example H n orm = 0.42; Generate the environmental change index vector E = [ΔV, ΔP, THD, T n orm, H n orm] = [0.05, 0.12, 0.03, 0.65, 0.42], where E represents the environmental change index vector. According to the index vector, using the unsupervised clustering method, classify the current environmental state into predefined environmental categories to obtain the environmental category identifier C_env. In this example, use the K-means clustering algorithm to classify the current environment into category 2 (medium load, relatively high temperature).

[0190] 5.2. Adaptive update of phase network parameters. According to the environmental category identifier C_env and the identification credibility index Q, adjust the key parameters for constructing the phase network: the phase difference threshold λ p Adjust: λ pn ew = λ p _base·(1 + ρ1·(C_env - C_base) / C_base - ρ2·(1 - min(Q, Q_max) / Q_max)); where λ pn ew represents the updated phase difference threshold; λ p _base represents the base phase difference threshold, set to 15°; ρ1 represents the environmental adjustment coefficient, set to 0.2; C_env represents the current environmental category; C_base represents the reference environmental category, set to 1; ρ2 represents the credibility adjustment coefficient, set to 0.15; Q represents the identification credibility index; Q_max represents the maximum credibility threshold, set to 30. Calculate the updated phase difference threshold: λ pnew = 15·(1 + 0.2·(2 - 1) / 1 - 0.15·(1 - min(21.55, 30) / 30)) = 15·(1 + 0.2 - 0.15·(1 - 21.55 / 30)) = 15·(1.2 - 0.15·(1 - 0.72)) = 15·(1.2 - 0.15·0.28) = 15·(1.2 - 0.042) = 15·1.158 = 17.37°. Connection strength weight adjustment: w_scale = 1 + μ·log(Q / Q_base); where w_scale represents the connection strength weight adjustment coefficient; μ represents the weight adjustment coefficient, set to 0.1; Q represents the identification credibility index; Q_base represents the reference credibility, set to 10. Calculate the connection strength weight adjustment coefficient: w_scale = 1 + 0.1·log(21.55 / 10) = 1 + 0.1·log(2.155) = 1 + 0.1·0.77 = 1 + 0.077 = 1.077. Based on the above calculations, generate the updated parameter set P n et, including the updated phase difference threshold λ pn ew = 17.37° and the connection strength weight adjustment coefficient w_scale = 1.077. Under extreme environmental conditions (e.g., Q < 5), trigger the network reconstruction mechanism to ensure the stability of the topological structure.

[0191] 5.3. Dynamic optimization of the decision threshold. Based on the statistics of decision results under different environmental conditions in historical data, optimize the decision threshold function η(σ): Collect environment-threshold-performance data: (σ i , η i , Perf i ). Use the Bayesian optimization method to optimize the threshold function parameters (η0, α, β, γ) to balance the false alarm rate and the miss detection rate. The objective function is defined as: F(η0, α, β, γ) = ω1·(1 - TPR) + ω2·FPR; where F represents the optimization objective function; ω1 represents the miss detection penalty weight, set to 0.7; TPR represents the true positive rate; ω2 represents the false alarm penalty weight, set to 0.3; FPR represents the false positive rate. Through the Bayesian optimization algorithm, obtain the optimized parameters: η0_opt = 0.55, α_opt = 0.09, β_opt = -0.48, γ_opt = 0.32. Generate the optimized threshold parameter set P_th = [η0_opt, α_opt, β_opt, γ_opt] = [0.55, 0.09, -0.48, 0.32], where P_th represents the optimized threshold parameter set.

[0192] 5.4 Incremental update of the template library. After obtaining high-confidence recognition results under specific environmental conditions, extract the corresponding topological invariant feature vector V t , and evaluate its difference from the existing templates: Calculate the similarity between the current feature vector and all templates: sim(k) = cos(V t , V s k ); where sim(k) represents the similarity with the k-th template; cos represents the cosine similarity; V t represents the current feature vector; V s k represents the k-th standard feature template. For example, calculate the similarity with 5 existing templates: sim(1) = 0.85; sim(2) = 0.62; sim(3) = 0.73; sim(4) = 0.58; sim(5) = 0.67. Determine whether the maximum similarity is lower than the new template threshold: max_sim = max{sim(k)}; where max_sim represents the maximum similarity; max represents the maximum function; sim(k) represents the similarity with the k-th template. In this example, max_sim = 0.85. If max_sim < θ n ew (where θ n ew represents the new template threshold, set to 0.9) and Q > Q_min (where Q_min represents the minimum confidence threshold, set to 15), then add the current feature vector as a new template to the template library. In this example, max_sim = 0.85 < θ n ew = 0.9, and Q = 21.55 > Q_min = 15, so add the current feature vector as a new template (index 6) to the template library. In this way, the incremental update of the template library is realized, improving the adaptability of the system to environmental changes.

[0193] 5.5. System Performance Evaluation and Feedback. Regularly evaluate the overall system performance, including indicators such as recognition accuracy, response time, and environmental adaptability: Recognition accuracy: Acc = (TP + TN) / (TP + TN + FP + FN); where Acc represents the recognition accuracy; TP represents the number of true positives; TN represents the number of true negatives; FP represents the number of false positives; FN represents the number of false negatives. For example, based on the last 100 judgments, TP = 78, TN = 15, FP = 5, FN = 2, then: Acc = (78 + 15) / (78 + 15 + 5 + 2) = 93 / 100 = 93%. Sensitivity (recall rate): Sen = TP / (TP + FN); where Sen represents the sensitivity; TP represents the number of true positives; FN represents the number of false negatives. For example, Sen = 78 / (78 + 2) = 78 / 80 = 97.5%. Specificity: Spe = TN / (TN + FP); where Spe represents the specificity; TN represents the number of true negatives; FP represents the number of false positives. For example, Spe = 15 / (15 + 5) = 15 / 20 = 75%. Average response time: RT_avg = mean(RT i ); where RT_avg represents the average response time; mean represents the mean function; RT i represents the response time of the i-th judgment. For example, the average response time of the last 100 judgments is 215 ms. Environmental adaptability: Adapt = 1 - std(Acc j ) / mean(Acc j ); where Adapt represents the environmental adaptability; std represents the standard deviation function; mean represents the mean function; Acc j represents the recognition accuracy under different environmental conditions j. For example, the recognition accuracies under 5 different environments are [93%, 91%, 89%, 92%, 90%], then: Adapt = 1 - std([0.93, 0.91, 0.89, 0.92, 0.90]) / mean([0.93, 0.91, 0.89, 0.92, 0.90]) = 1 - 0.015 / 0.91 = 1 - 0.0165 = 0.9835. Generate a performance evaluation report R_perf, including all the above performance indicators. Based on the evaluation results, determine the next optimization focus, form a closed-loop feedback mechanism, and continuously improve the system performance. In this embodiment, the specificity index is relatively low (75%), and it is necessary to focus on optimizing the decision threshold to reduce the false alarm rate.

[0194] To verify the effectiveness of this embodiment, the performance of three methods for identifying the current characteristics of the substation area was compared: the traditional frequency analysis method (Method A), the method based on phase encoding but without using topological features (Method B), and the method based on the correlation peak decision method proposed in this embodiment (Method C). These three methods were tested in different power grid environments, and the following results were obtained: In the standard load (residential area during the day) environment, the accuracy rate of Method A was 92%, that of Method B was 94%, and that of Method C was 95%; in the heavy load (commercial area during the day) environment, the accuracy rate of Method A was 85%, that of Method B was 89%, and that of Method C was 93%; in the harmonic interference (industrial area) environment, the accuracy rate of Method A was 65%, that of Method B was 78%, and that of Method C was 91%; in the distributed energy access (photovoltaic area) environment, the accuracy rate of Method A was 68%, that of Method B was 82%, and that of Method C was 90%; in the season transition period (obvious temperature and humidity changes) environment, the accuracy rate of Method A was 73%, that of Method B was 76%, and that of Method C was 89%; in the extreme weather (thunderstorm weather) environment, the accuracy rate of Method A was 58%, that of Method B was 71%, and that of Method C was 87%; the average accuracy rate of Method A was 73.5%; the average accuracy rate of Method B was 81.7%; the average accuracy rate of Method C was 90.8%.

[0195] It can be seen from this that this embodiment shows a high recognition accuracy rate in various test environments. Especially in complex environments such as harmonic interference, distributed energy access, and extreme weather, it shows obvious advantages. Compared with the traditional frequency analysis method, the average accuracy rate has increased by 17.3 percentage points; compared with the method that only uses phase encoding but does not use topological features, the average accuracy rate has increased by 9.1 percentage points. In addition, this embodiment also shows excellent adaptability to environmental changes. When the environmental conditions suddenly change from the standard load to the harmonic interference environment, the accuracy rates of Method A and Method B decrease by 27 and 16 percentage points respectively, while this embodiment only decreases by 4 percentage points, demonstrating excellent environmental adaptability.

[0196] According to another aspect of this application, a system for identifying the current characteristics of the substation area based on the correlation peak decision method is composed of a power frequency interference filter, an FIR band-pass filter, an IIR digital notch filter, a characteristic frequency point extraction module, and a correlation peak decision module. Among them, the power frequency interference filter is an algorithm module for filtering the power frequency interference of the alternating current in the substation area; the FIR band-pass filter is an algorithm module for filtering the power frequency harmonic interference; the IIR digital notch filter is an algorithm module for filtering the harmonic interference within the passband of the band-pass filter; the characteristic frequency point extraction module is an algorithm module for extracting specific frequency points from the filtered signal; the correlation peak decision algorithm is an algorithm module for finally identifying the current characteristics.

[0197] In a specific embodiment, taking the injection of a current characteristic signal of 833.3 Hz by the resistance switching method as an example. The sampling frequency at the sending end is set to Fs = 5000 Hz, the load on / off period is set to 1200 μs, and the load on / off duty cycle is set to 1:3, that is, in each load on / off cycle, the load is off for 800 μs and on for 400 μs. The sampling interval is dt = 1 / Fs, and the sampling time of the k-th sampling point is kdt. At this time, the load current can be expressed as I(k) = {(311 / R sin(2πf_0 kdt + q), mod(k, 6) = 0, 1; 0, mod(k, 6) = 2, 3, 4, 5)}; where mod(k, 6) represents the remainder of k divided by 6. The implementation of the method for identifying the current characteristics of the transformer area based on the correlation peak decision method includes the following steps:

[0198] Step 1: Design a power frequency interference filter. According to the generation method of the characteristic current at the sending end, it can be found that the least common period of the power frequency alternating current cycle and the load on / off cycle is 3 power frequency alternating current cycles. This characteristic can stagger the characteristic current signals of adjacent power frequency alternating current cycles. By subtracting the waveforms of adjacent power frequency alternating current cycles, the interference of 50.00 Hz power frequency alternating current can be effectively removed, and the characteristic current signal can be enhanced. Let the sampling frequency be Fs = 5400 Hz. At this time, 108 points can be sampled in one power frequency alternating current cycle. The characteristic current at this time is: I_1(k) = I(k + 108) - I(k).

[0199] Step 2: Design a FIR band-pass filter. The two frequency points of the characteristic current signal are 783.3 Hz and 883.3 Hz. Through FFT analysis of the information after filtering out the power frequency, it can be seen that in addition to the above two frequency points, there are still a large number of power frequency harmonic interferences. Therefore, a FIR band-pass filter is designed, and the two cut-off frequencies of the band-pass filter are 730 Hz and 930 Hz. The main parameters of the FIR band-pass filter are: sampling frequency: 5400 Hz; left and right sideband frequencies: 730 Hz, 930 Hz; tap coefficient length: 121; window function: hanming.

[0200] Step 3: Design an IIR digital notch filter. On the power grid, the 50 Hz harmonic interference is an important interference source. After the above FIR band-pass filter, there are still harmonic interferences at 750 Hz, 850 Hz, and 950 Hz. When performing Fourier transform, if the amplitudes of these harmonic components are relatively large, they will still affect the Fourier transform results of the characteristic frequency points. Therefore, they must be filtered out. For the harmonic interferences at 750 Hz, 850 Hz, and 950 Hz, a second-order IIR digital notch filter is used in this embodiment to filter them out. The transfer function of the digital notch filter with a notch angular frequency of ω0 in the z-plane is: H(z) = (z 2-2cos(ω0 / Fs)z + 1) / (z 2 -(1 - μ)2cos(ω0 / Fs)z + (1 - μ) 2 ); where z is the complex variable in the z-transform. The transfer function can be implemented using a second-order linear difference equation with constant coefficients: -b1y(n - 1) - b2y(n - 2). Where y(n - 1) is the previous output value of the filter; a0 = 1, a1 = -2cos(ω0 / Fs); a2 = 1; b1 = -(1 - μ)2cos(ω0 / Fs); b2 = (1 - μ) 2 . In this embodiment, the notch depth parameter μ = 0.01, Fs = 5400Hz. The 3 notch filters are in a series structure, corresponding to frequencies of 750Hz, 850Hz, and 950Hz.

[0201] Step 4: Obtain the signal amplitude through the DFT method. Intercept the signal after the above-mentioned specific frequency point IIR filter and perform periodic extension into a periodic signal, which can be expanded using the Fourier series as shown in the following formula: I(k) = a_0 / 2 + ∑ n=1 ∞ [a n cos(2πf n kdt) + b n sin(2πf n kdt) ]; where a n and b n represent the amplitude of the signal with frequency f n , a n = 2 / T ∑ k=t_0 t_0+T I(k) cos(2πf n kdt); b n = 2 / T ∑ k=t_0 t_0+T I(k) sin(2πf n kdt). Where dt represents the sampling period of the signal, I(k) represents the k-th sample value, T is the length of the intercepted signal, a_0 is the DC component of the signal, n is the harmonic order, and k is the discrete time index. In this embodiment, T is the length of 5 cycle periods, and the 2 characteristic frequency points of concern are: F_1 = 783.30Hz and F_2 = 883.30Hz. After obtaining the in-phase and quadrature amplitudes of the 2 characteristic frequency points, the amplitude of the characteristic current signal can be calculated according to the following expression: A(F n ) = sqrt(a n 2 + b n 2 ); C = [A(F_1) + A(F_2)] / 2. In order to obtain an , b n value, the cos(2πf corresponding to two frequency points n kdt), sin(2πf n kdt) values can be calculated and stored as coefficients in an array, and then through the multiply-accumulate formula, a n , b n value is obtained. In this embodiment, a n , b n value is calculated once for each cycle, and each cycle signal sample value is updated each time.

[0202] Step 5: The amplitude sequence and the standard sequence are subjected to a correlation operation according to the following formula, where the amplitude sequence is c(n), and it slides one point each time it is operated, and thus the operation result r cr (k)=∑ n=0 N-1 c(n + k)∙r(n). Where r cr (k) is the correlation value between the amplitude sequence c(n) and the standard sequence r(n)

[0203] Step 6: Make a decision based on the correlation peak. The correlation peak is greater than the values of other non-feature signal parts. In this embodiment, if the correlation peak value is greater than twice the correlation value in the case of pure noise, it is identified that the current signal has a correlation with the standard signal, that is, a feature signal in the current is identified, and the identification of the current characteristics in the substation area is completed.

[0204] The present invention elaborates in detail a method for identifying the current characteristics in a substation area based on the correlation peak decision method. Through steps such as generating and injecting a rotation phase-encoded feature signal, collecting and preprocessing multi-point phase signals, constructing a self-organizing phase network and extracting features, making a decision on the correlation peak of the phase difference and making an identification decision, and adaptive optimization and parameter update, high-precision identification in a complex power grid environment is achieved. Compared with traditional methods, the present invention has obvious advantages in the following aspects: by using phase encoding instead of frequency features, the problems of harmonic interference and frequency overlap are effectively avoided; introducing topological invariant features enables the system to have strong resistance to local data loss or distortion; adopting an adaptive correlation peak decision mechanism improves the reliability and adaptability of the decision; establishing a complete adaptive optimization and parameter update mechanism enables the system to continuously learn and adapt to environmental changes.

[0205] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all belong to the protection scope of the present invention.

Claims

1. A method for identifying the current characteristics of a power distribution area based on the relevant peak judgment method, characterized in that Including: Collecting the grid synchronous phase information, generating a characteristic signal with phase encoding and injecting it into the power line to form an injected current signal; Collecting the line current signals containing the injected current signal at N key nodes in the substation area, preprocessing and extracting the phase information to obtain a phase time series; accordingly, extracting the current characteristic phase difference pattern and generating a topology-invariant feature vector; where N is a natural number greater than 0; Performing a correlation analysis on the topology-invariant feature vector and a preset standard feature template, calculating the correlation peak value and establishing a decision criterion to determine the identification result and the identification credibility index; Accordingly, dynamically adjusting the system parameters and updating the phase network structure and the decision threshold; The steps of generating the topology-invariant feature vector include: For the phase time series, calculating the phase difference between nodes to form a phase difference matrix; accordingly, applying a self-organizing network construction algorithm to establish the phase network topology structure of the substation area; Performing a time series analysis on the phase network topology structures obtained in a continuous predetermined number of cycles, extracting the dynamic change characteristics to obtain a network evolution feature matrix; Extracting the topology-invariant feature vector with environmental invariance from the network evolution feature matrix; The steps of applying a self-organizing network construction algorithm to establish the phase network topology structure of the substation area include: Performing a time dimension stability evaluation on each element of the phase difference matrix, calculating the standard deviation to generate a phase difference stability matrix; Based on the phase difference stability matrix, generating an adaptive threshold matrix; combining it with the phase difference stability matrix to establish a connection strength with continuous weights and generating a connection weight matrix; Based on historical data, evaluating the time stability of each connection in the connection weight matrix, calculating the connection reliability index to generate a connection reliability matrix; combining it with the connection weight matrix to construct a comprehensive weight matrix to form the final phase network topology structure.

2. The method according to claim 1, wherein The steps of forming the injected current signal include: Reading the grid voltage signal, obtaining the grid synchronous phase information through zero-crossing detection and phase-locked loop technology, and generating a rotating phase encoding sequence; Based on a preset fixed carrier frequency, mapping the rotating phase encoding sequence onto the carrier signal to generate a phase-modulated carrier signal; segmenting it at pseudo-random time intervals to form an intermittent pulse sequence and injecting it into the power line to generate an injected current signal.

3. The method according to claim 2, wherein The steps of generating the rotating phase encoding sequence include: Reading the grid synchronous phase information, constructing a phase offset basic unit based on the Gray code principle to form a coding base sequence; Applying the Hamming code principle to the coding base sequence to increase the coding redundancy and generate a redundant coding sequence; analyzing the phase jump characteristics therein, and adjusting the phase value through a phase mapping function to generate an optimized redundant sequence; Converting the optimized redundant sequence into a differential phase coding form and combining it with a preset absolute phase reference point to construct a rotating phase encoding sequence.

4. The method according to claim 2, wherein The steps of forming the intermittent pulse sequence include: Obtaining a pseudo-random sequence and mapping it to the time interval domain to generate a pseudo-random time interval sequence; Analyzing the current characteristics and interference spectrum characteristics of the substation area, determining the optimal pulse duration based on the principle of maximizing the signal-to-noise ratio; accordingly, optimizing the pulse edge to generate an optimized window function; Obtain the line impedance characteristics and background noise level, dynamically adjust the pulse energy, and generate an adaptive amplitude sequence; Combine the phase-modulated carrier signal with the pseudo-random time interval sequence, the adaptive amplitude sequence, and the optimized window function to generate an intermittent pulse sequence.

5. The method according to claim 1, wherein The steps to obtain the phase time series include: Select a predetermined number of key nodes in the substation area, use high-precision current sensors to collect line current signals, and set the sampling frequency to a value higher than the carrier frequency; Segment the line current signal into synchronous signal segments of the same length, filter it through a band-pass filter with a center frequency equal to the carrier frequency to obtain a filtered signal; Apply the Hilbert transform to the filtered signal to construct an analytic signal, extract the instantaneous phase information, and obtain the phase time series of each node.

6. The method according to claim 5, wherein The steps to extract the instantaneous phase information and obtain the phase time series of each node include: Perform zero-mean processing on the filtered signal and remove abnormal spikes to generate a preprocessed signal; calculate its orthogonal component using the Hilbert transform to generate the Hilbert transform result; Based on the preprocessed signal and its Hilbert transform result, construct an analytic signal and perform instantaneous phase extraction and phase unwrapping processing to obtain an unwrapped phase sequence; Perform trend analysis on the unwrapped phase sequence, extract and remove the linear phase growth component, and perform phase normalization to obtain the final phase time series.

7. The method according to claim 1, wherein The steps to determine the identification result and the identification confidence index include: Through the injection experiment in the ideal environment, obtain the topological invariant feature vector under standard conditions as the standard feature template; perform multi-dimensional correlation analysis with the topological invariant feature vector, introduce the topological weighted correlation function, and obtain the correlation peak sequence; Based on the pre-stored environmental noise level and historical recognition results, construct an adaptive decision threshold function, compare the correlation peak sequence with the threshold, and output the decision result; Based on the adaptive decision threshold function, calculate the decision confidence index; combine it with the consistency of the continuous predetermined round of decision results to generate the final identification result and the identification confidence index.

8. The method according to claim 7, characterized in that, The steps to obtain the correlation peak sequence include: Calculate the preliminary similarity between the topological invariant feature vector and each standard feature template to generate a preliminary similarity vector; evaluate the discrimination ability of each feature dimension based on this to generate a feature weight vector; Based on the principle of topological data analysis, construct a topological importance weighted function, assign weights to the persistent homology feature, the Betti number feature, and the spectral feature to generate a topological importance function; Combine the feature weight vector with the topological importance function to construct the final weight, calculate the topological weighted correlation function, determine the maximum correlation value, and generate the correlation peak sequence.

Images