Outdoor distribution transformer state monitoring system based on edge computing

Abstract

The application discloses an outdoor power distribution transformer state monitoring system based on edge computing and belongs to the technical field of power distribution transformer state monitoring. The system is arranged on an edge computing node at the side of a power distribution transformer and comprises a data acquisition module for acquiring voltage waveforms, current waveforms and shell vibration waveforms. A time domain solution module is used for solving instantaneous power to obtain a power factor angle sequence and a differential excitation sudden change feature. A cycle segmentation module is used for segmenting each power frequency cycle of the voltage waveforms according to the zero-crossing time of the shell vibration waveforms to obtain voltage amplitude envelopes and current harmonic subbands. A correlation alignment module is used for calculating the cross-correlation delay of the voltage amplitude envelopes and the current harmonic subbands and aligning the excitation sudden change feature to generate a load dynamic response vector. A state judgment module is used for calculating the projection distance of the vector in a preset winding state feature space to determine whether the power distribution transformer is in a non-stable operation area.

Description

Technical Field

[0001] This invention relates to the field of distribution transformer condition monitoring technology, specifically to an outdoor distribution transformer condition monitoring system based on edge computing. Background Technology

[0002] Outdoor distribution transformers operate in environments with large temperature differences, are susceptible to wind and rain erosion, and are affected by load fluctuations. Therefore, online monitoring of winding conditions is crucial for power supply reliability. Existing monitoring methods generally combine electrical quantity analysis with independent vibration sensing. They calculate steady-state power factor and current harmonic content using RMS voltage and current values, or perform spectral envelope analysis on the tank surface vibration signal to determine mechanical loosening. These approaches treat electrical and mechanical responses as separate sources, failing to address their coupling characteristics and timing correlation within the power frequency cycle. Consequently, they struggle to reflect transient anomalies caused by minor winding deformations or inter-turn insulation degradation.

[0003] In monitoring devices deployed at the edge, data analysis often involves comparing preset thresholds or extracting a small number of statistical features before uploading, lacking in-depth fusion processing of the waveform's fine structure. When the load changes dynamically, there is an inherent delay between the voltage amplitude envelope and the current harmonic components caused by the winding's electromagnetic inertia and nonlinear magnetic circuit. Furthermore, excitation transients are often accompanied by sudden changes in the power factor angle. These implicit coupling information within the time series are discarded by conventional steady-state analysis methods, resulting in the inability to promptly identify early unstable states of the winding. Simultaneously, the judgment model relies on single or multi-dimensional threshold logic, lacking an overall similarity measurement of high-dimensional feature vectors, blurring the boundary between stable and unstable regions, and easily leading to false alarms or missed alarms during load abrupt changes.

[0004] Condition monitoring of distribution transformers requires effectively integrating synchronously acquired voltage, current, and vibration waveforms at edge computing nodes, extracting a high-dimensional response vector that takes into account both excitation mutation timing and load dynamic delay, and achieving adaptive and highly sensitive judgment of unstable operating areas based on the spatial projection deviation of this vector. Summary of the Invention

[0005] An outdoor distribution transformer condition monitoring system based on edge computing is provided. By using edge computing nodes to perform joint time-domain analysis and cross-correlation alignment of voltage waveform, current waveform and shell vibration waveform, load dynamic response vector is generated and the spatial projection distance of winding state characteristics is used to accurately determine the unstable operation zone of the distribution transformer.

[0006] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides an outdoor distribution transformer condition monitoring system based on edge computing. The system is directly deployed on the edge computing node on the distribution transformer side. Through an integrated signal acquisition and analysis architecture, it completes the entire process monitoring from raw waveform to condition determination locally, significantly reducing the dependence on cloud transmission bandwidth and centralized computing resources. At the same time, it has the ability to sensitively identify early unstable operating trends of winding condition and to subdivide faults.

[0007] The system's data acquisition module is used to collect initial monitoring data during the operation of the distribution transformer. This initial monitoring data includes voltage waveforms, current waveforms, and shell vibration waveforms. The acquired multi-source waveforms are strictly synchronized and aligned in the time domain, providing a time consistency basis for subsequent multi-physical quantity coupling analysis. This avoids timing errors introduced by asynchronous sampling, thereby improving the accuracy of the correlation extraction between excitation mutations and mechanical responses.

[0008] The time-domain calculation module performs instantaneous power calculations on the voltage and current waveforms in the time domain to obtain a power factor angle sequence. It then performs differential operations on this power factor angle sequence to obtain excitation abrupt change characteristics. By directly utilizing the instantaneous phase difference and power angle changes, it can capture subtle jumps in the core excitation state within the power frequency cycle. This effectively transforms the excitation angle disturbance caused by the small-energy arc discharge triggered in the initial stage of an inter-turn short circuit into a quantifiable abrupt change marker, without waiting for the fault to develop into a steady-state current distortion.

[0009] The period segmentation module divides the voltage waveform into each power frequency cycle based on the zero-crossing moment of the shell vibration waveform, obtaining the voltage amplitude envelope and current harmonic sub-band within each power frequency cycle. Using the zero-crossing characteristic of the vibration waveform itself as the segmentation benchmark ensures that the waveform truncation is naturally synchronized with the phase of the transformer's mechanical vibration, guaranteeing that each analysis segment strictly corresponds to a complete vibration cycle of the transformer's main structure. This ensures that the voltage amplitude fluctuations and changes in current harmonic components completely correspond to the load cycle and mechanical stress cycle, significantly suppressing spectral leakage and timing misalignment caused by asynchronous truncation.

[0010] The relevant alignment module calculates the cross-correlation delay between the voltage amplitude envelope and the current harmonic subband, and aligns this cross-correlation delay with the excitation mutation characteristics in time to generate a load dynamic response vector. By extracting the time lag relationship between voltage amplitude changes and harmonic energy conduction, and synchronously embedding the delayed step at the moment of excitation mutation into the same vector, the generated load dynamic response vector simultaneously carries the time delay characteristics of the electrical response, the harmonic-added lag, and the delay transient information brought about by the excitation mutation, providing a high-dimensional and highly sensitive joint feature expression for the winding mechanical state.

[0011] The status determination module determines whether the distribution transformer is in an unstable operating region based on the projection distance of the load dynamic response vector in a preset winding status feature space, and outputs the status determination result. It utilizes feature subspace projection to obtain the principal component components under normal conditions, and uses residual energy to measure the degree of deviation. This allows even small feature offsets caused by winding loosening, deformation, or insulation degradation to be significantly amplified in the projection distance, avoiding the insensitivity of fixed thresholds or single-parameter comparisons to small defects. It also possesses good anti-interference capabilities, reducing the probability of false alarms caused by normal load fluctuations.

[0012] As a preferred embodiment of the present invention, during the initial monitoring data acquisition, voltage and current transformers connected via edge computing nodes synchronously acquire the three-phase voltage and current signals from the high-voltage and low-voltage sides, respectively. A piezoelectric accelerometer attached to the surface of the distribution transformer tank acquires the shell vibration waveform. The resulting analog signals are then converted from analog to digital and stored as a time-series matrix with the same sampling frequency. This multi-channel synchronous acquisition and matrix storage method ensures precise time-base alignment of voltage, current, and vibration at the data source, providing time-delay-free raw data for instantaneous power calculation and vibration period segmentation.

[0013] Preferably, in the instantaneous power calculation, the instantaneous power sequence is obtained by multiplying the voltage and current of each phase point by point, and the active power envelope is extracted by low-pass filtering. The instantaneous phase angle of the voltage is obtained by Hilbert transform, and the power factor angle in each power frequency cycle is calculated by combining the current zero-crossing time, forming a power factor angle sequence. The sequence is subjected to first-order forward differential, and the differential points exceeding the preset angle threshold are marked as excitation mutation points, and their amplitude and direction are recorded as excitation mutation characteristics. This method directly generates the angle sequence from the three-phase original waveform in real time, without the need for the steady-state assumption of the phase-locked loop. It is extremely sensitive to excitation anomalies under transient disturbances, and the preset angle threshold is preferably dynamically and adaptively calculated based on the transformer's rated capacity and short-circuit impedance percentage, so that the detection sensitivity of excitation mutations is automatically adapted between transformers with different capacities and design parameters, avoiding missed detections or over-triggering caused by a fixed threshold.

[0014] In terms of period segmentation, preferably, the positive and negative zero-crossing moments of the shell vibration waveform are detected, and the interval between two adjacent zero-crossing moments in the same direction is taken as a vibration period interval. The voltage and current waveforms are truncated according to this interval. The voltage segment waveform is subjected to absolute value calculation and sliding peak hold processing to obtain the voltage amplitude envelope. The current segment waveform is subjected to multi-order bandpass filtering to separate the fundamental component and multiple higher harmonic components. The higher harmonic amplitude sequence is combined into a current harmonic sub-band. The center frequencies of the multi-order bandpass filters are integer multiples of the current fundamental frequency, and the quality factor decreases with the increase of the harmonic order. This design retains the effective components of higher harmonics while adaptively suppressing high-frequency noise, so that the harmonic sub-band contains more realistic winding nonlinear conduction information.

[0015] Preferably, when calculating the cross-correlation delay, cross-correlation is performed on the fundamental amplitude sequences of the voltage amplitude envelope and current harmonic subbands within the same vibration period interval. After mean removal, sliding dot product, and normalization, the sub-sampling precision base delay is obtained through three-point Gaussian interpolation fitting, and its confidence index is calculated. Simultaneously, the harmonic delay of each harmonic relative to the voltage envelope is calculated, and the difference between the base delay and each harmonic delay is used as an additional hysteresis sequence. Then, the excitation abrupt change moment is time-matched with the base delay to extract the change value of the base delay before and after the step. Finally, the base delay, the additional hysteresis sequence, and the change value of the base delay are concatenated to form the load dynamic response vector. The sub-sampling precision delay ensures that the response time difference caused by the micro-deformation of the winding can be accurately captured, while the additional hysteresis sequence reflects the asymmetric delay of current distortion in different frequency bands propagating along the winding. Together with the delay step caused by the excitation abrupt change, they constitute a comprehensive characterization of winding loosening, displacement, and insulation degradation.

[0016] During state determination, preferably, a pre-stored winding state reference matrix is ​​read, which consists of multiple load dynamic response reference vectors under normal conditions. Singular value decomposition (SVD) is performed on the winding state reference matrix, and the left singular vectors corresponding to the k largest singular values ​​are selected as the orthogonal basis matrix of the feature subspace based on a preset energy retention threshold. The orthogonal projection of the load dynamic response vector onto this feature subspace is calculated to obtain the projection component vector, and the residual vector and its Euclidean norm are calculated as the projection distance. The projection distance is compared with the state distance threshold; if it exceeds the threshold, it is determined to be in an unstable operating region and an unstable state label is output; otherwise, a stable state label is output. The subspace projection method based on singular value decomposition can retain the main characteristic patterns of normal state data, weaken the influence of noise and minor operating condition fluctuations, and therefore the residual is extremely sensitive to changes in abnormal structures, and the computational complexity is suitable for real-time inference at edge nodes.

[0017] Once an unstable operating zone is identified, this invention further provides a fault type confirmation method. The load dynamic response vector corresponding to the unstable state is marked as an abnormal response vector. The cosine value of the angle between this vector and each reference vector in the winding state reference matrix is ​​calculated. The fault type label corresponding to the reference vector with the largest cosine value is selected as the preliminary fault type determination. A pre-stored spectral analysis template in the edge computing node is then called to calculate the matching degree, and the fault type confirmation result is output. This process automatically completes the decision-making process from anomaly detection to type classification without human intervention, enabling edge devices to not only provide status alarms but also offer actionable maintenance suggestions.

[0018] Furthermore, preferably, voltage and current waveform segments within the time period corresponding to the abnormal response vector can be extracted from the time series matrix and subjected to Fast Fourier Transform (FFT) to extract supplementary feature components such as the fundamental amplitude decrease rate, unbalance coefficient, zero-sequence current amplitude, and negative-sequence current amplitude. These components are then concatenated with the abnormal response vector to obtain an extended abnormal response vector, and the projection distance is recalculated for secondary judgment. By introducing frequency domain energy distribution and symmetry component information, the illusion of transient disturbances caused by purely mechanical transients or external load impacts is further eliminated, significantly improving the accuracy of confirming electrically induced unstable states such as inter-turn short circuits and phase-to-phase imbalances.

[0019] Preferably, after collecting the initial monitoring data, the voltage variance sequence and current variance sequence are calculated using a sliding window variance calculation. When both the voltage and current variances exceed the corresponding thresholds and the time difference is within the preset synchronization tolerance range, it is identified as an impact disturbance moment. Data segments of a preset length before and after this moment are deleted, and the remaining initial monitoring data is used for subsequent processing. This preprocessing eliminates the end-to-end interference of external instantaneous impacts on the excitation mutation characteristics and load dynamic response vector calculation, enabling the condition monitoring system to maintain stable judgment consistency under disturbance environments such as lightning strikes and capacitor switching.

[0020] The technical effects and advantages provided by the present invention in the above technical solution are as follows: Using the zero-crossing moment of the shell vibration waveform as a time base, the voltage waveform is segmented into each power frequency cycle, forming voltage and current segments synchronized with the mechanical vibration. Absolute value calculation and sliding peak hold are performed on the voltage segments to obtain the voltage amplitude envelope reflecting the load fluctuation amplitude. Multi-order bandpass filtering is applied to the current segments to separate the fundamental wave and various higher harmonic components, forming the current harmonic sub-band. This vibration zero-crossing driven period segmentation method ensures that the electrical feature extraction is strictly aligned with the actual excitation cycle of the transformer, eliminating waveform truncation distortion caused by fixed time windows and guaranteeing the physical consistency of the envelope and harmonic sub-band on the time axis. Within the same vibration cycle interval, cross-correlation is performed on the voltage amplitude envelope and the fundamental amplitude sequence in the current harmonic sub-band to obtain the basic delay. Then, cross-correlation is performed with each harmonic amplitude sequence to obtain the harmonic delay and construct an additional hysteresis sequence. This basic delay, the additional hysteresis sequence, and the excitation mutation characteristics obtained from the power factor angle difference are time-aligned and concatenated into a load dynamic response vector. The cross-correlation delay between voltage and current components objectively quantifies the electromagnetic response inertia of the winding under dynamic load and the phase lag of harmonic distortion transmission. The excitation mutation feature directly embeds the abrupt change of the discontinuous excitation state of the iron core into the vector, so that the final dynamic response vector carries both electromechanical coupling delay information and electromagnetic transient markers. This can significantly improve the sensitivity and feature resolution of characterizing early unstable states such as inter-turn micro-discharge and local contact loosening of the winding.

[0021] A winding state reference matrix, composed of multiple load dynamic response reference vectors under normal operating conditions, is pre-stored in the edge computing node. A feature subspace reflecting the main variation pattern of the healthy winding response is extracted through singular value decomposition. For the load dynamic response vector to be tested, its orthogonal projection onto the feature subspace is calculated to obtain the projected component vector, and the Euclidean norm of the projected residual vector is obtained as the projection distance. This projection distance measures the statistical distance of the current response vector from the normal feature set, reflecting the degree of multidimensional anomalies with a single continuous scalar value, avoiding the tuning difficulties and insufficient adaptability caused by setting thresholds for multiple feature dimensions separately. The projection operation involves only matrix multiplication and vector subtraction, resulting in extremely low computational complexity, making it perfectly suitable for real-time processing by edge computing nodes. The projection distance exhibits good robustness to normal disturbances caused by load changes, while demonstrating high sensitivity to feature space deviations caused by winding state anomalies, effectively achieving adaptive partitioning between stable and unstable operating regions. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0023] Figure 1 This is a flowchart of an outdoor power distribution transformer condition monitoring system based on edge computing; Figure 2 This is a flowchart of the impact disturbance data identification and deletion process; Figure 3 This is a flowchart of the time-domain solution module extracting excitation mutation features; Figure 4 This is the process flowchart for the periodic segmentation module; Figure 5 This is a schematic diagram of the shell vibration waveform, voltage amplitude envelope, and positive zero-crossing moment; Figure 6 It is the curve of the load dynamic response vector projection distance and the result of the unstable state determination. Detailed Implementation

[0024] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0025] See Figure 1This invention provides an outdoor distribution transformer condition monitoring system based on edge computing. The system is deployed on an edge computing node on the distribution transformer side and includes a data acquisition module, a time-domain calculation module, a period segmentation module, a correlation alignment module, and a condition determination module. The data acquisition module collects initial monitoring data during the operation of the distribution transformer, including voltage waveforms, current waveforms, and shell vibration waveforms. The time-domain calculation module performs instantaneous power calculations on the voltage and current waveforms in the time domain to obtain a power factor angle sequence, and performs differential operations on the power factor angle sequence to obtain excitation mutation characteristics. The period segmentation module segments the voltage waveform into each power frequency cycle based on the zero-crossing time of the shell vibration waveform, obtaining the voltage amplitude envelope and current harmonic sub-band within each power frequency cycle. The correlation alignment module calculates the cross-correlation delay between the voltage amplitude envelope and the current harmonic sub-band, and aligns the cross-correlation delay with the excitation mutation characteristics in time to generate a load dynamic response vector. The status determination module is used to determine whether the distribution transformer is in an unstable operating area based on the projection distance of the load dynamic response vector in the preset winding status characteristic space, and outputs the status determination result.

[0026] Example 1: In specific implementation, please refer to Figure 2 The data acquisition module synchronously acquires three-phase voltage and current signals from the high-voltage and low-voltage sides of the distribution transformer via voltage and current transformers connected to the edge computing node. The primary side of the voltage transformer is connected in parallel to the three-phase terminals on both the high-voltage and low-voltage sides of the distribution transformer, while the secondary side output is connected to the analog input channel of the edge computing node. The primary side of the current transformer is connected in series to the three-phase lines on both the high-voltage and low-voltage sides of the distribution transformer, while the secondary side output is connected to the analog input channel of the edge computing node. The analog-to-digital conversion unit within the edge computing node synchronously samples the analog signals output from the voltage and current transformers at the same sampling frequency, set to 128 points per power frequency cycle, corresponding to 6400 samples per second at a 50Hz power system frequency. The synchronously acquired three-phase voltage and current signals are used as voltage and current waveforms, respectively.

[0027] Vibration waveforms of the distribution transformer tank are acquired using piezoelectric accelerometers attached to its surface. The piezoelectric accelerometers are installed at three locations: the center of the front, the center of the back, and the center of the side of the distribution transformer tank, with one sensor installed at each location. The piezoelectric accelerometers are attached to the tank surface using magnetic bases, with coupling agent applied between the magnetic bases and the tank surface to eliminate air gaps. The output signal of the piezoelectric accelerometers is connected to the input port of the charge amplifier on the edge computing node via a shielded cable. The charge amplifier converts the charge signal into a voltage signal, which is then converted to digital signal by the analog-to-digital conversion unit of the edge computing node at a sampling frequency of 128 points per power frequency cycle.

[0028] The edge computing node performs analog-to-digital conversion on the acquired three-phase voltage signals, three-phase current signals, and shell vibration waveforms output from three piezoelectric accelerometers. All converted signals are then aligned and stored at the same sampling frequency as a time-series matrix of initial monitoring data. Each row of the time-series matrix corresponds to a sampling time, and the columns sequentially contain six channels for the three-phase voltage signals, six channels for the three-phase current signals, and three channels for the shell vibration waveforms, totaling fifteen columns. During the alignment and storage process, the edge computing node uses a global timestamp counter to timestamp each sampling point with a timestamp accuracy of 1 microsecond, ensuring timing synchronization between the three-phase voltage signals, three-phase current signals, and shell vibration waveforms.

[0029] After collecting initial monitoring data, the edge computing nodes perform sliding window variance calculations on the voltage and current waveforms from the initial monitoring data. When calculating the voltage waveform variance, the sliding window length is set to the number of sampling points corresponding to 10 power frequency cycles, i.e., 1280 sampling points. The sliding window slides point-by-point along the time series of the voltage waveform, moving one sampling point at a time. At each sliding window position, the variance value of all voltage waveform data points within the sliding window is calculated, and the calculated variance value is recorded as the voltage variance value at the start of the corresponding window in the voltage variance sequence. When calculating the current waveform variance, the sliding window length is also set to 1280 sampling points corresponding to 10 power frequency cycles, and the sliding method is the same as for the voltage waveform variance calculation, obtaining the current variance sequence.

[0030] In practice, the first variance threshold is determined as follows: The sliding window variance of historical voltage waveform data of the distribution transformer under normal operating conditions is statistically analyzed, and the mean and standard deviation of the historical variance are calculated. The first variance threshold is set as the historical variance mean plus four times the standard deviation. The second variance threshold is determined as follows: The sliding window variance of historical current waveform data of the distribution transformer under normal operating conditions is statistically analyzed, and the historical variance mean and standard deviation are calculated. The second variance threshold is set as the historical variance mean plus four times the standard deviation. The statistical base for both the first and second variance thresholds is no less than 100 sets of historical data, with each set of historical data lasting 60 seconds.

[0031] The edge computing node detects the first moment when the voltage variance sequence exceeds a first variance threshold. Specifically, it compares each variance value in the voltage variance sequence with the first variance threshold point by point. When a voltage variance value greater than the first variance threshold is detected, the sampling time corresponding to the voltage variance value exceeding the first variance threshold is recorded, and the earliest sampling time among all sampling times that meet the condition is taken as the first moment. Simultaneously, it detects the second moment when the current variance sequence exceeds a second variance threshold. This is done by comparing each variance value in the current variance sequence with the second variance threshold point by point. When a current variance value greater than the second variance threshold is detected, the sampling time corresponding to the current variance value exceeding the second variance threshold is recorded, and the earliest sampling time among all sampling times that meet the condition is taken as the second moment.

[0032] The edge computing node calculates the time difference between the first and second moments, which is the absolute value of the difference between the timestamps of the first and second moments. This time difference is compared to a preset synchronization tolerance, which is set to the time length corresponding to two sampling points. Under the condition of a sampling frequency of 128 points per power frequency cycle, the preset synchronization tolerance is 2 divided by 6400 seconds, or 312.5 microseconds. When the time difference between the first and second moments is less than the preset synchronization tolerance, the edge computing node marks both moments as the impact disturbance moment. The impact disturbance moment represents the starting point when both the voltage and current waveforms are simultaneously subjected to external impact interference.

[0033] Edge computing nodes remove data segments of a preset length from the initial monitoring data before and after the impact disturbance time. The preset length is set to the number of sampling points corresponding to 5 power frequency cycles before and after the impact disturbance time, i.e., 640 sampling points before and after. The deletion operation removes all rows in the time series matrix of the initial monitoring data whose timestamps fall within the range of the impact disturbance time minus 640 sampling points to the impact disturbance time plus 640 sampling points from the time series matrix. If there are multiple impact disturbance times, the deletion operation is performed separately for each impact disturbance time, and the deletion intervals are merged when the deletion ranges overlap.

[0034] The remaining initial monitoring data after deleting the impact disturbance data segment is used as the updated initial monitoring data, which is then used in subsequent instantaneous power calculation steps. The updated initial monitoring data is saved as a new time series matrix with the same column structure as the original time series matrix. The sampling frequency of the new time series matrix remains unchanged, the temporal continuity is broken at the deletion interval, and the data segments are stored contiguously in their original chronological order.

[0035] Example 2: In specific implementation, please refer to Figure 3 The time-domain solution module reads the time series matrix after impact disturbance segment deletion processing from the storage unit of the edge computing node. Each column of the time series matrix corresponds to the data sequence of one acquisition channel. The time-domain solution module extracts the A-phase voltage signal column, B-phase voltage signal column, C-phase voltage signal column, A-phase current signal column, B-phase current signal column, and C-phase current signal column from the time series matrix respectively.

[0036] The time-domain solution module multiplies the voltage signal of each phase in the time series matrix with the current signal of the corresponding phase point by point. For phase A, it multiplies the voltage value at the m-th sampling point of the phase A voltage signal sequence. The time-domain solution module for the current value at the m-th sampling point of the phase A current signal sequence multiplies the voltage signal of each phase in the time series matrix with the current signal of the corresponding phase point by point. For phase A, the voltage value at the m-th sampling point of the phase A voltage signal sequence is... The current value at the m-th sampling point of the A-phase current signal sequence Multiplying them together yields the m-th instantaneous power value of the instantaneous power sequence of phase A. Here, the subscript A represents phase A, the subscript m represents the sampling point number, where m ranges from 1 to N, and N is the total number of rows in the time series matrix. The same operation is performed on phases B and C to obtain the instantaneous power sequences of phase B and phase C. The mathematical expression of the three-phase instantaneous power sequences is: multiplying them yields the m-th instantaneous power value of the phase A instantaneous power sequence. In this context, subscript A represents phase A, subscript m represents the sampling point number, where m ranges from 1 to N, and N is the total number of rows in the time series matrix. The same operation is performed on phases B and C to obtain the instantaneous power sequences of phases B and C. The mathematical expression of the three-phase instantaneous power sequences is:

[0037] in, Indicates phase identification, The value range of is {A,B,C}, which represent phase A, phase B and phase C respectively; express The instantaneous voltage value of the phase voltage signal at sampling point m, in volts; express The instantaneous current value of the phase current signal at sampling point m, in amperes; express The instantaneous power value at the sampling point with serial number m, in watts; m represents the sampling point number, and the value of m ranges from 1 to N, where N is the total number of rows in the time series matrix.

[0038] The time-domain solution module performs low-pass filtering on the instantaneous power sequence of each phase to extract the active power envelope. The low-pass filter employs a second-order Butterworth infinite impulse response filter with a cutoff frequency set to 5 Hz. The filter coefficients are obtained from the s-domain to the z-domain using the bilinear transform method. The low-pass filter performs point-by-point filtering on the instantaneous power sequence, removing fluctuation components above 5 Hz and outputting a smooth active power envelope sequence.

[0039] The time-domain solution module performs Hilbert transforms on the three-phase voltage signals. The Hilbert transform is implemented using a finite impulse response (FIR) Hilbert transform filter of order 64, with filter coefficients optimized using the Parks-McClellan algorithm. The A-phase voltage signal sequence is input into the Hilbert transform filter, outputting an A-phase quadrature voltage component sequence. The voltage value at each sampling point in the A-phase quadrature voltage component sequence forms a pair of quadrature components with the corresponding sampling point in the original A-phase voltage signal sequence. The same Hilbert transform process is performed on the B-phase and C-phase voltage signal sequences to obtain the B-phase and C-phase quadrature voltage component sequences, respectively.

[0040] The time-domain solution module calculates the instantaneous phase angle of each phase voltage based on the three-phase voltage signals and their corresponding orthogonal voltage components. For the m-th sampling point in the A-phase voltage signal sequence, the instantaneous phase angle of the A-phase voltage is... By calculating the quadrature voltage components of phase A With phase A voltage signal The arctangent value is obtained by the following formula: At the same time, according to and The sign of the voltage is determined by quadrant selection, and the arctangent result is expanded to the range of 0 to 2π. The instantaneous phase angles of phase B and phase C voltages are calculated using the same method.

[0041] The time-domain calculation module calculates the power factor angle for each power frequency cycle based on the time difference between the zero-crossing moment of each phase current signal and the instantaneous phase angle of the corresponding phase voltage. Zero-crossing detection employs a positive slope zero-crossing detection method. Specifically, it searches for two adjacent sampling points in the current signal sequence that satisfy the condition that the current value of the previous sampling point is negative and the current value of the next sampling point is positive. The precise zero-crossing moment is determined through linear interpolation between these two sampling points. Each positive slope zero-crossing moment of the A-phase current signal is used as a reference moment to begin a power frequency cycle, and the instantaneous phase angle of the A-phase voltage is read at this reference moment. The instantaneous phase angle of phase A voltage under ideal sinusoidal conditions corresponds to the zero-crossing phase angle of the phase A voltage sine wave, and the phase A power factor angle. The instantaneous phase angle of phase A voltage Taking the negative value of the difference from 2π, that is... This ensures that the power factor angle is 0 for purely resistive loads and positive for inductive loads. A power factor angle is calculated for each power frequency cycle, and these are arranged in chronological order to form the A-phase power factor angle sequence. The B-phase and C-phase power factor angle sequences are formed using the same method.

[0042] The time-domain solution module performs a first-order forward difference operation on the power factor angle sequence of each phase. For the q-th power factor angle in the power factor angle sequence of phase A... and the (q+1)th power factor angle Perform a first-order forward difference and calculate the difference angle. Where q ranges from 1 to Q-1, and Q is the length of the power factor angle sequence. First-order forward differential operations are performed simultaneously on the B-phase power factor angle sequences and the C-phase power factor angle sequences to obtain the A-phase differential angle sequences, B-phase differential angle sequences, and C-phase differential angle sequences.

[0043] The preset angle threshold is an adaptive threshold dynamically calculated based on the rated capacity and short-circuit impedance percentage of the distribution transformer. The adaptive threshold is calculated as follows: Let the rated capacity of the distribution transformer be... The unit is megavolt-amperes, and the short-circuit impedance percentage is... Then the adaptive threshold for: .

[0044] in, This indicates the rated capacity of the distribution transformer, in megavolt-amperes (MVA), and the value is taken from the nameplate parameters of the distribution transformer. This indicates the percentage of short-circuit impedance of the distribution transformer, expressed as a percentage, and the value is taken from the factory test report of the distribution transformer. The preset angle threshold is expressed in radians; the constant 0.5 is the reference angle threshold coefficient, which is determined by statistical analysis of the slope of the linear segment of the typical core excitation characteristic curve of the distribution transformer. This coefficient is obtained by offline testing based on the variation law of excitation inductance of distribution transformers of different capacities. The calibration method is to measure the steady-state fluctuation range of the phase difference between the excitation current and voltage when the distribution transformer is closed under no-load conditions, and take the phase difference change corresponding to the half peak value of the fluctuation range as the reference angle threshold coefficient.

[0045] In practical implementation, the time-domain solution module compares the absolute value of each difference point in the A phase difference angle sequence, B phase difference angle sequence, and C phase difference angle sequence with a preset angle threshold. Perform a comparison one by one. For the q-th difference angle in the A phase difference angle sequence... ,when Greater than At that time, the difference angle The corresponding location is marked as the A-phase excitation abrupt change point. The amplitude of the A-phase excitation abrupt change point is recorded as follows: The value, and the direction of the excitation mutation point of phase A are recorded as follows. The symbols are positive for a positive step in the power factor angle and negative for a negative step. The same comparison and marking method is used for the B-phase and C-phase differential angle sequences to obtain the B-phase and C-phase excitation abrupt change points. The amplitude and direction of the three-phase excitation abrupt change points together constitute the excitation abrupt change characteristics. The time-domain solution module stores the excitation abrupt change characteristics in the form of a structure array. Each element in the structure array contains four fields: the phase of the excitation abrupt change point, the sampling point number, the amplitude, and the direction. The structure array is arranged in ascending order of the sampling point number and written to the shared storage area of ​​the edge computing node for subsequent modules to read.

[0046] Example 3: In specific implementation, please refer to Figure 4 The period segmentation module reads the time series matrix and the shell vibration waveforms output by three piezoelectric accelerometers from the storage unit of the edge computing node. The period segmentation module selects the channel with the largest signal amplitude among the three shell vibration waveforms as the main vibration signal channel. When detecting the main vibration signal channel, the period segmentation module scans each sampling point in the shell vibration waveform sequence point by point, identifying the positive and negative zero-crossing moments. The positive zero-crossing moment is defined as the moment when the amplitude of the previous sampling point in the shell vibration waveform sequence is negative and the amplitude of the next sampling point is positive. The precise positive zero-crossing moment is determined by linear interpolation between two adjacent sampling points. The linear interpolation method for the positive zero-crossing moment is as follows: assuming the amplitude of the previous sampling point is... The amplitude of the next sampling point is The timestamp of the previous sampling point is The timestamp of the next sampling point is Precise positive zero crossing time The calculation formula is expressed as follows:

[0047] in, This represents the amplitude of the previous sampling point in the positive zero-crossing interval of the shell vibration waveform sequence, in millivolts. This represents the amplitude of the next sampling point in the positive zero-crossing interval of the shell vibration waveform sequence, in millivolts. The timestamp of the previous sampling point in the shell vibration waveform sequence, in seconds; The timestamp, in seconds, represents the last sampling point in the shell vibration waveform sequence. This indicates the exact positive zero-crossing time calculated, in seconds.

[0048] The negative zero-crossing moment is defined as the moment when the amplitude of the previous sampling point in the shell vibration waveform sequence is positive and the amplitude of the next sampling point is negative. The precise negative zero-crossing moment is determined by linear interpolation between two adjacent sampling points, and the calculation method is the same as that for the positive zero-crossing moment.

[0049] The period segmentation module uses the interval between two adjacent zero-crossing moments in the same direction as a single oscillation period interval. For positive zero-crossing moments, the time interval between the r-th positive zero-crossing moment and the (r+1)-th positive zero-crossing moment is defined as a single oscillation period interval, where r is the sequence number of the positive zero-crossing moment. For negative zero-crossing moments, the time interval between the r-th negative zero-crossing moment and the (r+1)-th negative zero-crossing moment is defined as a single oscillation period interval. The period segmentation module records the start and end times of each oscillation period interval, with the start time being the r-th zero-crossing moment in the same direction and the end time being the (r+1)-th zero-crossing moment in the same direction. The start and end times of the oscillation period interval are in seconds and accurate to the microsecond level.

[0050] The period segmentation module truncates the voltage waveforms in the time series matrix according to the start and end times of the vibration period intervals. Specifically, for the A-phase, B-phase, and C-phase voltage signal columns in the time series matrix, all sampling points whose timestamps fall within the start and end times of each vibration period interval are extracted, forming the corresponding A-phase, B-phase, and C-phase voltage segmented waveforms for that vibration period interval. For each vibration period interval, the three-phase voltage signals are truncated independently, and the number of sampling points for each voltage segmented waveform after truncation varies, depending on the actual length of the vibration period interval.

[0051] The period segmentation module performs absolute value calculations on each voltage segment waveform. This calculation takes the absolute value of the amplitude at each sampling point in the voltage segment waveform, converting the voltage value of the negative half-cycle to a positive value. The voltage segment waveform after the absolute value calculation is then processed using a sliding peak hold method. This method employs a sliding window, with the window length set to one-quarter of the number of sampling points corresponding to one power frequency cycle of the voltage segment waveform sampling rate, i.e., 32 sampling points. The sliding window slides point-by-point along the absolute value sequence, outputting the maximum absolute value of all sampling points within the sliding window at each window position as the peak hold output value for that window position. Adjacent identical values ​​in the sliding peak hold output value sequence are merged, retaining only the peak hold value at the moment of change to obtain the voltage amplitude envelope.

[0052] In practical implementation, the period segmentation module truncates the current waveforms in the time series matrix according to the start and end times of the vibration period intervals. The truncation operation is performed separately for the A-phase, B-phase, and C-phase current signal trains. For each vibration period interval, all sampling points whose timestamps fall within the start and end times are extracted, forming the corresponding A-phase, B-phase, and C-phase current segmented waveforms for that vibration period interval. The voltage segmented waveforms and current segmented waveforms share the same set of start and end times for the vibration period intervals, ensuring the timing correspondence between them.

[0053] The period segmentation module performs multi-order bandpass filtering on each current segment waveform. This multi-order bandpass filtering is implemented using a cascaded finite impulse response bandpass filter bank, which includes a fundamental bandpass filter and multiple higher harmonic bandpass filters. The center frequency of the fundamental bandpass filter is set to 50 Hz, corresponding to the power system's operating frequency. The center frequencies of each higher harmonic bandpass filter are integer multiples of the fundamental frequency of the current segment waveform, specifically 100 Hz for the 2nd harmonic, 150 Hz for the 3rd harmonic, 200 Hz for the 4th harmonic, 250 Hz for the 5th harmonic, 300 Hz for the 6th harmonic, and 350 Hz for the 7th harmonic. The passband width of each bandpass filter is determined by the quality factor, defined as the ratio of the bandpass filter's center frequency to its bandwidth.

[0054] The quality factor of each bandpass filter decreases as the harmonic order increases. The quality factor of the fundamental bandpass filter is set to 20, the quality factor of the second harmonic bandpass filter is set to 16, the quality factor of the third harmonic bandpass filter is set to 14, the quality factor of the fourth harmonic bandpass filter is set to 12, the quality factor of the fifth harmonic bandpass filter is set to 10, the quality factor of the sixth harmonic bandpass filter is set to 8, and the quality factor of the seventh harmonic bandpass filter is set to 6. The specific values ​​of the aforementioned quality factor are determined based on the following principles: First, using the theoretical amplitude ratio of each harmonic in the current waveform under rated load of the distribution transformer as a benchmark, when the fundamental amplitude is normalized to 1, the typical amplitude ranges of the 2nd to 7th harmonics are 0.05~0.15, 0.03~0.10, 0.02~0.08, 0.01~0.06, 0.005~0.04, and 0.005~0.03, respectively; Second, the passband width is determined based on the signal-to-noise ratio of each harmonic. The passband width is inversely proportional to the harmonic amplitude; that is, the lower the harmonic amplitude, the wider the passband to obtain sufficient energy accumulation. Secondly, the theoretical values ​​are corrected through offline experimental calibration. The calibration method involves superimposing harmonic components with known amplitudes and frequencies onto a standard signal source, adjusting the quality factor of each order bandpass filter to minimize the harmonic amplitude extraction error. During calibration, a fundamental frequency of 50Hz and a sampling rate of 6400Hz are used as fixed conditions. Finally, the calibration results are rounded to even numbers to facilitate the symmetrical design of the digital filter coefficients. The reason for the decreasing quality factor is that the frequency distribution range of higher harmonics is larger relative to their center frequency, and the amplitude of higher harmonics is usually lower than that of the fundamental frequency. A wider passband can compensate for the insufficient signal-to-noise ratio caused by the low energy proportion of higher harmonics in the current waveform.

[0055] For distribution transformers of different models and capacities, the applicability of the above quality factor values ​​follows these adjustment rules: When the rated capacity of the distribution transformer is between 100kVA and 10000kVA and the short-circuit impedance percentage is between 4% and 10%, the above quality factor settings can be directly applied without adjustment. This range covers the vast majority of 10kV and 35kV outdoor distribution transformers. When the rated capacity of the distribution transformer exceeds the above range, the adjustment amount of the quality factor is determined by the formula...

[0056] Calculation, where The baseline quality factor (i.e., the set value corresponding to each of the aforementioned harmonics). This is the actual rated capacity (in kVA). With a base capacity of 1000kVA, the index of 0.15 is determined by statistical regression analysis of the harmonic spectra of transformers of different capacities; when the short-circuit impedance percentage of the distribution transformer exceeds the range of 4% to 10%, the adjustment amount of the quality factor is determined by the formula.

[0057] Calculation, where This represents the percentage of the actual short-circuit impedance. The baseline short-circuit impedance percentage is 6%, and the exponent 0.1 is determined by the influence of short-circuit impedance on harmonic conduction characteristics. After the above adjustments, the quality factors of each harmonic still need to be rounded to the nearest even number to ensure the symmetry of the digital filter coefficients.

[0058] Each bandpass filter is implemented using a finite impulse response (FIR) structure of the same order, with the FIR order set to 256. The coefficients of the FIR bandpass filter are designed using the window function method, with a Hamming window selected. The design process is as follows: The amplitude-frequency response of the ideal bandpass filter is calculated based on the target center frequency and quality factor. The amplitude-frequency response of the ideal bandpass filter has an amplitude of 1 in the target passband and 0 in the stopband. The upper and lower cutoff frequencies in the passband are the center frequency minus and plus half the center frequency divided by the quality factor, respectively. An inverse discrete Fourier transform is performed on the ideal amplitude-frequency response to obtain an infinitely long impulse response sequence. The infinitely long impulse response sequence is truncated using a Hamming window with 257 coefficients to obtain the coefficients of the 256th order FIR bandpass filter. Each segmented current waveform passes sequentially through a fundamental bandpass filter and various higher harmonic bandpass filters. The fundamental bandpass filter outputs the fundamental component amplitude sequence, and each higher harmonic bandpass filter outputs the corresponding harmonic component amplitude sequence.

[0059] The period segmentation module arranges the amplitude sequences of each higher harmonic component in ascending order of harmonic order, and combines the amplitude sequences of the 2nd, 3rd, 4th, 5th, 6th, and 7th harmonic components column-wise to form current harmonic sub-bands. The current harmonic sub-band is a multi-column data structure, where each row corresponds to a sampling time, and each column corresponds to the amplitude of a specific harmonic component. For each vibration period interval, the period segmentation module outputs the A-phase, B-phase, and C-phase current harmonic sub-bands respectively, and stores them in association with the voltage amplitude envelope of the corresponding phase according to the vibration period interval number. These sub-bands are then written to the shared storage area of ​​the edge computing nodes for relevant alignment modules to read.

[0060] See Figure 5 In the graph, the horizontal axis represents time in seconds, with a display duration of 0 to 0.1 seconds. The left side of the vertical axis corresponds to the amplitude of the shell vibration waveform in millivolts (mV), and the right side corresponds to the voltage amplitude envelope in volts (V). The solid curve in the graph represents the shell vibration waveform, the dashed curve represents the voltage amplitude envelope, and the dotted vertical lines represent the positive zero-crossing moments of the shell vibration waveform.

[0061] As can be seen from the solid line curve of the shell vibration waveform, it exhibits periodic oscillations with amplitude fluctuating within a range of approximately ±1.5 millivolts. The waveform also shows some random disturbances, reflecting the actual situation of the distribution transformer tank shell being affected by mechanical vibration. The dashed line curve of the voltage amplitude envelope shows periodic changes, with a peak amplitude of approximately 320 volts and a trough of approximately 70 volts, indicating the power envelope characteristics of the voltage signal.

[0062] The vertical dotted lines corresponding to the positive zero-crossing moments are evenly distributed on the time axis in the figure, with approximately equal intervals, consistent with the characteristic of a 50 Hz power frequency cycle of about 20 milliseconds. The cycle segmentation module divides the voltage waveform into multiple oscillation cycle intervals based on these positive zero-crossing moments, providing a time reference for subsequent time-series analysis of the voltage amplitude envelope and current harmonic sub-bands.

[0063] The voltage amplitude envelope in the figure is formed after processing the absolute value and sliding peak value of the segmented voltage waveform. The dashed envelope curve closely follows the local peak value of the voltage waveform, highlighting the voltage amplitude change trend within each vibration cycle. The positive zero-crossing moment of the shell vibration waveform corresponds closely to the voltage amplitude envelope change cycle, ensuring that the vibration cycle interval intercepted by the cycle segmentation module has accurate time boundaries.

[0064] Example 4: In practice, the alignment module reads the voltage amplitude envelope and current harmonic sub-bands, which are associated and stored according to the vibration period interval number, from the shared storage area of ​​the edge computing node's period segmentation module. Simultaneously, it reads the excitation mutation characteristic structure array output by the time-domain solution module from the shared storage area. The alignment module processes each vibration period interval sequentially, following the chronological order of the vibration period intervals.

[0065] The relevant alignment module performs mean-removal processing on the read A-phase voltage amplitude envelope sequence and the fundamental amplitude sequence of the A-phase current harmonic sub-band. The arithmetic mean of the amplitudes of all sampled points in the A-phase voltage amplitude envelope sequence is calculated, and this arithmetic mean is subtracted from the amplitude of each sampled point in the A-phase voltage amplitude envelope sequence to obtain a zero-mean voltage envelope sequence. Similarly, the arithmetic mean of the amplitudes of all sampled points in the fundamental amplitude sequence of the A-phase current harmonic sub-band is calculated, and this arithmetic mean is subtracted from the amplitude of each sampled point in the fundamental amplitude sequence to obtain a zero-mean fundamental amplitude sequence.

[0066] The alignment module keeps the zero-mean fundamental amplitude sequence fixed on the time axis and slides the zero-mean voltage envelope sequence point by point within a preset delay search range. The preset delay search range is set to the number of sampling points corresponding to the period from -4 / 5 to 4 / 5 of the power frequency cycle. Using 128 sampling points per power frequency cycle, the preset delay search range is from -32 to 32 sampling points, totaling 65 sliding displacement positions. At each sliding displacement, the sum of the products of the corresponding points in the zero-mean voltage envelope sequence and the zero-mean fundamental amplitude sequence is calculated. This sum of products is then divided by the product of the standard deviations of the zero-mean voltage envelope sequence and the zero-mean fundamental amplitude sequence to obtain the normalized cross-correlation function value sequence.

[0067] Let the sliding displacement d represent the number of sliding sampling points of the zero-mean voltage envelope sequence relative to the zero-mean fundamental amplitude sequence. The value of d is an integer ranging from -32 to +32. The normalized cross-correlation function value The calculation method for d is as follows: multiply the amplitude of the sampling point with index n in the zero-mean voltage envelope sequence with the amplitude of the sampling point with index n minus d in the zero-mean fundamental amplitude sequence, accumulate all product values ​​within the overlapping sampling point range of the two sequences, divide the sum by the number of overlapping sampling points, and then divide by the product of the standard deviation of the zero-mean voltage envelope sequence and the standard deviation of the zero-mean fundamental amplitude sequence. The standard deviation of the zero-mean voltage envelope sequence is the square root of the mean of the squares of the amplitudes of all sampling points in the zero-mean voltage envelope sequence, and the standard deviation of the zero-mean fundamental amplitude sequence is the square root of the mean of the squares of the amplitudes of all sampling points in the zero-mean fundamental amplitude sequence.

[0068] In practice, the alignment module iterates through all 65 normalized cross-correlation function values ​​in the sequence and finds the sliding displacement corresponding to the maximum value. The found sliding displacement is recorded as a candidate delay. .

[0069] The correlation alignment module performs three-point Gaussian interpolation fitting on the normalized cross-correlation function values ​​of the two discrete points before and after the candidate delay value. Candidate delay values ​​are then selected. Normalized cross-correlation function at location Candidate delay Normalized cross-correlation function value at the minus 1 position -1 and candidate delay Normalized cross-correlation function value at position plus 1 +1, for a total of three discrete points. The calculation method for three-point Gaussian interpolation fitting is as follows: Let the Gaussian function be of the form... Where parameter a represents the peak height of the Gaussian function, which is dimensionless; parameter b represents the width factor of the Gaussian function, in units of the negative square of the number of sampling points; parameter The delay corresponding to the peak value of the Gaussian function is expressed in units of sampling points; d represents the sliding displacement, also expressed in units of sampling points. Taking the natural logarithm over the three discrete points transforms the Gaussian function into a quadratic parabola. Then, by substituting the logarithmic form of the three discrete points and solving the system of equations consisting of the three equations, the delay corresponding to the peak value of the Gaussian function can be obtained. The calculation expression:

[0070] in, This represents the candidate delay amount, in units of sampling points, and its value is the sliding displacement that maximizes the value of the normalized cross-correlation function. The sliding displacement is The normalized cross-correlation function value when reduced by 1, dimensionless; The sliding displacement is The normalized cross-correlation function value at time is dimensionless; The sliding displacement is The normalized cross-correlation function value when 1 is added is dimensionless. This represents the base delay obtained after fitting with three-point Gaussian interpolation, with units of the number of sampling points. It is a real number.

[0071] In practical implementation, the basic delay amount of subsampling precision is calculated. Then, the relevant alignment module will Convert to time units, the conversion method is as follows Divide by the sampling frequency, where the sampling frequency is 6400 samples per second, to obtain the base delay in seconds. The relevant alignment module uses this base delay. The corresponding normalized cross-correlation function value is used as the confidence index for the basic delay. The confidence index is set to... Substitute into Gaussian function The obtained peak value peak It is obtained by taking the exponent of the constant term in the quadratic parabola fitting. The confidence index ranges from 0 to 1, with the value closer to 1 indicating a higher degree of confidence in the match.

[0072] The relevant alignment module repeats all the above steps for the voltage amplitude envelopes of phase B and phase C and the fundamental amplitude sequences of the corresponding phase current harmonic subbands to obtain the basic delay of phase B and the confidence index of the basic delay of phase B, as well as the basic delay of phase C and the confidence index of the basic delay of phase C.

[0073] In practical implementation, the alignment module performs cross-correlation calculations on the voltage amplitude envelope within the same vibration period interval with the amplitude sequences of each harmonic in the current harmonic sub-band. Cross-correlation calculations are performed on the voltage amplitude envelope of phase A with the 2nd, 3rd, 4th, 5th, 6th, and 7th harmonic amplitude sequences of the phase A current harmonic sub-band. The steps of the cross-correlation calculation are the same as those for the fundamental amplitude sequence, including mean removal, point-by-point sliding calculation of the normalized cross-correlation function value sequence, finding the candidate delay corresponding to the maximum value, and three-point Gaussian interpolation fitting to obtain the harmonic delay corresponding to each harmonic with sub-sampling precision. The 2nd harmonic delay is calculated for the 2nd harmonic amplitude sequence, the 3rd harmonic delay is calculated for the 3rd harmonic amplitude sequence, and so on, resulting in a total of 6 harmonic delays. The same processing is performed on phases B and C, yielding their respective 6 harmonic delays.

[0074] The relevant alignment module uses the difference between the basic delay of phase A and the delay of each harmonic of phase A within the same vibration period interval as the additional hysteresis sequence of phase A. The additional hysteresis sequence contains six elements: the additional hysteresis of the 2nd, 3rd, 4th, 5th, 6th, and 7th harmonics, respectively. Each element is calculated by subtracting the corresponding harmonic delay from the basic delay. The corresponding additional hysteresis sequences are calculated for phases B and C respectively.

[0075] In practical implementation, the alignment module performs time-series matching between the abrupt change time of the excitation mutation feature and the basic delay amount. The specific time-series matching method is as follows: the time point of the sampling point number of each excitation mutation point in the excitation mutation feature structure array is matched with the time point of the basic delay amount within each vibration cycle interval. If the sampling point number of an excitation mutation point falls within the range of a vibration cycle interval, a matching relationship is established between the sampling point number of that excitation mutation point and the basic delay amount of that vibration cycle interval. Within the vibration cycle interval where the basic delay amount undergoes a step change, the alignment module extracts the change value of the basic delay amount before and after the step change. The method for identifying the step change in the basic delay amount is as follows: the basic delay amount of the w-th vibration cycle interval is compared with the basic delay amount of the w+1-th vibration cycle interval. If the absolute difference between the two basic delay amounts is greater than a preset step threshold, it is determined that a step change has occurred in the basic delay amount. The preset step threshold is set to the time length corresponding to two sampling points, i.e., 312.5 microseconds. The value before the step is the basic delay of the w-th oscillation period interval, and the value after the step is the basic delay of the w+1-th oscillation period interval. The change in the basic delay is the value after the step minus the value before the step.

[0076] The alignment module concatenates the basic delay, the additional hysteresis sequence, and the change in basic delay into a load dynamic response vector. For phase A, the load dynamic response vector is a 13-dimensional vector, with the following components in order: phase A basic delay, 2nd harmonic additional hysteresis, 3rd harmonic additional hysteresis, 4th harmonic additional hysteresis, 5th harmonic additional hysteresis, 6th harmonic additional hysteresis, 7th harmonic additional hysteresis, phase B basic delay, phase C basic delay, the average of the six elements of the phase B additional hysteresis sequence, the average of the six elements of the phase C additional hysteresis sequence, the change in phase A basic delay, and the average of the changes in phase B and phase C basic delay. Each load dynamic response vector corresponds to a specific vibration period interval and is stored in association with the start and end times of that interval. The alignment module assembles the load dynamic response vectors of all vibration period intervals into a load dynamic response vector sequence in chronological order and writes it to the shared storage area of ​​the edge computing node for the state determination module to read.

[0077] Example 5: In practical implementation, the state determination module reads the load dynamic response vector sequence output by the relevant alignment module from the shared storage area of ​​the edge computing node. Each load dynamic response vector is a 13-dimensional vector and is associated with a specific vibration period interval. The state determination module reads the winding state reference matrix pre-stored in the non-volatile storage unit of the edge computing node. Each row of the winding state reference matrix corresponds to a load dynamic response reference vector under normal operating conditions, and each column corresponds to a component of the load dynamic response vector. The construction process of the winding state reference matrix is ​​completed during the on-site calibration stage before the distribution transformer is put into operation: under multiple operating conditions of no-load operation and rated load operation of the distribution transformer, no less than 500 load dynamic response vectors are continuously collected and calculated. After manual confirmation that the distribution transformer state corresponding to these vectors is in normal operating condition, these vectors are arranged by row to form a winding state reference matrix with 500 rows and 13 columns. The winding state reference matrix is ​​stored in double-precision floating-point format and write-protected to prevent accidental modification.

[0078] The state determination module calculates the orthogonal projection of the load dynamic response vector onto the eigenspace spanned by the winding state reference matrix. The process of calculating the orthogonal projection involves performing singular value decomposition on the winding state reference matrix to obtain the left singular matrix, the singular value diagonal matrix, and the right singular matrix. The winding state reference matrix is ​​denoted as matrix [missing information]. ,matrix The matrix has 500 rows by 13 columns. Singular value decomposition transforms the matrix... Decomposed into , where the matrix A left singular matrix with 500 rows by 500 columns. It is a singular value diagonal matrix with 500 rows and 13 columns. It is a right singular matrix with 13 rows by 13 columns. Representation matrix The transpose of .

[0079] Singular value diagonal matrix The number of non-zero singular values ​​is equal to the matrix The rank of the singular values ​​is denoted as R. The R singular values ​​are arranged in descending order in the singular value diagonal matrix. On the diagonal, the first The singular values ​​are denoted as ,in The value range is from 1 to R. The state determination module determines the state based on a preset energy retention threshold. From the singular value diagonal matrix Before being selected The largest singular value, The determination of satisfies the following conditions:

[0080] in, Represents a singular valued diagonal matrix The Middle The singular value, that is, the nth singular value after sorting in descending order. One non-negative real number; Represents a singular valued diagonal matrix The total number of non-zero singular values ​​of the matrix, i.e. rank; This represents a preset energy retention threshold, with a value ranging from 0.95 to 0.99. In this embodiment... The value is set to 0.98, which means that 98% of the matrix energy is retained. This value is based on the premise that while ensuring that the feature subspace can cover most of the normal operation state change patterns, singular vectors corresponding to noise components with an energy ratio of less than 2% are removed. This represents the number of singular values ​​ultimately selected, and is also the number of columns in the orthogonal basis matrix of the eigenspace. Its value is determined by the smallest integer that satisfies the energy retention threshold condition.

[0081] The state determination module will use the left singular matrix Middle and front The first maximal singular value corresponding to the previous The singular vectors on the left column are extracted into orthogonal basis matrices for the eigenspace, denoted as matrix . ,matrix The dimension is 500 rows multiplied by Columns. Matrices Each column is an orthogonal basis vector, and all orthogonal basis vectors are pairwise orthogonal and have a magnitude of 1.

[0082] The state determination module calculates the orthogonal basis matrix. transpose matrix The product of the load dynamic response vector and the current vibration period interval. The load dynamic response vector is denoted as a column vector. ,vector The length of the matrix is ​​13. The dimension is Rows multiplied by 500 columns, and a vector of length 13. They cannot be directly multiplied, therefore vectors need to be processed. Perform dimensional expansion processing on the vector. Expanded into a 500-dimensional vector with the same dimensions as the row vectors of the winding state reference matrix, the expansion method is to... The 13 components are copied to the corresponding positions in the 13 columns of the winding state reference matrix, and the remaining 487 components are set to zero. The expanded vector is denoted as vector. ,vector The length of the matrix is ​​500. Calculate the matrix. with vector The product of these two factors yields the projection coefficient vector, which is denoted as [the product of the two factors]. ,vector The length is ,vector The first in Each component Representing vectors In the matrix The The projection length along the directions of the orthogonal basis vectors.

[0083] The state determination module will use orthogonal basis matrices With projection coefficient vector Multiplying them yields the projected component vector, which is denoted as . ,vector The length is 500, and the calculation method is as follows: Projection component vector This represents the best approximation of the load dynamic response vector in the eigenspace spanned by the winding state reference matrix.

[0084] The status determination module will use the extended vector corresponding to the load dynamic response vector. Subtract the projected component vector The residual vector is obtained, and the residual vector is denoted as... ,vector The length is 500. Calculate the residual vector. The residual energy value is obtained by summing the squares of all elements in the vector. Using this residual energy value as the basis for calculating the projected distance, the square root of the residual energy value is taken to obtain the Euclidean norm of the residual vector, which is the projected distance.

[0085] The state determination module compares the projected distance with the state distance threshold pre-stored in the non-volatile storage unit of the edge computing node. The state distance threshold is determined as follows: under the normal operation state of the distribution transformer, the projected distance of 100 load dynamic response vectors on the characteristic subspace spanned by the winding state reference matrix is ​​repeatedly calculated. The 100 projected distance values ​​are arranged in ascending order, and the 95th percentile is selected as the state distance threshold and stored in the edge computing node.

[0086] When the projected distance is greater than the state distance threshold, the state determination module determines that the distribution transformer is in an unstable operating zone and outputs an unstable state flag to the state output interface of the edge computing node. The unstable state flag is a Boolean variable with a value of true. When the projected distance is less than or equal to the state distance threshold, the state determination module determines that the distribution transformer is in a stable operating zone and outputs a stable state flag to the state output interface of the edge computing node. The stable state flag is a Boolean variable with a value of false.

[0087] In practical implementation, after the state determination module determines that the distribution transformer is in an unstable operating zone and outputs an unstable state flag, it marks the load dynamic response vector corresponding to the unstable state flag as an abnormal response vector. The abnormal response vector is the load dynamic response vector of the vibration period interval that caused the unstable determination, and its data type and dimension are the same as the original load dynamic response vector.

[0088] The state determination module calculates the cosine of the angle between the abnormal response vector and each reference vector in the winding state reference matrix. The winding state reference matrix contains 500 reference vectors. For the first reference vector in the winding state reference matrix... The cosine of the angle between the reference vector and the extended vector of the abnormal response vector is calculated as follows: Perform a dot product operation on the first reference vector, and divide the dot product result by the magnitude of the extended vector of the anomaly response vector and the first reference vector. The product of the magnitudes of the reference vectors, where the magnitude of the extended vector of the anomaly response vector is the square root of the sum of the squares of the elements in the extended vector of the anomaly response vector. The magnitude of the first reference vector is the first... The cosine of the included angle is calculated as the square root of the sum of the squares of the elements in each of the reference vectors. The range of the cosine value obtained by this method is from -1 to 1.

[0089] The state determination module iterates through all 500 cosine values ​​of the included angle and selects the fault type label corresponding to the reference vector with the largest included angle cosine value as the preliminary fault type determination result. The fault type labels are pre-stored in the fault type mapping table of the edge computing nodes. The fault type mapping table maps each row of the winding state reference matrix to a fault type label. The fault type labels include "inter-turn short circuit," "phase-to-phase short circuit," "multi-point grounding of the core," "winding deformation," and "normal." The mapping relationship is determined by the factory test data and fault simulation test data of the distribution transformer. When the fault type label corresponding to the reference vector with the largest included angle cosine value is "normal," the preliminary fault type determination result is still retained as an unstable state, but no specific fault type is specified.

[0090] Based on the preliminary fault type determination result, the status determination module retrieves the spectral analysis templates pre-stored in the non-volatile storage unit of the edge computing node, which are associated with the fault type labels. The spectral analysis template is a set of typical response vector templates corresponding to each fault type. Each set contains 5 template vectors, calculated from historical data under typical fault conditions of the distribution transformer. Each template vector has the same 13-dimensional dimension as the load dynamic response vector. The matching degree is calculated as follows: the cosine similarity between the abnormal response vector and the 5 template vectors corresponding to the fault type label in the spectral analysis template is calculated, and the maximum value of the 5 cosine similarities is taken as the matching degree. When the matching degree is greater than the preset matching degree threshold of 0.85, the fault type label is output as a fault type confirmation result; when the matching degree is less than or equal to the preset matching degree threshold of 0.85, an "unconfirmed fault" label is output.

[0091] In practical implementation, after marking the abnormal response vector, the state determination module extracts voltage and current waveform segments within the abnormal time period corresponding to the abnormal response vector from the time series matrix. The abnormal time period is determined by the start and end times of the vibration period interval associated with the abnormal response vector. All rows in the time series matrix whose timestamps fall within the range of the start and end times are extracted, and the A-phase voltage signal column, B-phase voltage signal column, C-phase voltage signal column, A-phase current signal column, B-phase current signal column, and C-phase current signal column are extracted from them to form voltage and current waveform segments, respectively.

[0092] The state determination module performs Fast Fourier Transform (FFT) on both the voltage and current waveform segments. The FFT uses a radix-2 decimation-time algorithm, and the number of transformation points is selected as a power of 2 that is not less than the minimum number of sampling points for both the voltage and current waveform segments. Before the transformation, the voltage and current waveform segments are weighted using a Hanning window to suppress spectral leakage. The FFT outputs abnormal voltage and current spectra. The abnormal voltage spectrum includes the spectral results of the three-phase voltage, and the abnormal current spectrum includes the spectral results of the three-phase current.

[0093] The status determination module extracts the fundamental amplitude descent rate and unbalance coefficient from the abnormal voltage spectrum. The fundamental amplitude descent rate is calculated as follows: Obtain the fundamental amplitudes of phase A, phase B, and phase C from the abnormal voltage spectrum; take the arithmetic mean of the three-phase fundamental amplitudes; divide the difference between the rated phase voltage amplitude and the average of the three-phase fundamental amplitudes by the rated phase voltage amplitude to obtain the fundamental amplitude descent rate. The fundamental amplitude descent rate is a real number between 0 and 1. The unbalance coefficient is calculated as follows: Calculate the difference between the maximum and minimum values ​​of the fundamental amplitudes of phase A, phase B, and phase C; divide this difference by the arithmetic mean of the three-phase fundamental amplitudes to obtain the unbalance coefficient. The unbalance coefficient is a non-negative real number.

[0094] The state determination module extracts the zero-sequence current amplitude and negative-sequence current amplitude from the abnormal current spectrum. The zero-sequence current amplitude is calculated by superimposing the complex amplitudes of each harmonic in the current spectra of phases A, B, and C. The superposition method involves adding the complex amplitudes at the same frequency points in the three-phase current spectra, and dividing the magnitude of the superposition result at the fundamental frequency point by 3 to obtain the zero-sequence current amplitude, in amperes. The negative-sequence current amplitude is calculated by using a symmetric component transformation matrix to convert the fundamental complex amplitudes of phases A, B, and C into positive-sequence, negative-sequence, and zero-sequence components, respectively. The magnitude of the negative-sequence component is taken as the negative-sequence current amplitude, in amperes.

[0095] The state determination module incorporates the fundamental amplitude decay rate, unbalance coefficient, zero-sequence current amplitude, and negative-sequence current amplitude as supplementary feature components into the abnormal response vector to obtain an extended abnormal response vector. The extended abnormal response vector has 13 plus 4 equals 17 dimensions. The first 13 components are the same as the original abnormal response vector, while the 14th to 17th components are the fundamental amplitude decay rate, unbalance coefficient, zero-sequence current amplitude, and negative-sequence current amplitude, respectively. The state determination module uses the extended abnormal response vector as a new load dynamic response vector, re-executes orthogonal projection and projection distance calculation, and compares the calculated projection distance with a state distance threshold for a secondary determination. The secondary determination steps are the same as the initial determination steps. If the projection distance in the secondary determination is still greater than the state distance threshold, an unstable state marker is finally output along with the extended abnormal response vector.

[0096] See Figure 6In the figure, the horizontal axis represents the vibration period interval number, ranging from 0 to 500, and the vertical axis represents the dimensionless value of the projected distance. The solid curve reflects the changing trend of the projected distance of the load dynamic response vector corresponding to each vibration period interval calculated by the state determination module on the characteristic subspace spanned by the winding state reference matrix. The curve fluctuates relatively smoothly overall, with the projected distance mostly between 0.0 and 0.2, indicating that the operating state of the distribution transformer is basically stable.

[0097] The dashed line in the diagram represents a preset state distance threshold, approximately 0.285. This threshold is determined based on the statistical distribution of projected distances under normal operating conditions and is used to distinguish between stable and unstable operating zones. When the projected distance exceeds this threshold, the state determination module determines that the distribution transformer in the corresponding vibration cycle interval is in an unstable operating zone.

[0098] In the figure, the points marked with hollow squares represent vibration period intervals determined to be in an unstable state, and the corresponding projected distances all exceed the threshold of the dashed line. There are three obvious groups of unstable state markers in the figure, which are concentrated at vibration period intervals of approximately 50, 200, and 350. The peak projected distances at these locations significantly exceed the threshold, with the highest projected distance reaching approximately 0.42, which is significantly higher than the state distance threshold, indicating that the state of the distribution transformer windings is abnormal during this period.

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

Claims

1. An outdoor distribution transformer condition monitoring system based on edge computing, characterized in that, The system is deployed on the edge computing nodes on the distribution transformer side, including: The data acquisition module is used to collect initial monitoring data during the operation of the distribution transformer, including voltage waveform, current waveform, and shell vibration waveform. The time-domain solution module is used to perform instantaneous power solution on the voltage waveform and the current waveform in the time domain to obtain the power factor angle sequence, and to perform differential operation on the power factor angle sequence to obtain the excitation mutation characteristics. The period segmentation module is used to segment the voltage waveform into each power frequency period according to the zero-crossing time of the shell vibration waveform, so as to obtain the voltage amplitude envelope and current harmonic sub-band in each power frequency period; The correlation alignment module is used to calculate the cross-correlation delay between the voltage amplitude envelope and the current harmonic sub-band, and to align the cross-correlation delay with the excitation mutation characteristics in time to generate a load dynamic response vector. The status determination module is used to determine whether the distribution transformer is in an unstable operating area based on the projection distance of the load dynamic response vector in a preset winding status feature space, and output the status determination result.

2. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 1, characterized in that, The initial monitoring data collected during the operation of the distribution transformer includes: The voltage transformers and current transformers connected through the edge computing nodes synchronously collect the three-phase voltage signals and three-phase current signals of the high-voltage side and low-voltage side of the distribution transformer, which are used as the voltage waveform and the current waveform. The vibration waveform of the housing is collected by a piezoelectric accelerometer attached to the surface of the distribution transformer tank. The three-phase voltage signal, the three-phase current signal, and the shell vibration waveform are subjected to analog-to-digital conversion and stored as a time series matrix of the initial monitoring data with the same sampling frequency.

3. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 2, characterized in that, The instantaneous power calculation of the voltage and current waveforms in the time domain to obtain a power factor angle sequence, and the differential operation of the power factor angle sequence to obtain excitation mutation characteristics, includes: The voltage signal of each phase in the time series matrix is ​​multiplied point by point with the current signal of the corresponding phase to obtain the instantaneous power sequence, and the instantaneous power sequence is low-pass filtered to extract the active power envelope. The three-phase voltage signals are subjected to Hilbert transform to obtain the corresponding quadrature voltage components. The instantaneous phase angle of each phase voltage is calculated based on the three-phase voltage signals and the quadrature voltage components. Based on the time difference between the zero-crossing moment of each phase current signal and the instantaneous phase angle of the corresponding phase voltage, the power factor angle in each power frequency cycle is calculated and arranged in chronological order to form the power factor angle sequence. The power factor angle sequence is subjected to first-order forward difference to obtain a difference angle sequence. The difference points in the difference angle sequence whose absolute values ​​exceed a preset angle threshold are marked as excitation mutation points, and the amplitude and direction of the excitation mutation points are recorded as the excitation mutation features.

4. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 3, characterized in that, The preset angle threshold is an adaptive threshold dynamically calculated based on the rated capacity and short-circuit impedance percentage of the distribution transformer.

5. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 3, characterized in that, The step of dividing the voltage waveform into each power frequency cycle based on the zero-crossing time of the shell vibration waveform to obtain the voltage amplitude envelope and current harmonic sub-band within each power frequency cycle includes: The positive and negative zero-crossing times of the shell vibration waveform are detected, and the interval between two adjacent zero-crossing times in the same direction is taken as a vibration period interval. The start and end times of each vibration period interval are recorded. The voltage waveform in the time series matrix is ​​truncated according to the start and end times of the vibration period interval to obtain the voltage segment waveform corresponding to each vibration period interval. The absolute value calculation and sliding peak hold processing are performed on each of the voltage segment waveforms to obtain the voltage amplitude envelope; The current waveform in the time series matrix is ​​truncated according to the start and end times of the vibration period interval to obtain the current segment waveform corresponding to each vibration period interval. Then, multi-order bandpass filtering is performed on each current segment waveform to separate the fundamental component and multiple higher harmonic components. The amplitude sequence of the higher harmonic components is combined into the current harmonic subband.

6. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 5, characterized in that, When performing multi-order bandpass filtering on each of the current segment waveforms, the center frequency of each order bandpass filter is an integer multiple of the fundamental frequency of the current segment waveform, and the quality factor of each order bandpass filter decreases as the harmonic order increases.

7. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 5, characterized in that, The step of calculating the cross-correlation delay between the voltage amplitude envelope and the current harmonic sub-band, and aligning the cross-correlation delay with the excitation abrupt change characteristics in time to generate a load dynamic response vector includes: The voltage amplitude envelope and the fundamental amplitude sequence of the current harmonic subband within the same vibration period interval are cross-correlated to obtain the delay time corresponding to the peak value of the cross-correlation function, which is used as the basic delay amount within the vibration period interval. The voltage amplitude envelope within the same vibration period interval is cross-correlated with the harmonic amplitude sequence of each harmonic in the current harmonic subband to obtain the harmonic delay corresponding to each harmonic, and the difference between the basic delay and the delay of each harmonic is used as the additional hysteresis sequence. The abrupt change time of the excitation mutation feature is matched with the basic delay amount in time sequence. Within the vibration period interval where the basic delay amount undergoes a step change, the change value of the basic delay amount before and after the step change is extracted. The basic delay, the additional hysteresis sequence, and the change value of the basic delay are concatenated to form the load dynamic response vector.

8. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 7, characterized in that, The step of performing a cross-correlation operation between the voltage amplitude envelope and the fundamental amplitude sequence of the current harmonic sub-band within the same vibration period interval, and obtaining the delay time corresponding to the peak value of the cross-correlation function as the basic delay amount within the vibration period interval, specifically includes: The fundamental amplitude sequences of the voltage amplitude envelope and the current harmonic subband are respectively subjected to mean-removal processing to obtain a zero-mean voltage envelope sequence and a zero-mean fundamental amplitude sequence. The zero-mean fundamental amplitude sequence is kept fixed on the time axis, and the zero-mean voltage envelope sequence is slid point by point within a preset delay search range. At each sliding displacement, the sum of the products of the corresponding points of the zero-mean voltage envelope sequence and the zero-mean fundamental amplitude sequence is calculated and divided by the product of the standard deviations of the two sequences to obtain the normalized cross-correlation function value sequence. Traverse the normalized cross-correlation function value sequence, find the sliding displacement corresponding to the maximum value, and record the sliding displacement as a candidate delay value; Three-point Gaussian interpolation fitting is performed on the normalized cross-correlation function values ​​of the two discrete points before and after the location of the candidate delay. The displacement corresponding to the extreme point after interpolation is used as the basic delay with subsampling precision, and the normalized cross-correlation function value corresponding to the basic delay is used as the confidence index of the basic delay.

9. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 7, characterized in that, The determination of whether the distribution transformer is in an unstable operating region based on the projection distance of the load dynamic response vector in the preset winding state characteristic space, and the output of the state determination result, includes: Read the winding state reference matrix pre-stored in the edge computing node. The row vector of the winding state reference matrix corresponds to multiple load dynamic response reference vectors under normal operating conditions, and the column vector of the winding state reference matrix corresponds to each component of the load dynamic response vector. Calculate the orthogonal projection of the load dynamic response vector onto the eigenspace spanned by the winding state reference matrix to obtain the projection component vector, and calculate the residual vector between the projection component vector and the load dynamic response vector; Calculate the Euclidean norm of the residual vector and use it as the projection distance; The projection distance is compared with the state distance threshold pre-stored in the edge computing node. When the projection distance is greater than the state distance threshold, the distribution transformer is determined to be in an unstable operating area, and an unstable state flag is output. When the projection distance is less than or equal to the state distance threshold, the distribution transformer is determined to be in a stable operating area, and a stable state flag is output.

10. The outdoor distribution transformer condition monitoring system based on edge computing according to claim 9, characterized in that, The calculation of the orthogonal projection of the load dynamic response vector onto the eigenspace spanned by the winding state reference matrix to obtain the projection component vector, and the calculation of the residual vector between the projection component vector and the load dynamic response vector, specifically includes: The winding state reference matrix is ​​subjected to singular value decomposition to obtain a left singular matrix, a singular value diagonal matrix and a right singular matrix. The first k largest singular values ​​are selected from the singular value diagonal matrix according to a preset energy retention threshold. The first k columns of left singular vectors corresponding to the first k largest singular values ​​in the left singular matrix are extracted as orthogonal basis matrices of the feature subspace. Calculate the product of the transpose of the orthogonal basis matrix and the load dynamic response vector to obtain the projection coefficient vector of the load dynamic response vector on the feature subspace; Multiply the orthogonal basis matrix by the projection coefficient vector to obtain the projection component vector; The residual vector is obtained by subtracting the projection component vector from the load dynamic response vector, and the sum of squares of each element in the residual vector is calculated to obtain the residual energy value. The residual energy value is used as the basis for subsequent projection distance calculation.

Images