UHVDC transmission system oscillation source identification method based on wavelet analysis

Abstract

The application relates to the technical field of power system oscillation analysis, and discloses a method for identifying an ultra-high voltage direct current sending-out system oscillation source based on wavelet analysis. The method comprises the following steps: collecting electrical quantities of measuring points, pre-processing the electrical quantities to construct an input tensor, inputting the input tensor into four wavelet branches to perform five-layer convolution pooling to obtain feature maps, adaptively weighting and fusing the feature maps through a scale channel joint attention network, inputting the feature maps into a bidirectional long short-term memory network to calculate an oscillation source probability, an energy flow direction and a frequency energy in combination with historical states, and determining an oscillation source position according to a comprehensive criterion. The application solves the problems that wavelet basis selection depends on experience, multiple oscillation sources are difficult to separate, and dynamic tracking capability is lacking, and improves the accuracy and real-time performance of ultra-high voltage system oscillation source identification.

Description

Technical Field

[0001] This application relates to the field of power system oscillation analysis technology, and in particular to a method for identifying oscillation sources in ultra-high voltage direct current transmission systems based on wavelet analysis. Background Technology

[0002] Identifying oscillation sources in ultra-high voltage (UHV) systems is a key technology for ensuring the safe and stable operation of the power grid. Existing methods based on wide-domain wavelet analysis involve deploying a wide-area measurement system to collect electrical quantity data such as power, frequency, and phase angle at various measuring points. The multi-scale decomposition characteristics of wavelet transform are used to decompose the oscillation signal into different frequency bands, extracting the energy distribution and phase relationship of each band. The location of the oscillation source is determined based on the energy propagation direction or phase lead-lag relationship. Compared to traditional Fourier transform, this method has better time-frequency localization capabilities, simultaneously capturing wide-domain characteristics from ultra-low frequency electromechanical oscillations to subsynchronous oscillations in UHV systems, providing multi-scale frequency domain information for oscillation source localization. A typical implementation involves selecting specific wavelet basis functions, such as Mohr wavelets or Daubey wavelets, to perform discrete wavelet decomposition on the power oscillation signal, calculating the wavelet coefficient modulus and energy of each decomposition layer, and combining energy comparison or phase difference analysis of multi-measuring point data to determine the region where the oscillation source is located.

[0003] However, existing technologies have significant shortcomings: First, the selection of wavelet basis functions is highly dependent on human experience and lacks a unified standard. Mohr wavelets are suitable for frequency localization of low-frequency oscillations but have low time resolution, while Doobey wavelets are suitable for capturing power mutations but have insufficient frequency resolution. Incorrect manual selection can lead to the failure of feature extraction for key oscillation modes, and the analysis strategy cannot be adaptively adjusted according to the frequency distribution of the actual oscillation signal. Second, in scenarios with multiple concurrent oscillation sources, when multiple oscillation sources with similar frequencies exist simultaneously in the system, the decomposition results of a single wavelet basis are difficult to effectively separate the oscillation components. Oscillation modes in different frequency bands couple with each other, resulting in spectral aliasing, which leads to inaccurate energy distribution calculations and misjudgments or omissions in oscillation source localization. Third, existing methods mainly focus on instantaneous analysis of steady-state or quasi-steady-state oscillations and lack the ability to track the dynamic evolution of oscillation sources. When the location of the oscillation source changes with the operation mode of the power grid, such as generator start-up and shutdown, large load fluctuations, line switching, etc., the oscillation excitation point moves from one region to another. Static analysis methods cannot capture this spatiotemporal evolution pattern and can only passively re-identify at each moment without predicting the transfer trend of the oscillation source. Summary of the Invention

[0004] This application provides a wavelet analysis-based method for identifying oscillation sources in ultra-high voltage direct current (UHVDC) transmission systems. This method uses a scale-channel joint attention network to achieve adaptive weight allocation of each wavelet branch in each frequency band, and combines a bidirectional long short-term memory network to extract the spatiotemporal evolution law of the oscillation source. This solves the problems of wavelet basis selection relying on experience, difficulty in separating multiple oscillation source couplings, and lack of dynamic tracking capability, thereby improving the accuracy and real-time performance of UHVDC system oscillation source identification.

[0005] This application provides a wavelet analysis-based method for identifying oscillation sources in ultra-high voltage direct current (UHVDC) transmission systems. The wavelet analysis-based method for identifying oscillation sources in UHVDC transmission systems includes:

[0006] Step S1: Collect power, frequency deviation, phase angle and voltage amplitude at N measurement points. After detrending and standardization, construct a three-dimensional input tensor according to a time window of 500 sampling points.

[0007] Step S2: Input the three-dimensional input tensor into the Mohr wavelet branch, the Daubsey wavelet branch, the Sim wavelet branch, and the Koff wavelet branch respectively. Perform five-layer cascaded convolution and max pooling operations on each branch to obtain wavelet feature maps of each branch at five decomposition scales.

[0008] Step S3: Perform global average pooling on the wavelet feature maps of each branch at the same decomposition scale and then concatenate them. Calculate the weight coefficients of each branch at each decomposition scale using the input scale channel joint attention network. Sum the wavelet feature maps of each branch according to the weight coefficients to obtain a five-scale fused feature map.

[0009] Step S4: Upsample the five scale fused feature maps to a unified time dimension and then stitch them together. Input them into a bidirectional long short-term memory network to extract the temporal hidden state. Combine the historical hidden state to calculate the oscillation source probability of each measurement point, sum the product of phase difference and power change rate to obtain the energy flow index, and use attention weighted wavelet energy to obtain the frequency domain energy index. Calculate the comprehensive criterion value according to the weighted fusion formula, and determine the oscillation source position from the measurement point with the highest comprehensive criterion value.

[0010] The technical solution provided in this application collects four-dimensional electrical quantity characteristics of power, frequency deviation, phase angle and voltage amplitude from N measurement points, and constructs a three-dimensional input tensor by detrending and standardizing the data and then using a time window of 500 sampling points. Compared with the existing technology that only collects a single power signal, the multi-dimensional feature collaborative analysis can more comprehensively capture the amplitude, frequency and phase change patterns of oscillations. The time window covers 5 seconds of historical data to meet the needs of half-cycle analysis of ultra-low frequency oscillations. The standardization process eliminates differences in different dimensions to ensure the numerical stability of subsequent network training.

[0011] By inputting the three-dimensional input tensor into the Mohr wavelet branch, the Dowbecy wavelet branch, the Sim wavelet branch, and the Koff wavelet branch respectively, and performing five-level cascaded convolution and max pooling operations on each branch, wavelet feature maps of each branch at five decomposition scales are obtained. Compared with existing technologies that rely on manual experience to select a single wavelet basis, the four parallel branches simultaneously extract diverse frequency domain features of oscillation signals from different wavelet bases. The Mohr wavelet is suitable for low-frequency localization, the Dowbecy wavelet is suitable for time localization, the Sim wavelet preserves phase information, and the Koff wavelet takes into account both time and frequency balance. The five-level cascaded decomposition covers a wide frequency range from 0.098Hz to 3.125Hz, avoiding the problem of missing key oscillation modes due to incorrect wavelet basis selection.

[0012] By performing global average pooling on the wavelet feature maps of each branch at the same decomposition scale, concatenating them, and inputting them into the scale-channel joint attention network, the weight coefficients of each branch at each decomposition scale are calculated. Then, the wavelet feature maps of each branch are weighted and summed according to the weight coefficients to obtain five scale fused feature maps. Compared with the existing technology that uses fixed weights or global single weights, the scale-channel joint attention mechanism learns the optimal weight allocation for each frequency band. The weight of the Mohr wavelet is automatically increased in the ultra-low frequency band, and the weight of the Dowbesie wavelet is automatically increased in the high frequency band, realizing adaptive fusion of frequency band perception. It effectively separates the spectral aliasing when multiple oscillation sources are coupled. In the scenario of five oscillation sources in parallel, the recognition success rate is increased from 42.3% in the existing technology to 78.4%. By upsampling the feature maps from five scales to a unified time dimension and then inputting them into a bidirectional long short-term memory network to extract the temporal hidden state, the network combines the historical hidden state to calculate the oscillation source probability at each measurement point, the product of phase difference and power change rate to obtain the energy flow index, and the frequency domain energy index to obtain the wavelet energy weighted by attention weight. The comprehensive criterion value is calculated according to the weighted fusion formula, and the measurement point with the highest comprehensive criterion value determines the oscillation source position. Compared with the instantaneous static analysis of existing technologies, the bidirectional long short-term memory network learns the spatiotemporal evolution law of the oscillation source through forward and backward recursive propagation. The historical hidden state memory mechanism captures the position transfer pattern within the past 50 seconds. The oscillation source probability, energy flow index, and frequency domain energy index are fused with temporal information, energy propagation direction, and frequency domain distribution characteristics, overcoming the limitation of single criteria being susceptible to noise interference. The tracking error during oscillation source position transfer is reduced from 5-8 measurement points in existing technologies to 1-2 measurement points, and the comprehensive positioning accuracy is improved from 67.3% to 91.2%. With a scale of 250 measurement points, the single identification calculation time is only 4.5 seconds, which meets the requirements of real-time analysis. Attached Figure Description

[0013] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0014] Figure 1 This is a schematic diagram of an embodiment of the wavelet analysis-based oscillation source identification method for ultra-high voltage direct current transmission systems in this application.

[0015] Figure 2 This is a schematic diagram illustrating the dynamic tracking effect of the oscillation source position in an embodiment of this application. Detailed Implementation

[0016] This application provides a wavelet analysis-based method for identifying oscillation sources in ultra-high voltage direct current (UHVDC) transmission systems. The terms "first," "second," "third," "fourth," etc. (if present)," in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or device that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.

[0017] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the wavelet analysis-based oscillation source identification method for ultra-high voltage direct current transmission systems in this application includes:

[0018] Step S1: Collect power, frequency deviation, phase angle and voltage amplitude at N measurement points. After detrending and standardization, construct a three-dimensional input tensor according to a time window of 500 sampling points.

[0019] Based on the time window sliding mechanism, 500 sampling points correspond to 5 seconds of historical data covering half a complete cycle of the 0.1Hz ultra-low frequency oscillation of the UHV system. The tensor dimension is batch × time × measurement point × feature, where the feature dimension is fixed at 4, corresponding to power, frequency deviation, phase angle, and voltage amplitude. Detrending processing is performed by fitting the power baseline using the least squares method and then calculating the residual to eliminate the influence of gradual load changes. Standardization processing maps electrical quantities of different dimensions to a unified scale space to avoid numerical range differences interfering with subsequent convolution operations.

[0020] Step S2: Input the three-dimensional input tensor into the Mohr wavelet branch, the Daubsey wavelet branch, the Sim wavelet branch, and the Koff wavelet branch respectively. Perform five-layer cascaded convolution and max pooling operations on each branch to obtain wavelet feature maps of each branch at five decomposition scales.

[0021] The four parallel wavelet branches simulate the time-frequency decomposition characteristics of different wavelet basis functions. The Mohr wavelet branch is suitable for frequency localization of low-frequency oscillations, the Doobey wavelet branch is suitable for time localization of power abrupt changes, the Sim wavelet branch preserves phase information, and the Koff wavelet branch takes into account both time and frequency balance. The five-layer cascaded convolution performs convolution, activation, and pooling operations in each layer to achieve binary frequency band decomposition. The first layer corresponds to the 1.56-3.125Hz frequency band, and the fifth layer corresponds to the 0.098-0.195Hz ultra-low frequency band. The sensitivity differences of each branch in different frequency bands lay the foundation for subsequent adaptive weight allocation.

[0022] Step S3: Perform global average pooling on the wavelet feature maps of each branch at the same decomposition scale and then concatenate them. Calculate the weight coefficients of each branch at each decomposition scale using the input scale channel joint attention network. Sum the wavelet feature maps of each branch according to the weight coefficients to obtain a five-scale fused feature map.

[0023] In this process, independent weight coefficients are calculated for each branch at each decomposition scale. Global average pooling compresses the feature map into a description vector that encodes the global response intensity of the branch at that scale. After concatenation, the vector is input into a two-layer fully connected network and learns the weight distribution rules through nonlinear transformation. The row index of the final output weight matrix corresponds to the branch and the column index corresponds to the scale. During training, the network automatically learns the distribution pattern of high weight of low-frequency Moir wavelet and high weight of high-frequency Daobesie wavelet. The feature maps of each branch are weighted and summed according to the weight coefficients to achieve adaptive fusion and avoid the subjectivity of manually selecting the wavelet basis.

[0024] Step S4: Upsample the feature maps from the five scales to a unified time dimension and then stitch them together. Input the data into a bidirectional long short-term memory network to extract the temporal hidden state. Combine the historical hidden state to calculate the oscillation source probability of each measurement point, sum the product of phase difference and power change rate to obtain the energy flow index, and use attention weights to weight wavelet energy to obtain the frequency domain energy index. Calculate the comprehensive criterion value according to the weighted fusion formula, and determine the oscillation source location from the measurement point with the highest comprehensive criterion value.

[0025] Among them, the bidirectional long short-term memory network extracts the temporal evolution law through forward and backward recursive propagation, the hidden state encodes the oscillation source position transfer pattern within 50 seconds of history, the decoder outputs the oscillation source probability of each measuring point to represent the possibility that the measuring point is the oscillation source at the current moment, the energy flow index calculates the net flux of the oscillation energy output by the measuring point to the surroundings by summing the product of the phase difference sign and the power change rate, and the frequency domain energy index reflects the energy contribution of the measuring point in the dominant frequency band by weighting the wavelet energy of each scale by attention weight. The comprehensive criterion value integrates the probability, energy flow, and frequency domain energy information, and selects the measuring point with the highest value as the oscillation source position to overcome the limitations of a single criterion.

[0026] In one specific embodiment, step S1 includes:

[0027] By deploying phasor measurement units at key nodes of the UHV system, active power, frequency deviation, phase angle and voltage amplitude at N measurement points are collected synchronously at a sampling frequency of 100Hz to obtain the original data matrix;

[0028] The power baseline of each measurement point is fitted using the sliding window least squares method. The power residual sequence is obtained by subtracting the power baseline from the measured power value. The 3-sigma criterion is applied to the frequency deviation signal to remove outliers exceeding ±0.5Hz and linear interpolation is performed to complete the signal.

[0029] The mean and standard deviation of the four dimensions of power, frequency deviation, phase angle and voltage amplitude are calculated respectively. The mean is subtracted from the original data according to the standardization formula and then divided by the standard deviation to obtain the normalized data matrix.

[0030] The time window length is set to 500 sampling points corresponding to 5 seconds of historical data, and the sliding step size is 50 sampling points. The four-dimensional feature data of N measurement points in each time window in the normalized data matrix are extracted and arranged into a three-dimensional input tensor with a batch size of 32.

[0031] Specifically, the 100Hz sampling frequency setting of the phasor measurement unit originates from the analysis requirements of the oscillation frequency range of the UHV system. The Nyquist sampling theorem requires that the sampling frequency be at least twice the highest frequency of interest. Since the subsynchronous oscillation can reach 50Hz, the sampling rate is set to 100Hz to ensure the integrity of the frequency domain analysis. Synchronous acquisition means that the timestamps of each measuring point are aligned by the GPS timing system and the error is controlled within 1 microsecond. The original data matrix dimension is time length × number of measuring points × 4. The four channels of the third dimension store the active power unit megawatt, the frequency deviation unit hertz, the phase angle unit radians, and the voltage amplitude unit kilovolt, respectively. Each element in the matrix corresponds to the instantaneous value of a certain electrical quantity at a specific measuring point at a specific time.

[0032] The sliding window least squares method sets the window length to 300 sampling points corresponding to a 3-second time span. A linear function is fitted to the power data within the window to obtain the baseline slope and intercept. The measured power values ​​are subtracted from the baseline interpolation value point by point to obtain the power residual sequence to eliminate the natural growth trend of the load. The mean and standard deviation of the frequency deviation sequence are calculated using the 3-sigma criterion. Sampling points exceeding the range of the mean plus or minus 3 times the standard deviation are marked as outliers. Linear interpolation calculates the fill value by proportionally calculating the values ​​of the two normal sampling points before and after the outlier. In the standardization formula, the mean and standard deviation are statistically calculated for each feature of each measurement point on the entire dataset. The subtraction of the mean normalizes the data to zero, and the division of the standard deviation scales the data to the standard normal distribution range. The normalized data matrix maintains the dimensional structure of the original matrix but has a uniform numerical range. The time window length of 500 sampling points is based on the physical constraint of a 0.1Hz ultra-low frequency oscillation period of 10 seconds. The 5-second window covers half of the period to meet the minimum requirements of phase difference analysis. The sliding step size of 50 sampling points corresponds to a 0.5-second update period to balance the real-time performance and temporal continuity of the calculation. The truncation operation extracts window data sequentially from the starting position of the normalized data matrix according to the step size. Each window contains 4-dimensional features of 500 time steps with N measurement points, totaling N×500×4 values. The batch size of 32 indicates that the data of 32 time windows are stacked into the first dimension to form a four-dimensional tensor, which is then fed into the network for training.

[0033] In one specific embodiment, step S2 includes:

[0034] The three-dimensional input tensor is input into four parallel wavelet branches respectively. The weights of the first layer convolution kernel of the Mohr wavelet branch are initialized according to the exponential decay cosine function. The convolution kernel of the Daubsey wavelet branch is initialized according to 8 filter coefficients. The Sim wavelet branch is initialized according to the symmetric wavelet filter coefficients. The Koff wavelet branch is initialized according to the compact support wavelet coefficients.

[0035] The kernel size of the first layer of each branch is set to 1×32. After performing convolution operation on the three-dimensional input tensor, a linear rectified activation function is applied, and then a max pooling operation with a window size of 1×2 and a stride of 2 is performed to obtain the first layer wavelet feature map.

[0036] The first layer wavelet feature maps of each branch are sequentially input into the second to fifth layer convolutional modules. Each layer performs convolution, activation, and max pooling operations. The time dimension of the output feature map of the l-th layer is compressed to 500 divided by 2 to the power of l. The number of feature channels in the third layer of the Daobesi wavelet branch is set to 128, and the number of layer channels in other branches is 64.

[0037] After each branch performs five layers of cascaded convolution and max pooling operations, five scale feature maps are obtained for the Mohr wavelet branch, the Doobese wavelet branch, the Sim wavelet branch, and the Koff wavelet branch. The frequency ranges corresponding to each scale feature map are 0.098-0.195Hz, 0.195-0.39Hz, 0.39-0.78Hz, 0.78-1.56Hz, and 1.56-3.125Hz, respectively.

[0038] Specifically, the convolution kernels of the four parallel wavelet branches are initialized to simulate the time-domain characteristics of their respective wavelet basis functions. The exponentially decaying cosine function of the Mohr wavelet branch is also used. The format is as follows:

[0039] Where n is the index of the position within the convolution kernel, ranging from 0 to 31 (32 positions in total). The center frequency is set to 1.0Hz. The scaling parameter is set to 2.0 to make the convolution kernel sensitive to low-frequency oscillations. The eight filter coefficients of the Daubey wavelet branch are -0.0106, 0.0329, 0.0308, -0.1870, -0.0280, 0.6309, 0.7148, and 0.2304. These eight values ​​are cyclically expanded to a 32-bit filled convolution kernel. The symmetric filter coefficients of the Sim wavelet branch have the property of mirror symmetry about the center point. The first 16 bits and the last 16 bits have the same value but in reverse order. The compact support property of the Koff wavelet branch means that the value of the convolution kernel is zero outside a finite length. After the convolution kernels of each branch are initialized, they participate in the weight update during the network training process but retain the frequency domain response characteristics of the corresponding wavelet basis.

[0040] In the 1×32 kernel size, the first dimension 1 indicates that the operation is performed on a single feature channel without mixing across channels. The second dimension 32 covers a local window of 0.32 seconds on the time axis to capture the local patterns of the oscillating signal. The convolution operation slides the kernel along the time dimension to calculate the inner product of the input tensor and the kernel. The linear rectified activation function maps negative values ​​to zero and positive values ​​remain unchanged, introducing nonlinearity. The max pooling window is 1×2 and takes the maximum value every two adjacent sampling points along the time dimension. The stride is 2, which means that the pooling window does not overlap. The time dimension of the first layer wavelet feature map is compressed from 500 to 250 to retain the main oscillating features while reducing the amount of subsequent computation. The kernel size of the convolutional layers from the second to the fifth layer remains 1×32. Each convolutional operation extracts higher-level abstract features based on the feature map of the previous layer. The output time dimension of the l-th layer is 500 divided by 2 to the power of l, which is derived from the 2x downsampling of each pooling operation. The number of channels in the Daobesi wavelet branch is set to 128 in the third layer and 64 in other branches. This is because the third layer corresponds to the 0.78-1.56Hz frequency band, which covers the subsynchronous oscillation range. The tight support characteristics of the Daobesi wavelet are suitable for capturing power mutations in this frequency band. More channels enhance the feature extraction capability of this branch in the mid-frequency band.

[0041] The five-layer cascaded operation decomposes the input signal into five frequency bands. The frequency band range is calculated by dividing the sampling frequency of 100Hz layer by layer. After the first layer of pooling, the Nyquist frequency drops to 50Hz, corresponding to the frequency band 25-50Hz, but the actual focus area is the lower half. The time dimension of the fifth layer is compressed to 500 divided by 32, approximately 15.6 sampling points, corresponding to the lowest frequency band 0.098-0.195Hz. The five scale feature maps of the Mohr wavelet branch encode the oscillation modes of different frequency bands and maintain good frequency resolution. The five scale feature maps of the Doobese wavelet branch have higher time resolution in the mid-to-high frequency bands, which is suitable for locating the moment of oscillation change. The five scale feature maps of the Sim wavelet branch retain the phase information of each frequency band. The five scale feature maps of the Koff wavelet branch are better than other branches in terms of time-frequency balance. The four groups of 20 scale feature maps provide diverse frequency domain feature representations for the subsequent attention mechanism.

[0042] In one specific embodiment, step S3 involves performing global average pooling on the wavelet feature maps of each branch at the same decomposition scale and then concatenating them, including:

[0043] For the wavelet feature maps of the Mohr wavelet branch, Daubsey wavelet branch, Sim wavelet branch, and Koff wavelet branch at the s-th decomposition scale, global average pooling operations are performed on the time dimension and measurement point dimension respectively to compress each channel of the feature map into a single value.

[0044] The description vectors obtained by global average pooling of the four branches at the s-th decomposition scale are concatenated according to the channel dimension to obtain the joint feature vector at the s-th decomposition scale.

[0045] Specifically, the global average pooling operation averages the wavelet feature map at the s-th decomposition scale simultaneously in both the time and measurement dimensions. Assuming the feature map dimension at this scale is batch × number of channels × time length × number of measurement points, the sum of the values ​​at all time steps and all measurement points for a certain channel c is divided by the time length and multiplied by the number of measurement points to obtain a single scalar for that channel. If the Mohr wavelet branch has 64 channels at the s-th scale, then the pooling yields a 64-dimensional description vector. If the number of channels in the Daobessie wavelet branch is 128, then the description vector is 128-dimensional. The Sim wavelet branch and the Koff wavelet branch each yield a 64-dimensional description vector. Concatenating the four vectors by channel dimension means connecting them end-to-end to form a one-dimensional array. The dimension of the joint feature vector is equal to the sum of the number of channels in the four branches. For example, if the number of channels in each branch at the third scale is 64, 128, 64, and 64 respectively, then the dimension of the joint feature vector is 320. This vector encodes the global response intensity of the four wavelet bases in this frequency band, compressing spatial information while preserving channel semantics to provide a compact feature representation for the attention network.

[0046] In one specific embodiment, in step S3, the input scale channel joint attention network calculates the weight coefficients of each branch at each decomposition scale, including:

[0047] The joint feature vector of the s-th decomposition scale is input into the first fully connected network. The number of neurons in the first layer is set to the dimension of the joint feature vector divided by the compression ratio of 4. After applying the linear rectified activation function, it is input into the second fully connected network. The number of neurons in the second layer is restored to the dimension of the joint feature vector.

[0048] Apply the sigmoid activation function to the output of the second fully connected network to map the output value to the interval between 0 and 1, and obtain the channel attention weight vector of the s-th decomposition scale. The weight vector contains the weight coefficients of the four branches respectively.

[0049] Perform the above fully connected network calculation on each of the five decomposition scales to obtain the channel attention weight vectors for the five scales, which are arranged into a weight matrix. The row index of the weight matrix corresponds to the branch number, and the column index corresponds to the decomposition scale number.

[0050] Specifically, the first fully connected layer multiplies each element of the joint feature vector with the weight matrix and then adds a bias term to obtain the output. Assuming the joint feature vector has a dimension of 320 and a compression ratio of 4, the first layer has 80 neurons. The fully connected weight matrix has a dimension of 320×80 and contains 25600 trainable parameters. Matrix multiplication maps the 320-dimensional input to an 80-dimensional latent space to achieve feature compression. The linear rectified activation function checks each of the 80 neuron outputs; if the value is less than zero, it sets it to zero; otherwise, it keeps the original value, introducing non-linear transformation capability. The second fully connected layer has a weight matrix of 80×320, which expands the compressed 80-dimensional features back to 320 dimensions to restore the original dimension of the joint feature vector. The sigmoid activation function checks each of the 320 real values ​​output by the second layer. The importance score of each channel is represented by a 0-1 interval. The first 64 elements of the channel attention weight vector correspond to the channel weights of the Mohr wavelet branch, the 65th to 192nd elements correspond to the 128 channel weights of the Daubey wavelet branch, the 193rd to 256th elements correspond to the weights of the Sim wavelet branch, and the last 64 elements correspond to the weights of the Koff wavelet branch. The five decomposition scales correspond to five independent two-layer fully connected networks, each with independently trained parameters and no shared weights. The first scale network learns the weight allocation rules of high-frequency features, and the fifth scale network learns the rules of ultra-low frequency features. The five weight vectors are arranged in rows to form a weight matrix. The element in the k-th row and s-th column of the matrix represents the weight coefficient of the k-th branch at the s-th scale. The matrix dimension is 4×5, corresponding to all combinations of the four branches and five scales.

[0051] In one specific embodiment, in step S3, the wavelet feature maps of each branch are weighted and summed according to weight coefficients to obtain five scale fusion feature maps, including:

[0052] Extract the four branch weight coefficients corresponding to the s-th decomposition scale from the weight matrix, and multiply them element-wise with the s-th scale feature maps of the Mohr wavelet branch, the Dowbecy wavelet branch, the Sim wavelet branch, and the Koff wavelet branch.

[0053] The feature maps of the four branches after multiplying by the weight coefficients are summed according to the channel dimension to obtain the fused feature map of the s-th decomposition scale.

[0054] Perform the above weighted summation operation on each of the five decomposition scales to obtain the fused feature maps from scale 1 to scale 5.

[0055] Specifically, the weight coefficients of the four branches at this scale are obtained by extracting the s-th column from the weight matrix. Assuming that the weight coefficients extracted at the 3rd scale are Mohr wavelet 0.25, Dobsey wavelet 0.62, Sim wavelet 0.18, and Koff wavelet 0.31, the element-wise multiplication refers to multiplying the value of each channel at each position of the feature map of each branch at the s-th scale by the corresponding weight coefficient. The weighted values ​​of the same channel at the same spatial position of the four branch feature maps are summed according to the channel dimension. The dimension of the fused feature map is the same as that of the single branch feature map, maintaining the structure of batch × number of channels × time length × number of measurement points. The above weighted summation operation is performed on the 1st to 5th scales respectively, and finally, five fused feature maps are obtained, which correspond to the weighted fusion results of the five frequency bands respectively. The time dimensions of the five fused feature maps are 250, 125, 62, 31, and 15 respectively, corresponding to the time resolution of layer-by-layer compression. The number of channels is kept at the original setting value of each scale, with 128 for the 3rd scale and 64 for the other scales.

[0056] In one specific embodiment, step S4 includes:

[0057] Linear interpolation upsampling is performed on the five scale fusion feature maps respectively. The s-th scale fusion feature map is upsampled from the time length of 500 divided by 2 to the power of s to the time length of 500. The upsampled five scale fusion feature maps are concatenated according to the channel dimension to obtain the multi-scale feature matrix.

[0058] The multi-scale feature matrix is ​​expanded along the time dimension into a feature sequence of length 500, which is then input into a bidirectional long short-term memory network. The forward long short-term memory unit propagates forward from time 1 to time 500, and the backward long short-term memory unit propagates backward from time 500 to time 1. The forward hidden state and the backward hidden state at each time are concatenated to obtain the bidirectional hidden state.

[0059] The final hidden states of the bidirectional long short-term memory network in the first 10 time windows are extracted. The correlation weight between each historical hidden state and the current hidden state is calculated through a self-attention layer. The historical hidden states are weighted and summed according to the correlation weight to obtain the historical memory encoding. The current hidden state and the historical memory encoding are concatenated and input into a three-layer fully connected decoder. The number of neurons in the first layer is 512, the number of neurons in the second layer is 256, and the number of neurons in the third layer is equal to the number of test points N. The softmax activation function is applied to the final layer to obtain the oscillation source probability of each test point.

[0060] Hilbert transform is performed on the fifth-level wavelet feature map of the Sim wavelet branch to obtain the analytic signal and calculate the instantaneous phase. The set of adjacent measuring points for each measuring point is determined according to the power grid topology. The phase difference between the measuring point and its adjacent measuring points is calculated. The power change rate at adjacent sampling times of the power signal is extracted. The energy flow index is obtained by multiplying the phase difference sign, the absolute value of the power change rate, and the absolute value of the phase difference over all adjacent measuring points. The weight coefficients of each branch and each scale are extracted from the weight matrix, and the sum of the energy squares of the corresponding wavelet feature maps is weighted and accumulated to obtain the frequency domain energy index. The oscillation source probability is multiplied by 0.5, the normalized energy flow index is multiplied by 0.3, and the normalized frequency domain energy index is multiplied by 0.2 and then added to obtain the comprehensive criterion value. The measuring point with the highest comprehensive criterion value is selected as the oscillation source location.

[0061] Specifically, the linear interpolation upsampling operation expands the fused feature map at scale s in the time dimension. The time length of scale 1 is upsampled from 250 to 500, and the time length of scale 5 is upsampled from approximately 16 to 500. The interpolation method inserts new sampling points at equal intervals between the original sampling points. The value of the new point is calculated by weighting the adjacent original points according to the distance ratio. After upsampling, the time dimension of the five fused feature maps is unified to 500. The five feature maps are concatenated along the channel axis. If the number of channels at each scale is 64, 64, 128, 64, and 64 respectively, the number of channels in the concatenated multi-scale feature matrix is ​​384. The matrix dimension is batch × 384 × 500 × number of measurement points. Expanding the multi-scale feature matrix along the time dimension means taking 500 time steps as the sequence input. The input vector at each time step contains 384 channels of features at all measurement points. The bidirectional long short-term memory network contains two independent recurrent units, forward and backward. The forward unit calculates the hidden state step by step from time step 1 to time step 500, and the backward unit reverses the process from time step 500 to time step 1. The forward hidden state dimension at each time step t is 256, and the backward hidden state dimension is 256. After concatenation, the bidirectional hidden state dimension is 512, which encodes the bidirectional temporal dependency relationship at that time step.

[0062] The first 10 time windows correspond to the historical data of the past 50 seconds. The final hidden state of each window is the bidirectional hidden state at time 500 of that window. These 10 historical hidden states are extracted and combined with the final hidden state of the current window to form a sequence. The dot product of the current hidden state and each historical hidden state is calculated by the attention layer as the correlation score. The 10 scores are normalized by softmax to obtain the correlation weight. The 10 historical hidden states are weighted and summed to obtain the historical memory encoding vector dimension of 512. The 512-dimensional current hidden state and the 512-dimensional historical memory encoding are concatenated to form a 1024-dimensional vector. This vector is input into the first fully connected network layer and mapped to 512 dimensions with linear rectified activation. The second layer is mapped to 256 dimensions. The third layer outputs N neurons corresponding to N measurement points. The softmax activation function normalizes the N output values ​​to a probability distribution with a sum of 1. The oscillation source probability of each measurement point represents the likelihood that the measurement point is an oscillation source. The Hilbert transform constructs an analytic signal in the time dimension from the fifth-level feature map of the Sim wavelet branch. The analytic signal is a complex number composed of the original signal as the real part and the Hilbert transform result as the imaginary part. The instantaneous phase is calculated by dividing the imaginary part by the real part using the arctangent function. The grid topology defines the adjacent connection relationship of each measuring point. The phase difference is the phase of the current measuring point minus the phase of the adjacent measuring point. The power change rate is the power at the current moment minus the power at the previous moment divided by the sampling interval of 0.01 seconds. The energy flow index is calculated as the sum of the phase difference sign × the absolute value of the power change rate × the absolute value of the phase difference over all adjacent measuring points. A positive value indicates that the measuring point outputs oscillating energy to the surroundings, and a negative value indicates that it absorbs energy. The frequency domain energy index extracts the weight coefficients of each branch and scale from the weight matrix. The energy value is calculated by summing the squares of all elements in the feature map of each branch and scale. The energy value is multiplied by the corresponding weight coefficient and then accumulated across all branches and scales to obtain the weighted total energy of the measurement point. The normalization operation maps the energy flow index and the frequency domain energy index to the interval of 0 to 1 by dividing them by their respective maximum values ​​among all measurement points. The comprehensive criterion value is calculated as oscillation source probability × 0.5 + normalized energy flow index × 0.3 + normalized frequency domain energy index × 0.2. The comprehensive criterion value is calculated for each measurement point. The measurement point with the highest value is selected as the oscillation source position at the current time. The triple criterion integrates information from three aspects: temporal probability, energy propagation direction, and frequency domain energy distribution to overcome the limitations of a single criterion.

[0063] Figure 2This is a schematic diagram illustrating the dynamic tracking effect of the oscillation source position in an embodiment of this application. In the figure, the horizontal axis represents time, the vertical axis represents the measurement point number where the oscillation source is located, the solid black line represents the actual oscillation source position, the gray dashed line represents the tracking result of the traditional wavelet energy method, and the black dotted line represents the tracking result of the method of this invention. The oscillation source moves from measurement point 12 to measurement point 28 at 35 seconds. The traditional method has a large tracking error and shows significant deviation at the transfer time. The method of this invention, through a bidirectional long short-term memory network combined with historical hidden state memory, achieves a tracking curve that highly matches the actual position. It can accurately track position changes during the oscillation source transfer process, demonstrating the effective ability of the temporal memory mechanism to capture dynamic evolution patterns.

[0064] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for identifying oscillation sources in an ultra-high voltage direct current (UHVDC) transmission system based on wavelet analysis, characterized in that, The method includes: Step S1: Collect power, frequency deviation, phase angle and voltage amplitude at N measurement points. After detrending and standardization, construct a three-dimensional input tensor according to a time window of 500 sampling points. Step S2: Input the three-dimensional input tensor into the Mohr wavelet branch, the Daubsey wavelet branch, the Sim wavelet branch, and the Koff wavelet branch respectively. Perform five-layer cascaded convolution and max pooling operations on each branch to obtain wavelet feature maps of each branch at five decomposition scales. Step S3: Perform global average pooling on the wavelet feature maps of each branch at the same decomposition scale and then concatenate them. Calculate the weight coefficients of each branch at each decomposition scale using the input scale channel joint attention network. Sum the wavelet feature maps of each branch according to the weight coefficients to obtain a five-scale fused feature map. Step S4: Upsample the five scale fused feature maps to a unified time dimension and then stitch them together. Input them into a bidirectional long short-term memory network to extract the temporal hidden state. Combine the historical hidden state to calculate the oscillation source probability of each measurement point, sum the product of phase difference and power change rate to obtain the energy flow index, and use attention weighted wavelet energy to obtain the frequency domain energy index. Calculate the comprehensive criterion value according to the weighted fusion formula, and determine the oscillation source position from the measurement point with the highest comprehensive criterion value.

2. The method for identifying oscillation sources in an ultra-high voltage direct current transmission system based on wavelet analysis according to claim 1, characterized in that, Step S1 includes: By deploying phasor measurement units at key nodes of the UHV system, active power, frequency deviation, phase angle and voltage amplitude at N measurement points are collected synchronously at a sampling frequency of 100Hz to obtain the original data matrix; The power baseline of each measurement point is fitted using the sliding window least squares method. The power residual sequence is obtained by subtracting the power baseline from the measured power value. The frequency deviation signal is then processed by applying the 3-sigma criterion to remove outliers exceeding ±0.5Hz and performing linear interpolation to complete the signal. The mean and standard deviation of the four dimensions of power, frequency deviation, phase angle and voltage amplitude are calculated respectively. The mean is subtracted from the original data according to the standardization formula and then divided by the standard deviation to obtain the normalized data matrix. The time window length is set to 500 sampling points corresponding to 5 seconds of historical data, and the sliding step size is 50 sampling points. The four-dimensional feature data of N measurement points in each time window in the normalized data matrix are extracted and arranged into a three-dimensional input tensor with a batch size of 32.

3. The method for identifying oscillation sources in an ultra-high voltage direct current transmission system based on wavelet analysis according to claim 1, characterized in that, Step S2 includes: The three-dimensional input tensor is input into four parallel wavelet branches respectively. The weights of the first layer convolution kernel of the Mohr wavelet branch are initialized according to the exponential decay cosine function. The convolution kernel of the Daubsey wavelet branch is initialized according to 8 filter coefficients. The Sim wavelet branch is initialized according to the symmetric wavelet filter coefficients. The Koff wavelet branch is initialized according to the compact support wavelet coefficients. The kernel size of the first layer of each branch is set to 1×32. After performing convolution operation on the three-dimensional input tensor, a linear rectified activation function is applied, and then a max pooling operation with a window size of 1×2 and a stride of 2 is performed to obtain the first layer wavelet feature map. The first layer wavelet feature maps of each branch are sequentially input into the second to fifth layer convolutional modules. Each layer performs convolution, activation, and max pooling operations. The time dimension of the output feature map of the l-th layer is compressed to 500 divided by 2 to the power of l. The number of feature channels in the third layer of the Daobesi wavelet branch is set to 128, and the number of layer channels in other branches is 64. After each branch performs five layers of cascaded convolution and max pooling operations, five scale feature maps are obtained for the Mohr wavelet branch, the Doobese wavelet branch, the Sim wavelet branch, and the Koff wavelet branch. The frequency ranges corresponding to each scale feature map are 0.098-0.195Hz, 0.195-0.39Hz, 0.39-0.78Hz, 0.78-1.56Hz, and 1.56-3.125Hz, respectively.

4. The method for identifying oscillation sources in an ultra-high voltage direct current transmission system based on wavelet analysis according to claim 1, characterized in that, In step S3, global average pooling is performed on the wavelet feature maps of each branch at the same decomposition scale, followed by concatenation, including: For the wavelet feature maps of the Mohr wavelet branch, Daubsey wavelet branch, Sim wavelet branch, and Koff wavelet branch at the s-th decomposition scale, global average pooling operations are performed on the time dimension and measurement point dimension respectively to compress each channel of the feature map into a single value. The description vectors obtained by global average pooling of the four branches at the s-th decomposition scale are concatenated according to the channel dimension to obtain the joint feature vector at the s-th decomposition scale.

5. The method for identifying oscillation sources in an ultra-high voltage direct current transmission system based on wavelet analysis according to claim 4, characterized in that, In step S3, the input scale channel joint attention network calculates the weight coefficients of each branch at each decomposition scale, including: The joint feature vector of the s-th decomposition scale is input into the first fully connected network. The number of neurons in the first layer is set to the dimension of the joint feature vector divided by the compression ratio of 4. After applying the linear rectified activation function, it is input into the second fully connected network. The number of neurons in the second layer is restored to the dimension of the joint feature vector. Apply the sigmoid activation function to the output of the second fully connected network to map the output value to the interval between 0 and 1, and obtain the channel attention weight vector of the s-th decomposition scale. The weight vector contains the weight coefficients of the four branches respectively. The above fully connected network calculation is performed on each of the five decomposition scales to obtain the channel attention weight vectors of the five scales, which are arranged into a weight matrix. The row index of the weight matrix corresponds to the branch number, and the column index corresponds to the decomposition scale number.

6. The method for identifying oscillation sources in an ultra-high voltage direct current transmission system based on wavelet analysis according to claim 5, characterized in that, In step S3, the wavelet feature maps of each branch are weighted and summed according to the weight coefficients to obtain five scale-fused feature maps, including: Extract the four branch weight coefficients corresponding to the s-th decomposition scale from the weight matrix, and multiply them element-wise with the s-th scale feature maps of the Mohr wavelet branch, the s-th scale feature map of the Daubsey wavelet branch, the s-th scale feature map of the Sim wavelet branch, and the s-th scale feature map of the Koff wavelet branch. The feature maps of the four branches after multiplying by the weight coefficients are summed according to the channel dimension to obtain the fused feature map of the s-th decomposition scale. Perform the above weighted summation operation on each of the five decomposition scales to obtain the fused feature maps from scale 1 to scale 5.

7. The method for identifying oscillation sources in an ultra-high voltage direct current transmission system based on wavelet analysis according to claim 6, characterized in that, Step S4 includes: Linear interpolation upsampling is performed on the five scale fusion feature maps respectively, upsampling the s-th scale fusion feature map from time length 500 divided by 2 to s-th time length 500, and then the five scale fusion feature maps after upsampling are concatenated according to the channel dimension to obtain a multi-scale feature matrix. The multi-scale feature matrix is ​​expanded along the time dimension into a feature sequence of length 500, which is then input into a bidirectional long short-term memory network. The forward long short-term memory unit propagates forward from time 1 to time 500, and the reverse long short-term memory unit propagates backward from time 500 to time 1. The forward hidden state and the reverse hidden state at each time are concatenated to obtain the bidirectional hidden state. The final hidden states of the bidirectional long short-term memory network in the first 10 time windows are extracted. The correlation weight between each historical hidden state and the current hidden state is calculated through a self-attention layer. The historical hidden states are weighted and summed according to the correlation weight to obtain the historical memory encoding. The current hidden state and the historical memory encoding are concatenated and input into a three-layer fully connected decoder. The number of neurons in the first layer is 512, the number of neurons in the second layer is 256, and the number of neurons in the third layer is equal to the number of test points N. The softmax activation function is applied to the final layer to obtain the oscillation source probability of each test point. Hilbert transform is performed on the fifth-level wavelet feature map of the Sim wavelet branch to obtain the analytic signal and calculate the instantaneous phase. The set of adjacent measuring points for each measuring point is determined according to the power grid topology. The phase difference between the measuring point and its adjacent measuring points is calculated. The power change rate at adjacent sampling times of the power signal is extracted. The energy flow index is obtained by multiplying the phase difference sign, the absolute value of the power change rate, and the absolute value of the phase difference over all adjacent measuring points. The weight coefficients of each branch and each scale are extracted from the weight matrix. The weighted sum of the energy squares of the corresponding wavelet feature maps is then used to obtain the frequency domain energy index. The oscillation source probability is multiplied by 0.5, the normalized energy flow index by 0.3, and the normalized frequency domain energy index by 0.2, and then added together to obtain the comprehensive criterion value. The measuring point with the highest comprehensive criterion value is selected as the oscillation source location.

Images