Transformer area current characteristic identification method based on correlation peak judgment method

Abstract

The invention discloses a transformer area current characteristic identification method based on a correlation peak judgment method, and the method comprises the steps: collecting the synchronous phase information of a power grid, and generating and injecting a unique phase coding characteristic signal; current signals are collected at key nodes of the transformer area, and phase information is extracted; constructing a dynamic self-adaptive phase relation network, and extracting a feature vector with topology invariance; correlation analysis is carried out on the extracted feature vectors and a standard template, and an identification result is determined through adaptive correlation peak judgment; and dynamically adjusting system parameters according to an identification result. The method of combining phase coding and topology invariant features is adopted, and the method has the advantages of being high in anti-interference capacity, good in environment adaptability and high in recognition precision.

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] The identification of 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 early warning, load identification, unauthorized power consumption detection, and power quality assessment, which is of great 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 feature signals in specific frequency bands 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 significantly decreases 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 render the preset models ineffective; 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 current characteristic identification technology of distribution transformer areas 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 feature signal with phase encoding and inject it into the power line to form an injected current signal;

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

[0009] Perform correlation analysis on the topological invariant feature vector and the 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 and update the phase network structure and decision threshold.

[0010] Beneficial effects: The present invention can effectively avoid the problems of harmonic interference and frequency overlap, 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 flowchart of the steps of a method for identifying the current characteristics in the substation area based on the correlation peak decision method provided by the embodiment of the present application.

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

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

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

[0015] Figure 5 It is a flowchart of the steps for determining the identification result and the identification credibility index provided by the embodiment of the present application. Detailed Embodiment

[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 of 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 the purpose of clearly demonstrating the step - by - step process 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 relevant - 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 power 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 power 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. At a predetermined key node within the power distribution area, collect the line current signal containing the injected current signal, pre - process the line current signal 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 carried out to extract clearer and more meaningful content.

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

[0024] Specifically, the dynamic and adaptive phase - relationship network will automatically optimize according to the changes in 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 current characteristics. The topologically invariant feature vector is stable and can remain consistent even when the network structure changes, thus being used for subsequent analysis or applications.

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

[0026] Specifically, the relevant peak value represents the level of similarity. The peak value represents the matching degree between the feature vector and the standard template. The higher the peak value, the more similar the two are. The decision criterion can ensure the scientific nature 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 identification of current characteristics, so that the system can better adapt to the current conditions. Update the phase network structure so that it can capture the characteristics of the current more accurately. 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 feature vector 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 the phase-locked loop technology to track the changes of the grid frequency and phase 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 rotational phase encoding sequence Φ e . Select a fixed carrier frequency f k (usually 700 Hz), map the rotational 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, which 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 , and the duration of each segment is 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, generating a line injection current signal I k (t).

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

[0037] Read the grid synchronization phase information, construct a basic unit of phase offset 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, and 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 rotational phase encoding sequence, ensuring that the rotational phase encoding sequence has the characteristics of forward error correction ability and a single sharp peak autocorrelation function.

[0041] Specifically, read the grid synchronization phase reference φ 0 , design a basic unit of phase offset based on the Gray code principle, including four basic phase states of 0°, 90°, 180°, and 270°, and 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: P 1 = B 1 Θ B 2 Θ B 4 ; P 2 = B 1 Θ B 3 Θ B 4 ; P 3 = B 2 Θ B 3 ΘB 4 ; where B 1 to B 4 are the 4-bit values of the base sequence, and P 1 to P 3 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 capability; the periodic autocorrelation function has a single sharp peak; the non-periodic autocorrelation sidelobe peak does not exceed 0.3 times 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, and 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 for controlling the duration of each pulse;

[0046] Based on the optimal pulse duration, optimize the pulse edge to generate an optimized window function for smoothing the pulse edge;

[0047] Obtain the line impedance characteristics and the background noise level, dynamically adjust the pulse energy, and 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 , with a value range of [0.5s, 2.5s]. Analyze the current characteristics of the substation area and the interference spectrum characteristics, and based on the principle of maximizing the signal-to-noise ratio, 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 the background noise level, dynamically adjust the pulse energy to balance the detection performance and the interference to the power grid: A(t) = A 0 ·[1 + α·sin(2πf 0t)]·[1 +β·N(t)], where A 0 is the reference amplitude, α is the power frequency modulation coefficient (typical value 0.1), f 0 is the power frequency, β is the noise response coefficient, N(t) is the normalized noise index, to generate the adaptive amplitude sequence A(t). To reduce spectral leakage and adjacent frequency point interference, optimize the design of the pulse edge, and apply the improved Hamming window function: w(t) = 0.54 - 0.46·cos(2πt / T p ) +γ·t / T p ·(1 - t / T p ), where γ is the adaptive parameter, dynamically adjusted according to the power grid noise characteristics, to generate the optimized window function w(t). Combining the phase modulation 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), construct the complete discontinuous pulse sequence P k (t): 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 .

[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 distribution area, and use high-precision current sensors to collect line current signals, and set the sampling frequency to a value higher than the carrier frequency;

[0052] S22. Segment the line current signals 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 the 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 power distribution area (such as the low-voltage side of the transformer, important branch lines, end-user access points, etc.), and 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. 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) of the same length, with each segment corresponding to a complete injection pulse period to ensure the time alignment of signals at different nodes. For the synchronous signal segment I i，j (t), design a band-pass filter with a center frequency of f k 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). Apply the Hilbert transform to the filtered signal I i，j '(t) to construct an analytic signal, extract the instantaneous phase information, and obtain 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 and obtaining the phase time series of each node include:

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

[0058] Calculate the orthogonal component of the preprocessed signal using the Hilbert transform to generate the Hilbert transform result;

[0059] Based on the preprocessed 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 and obtain 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, zero-mean 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 to generate the preprocessed signal I i，j ''(t). Calculate the orthogonal component of the preprocessed signal I i，j ''(t) using the improved Hilbert transform to generate 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 and τ is the integration variable; to reduce the computational complexity, implement it indirectly using FFT: calculate the FFT of the preprocessed signal I i，j ''(t) to obtain F(ω); construct the frequency response of the Hilbert transformer H(ω): 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 raw phase sequence θ i，j _raw(t). Due to the range limitation of the atan2 function in [-π, π], there are phase jumps in the raw 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 output is the unwrapped phase sequence θ i，j _unwrap(t). For the unwrapped phase sequence θ i，j _unwrap(t), perform trend analysis, 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 k is 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, map it to the interval [0, 2π), and 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 structure of the substation area;

[0067] S33. Perform time series analysis on the phase network topology structures obtained in consecutive predetermined 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 , a new self-organizing network construction algorithm is applied to establish the phase network topology G of the substation area j . 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λ p , no connection is established between nodes i and k; based on the connection strength, a phase network topology G represented by a weighted adjacency matrix is constructed j .

[0070] Perform time series analysis on the phase network topologies G j (j = 1, 2,..., M) obtained in consecutive multiple periods, extract their dynamic change characteristics, and obtain the 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: extract the basic topological parameters for the phase network topology G at each time point j , 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 the 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 the parameter stability vector S_P. Apply the improved Louvain algorithm to perform community division on the phase network topology G j : In the initial state, each node is used as 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, and ci Let \(C_i\) be the community to which node \(i\) belongs, \(\delta\) be the Kronecker function, and \(W'(i, j)\) be the edge weight between nodes \(i\) and \(j\). The evolutionary relationship of communities at different time points is tracked through a matching algorithm, and a community structure evolution sequence \(C\_comm(j)\) is generated. Based on a time sliding window, dynamic characteristic indicators of the network structure are calculated: Structure change rate: a normalized measure of the topological difference between adjacent time points; Community stability: the proportion of community members remaining unchanged; Key node migration: the proportion of nodes with a significant change in centrality ranking; Link volatility: the proportion of newly added / removed links, which are combined to form a network dynamic characteristic vector \(D\). n et. 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 a 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 (an algebraic representation of the network cycle structure); generating topological spectral features (a spectral representation of the network topological structure). Combine the above features to form the topological invariant feature vector \(V\). t as the core basis for feature identification. t

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

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

[0074] Based on the phase difference stability matrix, generate an adaptive threshold matrix;

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

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

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

[0078] Specifically, perform a time - dimension stability assessment on each element of the phase - difference matrix D j and calculate the standard deviation: σ_D(i ， k) = std(D j (i, k, t)), t ∈ [t 1 , t 2 , where [t 1 , t 2 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 α, β, γ are parameters automatically adjusted according to the current environmental noise, and generate the 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, to generate the connection weight matrix W. Based on historical data, evaluate the temporal 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 W_hist is the historical connection weight record, and var and mean represent variance and mean respectively, to generate the connection reliability matrix R. 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); Apply the minimum spanning tree algorithm in graph theory to ensure network connectivity while removing possible noisy connections. Specific steps: Construct a fully connected 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 , which is 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 multi-scale simplicial complexes based on the network evolution feature matrix, and form 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 to generate a Betti number sequence vector;

[0085] Calculate the topological spectral features based on the network evolution feature matrix, including feature spectral energy, spectral gap, and spectral moment, to 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 topologically 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}, for each threshold ε i : Calculate the pairwise distance matrix D_E of the feature vectors; When d(x, y) ≤ ε iWhen a connection is established between x and y; when all edges (x, y, z) exist, a 2-simplex (triangle) is formed; higher-order simplices are constructed in this way; a sequence of simplicial complexes K(ε i ) is obtained. Apply the persistent homology algorithm to the sequence of simplicial complexes K(ε i ): For each threshold ε i , calculate the homology groups H 0 (K(ε i )), H 1 (K(ε i )), H 2 (K(ε i )), where H 0 corresponds to the number of connected components; H 1 corresponds to the number of one-dimensional holes; H 2 corresponds to the number of two-dimensional cavities. Track the "birth" and "death" of homology features to generate the persistence intervals 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: Represent each persistence interval in the two-dimensional plane with (ε βi rth , ε death ) as coordinates; Calculate the statistical features of the persistence diagram: the number of persistence intervals, the maximum persistence duration, the persistence 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 persistent 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, and H* represents the homology group, to generate the sequence vector β of Betti numbers. Calculate the topological spectral features based on the network evolution feature matrix E t : 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 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 = [w 1 ·PD; w 2 ·β; w 3 ·S_topo], where w 1 、w 2 、w 3 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 characteristics 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 of determining the identification result and the identification credibility index include:

[0090] S41. Through the injection experiment in the 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 best matching template index;

[0093] S44. Based on the adaptive decision threshold function, calculate the decision confidence index; 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 s , 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 lies in introducing the topological weighted correlation function: R(k, τ) = Σ i w(i)·V t (i + τ)·Vs k (i), where w(i) is the topological importance weight, τ is the relative time shift, and the maximum correlation value is calculated for each template k: R max (k) = max{R(k, τ)}, obtaining the correlation peak sequence R max . Perform adaptive correlation peak decision. Based on the environmental noise level and historical recognition results, design an adaptive decision threshold function η(σ), where σ represents the current environmental noise index. When the condition: max{R max (k)} > η(σ)·R 0 is satisfied, it is determined that the characteristic signal exists, where R 0 is the reference correlation value in a 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)} - η(σ)·R 0 / [η(σ)·R 0 ; Based on the confidence C and the consistency of consecutive multi-round decisions, generate the final identification result Y and the identification credibility index Q.

[0095] Among them, the specific process of performing 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)|); Degree of harmonic interference: H_dist = ∑ i P(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 environmental 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 an environment-threshold-performance dataset {(σ_hist, η_hist, Perf_hist)}; Use the K-nearest neighbor algorithm to find the K historical environmental conditions most similar to the current environment σ; Extract the average optimal threshold of these K environments: η_KNN = ∑ i=1 K w i ·η_hist, i / ∑ i=1 Kw i , where w i is the weight based on environmental similarity; generate the historical reference threshold η_KNN. Based on the environmental noise vector σ and the historical reference threshold η_KNN, construct an adaptive threshold function: η(σ) = η_0 + α·exp(β·||σ||) + γ·η_KNN, where η_0 is the base threshold, and α, β, and γ are system parameters obtained by offline optimization of historical data, and generate the adaptive threshold η_cur of the current environment. Based on the adaptive threshold η_cur, the correlation peak sequence R max and the peak significance vector SIG, construct a multi-dimensional decision criterion: D_Peak(k) = {1, if R max (k) > η_cur·R 0 && SNR(k) > SNR_min && Sharp(k) > Sharp_min; 0, otherwise}; where R 0 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 feature signal, and the decision result D = 0; if there exists D_Peak(k) = 1, it is determined that there is a feature 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 reference correlation value R 0 , calculate the decision confidence: C_single=[max{R max (k)}-η_cur·R 0 / [η_cur·R 0; When the judgment result D = 0, the confidence level is negative, indicating the degree of certainty tending to no feature signal, and a single-round confidence level C_single is generated. According to the single-round confidence level C_single and the peak significance vector SIG, the weighted confidence level is calculated: C_weight = C_single·[α·SNR(k_best)+β·Sharp(k_best)] / (α + β), where α and β are weight coefficients, and the weighted confidence level C_weight is generated. Store the judgment result D, the best matching template index k_best, and the weighted confidence level C_weight of the current round into the judgment history cache, retain the judgment results of the most recent M rounds, and form a judgment history sequence H_D: H_D ={(D 1 , k_best 1 , C_weight 1 ), (D 2 , k_best 2 , C_weight 2 ),..., (D m , k_best m , C_weight m )}. Analyze the consistency in the judgment history sequence H_D: Calculate the judgment result consistency rate: Con_D = count(D i = D 1 ) / M; Calculate the best template consistency rate: Con_k = count(k_best i = k_best 1 ) / count(D i = 1); Calculate the mean and variance of the confidence level: Mean_C = mean(C_weight i ); Var_C = var(C_weight i ); Generate a judgment consistency vector Con. Based on the judgment consistency vector Con and the current judgment result D, perform multi-round fusion decision-making: Calculate the fusion confidence level: C_fusion = Mean_C·(1 + γ·Con_D - Δ·sqrt(Var_C)), where γ and Δ are weighted coefficients; Set a fusion threshold λ_fusion; When C_fusion > λ_fusion, output the identification result Y = 1, otherwise output the identification result Y = 0; 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 initial similarity between the topological invariant feature vector and each standard feature template to generate an initial similarity vector;

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

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

[0101] Combine the feature weight vector and the topological importance function to construct the final weight, calculate the topological 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 initial 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, to generate the initial 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 -is the statistic of misidentified samples; normalize the discrimination ability into weights: w(i) = d i / ∑ j d j , and 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, calculate the maximum correlation value for each template k: R max (k) = max{R(k, τ)}, τ ∈ [-τ max , τ max , and generate the correlation peak sequence R max . Perform peak significance evaluation on the correlation peak sequence R max : 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 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 an 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-making under different environmental conditions in historical data to 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 realize the 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 substation area power grid, 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 to 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-making under different environmental conditions in historical data, optimize the decision threshold function η(σ) so that it can more accurately adapt to various noise environments. The specific optimization process uses the Bayesian optimization method to balance the false alarm rate and the missed detection rate to 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 , realize the incremental update of the template library, and improve the environmental adaptability of the system. Regularly evaluate the overall performance of the system, including indicators such as recognition accuracy, response time, and environmental adaptability, to 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 substation area based on the correlation peak decision method realizes high-precision identification of the current characteristics of the substation area in a complex power grid environment through technical means combining phase encoding and topological invariants. The implementation process of this method will be described in detail below with specific numerical examples.

[0112] Step 1. Generation and injection of the rotated phase-encoded 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 technology, and the phase-locked loop technology 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 technology, 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. Design of phase coding redundancy. 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: P 1 = B 1 Θ B 2 Θ B 4 ; where P 1 is the first parity bit; B 1 , B 2 , B 4 are the values of the 1st, 2nd, and 4th bits of the base sequence respectively; Θ represents the exclusive OR operation. P 2 = B 1 Θ B 3 Θ B 4 ; where P 2 is the second parity bit; B 1 , B 3 , B 4 are the values of the 1st, 3rd, and 4th bits of the base sequence respectively; Θ represents the exclusive OR operation. P 3 =B 2 Θ B 3 Θ B 4 ; where P 3 is the third parity bit; B 2 , B 3 , B 4 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, map the phase values to 2-bit encoding: 0° = 00, 90° = 01, 180° = 11, 270° = 10, then: B 1 = 00, B 2 = 01, B 3 = 10, B 4 = 11; Calculate the parity bit: P 1 = 00 Θ 01 Θ 11 = 10 = 270°; P 2 = 00 Θ 10 Θ 11 = 01 = 90°; P 3 = 01 Θ 10 Θ 11 = 00 = 0°; Generate the redundant coding sequence Φ r = [0°, 90°, 270°, 180°, 270°, 90°, 0°], where Φ r represents the phase coding sequence including the redundant parity bit.

[0118] 1.2.3. Phase transition optimization. Analyze the phase transition characteristics in the redundant coding 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 coding sequence; i represents the sequence index, and the value range is [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 encoding sequence after jump 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 at adjacent time points; Φ r '(i) represents the i-th phase value of the optimized redundant sequence; i represents the sequence index, and the value range is [1, n - 1], where 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 encoding sequence.

[0122] 1.2.5. Final encoding sequence generation. Combine the differential phase sequence with the absolute phase reference point to construct the final rotational phase encoding sequence: Φ e (1) = Φ r '(1) = 0°; where Φ e (1) represents the first phase value of the final encoding 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 encoding sequence; Φ e (i) represents the i-th phase value of the final encoding 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 encoded 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 rotation phase encoding sequence Φ e = [0°, 90°, 270°, 180°, 270°, 90°, 0°], where Φ e represents the final rotation phase encoding sequence. Examine the characteristics of the final encoded sequence: the minimum phase difference is not less than 45°: the minimum phase difference is 90°, meeting the condition; having forward error correction ability: achieved through Hamming code check bits, capable of correcting 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, set to 1.5A in this embodiment; f k represents the carrier frequency; t represents time; φ 0 represents the power grid synchronous phase reference; Φ e (i) represents the i-th phase value of the final encoded sequence; t ∈ [t i , t i+1 indicates that t is within the i-th time interval. In this embodiment, each phase duration is 20 ms. For t = 0 - 20 ms, the phase is 0°, and the signal expression is: Sk (t) = 1.5·sin(2π·700·t + 32.7° + 0°) = 1.5·sin(2π·700·t + 32.7°); For t = 20 - 40 ms, 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. Realization of discontinuous pulse injection.

[0126] 1.4.1. Generation of pseudo - random time intervals. 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 i represents 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 X 0 = 12345, and calculate the pseudo - random sequence: X 1 = (1664525·12345 + 1013904223) mod 2 32 = 1666278168; X 2 = (1664525·1666278168 + 1013904223) mod 2 32 = 1253251530; X 3 = (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.276 s; 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...。Get 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, measure the signal-to-noise ratio at different pulse durations, and select the optimal pulse duration: 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, determine the optimal pulse duration T p = 100ms, where T p represents the pulse duration.

[0128] 1.4.3. Adaptive adjustment of pulse energy. Dynamically adjust the pulse energy according to the line impedance characteristics and background noise level. The adjustment formula is: A(t) = A 0 ·[1 + α·sin(2πf 0 t)]·[1 + β·N(t)]; where A(t) represents the adaptive amplitude; A 0 represents the reference amplitude, set to 1.5A; α represents the power frequency modulation coefficient, set to 0.1; f 0It 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. Assume that 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 an adaptive amplitude sequence A(t), where A(t) represents the adaptive amplitude sequence.

[0129] 1.4.4. Pulse waveform optimization design. 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 its value range is [0, T p ; T p represents 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, a complete optimized window function w(t) is generated, where w(t) represents the optimized window function.

[0130] 1.4.5. Discontinuous pulse sequence generation. 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 is constructed: P k (t) = A(t)·w(t - tj )·sin(2πf k t + φ 0 + Φ e (j)); where P k (t) represents an intermittent pulse sequence; t ∈ [t j , t j +T p represents that t is within the j-th pulse time; A(t) represents an adaptive amplitude sequence; w(t - t j ) represents a 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 the time interval between the end of the j-th pulse and the start of the j+1-th pulse.

[0131] Pulse starting time calculation: t j+1 = t j + T p + T r (j); where t j+1 represents the starting time of the j+1-th pulse; t j represents the starting time of the j-th pulse; T p represents the pulse duration; T r (j) represents the j-th pseudo-random time interval. Assume t 1 = 0, then: t 2 = 0 + 0.1 + 1.276 = 1.376s; t 3 = 1.376 + 0.1 + 1.084 = 2.56s; t 4 = 2.56 + 0.1 + 2.258 = 4.918s... The intermittent 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. Acquisition of key node signals. Select N = 5 key nodes in the power distribution 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. Segmentation and synchronization of current signals. Based on the zero-phase mark of the power grid, segment the collected line current signals I i (t) into synchronized 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 the 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 synchronized 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 as f k / BW = 700 / 20 = 35. Apply this filter to each synchronized signal segment to obtain the filtered signal.

[0136] 2.4. Phase information extraction.

[0137] 2.4.1. Preprocessing of sampling signals. Perform zero-mean processing on the filtered signal I i，j '(t) 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 value of the filtered signal. At the same time, the local polynomial fitting method is applied to remove abnormal spikes and generate the preprocessed signal I i，j ''(t).

[0138] 2.4.2. Analytical signal construction. The Hilbert transform is indirectly implemented using FFT to 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 of 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 analytical signal: z(t) = I i，j ''(t) + j·H{I i，j ''(t)}; where z(t) represents the analytical signal; j represents the imaginary unit; I i，j ''(t) represents the preprocessed signal; H{I i，j ''(t)} represents the Hilbert transform result of the preprocessed signal. Extract the instantaneous phase of the analytical 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 Hilbert transform result; I i，j ''(t) represents the preprocessed signal.

[0140] 2.4.4. Phase unwrapping. Due to the range limitation of the atan2 function within [-π, π], the original phase sequence θ i，j _raw(t) has phase jumps and needs to be phase unwrapped. 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 raw phase value of the k-th sampling point; θ i，j _raw(k + 1) represents the raw 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. For the unwrapped phase sequence θ i，j _unwrap(t), perform trend analysis, 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, perform phase normalization, mapping to the [0, 2π) interval: θ 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 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 represented 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 in the jth time segment; θ k，j (t) represents the phase time series of node k in the jth 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 phase network topology.

[0145] 3.2.1. Evaluation of phase difference stability. Conduct a stability evaluation of each element of the phase difference matrix D j in the time dimension, 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 ∈ [t 1 , t 2 represents the analysis time window. For example, for nodes 1 and 2, within the time window t ∈ [0, 100ms], take the phase differences at 10 time points to calculate the standard deviation: D j (1, 2, t 1 ) = 78.3°; D j (1, 2, t 2 ) = 77.9°;...; D j (1, 2, t 10 ) = 79.1°; σ_D(1, 2) = 0.5°. Similarly, calculate the phase difference stability of all node pairs 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. Multi-layer weighted connection establishment. 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 decay coefficient, set to 100. W(i, k) = 0; where |Dj |(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 for 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 Node 2 and Node 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 for all node pairs to generate the connection reliability matrix R.

[0150] 3.2.5. Topology optimization. Combining the connection weight matrix W and the connection reliability matrix R, 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 the W' value; select the edge with the largest weight to add to the spanning tree if no cycle 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 , which is represented by an adjacency matrix.

[0151] 3.3. Analysis of the time evolution of the phase network.

[0152] 3.3.1. Extraction of network structure parameters. For the phase network topology structure G at each time point j extract the basic topology 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 tightness of the connections between the neighbors of the node; 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 centrality of each node is: [0.2, 0.5, 0.1, 0.4, 0.1]. Spectral eigenvalue 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 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)); 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, which takes the value 1 when c i =c j , and 0 otherwise. 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. Track the evolution relationship of communities at different time points through a matching algorithm to generate the community structure evolution sequence C_comm(j).

[0155] 3.3.4. Network Dynamic Feature Extraction. Based on a time sliding window, calculate the dynamic feature 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 that remain 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 to form the network dynamic feature vector D n et.

[0156] 3.3.5. Network Evolution Feature Matrix Construction. 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 feature vector D n et to construct a complete network evolution feature matrix E t : E t = [P(j); S_P; C_comm(j); D n et]; where Et Denote the network evolution feature matrix; each column corresponds to a time point, and each row corresponds to a feature dimension.

[0157] 3.4. Topological invariant feature extraction.

[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 between feature vectors. For example, for feature vectors v 1 and v 2 , the Euclidean distance is 0.25. When d(x, y) ≤ ε i , establish a connection between x and y. For example, at the threshold ε 2 = 0.2, no connection is established between v 1 and v 3 because d(v 1 , v 3 ) = 0.32 > 0.2. When all edges (x, y, z) exist, form a 2-simplex (triangle); for example, at the threshold ε 3 = 0.3, a triangle is formed between v 1 , v 2 , v 4 because d(v 1 , v 2 ) = 0.25; d(v 1 , v 4 ) = 0.18; d(v 2 , v 4 ) = 0.22, all of which are less than 0.3. And so on to construct higher-order simplices, obtaining 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 ): Calculate the homology groups H i (K(ε 0 ))、H i (K(ε 1 ))、H i (K(ε 2 )) for each threshold ε i ; H 0 represents the number of connected components; H 1Represents the number of one-dimensional holes; H 2 Represents the number of two-dimensional cavities. For example, at the threshold ε 1 = 0.1, H 0 = 5 (5 independent connected components), H 1 = 0 (no holes), H 2 = 0 (no cavities); at the threshold ε 3 = 0.3, H 0 = 1 (all points are connected), H 1 = 2 (there are 2 holes), H 2 = 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. Persistence Diagram Construction and Feature Extraction. Construct a persistence diagram based on the persistent homology feature set PH: Represent each persistence interval on a two-dimensional plane with coordinates (ε βi rth , ε death ); for example, the point (0.2, 0.4) represents the persistence interval of the hole mentioned above. Calculate the statistical features of the persistence diagram: Number of persistence intervals: for example, there are 10 persistence intervals in total; Maximum persistence duration: for example, max(pers) = 0.3; Persistence duration distribution: for example, [0.1, 0.1, 0.2, 0.2, 0.3, 0.1, 0.2, 0.1, 0.1, 0.2]; Persistence interval density: for example, the persistence interval density 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. Betti Number Sequence Extraction. 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. Topological spectrum feature calculation. 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}; for example, the eigenvalue spectrum is {0, 0.8, 1.5, 2.3, 3.4}. Extract the spectrum features: Spectrum feature energy: ∑λ i 2 = 0 2 + 0.8 2 + 1.5 2 + 2.3 2 + 3.4 2 = 19.74; Spectral gap: λ 2 - λ 1 = 0.8 - 0 = 0.8; Spectral moment: ∑λ i k / n, for example, the first-order 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 = [w 1 ·PD; w 2 ·β; w 3 ·S_topo]; where w 1 、w 2 、w 3 represent feature weights, which are adaptively adjusted according to the current environment and historical recognition results. In this embodiment, according to the historical recognition performance, set w 1 = 0.5, w 2 = 0.3, w 3 = 0.2. The feature vector design ensures: being robust to topological structure deformation; being insensitive to noise and outliers; retaining the essential topological characteristics of the network; being able to effectively distinguish the topological patterns of different feature signals.

[0164] Step Four. Phase difference related 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 , which is used as the standard feature template for subsequent comparison. To adapt to different environmental conditions, establish multiple templates 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 Vt With template library V s k , 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 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 d i represents the discrimination ability of the i-th feature dimension; μ i + represents the mean value of correctly recognized samples on dimension i; μ i - represents the mean value of misrecognized samples on dimension i; σ i + represents the standard deviation of correctly recognized samples on dimension i; σ i - represents the standard deviation of misrecognized samples on dimension i. For example, for the 1st 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 the topological importance weighting function: for persistent homology features, the weight is proportional to the duration. For example, the weight of a feature with a duration of 0.3 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 β 1It remains unchanged within the threshold interval [0.3, 0.4], and 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, then the weight is 0.85. The topological importance function w_topo(i) is comprehensively generated.

[0171] 4.2.4. Calculation of multi-dimensional 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 correlation peak sequence R max = [0.92, 0.65, 0.78, 0.61, 0.70], where R max represents the correlation peak sequence.

[0172] 4.2.5. Peak significance evaluation. For the relevant peak sequence R max Perform 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, which includes 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(|Ii，j '(t)| 2 ) where \(P_{avg}\) represents the average power of the signal; \(mean\) represents the mean function; \(I\ i，j '(t)\) represents the filtered signal; \(|\cdot|\) represents the absolute value. For example, \(P_{avg}=1.25\). Signal volatility: \(F_{sig}=\frac{std(I i，j '(t))}{mean(|I i，j '(t)|)}\); where \(F_{sig}\) represents signal volatility; \(std\) represents the standard deviation function; \(mean\) represents the mean function; \(I\ i，j '(t)\) represents the filtered signal; \(|\cdot|\) represents the absolute value. For example, \(F_{sig}=0.35\). Degree of harmonic interference: \(H_{dist}=\sum_{}^{} i \frac{P(f_{harm},i)}{P_{total}}\); where \(H_{dist}\) represents the degree of harmonic interference; \(P(f_{harm},i)\) represents the power at the \(i\)-th harmonic frequency point; \(P_{total}\) represents the total power; \(\sum_{}^{} i \) represents the sum over all harmonic frequency points. For example, \(H_{dist}=0.15\). Phase stability: \(S_{phase}=std(\theta i，j (t))\); where \(S_{phase}\) represents phase stability; \(std\) represents the standard deviation function; \(\theta i，j (t)\) represents the phase time series. For example, \(S_{phase}=0.08\). Generate the environmental noise vector \(\sigma = [P_{avg},F_{sig},H_{dist},S_{phase}]=[1.25,0.35,0.15,0.08]\), where \(\sigma\) represents the environmental noise vector.

[0175] 4.3.2, Historical threshold data analysis. Read the historical decision data, extract the optimal thresholds and corresponding performance metrics under different environmental conditions: establish the environmental-threshold-performance dataset \(\{(\sigma_{hist},\eta_{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 \(\sigma\). For example, the three most similar historical environments and their corresponding optimal thresholds are: \((\sigma_{hist} 1 ,\eta_{hist} 1 ) = ([1.20,0.32,0.14,0.07],0.75)\); \((\sigma_{hist} 2 ,\eta_{hist} 2 ) = ([1.28,0.38,0.17,0.09],0.80)\); \((\sigma_{hist} 3 ,\eta_{hist} 3 ) = ([1.18,0.33,0.13,0.06],0.73)\).

[0176] Calculate the environmental similarity weight: 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 3 historical environments: d(σ, σ_hist 1 ) = 0.08, w 1 = 12.5; d(σ, σ_hist 2 ) = 0.05, w 2 = 20.0; d(σ, σ_hist 3 ) = 0.10, w 3 = 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 neighbors; 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 summation over 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 the 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 base 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 of the current environment = 0.792, where η_cur represents the adaptive threshold of the current environment.

[0179] 4.3.4. Peak decision criterion design. Based on the adaptive threshold η_cur, the correlation peak sequence R max and the peak significance vector SIG, construct a multi-dimensional decision criterion: D_Peak(k) = 1; where R max (k) > η_cur·R 0 && SNR(k) > SNR_min && Sharp(k) > Sharp_min indicates 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; R 0 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 indicates that the k-th template does not meet the decision condition.

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

[0181] 4.3.5. Comprehensive decision of multiple templates. 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}; 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 satisfies the decision condition. In this example, the template indices that satisfy the condition 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 decision result D = 1 and the best matching template index k_best = 1.

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

[0183] 4.4.1. Calculation of the confidence of a single-round decision. Based on the correlation peak sequence R max, adaptive threshold η_cur and pure noise reference correlation value R 0 , calculate the decision confidence: C_single = [max{R max (k)} - η_cur·R 0 ] / [η_cur·R 0 ]; 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; R 0 Indicates the baseline correlation value in a 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 a negative value, indicating the degree of confidence that there is no characteristic signal. Generate a 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 saliency vector SIG, the weighted confidence is calculated: C_weight = C_single · [α · SNR (k_best) + β · Sharp (k_best)] / (α + β); where C_weight represents weighted confidence; C_single represents single-round confidence; α represents the signal-to-noise ratio weight, which is set to 0.7; β represents the sharpness weight, which is set to 0.3; SNR (k_best) represents the peak signal-to-noise ratio of the best matching template; Sharp (k_best) represents the peak sharpness of the best matching 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 weighted confidence C_weight = 23.55.

[0185] 4.4.3. Multi-round decision result cache. The current round decision result D, the best matching template index k_best and the weighted confidence C_weight are stored in the decision history cache, and the decision results of the most recent M = 5 rounds are retained to form a decision history sequence H_D: H_D = {(D 1 , k_best 1 , C_weight 1 ) ,(D 2 , k_best 2 , C_weight 2),..., (D m , k_best m , C_weight m )}; where H_D represents the decision history sequence; D i represents the decision 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 of the i-th round. For example, the decision history sequence after adding the result of the current round 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. Decision Consistency Analysis. Analyze the consistency in the decision history sequence H_D: Calculate the decision result consistency rate: Con_D = count(D i = D 1 ) / M; where Con_D represents the decision result consistency rate; count(D i = D 1 ) represents the number of historical decisions that are consistent with the latest decision 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_best 1 ) / count(D i = 1); where Con_k represents the best template consistency rate; count(k_best i = k_best 1 ) represents the number of historical best templates that are consistent with the latest best template; count(D i = 1) represents the number of historical decisions with a decision 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 of 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 iRepresents the weighted confidence of the i-th round. For example, Var_C = 171.25. Generate the decision consistency vector Con = [Con_D, Con_k, Mean_C, Var_C] = [0.8, 0.75, 13.86, 171.25], where Con represents the decision 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-making: Calculate the fusion confidence: C_fusion = Mean_C·(1 + γ·Con_D - Δ·sqrt(Var_C)); where C_fusion represents the fusion confidence; Mean_C represents the mean of the weighted confidence; γ represents the decision consistency weight, set to 0.5; Con_D represents the decision result agreement rate; Δ 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·(1 + 0.5·0.8 - 0.01·sqrt(171.25)) = 13.86·(1 + 0.4 - 0.01·13.09) = 13.86·(1.4 - 0.1309) = 13.86·1.2691 = 17.59. Set the fusion threshold λ_fusion = 5; when C_fusion > λ_fusion, output the identification result Y = 1, otherwise output the identification result Y = 0. In this example, C_fusion = 17.59 > λ_fusion = 5, so the identification result Y = 1. Calculate the final identification credibility index Q: Q = |C_fusion|·(1 + ε·Con_k); where Q represents the identification credibility index; |C_fusion| represents the absolute value of the fusion confidence; ε represents the template consistency weight coefficient, set to 0.3; Con_k represents the best template agreement rate. Calculate the identification credibility index: Q = 17.59·(1 + 0.3·0.75) = 17.59·(1 + 0.225) = 17.59·1.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 substation power grid, and combine environmental data such as temperature and humidity to calculate the environmental change index vector E: Grid voltage volatility: ΔV = max(V) - min(V) / mean(V), for example, ΔV = 0.05; Load change rate: ΔP = std(P) / mean(P), for example, ΔP = 0.12; Total harmonic distortion rate: THD = sqrt(∑V_h 2 / V_1 2 ), for example, THD = 0.03; Temperature index: T n orm=(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, use the unsupervised clustering method to 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, 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: Phase difference threshold λ p Adjustment: λ 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 basic 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: λpn ew = 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 topological stability.

[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 a high-confidence identification result 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 system's adaptability 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 key points for the next optimization, 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 performances of three methods for identifying the current characteristics of the substation area were compared: the traditional frequency analysis method (Method A), the method based on phase coding 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 seasonal 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 coding without using 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. Based on this characteristic, the characteristic current signals of adjacent power frequency alternating current cycles can be staggered. By subtracting the waveforms of adjacent power frequency alternating current cycles, the interference of the 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 the power frequency, it can be seen that in addition to the above two frequency points, there are also 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, in this embodiment, a second-order IIR digital notch filter is used for filtering. The transfer function of the notch filter with angular frequency ω 0 in the z-plane is: H(z) = (z2 -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: -b 1 y(n - 1) - b 2 y(n - 2). Where y(n - 1) is the previous output value of the filter; a 0 = 1, a 1 = -2cos(ω 0 / Fs); a 2 = 1; b 1 = -(1 - μ)2cos(ω 0 / Fs); b 2 = (1 - μ) 2 . In this embodiment, the notch depth parameter μ = 0.01 and Fs = 5400 Hz. Three notch filters are in series, corresponding to frequencies of 750 Hz, 850 Hz, and 950 Hz.

[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 to form a periodic signal, which can be expanded using 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 nkdt). 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.30 Hz and F_2 = 883.30 Hz. 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 the a n , b n values of the 2 characteristic frequency points, the cos(2πf n kdt), sin(2πf n kdt) values corresponding to the 2 frequency points can be calculated and stored as coefficients in an array, and then the a n , b n values can be obtained through the multiply-accumulate formula. In this embodiment, the a n , b n values are calculated once for each cycle period, and each time a cycle signal sample value is updated.

[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 parts without characteristic signals. In this embodiment, if the correlation peak value is greater than 2 times 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, it is identified that there is a characteristic signal in the current, and the identification of the current characteristics in the power distribution area is completed.

[0204] The present invention elaborates in detail a method for identifying the current characteristics of a transformer substation area based on the correlation peak decision method. Through steps such as generating and injecting rotation phase-encoded characteristic signals, collecting and preprocessing multi-point phase signals, constructing a self-organizing phase network and extracting features, making a phase difference correlation peak decision and 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 characteristics, problems of harmonic interference and frequency overlap are effectively avoided; Topological invariant features are introduced, enabling the system to have strong resistance to local data loss or distortion; An adaptive correlation peak decision mechanism is adopted to improve the reliability and adaptability of the decision; A complete adaptive optimization and parameter update mechanism is established, enabling 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 current characteristics of a transformer substation based on a correlation peak judgment method, characterized in that: include: Collect the synchronous phase information of the power grid, generate a characteristic signal with phase coding and inject it into the power line to form an injection current signal; At predetermined key nodes in the substation area, line current signals including injected current signals are collected, phase information is preprocessed and extracted, and a phase time series is obtained; based on this, a current characteristic phase difference pattern is extracted to generate a topology-invariant feature vector; The topology invariant feature vector is correlated with the preset standard feature template, the correlation peak is calculated and the decision criterion is established to determine the identification result and the identification credibility index; accordingly, the system parameters are dynamically adjusted, and the phase network structure and decision threshold are updated.

2. The method according to claim 1, characterized in that The steps of forming the injection current signal include: Read the grid voltage signal, obtain the grid synchronization phase information through zero crossing detection and phase-locked loop technology, and generate a rotating phase coding sequence; Based on a preset fixed carrier frequency, the rotating phase coding sequence is mapped onto the carrier signal to generate a phase modulated carrier signal; it is segmented at pseudo-random time intervals to form an intermittent pulse sequence and injected into the power line to generate an injection current signal.

3. The method according to claim 2, characterized in that The steps of generating a rotating phase encoding sequence include: Read the grid synchronization phase information, build the phase shift basic unit based on the Gray code principle, and form the coding base sequence; Apply the Hamming code principle to the coding base sequence, increase the coding redundancy, and generate a redundant coding sequence; analyze the phase jump characteristics, adjust the phase value through the phase mapping function, and generate an optimized redundant sequence; The optimized redundant sequence is converted into a differential phase encoding form and combined with a preset absolute phase reference point to construct a rotational phase encoding sequence.

4. The method according to claim 2, characterized in that: The steps of forming a discontinuous pulse train include: Obtain a pseudo-random sequence and map it to a time interval domain to generate a pseudo-random time interval sequence; Analyze the current characteristics and interference spectrum characteristics of the substation area, and determine the optimal pulse duration based on the principle of maximizing the signal-to-noise ratio; based on this, optimize the pulse edge and generate an optimized window function; Obtain line impedance characteristics and background noise level, dynamically adjust pulse energy, and generate adaptive amplitude sequences; The phase modulated carrier signal is combined with a pseudo-random time interval sequence, an adaptive amplitude sequence and an optimized window function to generate a discontinuous pulse sequence.

5. The method according to claim 1, characterized in that The steps to obtain the phase time series include: Select a predetermined number of key nodes in the substation area, use a high-precision current sensor to collect line current signals, and set the sampling frequency to a value higher than the carrier frequency; The line current signal is segmented into synchronous signal segments with the same length, and filtered through a bandpass filter whose center frequency is the carrier frequency to obtain a filtered signal; The Hilbert transform is applied to the filtered signal to construct the analytical signal, extract the instantaneous phase information, and obtain the phase time series of each node.

6. The method according to claim 5, characterized in that The steps of extracting instantaneous phase information and obtaining the phase time series of each node include: The filtered signal is processed to zero mean and abnormal peaks are removed to generate a preprocessed signal; its orthogonal components are calculated using Hilbert transform to generate a Hilbert transform result; Based on the preprocessed signal and its Hilbert transform result, the analytical signal is constructed and the instantaneous phase is extracted and the phase unwrapping is performed to obtain the unwrapped phase sequence; The unwrapped phase sequence is subjected to trend analysis, the linear phase growth component is extracted and removed, and the phase is normalized to obtain the final phase time series.

7. The method according to claim 1, characterized in that The steps to generate topologically invariant eigenvectors include: For the phase time series, the phase difference between nodes is calculated to form a phase difference matrix; based on this, the self-organizing network construction algorithm is used to establish the phase network topology structure of the substation area; Perform time series analysis on the phase network topology structure obtained for a predetermined number of consecutive cycles, extract dynamic change characteristics, and obtain the network evolution feature matrix; Extract topologically invariant feature vectors with environmental invariance from the network evolution feature matrix.

8. The method according to claim 7, characterized in that The steps of applying the self-organizing network construction algorithm to establish the phase network topology structure of the substation area include: Perform time dimension stability evaluation on each element of the phase difference matrix, calculate the standard deviation, and generate a phase difference stability matrix; Based on the phase difference stability matrix, an adaptive threshold matrix is ​​generated; the adaptive threshold matrix is ​​combined with the phase difference stability matrix to establish the connection strength of the continuous weights and generate a connection weight matrix; Based on historical data, the time stability of each connection in the connection weight matrix is ​​evaluated, the connection reliability index is calculated, and the connection reliability matrix is ​​generated; it is combined with the connection weight matrix to construct a comprehensive weight matrix to form the final phase network topology.

9. The method according to claim 1, characterized in that: The steps for determining the identification results and identification confidence indicators include: Through injection experiments under ideal conditions, the topologically invariant feature vector under standard conditions is obtained as a standard feature template; multi-dimensional correlation analysis is performed between it and the topologically invariant feature vector, and a topological weighted correlation function is introduced to obtain a correlation peak sequence; Based on the pre-stored environmental noise level and historical recognition results, an adaptive decision threshold function is constructed, the relevant peak sequence is compared with the threshold, and the decision result is output; Based on the adaptive decision threshold function, the decision confidence index is calculated; it is combined with the consistency of the consecutive predetermined rounds of decision results to generate the final recognition result and recognition confidence index.

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

Images