A method and system for diagnosing transformer winding damage under short circuit impulse

Abstract

The present application relates to the technical field of transformer diagnosis, and discloses a kind of short-circuit impulse under transformer winding damage diagnosis method and system, the method includes obtaining original impulse waveform data of previous short-circuit event, and impulse time domain feature extraction is carried out, and waveform feature set is obtained;According to the waveform feature set, winding stress mapping solution is carried out, and the residual stress of each spatial node of winding is obtained;Residual stress and pre-stored historical stress record are combined to calculate stress impact frequency, and the distribution characteristics of impact frequency are obtained, and stress time sequence accumulation calculation is carried out, and the stress superposition result of each spatial node is obtained;According to stress superposition result, deduce cumulative deformation value, and identify fatigue fracture risk accordingly, and obtain the winding fracture risk identification of each spatial node;Based on winding fracture risk identification and cumulative deformation value, fit deformation evolution track, determine rate abnormal characteristics, extract early damage characteristics, and package and generate early damage warning instruction.The present application can accurately identify the local cumulative fatigue risk of winding in time.

Description

Technical Field

[0001] This invention relates to the field of transformer diagnostic technology, and in particular to a method and system for diagnosing transformer winding damage under short-circuit impact. Background Technology

[0002] Currently, power transformer windings are subjected to short-circuit impacts for extended periods. Residual stress accumulates with varying waveform peak values, durations, and loading sequences, potentially leading to localized fatigue conditions that precede overall deformation. On-site operations require continuous acquisition of impact time-domain characteristics, tracking of stress history, and estimation of material deformation. Therefore, it is necessary to combine intelligent sensors, multi-source stress modeling, and risk classification and location technologies to achieve early identification and equipment locking.

[0003] In one existing technology, transformer short-circuit impact current data is collected, and the corresponding short-circuit electrodynamic force is calculated based on the current amplitude. This electrodynamic force is then converted into equivalent winding stress. Damage is calculated for each stress range based on a preset stress range and by referencing material fatigue curves. The damage values ​​for each range are then summed to obtain the cumulative fatigue damage value of the winding. This damage value is then compared with a set threshold to output the remaining lifespan or risk level. However, the influence of the decay time constant of the DC component of the short-circuit current and the difference in the impact waveform envelope on the transient electrodynamic impact work of the winding can lead to equivalent stress distortion. The correlation and coupling of historical impact sequences with the degradation of winding clamping force and the accumulation of local residual strain can also ignore the nonlinear hysteresis effect of the material's mechanical response. Therefore, when the local plastic accumulation in the stress concentration area at the winding end enters the critical state before the macroscopic overall deformation, the existing technology is prone to masking local danger points, resulting in a serious lag in warning time and misjudgment of risk level.

[0004] In summary, existing technologies have the problem of difficulty in timely and accurately identifying the risk of localized cumulative fatigue in windings. Summary of the Invention

[0005] This invention provides a method and system for diagnosing transformer winding damage under short-circuit impact, in order to promptly and accurately identify the risk of localized cumulative fatigue in the winding.

[0006] In a first aspect, to solve the above-mentioned technical problems, the present invention provides a method for diagnosing transformer winding damage under short-circuit impact, comprising:

[0007] Obtain the original impact waveform data of each short circuit event, and extract the impact time domain features from the original impact waveform data to obtain a waveform feature set;

[0008] The winding stress is solved by performing a winding stress mapping based on the waveform feature set to obtain the residual stress at each spatial node of the winding.

[0009] Based on the residual stress and the pre-stored historical stress records, the frequency of stress impacts is statistically analyzed to obtain the impact frequency distribution characteristics. Then, based on the impact frequency distribution characteristics, the stress time series accumulation calculation is performed to obtain the stress superposition results of each spatial node.

[0010] Based on the stress superposition results, deformation time series deduction calculations are performed to obtain the cumulative deformation values ​​of each spatial node;

[0011] Based on the cumulative deformation values, fatigue fracture risk is identified to obtain the winding fracture risk identifier for each spatial node.

[0012] Based on the winding fracture risk identifier, the cumulative deformation value, and the corresponding node coordinates, a deformation sequence curve is fitted to obtain the initial deformation evolution trajectory. Based on the initial deformation evolution trajectory, rate anomaly characteristics are determined to obtain early damage characteristics. Based on the early damage characteristics, an early warning instruction is encapsulated to obtain an early damage early warning instruction containing the risk node coordinates.

[0013] Secondly, the present invention provides a transformer winding damage diagnosis system under short-circuit impact, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method described above.

[0014] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method described above.

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

[0016] (1) This invention obtains a set of waveform features by performing wavelet threshold denoising, peak amplitude extraction and over-limit duration calculation on the original impact waveform data of each short circuit event. Furthermore, it performs electromagnetic force conversion, finite element displacement solution and elastoplastic constitutive mapping to obtain the residual stress of each spatial node of the winding. Since the waveform features can reflect the short circuit impact intensity and the continuous action process, it solves the problem that the existing technology only calculates the equivalent stress based on the current amplitude and is difficult to reflect the impact time domain difference, thus improving the accuracy of local stress identification of the winding.

[0017] (2) This invention combines residual stress with historical stress records, uses stress interval marking, rain flow counting and impact time interval statistics to obtain the impact frequency distribution characteristics, and combines material fatigue curves, dynamic recovery modulus, stiffness degradation matrix and stress relaxation weight to accumulate stress over time. Since this process can comprehensively reflect historical impact frequency, loading sequence and material recovery effect, it solves the problem that the linear accumulation damage in the prior art is difficult to characterize the historical impact accumulation effect, and improves the reliability of fatigue damage assessment.

[0018] (3) This invention extracts multi-dimensional time series features such as stress peak, stress change gradient and high stress duration through a sliding window, and uses long short-term memory network and multilayer perceptron regression model to deduce cumulative deformation values. Then, it combines real-time operating parameters and Weibull probability distribution to identify fatigue fracture risk and generate early damage warning instructions. Since it can simultaneously characterize stress evolution, deformation trend and operating condition influence, it solves the problem of lagging identification of local cumulative fatigue risk of winding in the existing technology, and realizes timely identification and warning of winding damage risk under short circuit impact. Attached Figure Description

[0019] Figure 1 This is a schematic flowchart of the method for diagnosing transformer winding damage under short-circuit impact provided in the first embodiment of the present invention. Detailed Implementation

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

[0021] Reference Figure 1 The first embodiment of the present invention provides a method for diagnosing transformer winding damage under short-circuit impact, comprising the following steps:

[0022] S11, Obtain the original impact waveform data of each short circuit event, and extract the impact time domain features from the original impact waveform data to obtain a waveform feature set;

[0023] S12, Solve the winding stress mapping based on the waveform feature set to obtain the residual stress of each spatial node of the winding;

[0024] S13, Perform stress impact frequency statistics based on the residual stress and pre-stored historical stress records to obtain impact frequency distribution characteristics, and perform stress time sequence accumulation calculation based on the impact frequency distribution characteristics to obtain the stress superposition result of each spatial node;

[0025] S14, Perform deformation time series deduction calculation based on the stress superposition result to obtain the cumulative deformation value of each spatial node;

[0026] S15, fatigue fracture risk identification is performed based on the cumulative deformation value to obtain the winding fracture risk identifier of each spatial node;

[0027] S16. Based on the winding fracture risk identifier, the cumulative deformation value and the corresponding node coordinates, the deformation sequence curve is fitted to obtain the initial deformation evolution trajectory. Based on the initial deformation evolution trajectory, the rate anomaly characteristic is determined to obtain the early damage characteristic. Based on the early damage characteristic, the early warning instruction is encapsulated to obtain the early damage early warning instruction containing the risk node coordinates.

[0028] In step S11, the original impact waveform data of each short-circuit event is acquired, and the impact time-domain features of the original impact waveform data are extracted to obtain a waveform feature set, including:

[0029] Wavelet threshold denoising is performed on the original impact waveform data to obtain a smooth waveform sequence;

[0030] The signal amplitude of each sampling point in the smooth waveform sequence is traversed, the maximum signal amplitude is determined as the peak amplitude, and the peak time corresponding to the peak amplitude is recorded.

[0031] Search forward and backward from the peak time for the time points when the signal amplitude crosses the preset threshold for the over-limit signal, respectively, to obtain the over-limit start time and over-limit end time;

[0032] Subtracting the over-limit start time from the over-limit termination time yields the duration, and concatenating the peak amplitude and the duration in dimensional order yields the waveform feature set.

[0033] Specifically, the original impact waveform data comes from the waveform recordings saved after each short-circuit event is triggered. These recordings are written into the historical waveform storage area by the transformer short-circuit protection device or the winding end voltage acquisition unit when a short-circuit event occurs. The original impact waveform data includes the short-circuit event number, sampling point number, sampling time, and the signal amplitude corresponding to each sampling point. The signal amplitude is the voltage signal amplitude converted from the short-circuit current by the sampling resistor. During processing, the original impact waveform data is read one by one according to the short-circuit event number and arranged from earliest to latest according to the sampling time, forming a set of waveform sequences to be processed corresponding to each short-circuit event. Adjacent sampling points are arranged according to a fixed sampling interval, so that the peak amplitude, peak time, over-limit start time, and over-limit end time can all be determined by the sampling point position and sampling time.

[0034] When performing wavelet threshold denoising on the original impact waveform data, the signal amplitudes of two adjacent sampling points in the waveform sequence to be processed are first added together and divided by 2 to obtain the low-frequency approximation coefficients. Then, the signal amplitude of the previous sampling point is subtracted from the signal amplitude of the next sampling point and divided by 2 to obtain the high-frequency detail coefficients. The low-frequency approximation coefficients are further decomposed according to the above addition, subtraction, and division process. The number of decomposition layers is the base-2 logarithm of the number of sampling points, rounded down. If the value is less than 3, it is taken as 3; if it is greater than 8, it is taken as 8; otherwise, the value itself is taken. The absolute values ​​of the high-frequency detail coefficients of the highest decomposition layer are taken and sorted from smallest to largest. The median is read, and the median is divided by 0.6745 to obtain the noise standard deviation estimate. The noise standard deviation estimate is then multiplied by the square root of twice the natural logarithm of the number of sampling points to obtain the wavelet denoising threshold. For each high-frequency detail coefficient, a thresholding process is performed. High-frequency detail coefficients with absolute values ​​less than or equal to the wavelet denoising threshold are set to 0, while those with absolute values ​​greater than the wavelet denoising threshold retain their original signs and their absolute values ​​are subtracted from the wavelet denoising threshold. After thresholding, the low-frequency approximation coefficients and the processed high-frequency detail coefficients are synthesized layer by layer in reverse to obtain the smooth waveform sequence.

[0035] After obtaining the smoothed waveform sequence, the signal amplitude of each sampling point in the sequence is traversed from earliest to latest sampling time. During the traversal, a current maximum value storage bit and a current maximum value time storage bit are set. First, the signal amplitude of the first sampling point is written to the current maximum value storage bit, and the sampling time of the first sampling point is written to the current maximum value time storage bit. For each new sampling point, the signal amplitude of that sampling point is subtracted from the value in the current maximum value storage bit. If the difference is greater than 0, the signal amplitude of that sampling point replaces the value in the current maximum value storage bit, and the sampling time of that sampling point replaces the time in the current maximum value time storage bit. After the traversal is completed, the value in the current maximum value storage bit is determined as the peak amplitude, and the time in the current maximum value time storage bit is determined as the peak time.

[0036] The preset over-limit signal threshold is determined through historical data statistics. Historical short-circuit event records saved after the same transformer is put into operation are read. The signal amplitudes that did not trigger protection actions are extracted from the five power frequency cycles preceding each short-circuit event. The sampling rate of the waveform recording device is not less than 1kHz, and the five power frequency cycles correspond to no less than 100 sampling points. The absolute values ​​of the extracted signal amplitudes are taken and sorted from smallest to largest. The signal amplitude at the 95th percentile is read as the baseline fluctuation upper limit. The baseline fluctuation upper limit is multiplied by 1.2 to obtain the preset over-limit signal threshold, where 1.2 is the distinguishing margin coefficient between inrush current and short-circuit current for the same type of equipment under the system's maximum operating mode. The preset over-limit signal threshold is saved according to the transformer equipment number and remains consistent throughout the handling of short-circuit events for the same equipment.

[0037] After determining the peak time, a forward search is performed with the sampling point corresponding to the peak time as the center, comparing the signal amplitude of adjacent sampling points with the preset over-limit signal threshold point by point. When the signal amplitude of a sampling point closer to the peak time is greater than or equal to the preset over-limit signal threshold, and the signal amplitude of an adjacent sampling point farther from the peak time is less than the preset over-limit signal threshold, an over-limit start point is determined to exist in the adjacent sampling interval. The preset over-limit signal threshold is subtracted from the signal amplitude of the sampling point farther from the peak time to obtain the starting amplitude difference. The signal amplitude of the sampling point closer to the peak time is subtracted from the signal amplitude of the sampling point farther from the peak time to obtain the starting interval amplitude difference. The starting amplitude difference is divided by the starting interval amplitude difference to obtain the starting interval ratio. The starting interval ratio is multiplied by the time interval between the two sampling points, and the sampling time of the sampling point farther from the peak time is added to obtain the over-limit start time.

[0038] The system searches backward from the sampling point corresponding to the peak time, comparing the signal amplitude of adjacent sampling points with the preset over-limit signal threshold. When the signal amplitude of a sampling point closer to the peak time is greater than or equal to the preset over-limit signal threshold, and the signal amplitude of an adjacent sampling point farther from the peak time is less than the preset over-limit signal threshold, it is determined that there is an over-limit endpoint in the adjacent sampling interval. The preset over-limit signal threshold is subtracted from the signal amplitude of the sampling point closer to the peak time to obtain the endpoint amplitude difference. The signal amplitude of the sampling point farther from the peak time is subtracted from the signal amplitude of the sampling point closer to the peak time to obtain the endpoint interval amplitude difference. The endpoint amplitude difference is divided by the endpoint interval amplitude difference to obtain the endpoint interval ratio. The endpoint interval ratio is multiplied by the time interval between the two sampling points, and then the sampling time of the sampling point closer to the peak time is added to obtain the over-limit termination time.

[0039] After obtaining the over-limit start time and the over-limit end time, the over-limit end time is subtracted from the over-limit start time to obtain the duration. For each short-circuit event, a set of two-dimensional features consisting of the peak amplitude and the duration is generated, where the first dimension is the peak amplitude and the second dimension is the duration. Following the order of short-circuit event occurrence from earliest to latest, the two-dimensional features corresponding to each short-circuit event are sequentially written into the same set to obtain the waveform feature set. Each row in the waveform feature set corresponds to one short-circuit event, and the two dimensions in each row correspond to the peak amplitude and the duration, respectively.

[0040] In step S12, the winding stress mapping is solved based on the waveform feature set to obtain the residual stress at each spatial node of the winding, including:

[0041] The waveform feature set is input into a pre-constructed electromagnetic force conversion model to solve for the electromagnetic force distribution of the winding space leakage magnetic field, thereby obtaining the transient electromagnetic force distribution matrix.

[0042] Based on the transient electromagnetic force distribution matrix, the finite element displacement of each spatial node of the winding is solved to obtain the deformation displacement field.

[0043] The strain distribution tensor is obtained by dividing the displacement difference between adjacent spatial nodes in the deformation displacement field by the corresponding node spacing.

[0044] The strain distribution tensor is input into a pre-constructed elastoplastic bilinear kinematic reinforcement model, and constitutive relation mapping is performed on each spatial node of the winding to obtain the residual stress.

[0045] Specifically, the electromagnetic force conversion model is constructed before diagnostics are performed. Its data sources include transformer nameplate parameters, winding design dimension records, short-circuit impedance test records, and winding material parameter records. The input to the electromagnetic force conversion model is the waveform feature set, with the input dimensions set according to the order of peak amplitude and duration. The model is pre-programmed with the number of winding turns, average winding radius, winding axial height, winding short-circuit impedance, conductor cross-sectional area, effective conductor length, and leakage flux correction coefficient. The leakage flux correction coefficient is determined from the factory short-circuit impedance test records of transformers of the same model. The measured short-circuit impedance of each device is divided by the design short-circuit impedance to obtain the correction value for a single device. The correction values ​​for devices of the same model are then summed and divided by the number of devices to obtain the leakage flux correction coefficient.

[0046] When using the electromagnetic force conversion model, the peak amplitude of the waveform feature set is first read, and the peak amplitude is divided by the winding short-circuit impedance to obtain the peak value of the impulse current. Then, an impulse loading sequence is generated based on the duration. The time step of the impulse loading sequence is taken from the sampling interval of the original impulse waveform data. The loading sequence increases point by point from the start of the over-limit to the peak time, reaching the peak value of the impulse current at the peak time, and then decreases point by point to the end of the over-limit. The impulse current at each sampling point is multiplied by the number of winding turns to obtain the ampere-turn value. The ampere-turn value is multiplied by the leakage magnetic field correction coefficient to obtain the winding spatial leakage magnetic field strength. The winding spatial leakage magnetic field strength is multiplied by the impulse current at the corresponding sampling point, and then multiplied by the effective length of the conductor to obtain the electromagnetic force value at that sampling point. The electromagnetic force values ​​of each node are written in the order of the radial and axial nodes of the winding to obtain the transient electromagnetic force distribution matrix.

[0047] The transient electromagnetic force distribution matrix is ​​used to solve the finite element displacement of each spatial node of the winding. The finite element nodes are determined by the number of radial layers and axial segments of the winding. The radial direction is divided according to the insulation boundary between winding layers, and the axial direction is divided into 20 segments based on the average axial height of the winding. Each node is recorded with a node number, radial coordinates, axial coordinates, and the connection relationship with adjacent nodes. Before solving the finite element displacement, node stiffness values ​​are generated based on the conductor cross-sectional area, elastic modulus, and node spacing. The node stiffness value is obtained by multiplying the elastic modulus by the conductor cross-sectional area and then dividing by the node spacing. The displacement at the end of the winding, where the clamping position is pressed, is recorded as 0 as a boundary condition. The node electromagnetic force in the transient electromagnetic force distribution matrix is ​​divided by the corresponding node stiffness value to obtain the initial node displacement. Then, displacement balance iteration is performed according to the connection relationship between adjacent nodes. Iteration stops when the absolute value of the difference between two consecutive node displacements is less than one-thousandth of the absolute value of the previous displacement, thus obtaining the deformation displacement field.

[0048] After obtaining the deformation displacement field, the displacement values ​​of adjacent spatial nodes are read according to the node connection relationship. The displacement of the preceding spatial node is subtracted from the displacement of the subsequent spatial node to obtain the displacement difference. This displacement difference is then divided by the corresponding node spacing to obtain the strain value in the node connection direction. Radial adjacent nodes yield radial strain, axial adjacent nodes yield axial strain, and circumferential adjacent nodes yield circumferential strain. The radial strain, axial strain, and circumferential strain are written into the tensor storage area in node number order to obtain the strain distribution tensor. Each row in the strain distribution tensor corresponds to a spatial node, and the dimensions in each row correspond to radial strain, axial strain, and circumferential strain, respectively.

[0049] The elastoplastic bilinear kinematic hardening model is constructed before diagnostic testing. Its parameters include the elastic modulus, yield stress, tangent modulus, kinematic hardening modulus, and initial back stress of the winding conductor material. The elastic modulus and yield stress are read from the winding material parameter record; if the record is missing, it is read from the tensile test record of the same batch of winding conductors. The tangent modulus is obtained by dividing the elastic modulus by 100, a ratio determined by the statistical average of the slopes of the elastoplastic segment in the tensile tests of the same batch. The kinematic hardening modulus is the tangent modulus, and the initial back stress is set to 0. When using the elastoplastic bilinear kinematic hardening model, the strain at each node in the strain distribution tensor is first multiplied by the elastic modulus to obtain the trial stress. Then, the back stress is subtracted from the trial stress to obtain the yield discrimination stress. If the absolute value of the yield discrimination stress is not higher than the yield stress, the trial stress is written into the stress result for that node. When the absolute value of the yield criterion stress is higher than the yield stress, the strain portion exceeding the yield stress is multiplied by the tangent modulus and then added to the yield stress in the original strain direction to obtain the plastic correction stress of that node. Simultaneously, the plastic strain increment is multiplied by the kinematic hardening modulus and added to the back stress to obtain the updated back stress. After completing the mapping of all nodes, the node stress values ​​remaining after unloading are read to obtain the residual stress.

[0050] In step S13, stress impact frequency statistics are performed based on the residual stress and pre-stored historical stress records to obtain impact frequency distribution characteristics. Then, stress time-series accumulation calculations are performed based on these impact frequency distribution characteristics to obtain the stress superposition results for each spatial node, including:

[0051] The stress peak value and the time of occurrence of each short-circuit impact are read from the historical stress record and arranged in chronological order to obtain the historical stress event sequence.

[0052] The residual stress is added to the historical stress event sequence according to the current short-circuit impact time to obtain the updated stress event sequence;

[0053] According to the preset stress amplitude range, each stress peak in the updated stress event sequence is assigned a range to be marked, thus obtaining a stress range marking sequence;

[0054] According to the rainflow counting rule, the stress peak and valley changes in the stress interval marking sequence are identified by closed loop, and the number of stress cycles corresponding to each preset stress amplitude interval is obtained.

[0055] The frequency vector is arranged according to the stress cycle number and the time interval between adjacent impacts to obtain the impact frequency distribution characteristics.

[0056] Based on the impact frequency distribution characteristics, the fatigue damage increment is calculated by combining the pre-established fatigue curve of the winding material, and the fatigue damage evolution rate is obtained.

[0057] Based on the fatigue damage evolution rate, the elastic recovery stage data in the historical stress record are subjected to elastic recovery reduction correction to obtain the dynamic recovery modulus;

[0058] The stiffness degradation matrix is ​​obtained by superimposing the damage ratio based on the dynamic recovery modulus and the residual stress.

[0059] The cumulative amount of plastic deformation is obtained by implicit time integration based on the stiffness degradation matrix.

[0060] The cumulative amount of plastic deformation is substituted into a pre-constructed exponential decay function to obtain the stress relaxation weight. The residual stress and the historical stress record are then weighted and accumulated according to the stress relaxation weight to obtain the stress superposition result.

[0061] Specifically, the historical stress records are pre-stored in the transformer diagnostic database, and their data source is the records written after each short-circuit event has completed the winding stress mapping solution. The historical stress records include the short-circuit event number, impact time, peak stress, residual stress, elastic recovery stage start time, elastic recovery stage end time, and elastic recovery stage data; wherein the elastic recovery stage start time is the time when the current short-circuit impact exceeds its limit, and the elastic recovery stage end time is the time when the residual stress drops to 5% of the initial residual stress of the impact. After reading the peak stress and impact time corresponding to each short-circuit impact from the historical stress records, they are arranged in ascending order of impact time to obtain the historical stress event sequence. The residual stress calculated in this instance is added to the historical stress event sequence according to the current short-circuit impact time. If the current short-circuit impact time is later than the end of the sequence, it is written to the end of the sequence; if the current short-circuit impact time is between two historical times, it is written to the corresponding time position to obtain the updated stress event sequence.

[0062] The preset stress amplitude range is determined by the historical stress records and the current residual stress. First, all stress peak values ​​in the historical stress records are read, and the stress peak value corresponding to the current residual stress is added to obtain a set of stress peak values. The minimum value in the stress peak value set is rounded down to a multiple of 10 MPa as the lower limit of the range, and the maximum value in the stress peak value set is rounded up to a multiple of 10 MPa as the upper limit of the range. Starting from the lower limit, a new range is generated every 10 MPa increase until the upper limit is reached. When assigning ranges to each stress peak value in the updated stress event sequence, each stress peak value is compared sequentially with the lower and upper boundaries of each range. If the value is greater than or equal to the lower boundary and less than the upper boundary, the range number is written to the corresponding position of the stress peak value, resulting in the stress range marking sequence.

[0063] When identifying closed loops according to the rainflow counting rules, the stress peaks in the updated stress event sequence are first read in chronological order, and the corresponding interval numbers in the stress interval marking sequence are retained. Midpoints with continuous increases or decreases are deleted, while peaks and valleys with changed directions are retained. For every three consecutive peaks and valleys read, the stress amplitude difference between the last two peaks and valleys is calculated, and then the stress amplitude difference between the first two peaks and valleys is calculated. If the stress amplitude difference between the last two peaks and valleys is greater than or equal to the stress amplitude difference between the first two peaks and valleys, the first two peaks and valleys are identified as a closed loop, and the number of stress amplitude intervals corresponding to this closed loop is incremented by 1. If the above conditions are not met, the next peak and valley is read. After processing all peaks and valleys, unclosed peak and valley segments are written into the corresponding stress amplitude intervals in half-cycles, obtaining the stress cycle count corresponding to each preset stress amplitude interval.

[0064] After obtaining the stress cycle count, the times of adjacent impacts in the updated stress event sequence are read. The time interval between adjacent impact times is obtained by subtracting the time of the previous impact from the time of the subsequent impact. Each preset stress amplitude interval corresponds to a frequency vector. The first dimension of the frequency vector is the stress cycle count for that interval, the second dimension is the average time interval between adjacent impact times within that interval, and the third dimension is the shortest time interval within that interval. The frequency vectors of all intervals are arranged from low to high according to the interval amplitude to obtain the impact frequency distribution characteristics. When there are no impact records in a certain stress amplitude interval, the stress cycle count, average time interval, and shortest time interval for that interval are all written as 0.

[0065] The fatigue curve of the winding material is established before diagnosis, and its data comes from the fatigue test records of the winding conductor material and the factory material test records. During establishment, the stress amplitude and fracture cycle count in each set of test records are read. The stress amplitudes are arranged in ascending order, and the fracture cycle counts between adjacent stress amplitudes are filled using linear interpolation. When calculating the fatigue damage increment, the stress amplitude at the center of each interval in the impact frequency distribution characteristics is input into the winding material fatigue curve, and the corresponding fracture cycle count is read. The stress cycle count of the interval is then divided by the corresponding fracture cycle count to obtain the fatigue damage increment for that interval. The fatigue damage increments of each interval are accumulated sequentially and divided by the total time covered by the updated stress event sequence to obtain the fatigue damage evolution rate.

[0066] When performing elastic recovery reduction correction on the elastic recovery stage data in the historical stress record, first read the stress values ​​at the start and end of the elastic recovery stage, and subtract the end stress value from the start stress value to obtain the recovery stress difference. Then read the start and end deformation values ​​of the same stage, and subtract the end deformation value from the start deformation value to obtain the recovery deformation difference. Divide the recovery stress difference by the recovery deformation difference to obtain the basic recovery modulus. Subtract the fatigue damage evolution rate from 1 to obtain the reduction factor, and then multiply the basic recovery modulus by the reduction factor to obtain the dynamic recovery modulus. When the reduction factor is less than 0.05, write 0.05 to ensure that the dynamic recovery modulus does not have a negative value or an excessively small value.

[0067] In solving the stiffness degradation matrix, this embodiment uses a simplified discretization form with superimposed damage proportions. First, the initial stiffness matrix corresponding to the finite element nodes is read. The residual stress is divided by the yield stress of the winding material to obtain the node stress-damage ratio. The fatigue damage evolution rate is multiplied by the time interval between adjacent impact occurrences to obtain the node fatigue damage ratio. The node stress-damage ratio and the node fatigue damage ratio are added together to obtain the node damage ratio. A node damage ratio greater than 1 is written as 1. The node residual stiffness ratio is obtained by subtracting the node damage ratio from 1. Then, the stiffness value of the corresponding node in the initial stiffness matrix is ​​multiplied by the node residual stiffness ratio to obtain the degraded node stiffness value. All degraded node stiffness values ​​are written into the matrix according to the finite element node order to obtain the stiffness degradation matrix.

[0068] When performing implicit time integration based on the stiffness degradation matrix, the sampling interval of the original impact waveform data is used as the time step. Within each time step, the residual stress external force term corresponding to the current time step, the cumulative plastic deformation of the previous time step, and the stiffness degradation matrix are read first. The scalar norm of the external force term in the current time step is divided by the degraded stiffness value of the corresponding node in the stiffness degradation matrix to obtain the deformation increment of the current time step; where the scalar norm is the square root of the sum of squares of the force components in each direction. The deformation increment of the current time step is added to the cumulative plastic deformation of the previous time step to obtain the cumulative plastic deformation of the current time step. The calculation is completed for all time steps from earliest to latest sampling time to obtain the cumulative plastic deformation.

[0069] The exponential attenuation function is constructed before diagnostics are performed. Its input is the cumulative plastic deformation, and its output is the stress relaxation weight. The attenuation coefficient is determined by the stress relaxation records in the historical stress log. First, the data on the decrease in residual stress over time after each short-circuit impact is read. The cumulative plastic deformation corresponding to the point where the residual stress decreases to 36.8% of the initial residual stress after the impact is found. This cumulative plastic deformation is used as the attenuation benchmark value, and then 1 is divided by the attenuation benchmark value to obtain the attenuation coefficient. When the number of historical stress relaxation records is 0, the attenuation benchmark value is written into the median of the cumulative plastic deformation in the factory short-circuit withstand test records of the same model transformer. If the factory test records are missing, the attenuation coefficient uses the design value of the cumulative plastic deformation under rated short-circuit impact calculated during the transformer's design phase. In use, the cumulative plastic deformation is multiplied by the attenuation coefficient and the negative number is taken, then exponentially calculated to obtain the stress relaxation weight. When the stress relaxation weight is less than 0.05, it is written as 0.05. The residual stress and each historical residual stress in the historical stress record are multiplied by the corresponding stress relaxation weight, and then accumulated in chronological order to obtain the stress superposition result.

[0070] In step S14, deformation time series extrapolation calculations are performed based on the stress superposition results to obtain the cumulative deformation values ​​of each spatial node, including:

[0071] Based on the stress superposition result, the sequence is divided by a sliding window to obtain the stress window sequence;

[0072] Time-domain statistics are extracted from the stress window sequence to obtain multi-dimensional time series features including stress peak, stress change gradient, and high stress duration;

[0073] The multidimensional time series features are arranged in order of sampling time to obtain the input sequence;

[0074] The input sequence is fed into a pre-trained long short-term memory network for state updates to obtain the hidden layer state;

[0075] The hidden layer state is input into a pre-trained multilayer perceptron regression model and nonlinearly mapped to obtain the cumulative deformation value.

[0076] The duration of high stress is calculated by adding the durations of consecutive sampling intervals that are higher than the preset high stress discrimination threshold.

[0077] Specifically, when performing deformation time-series extrapolation calculations based on the stress superposition results, each spatial node is operated independently. The stress superposition results of that node are arranged according to the sampling time to form the stress time series of that node. All spatial nodes share the same LSTM and MLP weights, and the training labels for non-measurement point nodes are supplemented by finite element simulation. The length of the sliding window is determined by the duration of historical high stress. First, the duration of high stress for each short-circuit impact in the historical stress record is read, and the high stress durations are sorted from smallest to largest. The duration value corresponding to the 75th percentile is read, and the duration value is divided by the waveform recording sampling interval and rounded to the nearest integer to obtain the number of sampling points as the window length. If the window length is less than 10 sampling points, 10 sampling points are written. The sliding step size is half of the window length and rounded down. The first segment of the window is extracted from the starting point of the stress time series according to the window length, and then the window is moved backward according to the sliding step size until the end of the stress time series to obtain the stress window sequence.

[0078] When extracting time-domain statistics for the stress window sequence, the stress peak value, stress change gradient, and high stress duration are extracted for each stress window. The stress peak value is obtained by traversing the stress values ​​of each sampling point within the window to obtain the maximum stress value. The stress change gradient is obtained by subtracting the stress value of the previous sampling point from the stress value of the subsequent sampling point, then dividing by the time interval between the two sampling points, and summing all gradients within the window and dividing by the number of gradients to obtain the stress change gradient of that window. The preset high stress discrimination threshold is determined by the historical stress records. The stress peak values ​​of the historical stable operation phase are read and sorted from smallest to largest, and the 90th percentile is taken as the preset high stress discrimination threshold. When the historical stable operation phase data is insufficient, 1.2 times the stress peak value under rated operating conditions is taken as a substitute threshold. Sampling points within the window that are higher than the preset high stress discrimination threshold are marked as valid states. The duration of consecutive valid state sampling intervals is calculated and accumulated to obtain the high stress duration. The stress peak value, stress change gradient, and high stress duration are written in dimensional order to obtain the multidimensional time series features.

[0079] The input sequence is obtained by arranging the multidimensional time series features in the order of sampling time. Before inputting into the Long Short-Term Memory (LSTM) network, the input sequence is normalized. The normalization parameter is determined by the training samples. The minimum and maximum values ​​of each feature dimension in the training samples are read, the current feature value is subtracted from the minimum value, and then divided by the difference between the maximum and minimum values ​​to obtain the normalized feature value. When the normalized feature value is less than 0, it is taken as 0, and when it is greater than 1, it is taken as 1. The LSM network is trained before performing diagnosis. The training samples are obtained from the stress superposition results of historical short-circuit events and the winding deformation records obtained from maintenance measurements. The maintenance measurements are performed using a winding deformation frequency response tester, collecting deformation displacement data from 5 measuring points at the end and middle of the winding. The total number of training samples is not less than 200. If the actual measurements are insufficient, finite element simulation data is used to supplement the data to 200. The input of each training sample is the stress peak value, stress change gradient, and high stress duration arranged in time. The label is the cumulative deformation value obtained from maintenance measurements after the same short-circuit event.

[0080] The Long Short-Term Memory (LSTM) network employs a single-layer structure with an input dimension of 3 and 16 hidden layer nodes. Internally, the network includes an input gate, a forget gate, an output gate, and candidate states. At each time step, the input features are multiplied by the input weights, and the hidden layer state from the previous time step is multiplied by the cyclic weights. These two parts are summed, and a bias is added to obtain the gating value. This gating value is then input into a sigmoid compression function to obtain a gating ratio between 0 and 1. The sigmoid compression function is obtained by dividing 1 by the sum of 1 and the exponent value, where the exponent value is the negative of the gating value. Candidate states are compressed to between -1 and 1 using a hyperbolic tangent function. The forget gate controls the retention ratio of the previous memory state, the input gate controls the writing ratio of candidate states, and the output gate controls the output ratio of the current memory state. This process is repeated time-by-time to update the hidden layer state. During training, the average of squared differences is used as the error function, the weights are updated according to the backpropagation results, the learning rate is set to 0.001, the number of training epochs is set to 200, and the input weights, loop weights and biases after training are written to the model parameter file.

[0081] The multilayer perceptron regression model is jointly trained with the long short-term memory network before diagnostics are performed. The multilayer perceptron regression model has an input dimension of 16, 8 hidden layer nodes, and an output dimension of 1. The 16 values ​​in the hidden layer state are multiplied by the weights of the first layer, summed, and then the first layer bias is added to obtain the hidden layer input value. If the hidden layer input value is less than 0, it is written as 0; if the hidden layer input value is greater than or equal to 0, the original value is retained to obtain the hidden layer output value. The hidden layer output value is multiplied by the weights of the second layer, summed, and then the second layer bias is added to obtain the cumulative deformation value. After training, the input sequence is input into the pre-trained long short-term memory network to obtain the hidden layer state, and then the hidden layer state is input into the pre-trained multilayer perceptron regression model to obtain the cumulative deformation value.

[0082] In step S15, fatigue fracture risk is identified based on the accumulated deformation value to obtain the winding fracture risk identifier for each spatial node, including:

[0083] The accumulated deformation value and the pre-acquired initial physical property data of the winding are concatenated in the order of feature dimensions to obtain the deformation feature vector;

[0084] The deformation feature vector is sorted by Euclidean distance with each sample vector in the pre-established material fatigue attenuation mapping table to obtain the nearest neighbor sample sequence.

[0085] The sample attenuation coefficient is read from the nearest sample sequence according to the preset number of neighborhoods, the average value of the sample attenuation coefficient is calculated to obtain the fatigue attenuation coefficient, and the fatigue limit difference is calculated based on the fatigue attenuation coefficient to obtain the fatigue limit margin.

[0086] The top oil temperature, real-time load rate, and winding vibration fundamental frequency are read from the smart sensor. The operating condition threshold is mapped to the top oil temperature, the real-time load rate, and the winding vibration fundamental frequency to obtain the real-time margin threshold.

[0087] When the fatigue limit margin is less than the real-time margin threshold, the difference between the real-time margin threshold and the fatigue limit margin is mapped by a Weibull probability distribution to obtain a fracture probability value. Based on the fracture probability value, a preset probability interval is matched to obtain the winding fracture risk identifier.

[0088] Specifically, the initial physical property data of the windings is written into the equipment basic database before the transformer is put into operation. The data sources include transformer design records, factory inspection records, and winding material parameter records. The initial physical property data of the windings includes conductor cross-sectional area, winding axial height, average winding radius, number of winding turns, insulation layer thickness, initial clamping force, elastic modulus, yield stress, and initial fatigue limit. The cumulative deformation value is written into the first dimension, and then the initial physical property data of the windings is written in the above order to obtain the deformation feature vector. Before writing, each dimension is normalized, and the normalization parameters are determined by the minimum and maximum values ​​of the same dimension in the material fatigue attenuation mapping table.

[0089] The material fatigue attenuation mapping table is established before diagnostics are performed. Its data sources include fatigue test records of winding conductor materials, short-circuit impact simulation records, and maintenance-measured deformation records. The maintenance-measured deformation records are obtained by collecting deformation displacement data at the ends and middle of the winding using a winding deformation frequency response tester during periodic maintenance. Each record corresponds to a short-circuit event number and node number. Each sample record includes a sample vector and a sample attenuation coefficient. The dimensional order of the sample vector is consistent with the deformation feature vector. The sample attenuation coefficient is obtained by dividing the remaining fatigue limit by the initial fatigue limit, where the remaining fatigue limit is read from the withstandable stress amplitude under the current deformation state in the fatigue test, and the initial fatigue limit is read from the fatigue test of the same batch of winding conductor materials that have not experienced short-circuit impact. After construction, all sample records are written into the material fatigue attenuation mapping table according to their sample numbers.

[0090] When sorting the deformation feature vector with each sample vector in the material fatigue attenuation mapping table using Euclidean distance, the difference is calculated dimension by dimension for each sample vector. The difference in that dimension is obtained by subtracting the corresponding dimension value from the value in one dimension of the deformation feature vector. Each dimension difference is multiplied by itself to obtain the squared difference in that dimension. The squared differences of all dimensions are summed and the square root is taken to obtain the Euclidean distance between the deformation feature vector and the sample vector. The sample records are arranged in ascending order of Euclidean distance to obtain the nearest neighbor sample sequence. The preset number of neighbors is determined by the total number of samples in the material fatigue attenuation mapping table. The square root of the total number of samples is rounded down; a value less than 3 is written as 3, and a value greater than 9 is written as 9.

[0091] After reading the sample attenuation coefficients from the nearest neighbor sample sequence according to the preset number of neighborhoods, the read sample attenuation coefficients are added one by one and then divided by the number of reads to obtain the fatigue attenuation coefficient. When calculating the fatigue limit difference, the initial fatigue limit in the initial physical property data of the winding is first read, and the initial fatigue limit is multiplied by the fatigue attenuation coefficient to obtain the current fatigue limit. Then, the current maximum stress peak value in the stress superposition result is read, and the current fatigue limit is subtracted from the current maximum stress peak value to obtain the fatigue limit margin. When the fatigue limit margin is less than 0, a negative value is retained for subsequent fracture probability calculations.

[0092] The top oil temperature, real-time load rate, and winding vibration fundamental frequency are read by intelligent sensors at the same sampling time. The thresholds mapped by the operating condition thresholds are determined before diagnostics are performed. The top oil temperature threshold is read from the 95th percentile of the top oil temperature during historical stable operation. The real-time load rate threshold is read from the rated load rate on the transformer nameplate. The winding vibration fundamental frequency threshold is read from the 95th percentile of the winding vibration fundamental frequency offset during historical stable operation. The base value of the real-time margin threshold is 10% of the initial fatigue limit. If the top oil temperature is higher than the top oil temperature threshold, the excess is divided by the top oil temperature threshold and multiplied by the base value of the real-time margin threshold, then added to the real-time margin threshold. If the real-time load rate is higher than the real-time load rate threshold, the excess is divided by the real-time load rate threshold and multiplied by the base value of the real-time margin threshold, then added to the real-time margin threshold. If the winding vibration fundamental frequency offset is higher than the winding vibration fundamental frequency threshold, the excess is divided by the winding vibration fundamental frequency threshold and multiplied by the base value of the real-time margin threshold, then added to the real-time margin threshold, finally obtaining the real-time margin threshold.

[0093] When the fatigue limit margin is less than 0, the fracture probability value is written as 1, and the winding fracture risk indicator is written as a high-risk indicator; when the fatigue limit margin is greater than or equal to the real-time margin threshold, the fracture probability value is written as 0, and the winding fracture risk indicator is written as a low-risk indicator; when the fatigue limit margin is less than the real-time margin threshold, the fatigue limit margin is subtracted from the real-time margin threshold to obtain the margin difference. The shape parameter and scale parameter of the Weibull probability distribution are determined by the material fatigue attenuation mapping table before diagnosis. The shape parameter is written as 2, and the scale parameter reads the median of the margin difference corresponding to the fracture sample. The fracture sample is the sample with a sample attenuation coefficient less than 0.2 in the material fatigue attenuation mapping table. When the number of fracture samples is 0, the scale parameter reads the 75th percentile of the margin difference of all samples. The margin difference is divided by the scale parameter, squared, then the opposite number is taken and exponentially calculated. The result of the exponential calculation is subtracted from 1 to obtain the fracture probability value. When the fracture probability value is less than 0.3, a low-risk indicator is written; when it is greater than or equal to 0.3 but less than 0.7, a medium-risk indicator is written; and when it is greater than or equal to 0.7, a high-risk indicator is written, thus obtaining the winding fracture risk indicator.

[0094] In step S16, deformation sequence curve fitting is performed based on the winding fracture risk identifier, the cumulative deformation value, and the corresponding node coordinates to obtain an initial deformation evolution trajectory. Rate anomaly characteristics are then determined based on the initial deformation evolution trajectory to obtain early damage characteristics. Finally, early warning instructions are encapsulated based on these early damage characteristics to obtain an early damage early warning instruction containing the risk node coordinates, including:

[0095] Based on the early damage characteristics, industrial communication protocol fields are mapped to obtain early warning command messages;

[0096] According to a preset time extraction interval, the time coordinates and deformation coordinates are read from the initial deformation evolution trajectory to obtain the trajectory coordinate point set;

[0097] A synchronization word and a message length field are written into the frame header of the warning instruction message to obtain the message header field;

[0098] The early damage features and the trajectory coordinate point set are written into the data payload area of ​​the early warning instruction message to obtain the payload field;

[0099] Cyclic redundancy check (CRC) calculations are performed on the message header field and the payload field to obtain a check field. The check field is then encapsulated in the order of the message header field, the payload field, and the check field to obtain the early damage warning instruction.

[0100] Specifically, when fitting the deformation sequence curve based on the winding fracture risk identifier and the cumulative deformation value, the cumulative deformation value stored by time in the diagnostic database is first read, and then added to the currently calculated cumulative deformation value to obtain the deformation time series data. The fitting window length is determined by the winding fracture risk identifier: low risk identifier corresponds to the most recent 20 deformation time series data, medium risk identifier corresponds to the most recent 10 deformation time series data, and high risk identifier corresponds to the most recent 5 deformation time series data. If the number of deformation time series data is less than the corresponding window length, all existing deformation time series data are read. The relative time is obtained by subtracting the window start time from each sampling time within the window, and then the relative time, the squared value of the relative time, and the cumulative deformation value are written into the fitting data table.

[0101] When performing curve fitting, a quadratic curve is used to represent the deformation over time. Three coefficients are first established, corresponding to a constant term, a first-order time term, and a second-order time term, respectively. The relative time of each sampling point is substituted into the quadratic curve to obtain the calculated deformation value for that sampling point. The cumulative deformation value of that sampling point is subtracted from the calculated deformation value to obtain the fitting error. The fitting errors of all sampling points are multiplied by themselves and then summed to obtain the sum of squared errors. The three coefficients are obtained by minimizing the sum of squared errors, and each sampling time is substituted into the quadratic curve after obtaining the coefficients to obtain the fitted deformation value corresponding to each sampling time. When there is only one deformation time series data point, the cumulative deformation value of that point is used as the full-window fitted deformation value. When there are only two deformation time series data points, a first-order time term is used to complete the linear fitting. The initial deformation evolution trajectory is obtained by writing the correspondence between sampling time and fitted deformation value into the trajectory storage area.

[0102] When determining rate anomaly characteristics based on the initial deformation evolution trajectory, two adjacent trajectory points are first read. The deformation coordinates of the later trajectory point are subtracted from the deformation coordinates of the earlier trajectory point to obtain the deformation increment. The deformation increment is then divided by the time coordinate difference between the two trajectory points to obtain the deformation rate. The deformation rate is then subtracted from the earlier deformation rate and divided by the time difference between the two rates to obtain the deformation acceleration. The rate anomaly threshold is determined by the deformation rate during the historical stable operation phase. The deformation rates during the historical stable operation phase are sorted in ascending order of absolute value, the 95th percentile is read, and multiplied by 1.1 to obtain the rate anomaly threshold. The acceleration anomaly threshold is determined by the deformation acceleration during the historical stable operation phase. The deformation acceleration during the historical stable operation phase is sorted in ascending order of absolute value, the 95th percentile is read, and multiplied by 1.1 to obtain the acceleration anomaly threshold.

[0103] When the absolute value of any deformation rate in the initial deformation evolution trajectory is greater than the rate anomaly threshold, or the absolute value of any deformation acceleration is greater than the acceleration anomaly threshold, the rate anomaly flag is written as 1; if the thresholds are not exceeded, the rate anomaly flag is written as 0. The early damage feature consists of the winding fracture risk identifier, the current cumulative deformation value, the maximum deformation rate, the maximum deformation acceleration, and the rate anomaly flag. The writing order is: winding fracture risk identifier, current cumulative deformation value, maximum deformation rate, maximum deformation acceleration, and rate anomaly flag. The numerical fields in the early damage feature are retained to three decimal places, and the flag fields use integer encoding.

[0104] When mapping industrial communication protocol fields based on the early damage characteristics, field mapping rules are pre-constructed. In the field mapping rules, the synchronization word is fixed at hexadecimal AA55, the protocol version is written as 01, the device number is read from the transformer device number, and the message type is written as 10 in the early damage warning scenario. The low-risk identifier is mapped to 01, the medium-risk identifier is mapped to 02, and the high-risk identifier is mapped to 03. The current cumulative deformation value, maximum deformation rate, and maximum deformation acceleration are each multiplied by 1000, and the integer is written as a 4-byte signed integer. The rate anomaly flag is written as a 1-byte integer. After writing according to the above mapping rules, the basic fields of the warning command message are obtained.

[0105] The preset time extraction interval is determined by the total duration of the initial deformation evolution trajectory. The total trajectory duration is obtained by subtracting the starting time coordinate from the ending time coordinate of the initial deformation evolution trajectory. The total trajectory duration is then divided by 20 and rounded up to an integer multiple of the sampling interval to obtain the preset time extraction interval. Time coordinates and deformation coordinates are read from the initial deformation evolution trajectory according to the preset time extraction interval to obtain the trajectory coordinate point set. If the extracted data does not include an ending trajectory point, the ending trajectory point is appended to the trajectory coordinate point set. Each trajectory coordinate point is written in the order of time coordinate first, followed by deformation coordinate.

[0106] When encapsulating the early damage warning instruction, a synchronization word and a message length field are first written to the frame header of the warning instruction message to obtain the message header field. The message length field is calculated from the total number of bytes in the message header field, payload field, and check field. The early damage features and the trajectory coordinate point set are written to the data payload area of ​​the warning instruction message to obtain the payload field. Cyclic redundancy check calculation uses a 16-bit check register, with the initial value written to FFFF. During calculation, each byte in the message header field and the payload field is read sequentially, and the byte is XORed with the lower 8 bits of the check register, then shifted right 8 times consecutively. If the least significant bit is 1 before each right shift, it is XORed with A001 after the right shift; if the least significant bit is 0, it is directly right-shifted. After all bytes are processed, the check field is obtained. The early damage warning instruction is encapsulated in the order of the message header field, the payload field, and the check field to obtain the early damage warning instruction.

[0107] The second embodiment of the present invention provides a transformer winding damage diagnosis system under short-circuit impact, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the method described above.

[0108] It should be noted that the transformer winding damage diagnosis system provided in this embodiment of the invention is used to execute all the process steps of the transformer winding damage diagnosis method under short-circuit impact described in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0109] It should be noted that the system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0110] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for diagnosing transformer winding damage under short-circuit impact, characterized in that, include: Obtain the original impact waveform data of each short circuit event, and extract the impact time domain features from the original impact waveform data to obtain a waveform feature set; The winding stress is solved by performing a winding stress mapping based on the waveform feature set to obtain the residual stress at each spatial node of the winding. Based on the residual stress and the pre-stored historical stress records, the frequency of stress impacts is statistically analyzed to obtain the impact frequency distribution characteristics. Then, based on the impact frequency distribution characteristics, the stress time series accumulation calculation is performed to obtain the stress superposition results of each spatial node. Based on the stress superposition results, deformation time series deduction calculations are performed to obtain the cumulative deformation values ​​of each spatial node; Based on the cumulative deformation values, fatigue fracture risk is identified to obtain the winding fracture risk identifier for each spatial node. Based on the winding fracture risk identifier, the cumulative deformation value, and the corresponding node coordinates, a deformation sequence curve is fitted to obtain the initial deformation evolution trajectory. Based on the initial deformation evolution trajectory, rate anomaly characteristics are determined to obtain early damage characteristics. Based on the early damage characteristics, an early warning instruction is encapsulated to obtain an early damage early warning instruction containing the risk node coordinates.

2. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The step of extracting time-domain features from the original impact waveform data to obtain a waveform feature set includes: Wavelet threshold denoising is performed on the original impact waveform data to obtain a smooth waveform sequence; The signal amplitude of each sampling point in the smooth waveform sequence is traversed, the maximum signal amplitude is determined as the peak amplitude, and the peak time corresponding to the peak amplitude is recorded. Search forward and backward from the peak time for the time points when the signal amplitude crosses the preset threshold for the over-limit signal, respectively, to obtain the over-limit start time and over-limit end time; Subtracting the over-limit start time from the over-limit termination time yields the duration, and concatenating the peak amplitude and the duration in dimensional order yields the waveform feature set.

3. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The step of solving the winding stress mapping based on the waveform feature set to obtain the residual stress at each spatial node of the winding includes: The waveform feature set is input into a pre-constructed electromagnetic force conversion model to solve for the electromagnetic force distribution of the winding space leakage magnetic field, thereby obtaining the transient electromagnetic force distribution matrix. Based on the transient electromagnetic force distribution matrix, the finite element displacement of each spatial node of the winding is solved to obtain the deformation displacement field. The strain distribution tensor is obtained by dividing the displacement difference between adjacent spatial nodes in the deformation displacement field by the corresponding node spacing. The strain distribution tensor is input into a pre-constructed elastoplastic bilinear kinematic reinforcement model, and constitutive relation mapping is performed on each spatial node of the winding to obtain the residual stress.

4. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The step of statistically analyzing stress impact frequencies based on the residual stress and historical stress records to obtain impact frequency distribution characteristics includes: The stress peak value and the time of occurrence of each short-circuit impact are read from the historical stress record and arranged in chronological order to obtain the historical stress event sequence. The residual stress is added to the historical stress event sequence according to the current short-circuit impact time to obtain the updated stress event sequence; According to the preset stress amplitude range, each stress peak in the updated stress event sequence is assigned a range to be marked, thus obtaining a stress range marking sequence; According to the rainflow counting rule, the stress peak and valley changes in the stress interval marking sequence are identified by closed loop, and the number of stress cycles corresponding to each preset stress amplitude interval is obtained. The frequency vector is arranged according to the stress cycle number and the time interval between adjacent impacts to obtain the impact frequency distribution characteristics.

5. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The stress time-series accumulation calculation based on the impact frequency distribution characteristics to obtain the stress superposition result of each spatial node includes: Based on the impact frequency distribution characteristics, the fatigue damage increment is calculated by combining the pre-established fatigue curve of the winding material, and the fatigue damage evolution rate is obtained. Based on the fatigue damage evolution rate, the elastic recovery stage data in the historical stress record are subjected to elastic recovery reduction correction to obtain the dynamic recovery modulus; The stiffness degradation matrix is ​​obtained by superimposing the damage ratio based on the dynamic recovery modulus and the residual stress. The cumulative amount of plastic deformation is obtained by implicit time integration based on the stiffness degradation matrix. The cumulative amount of plastic deformation is substituted into a pre-constructed exponential decay function to obtain the stress relaxation weight. The residual stress and the historical stress record are then weighted and accumulated according to the stress relaxation weight to obtain the stress superposition result.

6. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The step of performing deformation time-series extrapolation calculations based on the stress superposition results to obtain the cumulative deformation values ​​of each spatial node includes: Based on the stress superposition result, the sequence is divided by a sliding window to obtain the stress window sequence; Time-domain statistics are extracted from the stress window sequence to obtain multi-dimensional time series features including stress peak, stress change gradient, and high stress duration; The multidimensional time series features are arranged in order of sampling time to obtain the input sequence; The input sequence is fed into a pre-trained long short-term memory network for state updates to obtain the hidden layer state; The hidden layer state is input into a pre-trained multilayer perceptron regression model and nonlinearly mapped to obtain the cumulative deformation value. The high stress duration is obtained by adding up the continuous sampling durations that are higher than the preset high stress discrimination threshold.

7. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The step of identifying fatigue fracture risk based on the accumulated deformation value to obtain the winding fracture risk identifier for each spatial node includes: The accumulated deformation value and the pre-acquired initial physical property data of the winding are concatenated in the order of feature dimensions to obtain the deformation feature vector; The deformation feature vector is sorted by Euclidean distance with each sample vector in the pre-established material fatigue attenuation mapping table to obtain the nearest neighbor sample sequence. The sample attenuation coefficient is read from the nearest sample sequence according to the preset number of neighborhoods, the average value of the sample attenuation coefficient is calculated to obtain the fatigue attenuation coefficient, and the fatigue limit difference is calculated based on the fatigue attenuation coefficient to obtain the fatigue limit margin. The top oil temperature, real-time load rate, and winding vibration fundamental frequency are read from the smart sensor. The operating condition threshold is mapped to the top oil temperature, the real-time load rate, and the winding vibration fundamental frequency to obtain the real-time margin threshold. When the fatigue limit margin is less than the real-time margin threshold, the difference between the real-time margin threshold and the fatigue limit margin is mapped by a Weibull probability distribution to obtain a fracture probability value. Based on the fracture probability value, a preset probability interval is matched to obtain the winding fracture risk identifier.

8. The method for diagnosing transformer winding damage under short-circuit impact according to claim 1, characterized in that, The process of encapsulating early warning instructions based on the early damage characteristics to obtain early damage early warning instructions containing the coordinates of risk nodes includes: Based on the early damage characteristics, industrial communication protocol fields are mapped to obtain early warning command messages; According to a preset time extraction interval, the time coordinates and deformation coordinates are read from the initial deformation evolution trajectory to obtain the trajectory coordinate point set; A synchronization word and a message length field are written into the frame header of the warning instruction message to obtain the message header field; The early damage features and the trajectory coordinate point set are written into the data payload area of ​​the early warning instruction message to obtain the payload field; Cyclic redundancy check (CRC) calculations are performed on the message header field and the payload field to obtain a check field. The check field is then encapsulated in the order of the message header field, the payload field, and the check field to obtain the early damage warning instruction.

9. A transformer winding damage diagnosis system under short-circuit impact, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements the method described in any one of claims 1 to 8.

Images