A railway mileage measurement method and system based on Fourier analysis

Abstract

The present invention discloses a method and system for measuring railway mileage based on Fourier analysis; the present invention relates to the field of railway signal technology; S100, Fourier transform: for discrete time domain signal x[n], discrete Fourier transform (DFT) is performed to obtain spectrum data X[k]:; wherein, N is the number of time domain signal samples; k is the index of the frequency domain; i is the imaginary unit; e is the base of the natural logarithm; S101, time-frequency analysis: short-time Fourier transform (STFT) operation is performed to provide information of the signal in the time domain and frequency domain, and the spectrum data X[k] is corrected to be new spectrum data X[k,t]; the present invention calculates the power spectrum density of the signal in real time, and dynamically adjusts the preset threshold value T according to its statistical characteristics, and the present invention can adaptively process various signal conditions. Avoid causing the system to erroneously judge the insulation node that the train is currently passing through, thereby affecting the subsequent positioning, mileage correction and error accumulation.

Description

Technical Field

[0001] The present invention relates to the field of railway signal technology, and in particular to a railway mileage measurement method and system based on Fourier analysis. Background Art

[0002] Railway signal infrastructure systems are complex and involve many interrelated factors. To ensure train operation safety and improve operational efficiency, an electrical inspection vehicle is generally used to conduct dynamic inspections of signal trackside equipment, providing data support for safety assessment and maintenance of the signal system (for details, see this document: Zhang Wei. Development and Application of Railway Electrical Inspection Vehicles [J]. Railway Computer Applications, 2009, 18 (3): 46-48).

[0003] The above data support the acquisition of positioning data including: vehicle real-time speed, longitude and latitude, etc. The on-board TAX box data channel is mainly integrated by the vehicle real-time speed through the automatic detection system to deduce the absolute mileage corresponding to the detection data (for the specific principle, please refer to this document: Xie Baofeng, Gao Limin, Meng Jinghui, et al. High-speed comprehensive detection train CTCS-3 level train control dynamic detection system [J]. Railway Technology Innovation, 2012 (1): 85-88). Therefore, the cumulative error is the main reason for the system positioning deviation; secondly, the detection vehicle’s route in the station line is uncertain. When the detection vehicle runs through the siding, the distance calculated by the system does not correspond to the absolute mileage of the line; finally, the data error and occasional errors of the positioning system itself are also one of the reasons for the positioning deviation.

[0004] In order to verify the output results of the positioning system and correct the accumulated errors of the system and improve the fault tolerance of the system, it is necessary to correct the estimated position of the system. The traditional technical solution is to use a preset insulation node as a signal relay detection method or as a signal node, calculate by counting the insulation nodes, determine the absolute value of the correction value and compare it with the preset threshold value, judge whether the positioning error is too large or too small, and use this as a criterion to alarm and notify ground personnel or on-board engineers to intervene. Or directly count the number of signals that trigger the insulation node, and record the corresponding relationship between the number of trigger signals and the mileage in the database. By looking up the corresponding mileage correction value in the database according to the number of trigger signals received and correcting it; the specific implementation principle can be found in the following four prior art documents:

[0005] (1) Meng Jinghui. Application of insulation joint positioning in railway professional inspection vehicles [J]. Railway Standard Design, 2017, 61 (3): 148-152;

[0006] (2) Yang Ji, Xu Qingyang, Duan Hehui, et al. Research on railway signal dynamic automatic detection system [J]. Railway Technology Innovation, 2021 (6): 98-103;

[0007] (3) Zhao Wenjun, Wang Ruifeng, She Chaofu, et al. A device for obtaining mileage correction value based on frequency shift signal of railway track insulation joint [P]. CN201220647494.9, 2012-11-30;

[0008] (4) Wang Xin, Wang Yuan, Wang Ping, et al. Evaluation and correction of mileage error in high-speed railway dynamic inspection vehicle inspection data [J]. Railway Standard Design, 2018 (7): 46-50.

[0009] However, this traditional technology has the following technical problems:

[0010] 1) When the preset threshold deviates from the objective threshold, the calculated distance value may deviate from the actual distance, causing the system to incorrectly judge the insulation section that the train is currently passing through, thus affecting the subsequent positioning and mileage correction.

[0011] 2) Inaccurate calculation of correction value: Due to the deviation of the distance value, the correction value calculated based on this value will also be inaccurate. This will directly affect the correction effect of the system on the absolute mileage of the train, resulting in a large error between the corrected mileage and the actual mileage.

[0012] 3) Increased cumulative error: As time goes by and the train continues to run, since there are deviations in each correction, these deviations will gradually accumulate, causing the positioning error of the system to become larger and larger. This will seriously affect the safety and accuracy of train operation.

[0013] Specifically, if Figure 2 As shown in the figure, when the midpoint of the insulating node is regarded as a signal detection node, the carrier frequency voltage signal of the current section near the midpoint of the insulating node crosses with the carrier frequency voltage signal of the next section, and the carrier frequency switching will occur at the same time, and the compensation capacitor detection system will also output a negative pulse; since the absolute position of the insulating node is known, the traditional technical method can obtain the deviation value of the current train position based on the automatic identification of the midpoint of the insulating node (which can be further understood in conjunction with this document: Meng Jinghui. Discussion on the application of insulating node positioning in railway professional inspection vehicles [J]. Railway Standard Design, 2017, 61 (3): 148-152).

[0014] If the above relationship is algorithmized, that is, the track circuit carrier voltage at time t is , the carrier voltage of the next section is ; In the next detection cycle (i.e., time t+1), the two carrier frequency voltages are , then the voltage relationship at the intersection is:

[0015] ;

[0016] The above formula is used to determine the intersection of the two carrier frequency voltages. If the carrier frequency switching or the compensation capacitor detection data has a negative pulse at the same time, it is determined that the location corresponds to the ground insulation node equipment. The location of the insulation point can be identified by combining the actual detection data with the simulation identification system and checking with the above formula. The location of the insulation node identified in the ground simulation identification system is shown in Figure 3 and Figure 4 The small dot mark in the middle;

[0017] When an insulation node is identified, there are two possible situations:

[0018] 1) Case 1: Absolute mileage of the train (best estimate) It is known that since the insulation nodes are distributed sequentially along the line, basic information such as the name, absolute mileage, longitude and latitude of the insulation node can be obtained from the database.

[0019] 2) Case 2: The absolute mileage of the train is unknown (such as data loss due to a fault such as a reset of the onboard host), then the basic information database of this route is searched according to the inspection plan (and the breakpoint information before the reset), and the distance between the longitude and latitude of the i-th insulation node after the starting position of the starting operation (or the reset mileage breakpoint) and the longitude and latitude of the train at time t is calculated in sequence If the distance value When it is less than the preset threshold T, it can be determined that the insulation node passed at time t is i, and its basic mileage is obtained through the database.

[0020] For case 2, the common approach is to set the system mileage coordinates at time t to , the system mileage coordinates at time t+1 are , the absolute mileage of the current insulation node obtained through the database is , and then get the system mileage correction value C:

[0021] ;

[0022] If the absolute value of C is greater than the preset threshold T, the positioning error is judged to be too large, and an alarm is issued to notify ground personnel or on-board engineers to intervene. If the absolute value of C is less than the preset threshold T, C is used to correct the current absolute mileage.

[0023] However, the practical effect of the above implementation is not good, because the setting of the preset threshold T is generally based on subjective assignment. If the setting of the preset threshold T deviates from the objectively reasonable threshold, the distance value The deviation of the judgment of the correction value C will further cause the system to incorrectly judge the insulation section that the train is currently passing, which will affect the subsequent positioning and mileage correction. If the deviation continues, since there is a deviation in each correction, these deviations will gradually accumulate, causing the positioning error of the system to become larger and larger. This will seriously affect the safety and accuracy of train operation.

[0024] That is, the technical problem to be solved by this application is: how to formulate an accurate preset threshold.

[0025] To this end, the present invention proposes a railway mileage measurement method and system based on Fourier analysis. Summary of the invention

[0026] In view of this, the embodiment of the present invention hopes to provide a method and system for measuring railway mileage based on Fourier analysis to solve or alleviate the technical problems existing in the prior art, namely, how to formulate an accurate preset threshold value, and at least provide a beneficial option for this;

[0027] The technical solution of the embodiment of the present invention is achieved as follows:

[0028] First, a method for measuring railway mileage based on Fourier analysis:

[0029] This technical solution uses the PSD algorithm (Power Spectral Density, PSD) to measure and assign a value to the preset threshold T. The PSD algorithm is a Fourier analysis technique that can understand the power distribution of a signal at different frequencies by calculating the power spectral density of the signal, thereby determining which frequency components are primary and which are secondary. This helps to pay more attention to the primary frequency components and ignore the secondary components when setting the preset threshold T, thereby improving the accuracy of the threshold setting. The PSD algorithm specifically includes three main steps, namely, steps S1 to S3.

[0030] (1) Preliminary steps:

[0031] First, the ground simulation recognition system obtains the time domain signal and air domain signal generated after the train passes through the insulation joint; the time domain signal should reflect the change of voltage parameters when the train passes through the insulation joint. The air domain signal should reflect the geographical coordinate position of the insulation joint, that is, the geographical coordinates of different insulation points along the railway line. Then the collected time domain signal is preprocessed, including conventional filtering, denoising and normalization steps, to reduce the interference components in the signal and highlight the main features.

[0032] (2) Step S1, Fourier analysis:

[0033] Apply Fourier transform, for the time domain signal, convert it to the frequency domain and obtain the spectrum data X[k]; for the spatial domain signal, associate it with the spectrum data X[k] to obtain a new data structure X'[k,t,s]; the above-mentioned "convert", "associate" and "obtain" steps are S100~S102 respectively.

[0034] (2.1) Step S100, Fourier transform:

[0035] For the discrete time domain signal x[n], perform discrete Fourier transform (DFT) to obtain the spectrum data X[k]:

[0036] ;

[0037] Where N is the number of time domain signal samples; k is the index of the frequency domain; i is the imaginary unit; e is the base of the natural logarithm (taken as 2.71828);

[0038] X[k] is the spectrum data, which represents the result of Fourier transform, that is, the magnitude (amplitude and phase) of the kth frequency component in the frequency domain;

[0039] n is the discrete time index of the time domain signal. In discrete time signal processing, the continuous time variable t is replaced by a discrete point n, where n ranges from 0 to N-1.

[0040] (2.2) Step S101, time-frequency analysis:

[0041] A short-time Fourier transform (STFT) operation is performed to provide information of the signal in the time domain and the frequency domain, and the spectrum data X[k] is corrected to be new spectrum data X[k, t]; specifically, it includes three sub-steps, namely, steps S1010 to S1012.

[0042] (2.2.1) Step S1010, define windowing and window functions:

[0043] Select a window function w[n] which determines the data segment to be analyzed at each time in the time domain.

[0044] For the discrete time domain signal x[n], multiply it with the window function to obtain the windowed signal . This step is to focus only on the part within the window when analyzing the signal:

[0045] ;

[0046] Where t is the starting position (or time index) of the window, which changes as the sliding window moves.

[0047] (2.2.2) Step S1011, executing Discrete Fourier Transform (DFT) algorithm:

[0048] ;

[0049] X[k,t] is the new spectrum data X[k,t], which is the result of short-time Fourier transform, where k is the index in the frequency domain and t is the starting position of the window (time index).

[0050] It should be noted that this formula is very similar to the expression in S100, but it is only applied to the windowed signal.

[0051] (2.2.3) Step S1012, iteration:

[0052] Move the sliding window to the next position, that is, by increasing the window starting position t by a step size (e.g., half the window size or less), and then repeating steps S1010 and S1011 until the entire signal has been analyzed.

[0053] The final result of the STFT is a series of new spectrum data X[k,t], where each spectrum corresponds to the frequency content of the signal in different time windows. These spectrum data can be viewed as a two-dimensional array or matrix, one dimension of which is the frequency index k and the other dimension is the time index t. These new spectrum data X[k,t] provide a time-frequency representation of the signal, correcting or supplementing the original spectrum data X[k] obtained only based on the DFT of the entire signal.

[0054] (2.3) Step S102, associating airspace information:

[0055] This is a data integration step. When executing S100 and S101, the spatial information S[t] of the spectrum data X[k] and the new spectrum data X[k, t] is obtained, including the geographical location of the isolation node signal source or the isolation node number at each time point t. Then data integration is performed:

[0056] S1020. For each time point t, record new spectrum data X[k, t].

[0057] S1021. Associate the spatial domain information S[t] corresponding to the time point t with the new spectrum data X[k,t], that is, create a data structure X'[k,t,s], which includes the original spectrum data and the corresponding spatial domain information.

[0058] In order to avoid confusion, the data structure X'[k,t,s] will be referred to as the "new data structure X'[k,t,s]" in the following text;

[0059] For S1021, a new spectrum data X[k,t] in the form of a two-dimensional array is used to store the spectrum data, and another one-dimensional array of spatial information S[t] or a more advanced data structure (such as a dictionary, which maps time points to spatial information) is used to store the spatial information. Then, when actually processing the data, the information in S[t] can be combined with the corresponding data in X[k,t] as needed.

[0060] Integrate the spatial domain information into a single data structure, and assume that each time point t has a corresponding spatial domain index s. The three-dimensional array is as follows:

[0061] ;

[0062] X'[k,t,s] is considered a three-dimensional array in this expression, where for each valid (k,t) pair, s is the spatial index obtained by mapping the time point t to the index s corresponding to the spatial information S[t]. If the [k,t,s] combination is invalid (for example, s is not obtained by a valid mapping), then X'[k,t,s] can be set to 0, NaN coefficient, or other placeholder representing invalid data.

[0063] (3) Step S2, analyze the power spectrum density:

[0064] Based on the spectrum data (ie, the new data structure X'[k, t, s] described above), the spatiotemporal power spectrum density, ie, the distribution of signal power in frequency, is calculated by the PSD algorithm. Specifically, the process includes three sub-steps, namely, steps S200 to S202.

[0065] (3.1) Step S200, calculate the estimated value P[k] of the power spectrum density:

[0066] For the discrete data structure X'[k,t,s], its power spectral density is estimated as:

[0067] ;

[0068] Where P[k] represents the estimated value of the power spectral density at frequency index k, is the square of the modulus of the frequency domain signal.

[0069] The above formula is a standard and classic way to calculate power spectrum density, but it may not be perfectly applicable to the spatial characteristics considered in this technical solution, because the expression focuses on the average power distribution of the signal within a certain frequency range, rather than the power at a single frequency point.

[0070] Therefore, the present technical solution further provides another better implementation method: each spatial point (insulation point) is considered as an independent signal source, then the power spectral density (PSD) is calculated independently for each spatial point, and finally these PSDs are combined to form a comprehensive PSD or kept separate for subsequent analysis.

[0071] (3.1.1) Step S2000, independent initialization:

[0072] A specific spatial point s (i.e., a specific spatial index s) is selected, corresponding to an isolation point. Since each isolation point is regarded as an independent signal source, the PSD will be calculated separately for each isolation point.

[0073] By fixing the spatial point s and traversing k and t, the signals of all frequencies and time points corresponding to the selected spatial point s are extracted from the multidimensional data structure X'[k,t,s].

[0074] (3.1.2) Step S2001, calculate the power spectrum density of a single spatial point:

[0075] ;

[0076] in, is the power spectral density estimate for a spatial point s at frequency index k, is the total number of time points used to calculate the PSD (i.e. the length of the signal in time or the length of the time window used for averaging), is the square of the signal modulus at a given frequency k (frequency domain index k), time t (window t), and spatial point s (spatial domain index s).

[0077] (3.1.3) Step S2002, combining the power spectrum densities of all spatial points:

[0078] There are three options to choose from:

[0079] 1) Subjective selection scheme: Based on the subjective initiative of the technician in the field, keep the PSD of each spatial point separate and analyze them independently. It is suitable for situations where the spectral characteristics of a specific insulation point are of interest. Then an estimate P[k] is independently formulated.

[0080] 2) Average power spectral density scheme: an overall PSD estimate can be obtained instead of being specific to a certain insulation point:

[0081] ;

[0082] in, is the total number of points in space, is the average power spectral density at each spatial point.

[0083] 3) Weighted Power Spectral Density Scheme: If some spatial points have different importance or uncertainty, the power spectral density can be combined using weighted averaging, where the contribution of each spatial point is adjusted according to its weight. The final output is a weighted estimated value P[k].

[0084] (3.2) Step S201, calculate the spatiotemporal power spectral density:

[0085] For the data structure X'[k,t,s] under the estimated value P[k], the power distribution of the time domain signal and the spatial domain signal in time and space is further calculated to form the time-space power distribution P[k]'. This helps to analyze the energy changes of the signal at different locations and times:

[0086] ;

[0087] Among them, the time-space power distribution P[k]' is the time-frequency power distribution at the frequency index k and window t (regarded as a time point). It gives the signal power at a specific time, frequency and spatial position.

[0088] However, it should be pointed out that the spatiotemporal power distribution P[k]' at this time still contains the spatial domain index s, so P[k]' at this time should be more appropriately regarded as a time-frequency power distribution in the Fourier correlation concept.

[0089] Since this technical solution needs to calculate the total spatiotemporal power distribution, if we want to further and more perfectly obtain spatial information, we need to consider how to process the spatial information, that is, to obtain an overall spatiotemporal power distribution P'[k,t] that does not depend on a specific spatial position, so we can average all spatial positions:

[0090] ;

[0091] in, is the total number of spatial points.

[0092] Finally, the overall spatiotemporal power distribution P'[k, t] can be further provided as a two-dimensional matrix, where each element represents the signal power at a given time point and frequency index. This matrix can be used to visualize the energy distribution of the signal in time and frequency.

[0093] (3.3) Step S202 (optional step), considering the propagation characteristics of the signal:

[0094] The signal of the isolation node may attenuate with the propagation distance, which should be considered when calculating the power spectrum density. The attenuation of the wireless signal can be expressed as a function of distance. This technical solution regards it as an exponential decay model, which means that the signal strength decreases exponentially with the increase of distance. In this way, the overall spatiotemporal power distribution P'[k, t] of step S201 is corrected.

[0095] Let the attenuation function A(d) represent the attenuation factor of the signal after propagation distance d, where d is the distance from the signal source to the receiving point. The exponential attenuation model is:

[0096] ;

[0097] Where α is the attenuation coefficient, which determines the rate at which the signal decays; e is the base of the natural logarithm (taken as 2.71828).

[0098] In order to apply this attenuation function to the overall spatiotemporal power distribution P'[k, t], it is necessary to know the distance corresponding to each spatial point s . Then use the decay function (decay factor) To correct the power value of each spatial point. However, since the overall spatiotemporal power distribution P'[k, t] has been averaged for all spatial points in step S201, it is necessary to combine this averaging process with the attenuation correction.

[0099] The operation mode of the technical solution is: in step S201, instead of directly performing spatial averaging, the weighted power of each spatial point is calculated first, and then the average is performed. In this way, the attenuation function can be directly applied to each spatial point. The corrected overall spatiotemporal power distribution P''[k, t] is calculated as:

[0100] ;

[0101] The power at each spatial point is The corresponding attenuation factor The weights are then averaged. Therefore, the contribution of farther spatial points to the overall power is reduced, reflecting the characteristic of signal attenuation with distance.

[0102] (4) Step S3, determine the preset threshold T:

[0103] Based on the statistical characteristics of the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P'[k,t], the preset threshold T is calculated for subsequent signal processing or decision making. Specifically, three schemes are included, which are represented by S300, S301 and S302 respectively.

[0104] (4.1) Step S300, single threshold solution:

[0105] A preset threshold T is set according to the distribution characteristics of the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P''[k, t]. The preset threshold T is set based on the statistical characteristics (mean or standard deviation) in the power spectrum density. Taking the overall spatiotemporal power distribution P''[k, t] as an example, it includes:

[0106] 1) Preset threshold based on mean: Assume that the background level of the power spectrum density is approximated by its mean. The preset threshold T can be set to the mean of the overall spatiotemporal power distribution P''[k, t] plus an offset, which can be fixed or based on a percentage of the mean:

[0107] T = μ + offset;

[0108] where μ is the mean of the overall spatiotemporal power distribution P''[k, t], calculated as:

[0109] ;

[0110] in, is the total number of frequency indices, is the total number of time points;

[0111] Here, offset is a positive offset, which is a fixed percentage of the mean (such as 10%):

[0112] offset=0.1μ;

[0113] Or it can be a fixed constant that can be adjusted according to the specific application.

[0114] 2) Preset threshold based on standard deviation:

[0115] This method takes into account the fluctuation of power spectral density. The preset threshold T can be set to the mean of the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P''[k, t] plus several times the standard deviation. This can capture more significant signal changes relative to the background level:

[0116] T = μ + σn;

[0117] where μ is the mean of the overall spatiotemporal power distribution P''[k, t], as described above; σ is the standard deviation of the overall spatiotemporal power distribution P''[k, t], calculated as:

[0118] ;

[0119] In this formula, n is a positive coefficient, for example, it can be 2 or 3, indicating that the preset threshold T is the mean plus 2 or 3 times the standard deviation.

[0120] (4.2) Step S301, dynamic threshold scheme:

[0121] The preset threshold T is dynamically adjusted according to the real-time statistical characteristics of the signal, that is, the mean value within the sliding window.

[0122] Suppose there is a sliding window with a width of W, which contains the power spectrum density data of the most recent W time points. For a given current spatiotemporal power distribution P[k]' or overall spatiotemporal power distribution P''[k, t], calculate the mean within the sliding window , and then use this mean to dynamically set the dynamic preset threshold T[t]; taking the overall spatiotemporal power distribution P'[k,t] as an example, the steps include:

[0123] S3010, calculating the mean value in the sliding window:

[0124] ;

[0125] in, is the total number of frequency indices k, and τ is the time point operation traversed within the sliding window.

[0126] S3011, set dynamic preset threshold:

[0127] Based on the mean within the sliding window , set the dynamic preset threshold T[t] as:

[0128] ;

[0129] or,

[0130] ;

[0131] in, is a positive offset, or is a positive percentage value. These two parameters can be adjusted according to the needs of specific applications.

[0132] (4.3) Step S302, multi-threshold scheme:

[0133] For different time domain and space domain conditions, multiple preset thresholds T can be set to form a threshold set. In practical applications, a suitable preset threshold T is selected according to the current state of the signal.

[0134] Second, a railway mileage measurement system based on Fourier analysis:

[0135] The system includes insulation joints arranged along the railway and at the back-end control end.

[0136] A ground simulation identification system, the ground simulation identification system is connected to the insulation node via wireless communication;

[0137] The ground simulation identification system is connected to a processor, and the processor is connected to a register;

[0138] The processor is used to analyze the electrical signal of the insulation node transmitted by the ground simulation identification system, wherein the electrical signal includes a time domain signal and a space domain signal;

[0139] The time domain signal and the air domain signal are as described in the preceding steps above, the time domain signal is the change of voltage parameters when the train passes through the insulation joint, and the air domain signal is the geographical coordinate position of the insulation joint, that is, the geographical coordinates of different insulation points along the railway line.

[0140] The register stores program instructions, and when the program instructions are executed by the processor, the processor executes the railway mileage measurement method based on Fourier analysis as described above.

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

[0142] 1. Objective signal processing: The present invention calculates the power spectrum density of the signal in real time and dynamically adjusts the preset threshold T according to its statistical characteristics. The present invention can adaptively process various signal conditions. That is, it can automatically adapt to different operating environments without manual intervention or frequent threshold adjustment. At the same time, by accurately setting the preset threshold T, the present invention can more accurately detect important features or events in the signal. This helps to reduce false positives and false negatives of insulation node mileage information, thereby improving the reliability and performance of the overall system.

[0143] 2. Enhanced robustness: Since the preset threshold T of the present invention is set based on the statistical characteristics of the real-time signal, the present solution has strong robustness to noise and signal fluctuations. This ensures that the system can maintain stable performance when facing various interferences and challenges. It avoids causing the system to mistakenly judge the insulation section that the train is currently passing through, thereby affecting the subsequent positioning, mileage correction and error accumulation.

[0144] 3. Flexibility and scalability: If the present invention selects a dynamic threshold solution or a multi-threshold solution, this flexibility enables the system to easily adapt to different application scenarios and requirements, while providing scalability to introduce more signal features and threshold conditions in the future. BRIEF DESCRIPTION OF THE DRAWINGS

[0145] Figure 1 It is a schematic diagram of the overall method flow of the present invention;

[0146] Figure 2 This is an example diagram of the test data at the insulation joint;

[0147] Figure 3 An example diagram of an insulation node identified in a ground simulation;

[0148] Figure 4 Another example of an insulation node identified in a ground simulation.

[0149] Figure 5 is a schematic diagram of the method flow of step S1 of the present invention;

[0150] Figure 6 This is a schematic diagram of the method flow of step S2 of the present invention;

[0151] Figure 7 It is a flowchart of the spatial domain index mapping mechanism in step S1021 of the present invention;

[0152] Figure 8 It is a schematic diagram of the system composition of the present invention. DETAILED DESCRIPTION

[0153] In order to make the above-mentioned purposes, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below;

[0154] Example 1: Please refer to Figure 1 When the time domain signal and air domain signal generated by the train passing through the insulation node are obtained based on the ground simulation recognition system, the collected time domain signal can be preprocessed first (because the geographical location of the insulation node is absolute, the air domain signal does not need to be preprocessed), including conventional filtering, denoising and normalization techniques to reduce the interference components in the signal and highlight the main features. Then execute steps S1~S4 to determine whether it is necessary to correct the absolute mileage of the train. Among them, steps S1~S3 will analyze the information of the correspondence between traditional insulation nodes and mileage through Fourier analysis technology, and assign a value to the preset threshold T. After step S4 calculates the train mileage correction value C, it determines whether it is necessary to correct the absolute mileage of the train based on the preset threshold T.

[0155] In this example, see Figure 5 Regarding step S1, Fourier analysis: for the time domain signal, convert it to the frequency domain and obtain the spectrum data X[k]; for the spatial domain signal, associate it with the spectrum data X[k] to obtain a new data structure X'[k,t,s]. It specifically includes three sub-steps, steps S100~S102:

[0156] Specifically, step S100, Fourier transform: for the discrete time domain signal x[n] captured from the insulation node, convert it into the frequency domain through discrete Fourier transform (DFT):

[0157] ;

[0158] in,

[0159] N represents the total number of time domain signal samples, that is, the number of data points collected from the insulation node in a specific time period.

[0160] k is an index in the frequency domain, used to identify different frequency components. In DFT, the value of k ranges from 0 to N-1.

[0161] i is an imaginary unit, satisfying: .

[0162] e is the base of natural logarithms, approximately equal to 2.71828.

[0163] X[k] is the spectrum data, which represents the size of the kth frequency component in the frequency domain after Fourier transform (including amplitude and phase information).

[0164] n is the discrete time index of the time domain signal, which is used to identify different time points. In discrete time signal processing, the continuous time variable t is replaced by a series of discrete points n, where n ranges from 0 to N-1.

[0165] It can be understood that DFT can convert the time domain signal into a set of spectrum data X[k], which reveals the intensity and phase information of different frequency components in the signal. For analyzing the specific frequency patterns generated when the train passes through the insulation section, these patterns are directly related to the passing event of the train and can be used for subsequent mileage determination.

[0166] Specifically, step S101, time-frequency analysis: perform short-time Fourier transform (STFT) operation to provide information about the signal in the time domain and frequency domain, and correct the spectrum data X[k] to the new spectrum data X[k, t]; this step can more accurately describe the characteristics of the signal in the time domain and frequency domain. STFT can reveal the frequency components of the signal in different time periods, which is particularly important for analyzing non-stationary signals (such as signals generated when a train passes through an insulating section). It further includes three sub-steps, steps S1010 to S1012.

[0167] Further, step S1010, define windowing and window function: first, select a window function w[n], which is used to define the data segment of each analysis in the time domain. The window function w[n] determines the trade-off between frequency resolution and time resolution.

[0168] The window function is a Hanning window or a Hamming window. The length of the window function is set as the window size or N. It determines the number of data points analyzed each time. Then, the discrete time domain signal x[n] is multiplied by the window function to obtain the windowed signal The purpose of this step is to focus only on the data within the window when analyzing the signal, so that the frequency components of the signal can be analyzed independently in different time periods.

[0169] Preferably, the window function w[n] is a Hamming Window:

[0170] ;

[0171] Where n is the discrete time index within the window, ranging from 0 to N-1; N is the window size, that is, the total number of data points within the window. This expression describes the shape of the Hamming window in the time domain.

[0172] Further, in step S1010, the discrete time domain signal x[n] is multiplied by the Hamming window function w[n] to obtain a windowed signal . This process is to focus only on the part of the signal within the window when analyzing it:

[0173] ;

[0174] Where t is the starting position (or time index) of the window, and the value of t changes as the sliding window moves across the signal. This expression ensures that only data points within the window are included in each analysis.

[0175] It is understandable that the purpose of windowing is to focus only on the data within the window when analyzing the signal, so that the frequency components of the signal can be analyzed independently in different time periods. By multiplying the signal with the window function, the spectrum leakage caused by the signal discontinuity at the edge of the window can be reduced, which helps to improve the accuracy of frequency analysis. The Hamming window is particularly suitable for application scenarios that require balancing the main lobe width and sidelobe attenuation. Through its specific time domain shape, the Hamming window provides an effective tool for time-frequency analysis of signals in the railway mileage measurement method based on Fourier analysis, ensuring accurate analysis of the signal frequency components in different time periods.

[0176] Furthermore, in step S1011, a discrete Fourier transform (DFT) is performed to obtain the frequency components of the signal in a specific time window:

[0177] ;

[0178] X[k,t] is the new spectrum data X[k,t], which is the result of short-time Fourier transform; k is the index in the frequency domain; t is the starting position of the window (time index), which changes as the sliding window moves.

[0179] It can be understood that the windowed time domain signal is transformed by DFT After conversion to the frequency domain, a set of spectrum data X[k, t] is obtained. These data reveal the frequency characteristics of the signal in a specific time window, including the amplitude and phase information of each frequency component. In the railway mileage measurement method, since specific frequency patterns are generated when a train passes through an insulating section, these patterns are directly related to the passing event of the train. Therefore, by analyzing these spectrum data, characteristic information related to the train passing event can be extracted, which can be used for subsequent mileage measurement. Therefore, step S1011 converts the windowed time domain signal to the frequency domain by performing DFT, and obtains a set of spectrum data with time-frequency characteristics. These data provide important characteristic information and analysis basis for subsequent railway mileage measurement.

[0180] It should be noted that the DFT here is performed on the windowed signal, so it can provide the frequency information of the signal in a local time period. Compared with the spectrum data obtained by performing DFT on the entire signal, this has a higher time resolution and can more accurately describe the time-varying characteristics of the signal.

[0181] It should be noted that the expression in step S1011 is very similar to the expression in step S100, but step S1011 is only applied to the windowed signal.

[0182] Further, step S1012, iteration: moving the sliding window to the next position, that is, by increasing the window start position t by a step size (for example, half of the window size or less), and then repeating steps S1010 and S1011 until the entire signal is analyzed.

[0183] The final result of the STFT is a series of new spectrum data X[k,t], where each spectrum corresponds to the frequency content of the signal in different time windows. These new spectrum data X[k,t] provide a time-frequency representation of the signal, correcting or supplementing the spectrum data X[k] originally obtained based only on the DFT of the entire signal.

[0184] It is understandable that these new spectrum data X[k, t] provide a time-frequency representation of the signal, that is, they describe the characteristics of the signal in both the time domain and the frequency domain. Compared with the spectrum data X[k] originally obtained based only on the DFT of the entire signal, these time-frequency representation data provide richer information and can reveal the frequency changes of the signal in different time periods. This is particularly important for analyzing non-stationary signals (such as the signals generated when a train passes through an insulating section), because the characteristics of these signals may change over time.

[0185] Preferably, these spectrum data are integrated into a two-dimensional array or matrix, one dimension of which is the frequency index k and the other dimension is the time index t. Taking a matrix as an example, its form is as follows:

[0186] ;

[0187] In this matrix, each element X[k, t] represents the spectrum value at time index t and frequency index k. The number of rows K of the matrix represents the number of frequency components, and the number of columns T represents the number of time windows, that is, the number of time periods in which the signal is analyzed.

[0188] In this matrix, k ranges from 0 to K-1, where K is the total number of frequency components. It represents the index of different frequencies and can be used to identify specific frequency components present in the signal.

[0189] In this matrix, t ranges from 0 to T-1, where T is the total number of time windows. It represents the time index, i.e., the different time periods during which the signal is analyzed.

[0190] As can be understood, this matrix provides a time-frequency representation of the signal, allowing personnel to observe the behavior of the signal at different times and frequencies simultaneously. In railway mileage determination, this time-frequency analysis is critical to accurately identify the specific frequency patterns generated when a train passes through an insulating joint, as these patterns may change over time and are an important basis for mileage determination.

[0191] In the railway mileage measurement method, these time-frequency representation data are crucial for accurately determining the location and time of the train passing through the insulation joint. By analyzing these data, staff can extract characteristic information related to the train passing event, such as the occurrence time and duration of specific frequency patterns. This information can be used in subsequent mileage measurement algorithms to improve the accuracy and reliability of the measurement.

[0192] Furthermore, the MATLAB execution program of the above steps S100-S101 is as follows:

[0193] % Let x be a discrete time domain signal and Fs be the sampling frequency

[0194] % N is the window size, step is the step size (half the window size or less)

[0195] % S1010: Windowed

[0196] windowSize = N; % window size

[0197] stepSize = step; % step length

[0198] nfft = 2^nextpow2(windowSize); % Number of points used for DFT, power of 2

[0199] window = hamming(windowSize); % Hamming window

[0200] % Initialize the result matrix

[0201] numWindows = 1 + floor((length(x) - windowSize) / stepSize); % Calculate the number of windows

[0202] X = zeros(nfft, numWindows); % Preallocate the result matrix

[0203] % S1011&S1012: Perform STFT (including DFT and iteration)

[0204] for t = 1:stepSize:numWindows*stepSize

[0205] % Ensure that the window does not exceed the signal boundary

[0206] windowIndex = min(max(1, t), length(x));

[0207] xw = x(windowIndex:min(windowIndex+windowSize-1, length(x))) .*window;

[0208] % Execute DFT

[0209] X(:, (t-1) / stepSize + 1) = fft(xw, nfft);

[0210] end

[0211] % X contains the STFT results of the signal, each column is the spectrum of a time window

[0212] In the above program, you need to fill in the following content to run:

[0213] 1) x: should be the input signal, a discrete time series;

[0214] 2) Fs: sampling frequency;

[0215] 3) N: Window size, which determines the number of data points analyzed each time;

[0216] 4) step: step size, which determines the interval of window movement;

[0217] 5) nfft: The number of points used when performing the DFT (the default in the above program is a power of 2), which can improve computational efficiency. Although nfft can be larger than windowSize, the extra points will not contain useful information.

[0218] The program will eventually output a result matrix X whose size is nfft * numWindows, where numWindows is calculated based on the signal length and step size.

[0219] In the above program, the signal x is divided into multiple overlapping segments, and the size of each segment is determined by windowSize. These segments are multiplied by a window function (here, a Hamming window) to reduce spectral leakage. The window function smoothly transitions to zero in the time domain, thereby reducing spectral interference caused by signal discontinuities. DFT (S1011) performs a fast Fourier transform (FFT) on each windowed signal segment to convert it from the time domain to the frequency domain. The fft function is a standard function in MATLAB for performing DFT. Iteration (S1012) moves the window and repeats the windowing and DFT steps until the entire signal is analyzed. The step size determines how fast the window moves and the time resolution of the result. A smaller step size will provide more time resolution, but the computational cost will also increase.

[0220] Specifically, step S102, associating spatial information: This is a data integration step. First, when executing S100 and S101, the spatial information S[t] of the spectrum data X[k] and the new spectrum data X[k, t] is obtained;

[0221] Specifically, when executing S100 and S101, not only the spectrum data is obtained, but also the geographical location or the number of the isolation node signal source at each time point t is recorded, and this information constitutes the spatial information S[t]. The spatial information S[t] includes the geographical location or the number of the isolation node signal source at each time point t. Then the data is integrated, and its steps include:

[0222] S1020. For each time point t, record new spectrum data X[k, t].

[0223] S1021. Associate the spatial domain information S[t] corresponding to the time point t with the new spectrum data X[k,t], that is, create a data structure X'[k,t,s], which includes the original spectrum data and the corresponding spatial domain information.

[0224] Furthermore, for S1021, a new spectrum data X[k,t] in the form of a two-dimensional array is used to store the spectrum data, and another spatial information S[t] in the form of a one-dimensional array or a more advanced data structure (such as a dictionary, which maps time points to spatial information) is used to store the spatial information. Then, when actually processing the data, the information in S[t] can be combined with the corresponding data in X[k,t] as needed.

[0225] It can be understood that when using a dictionary to map a time point to spatial information, a key-value pair data structure is actually created, where the key is the time point t and the value is the corresponding spatial information S[t].

[0226] For example, a dictionary time_to_space_map is set up in Python, which maps time points to spatial information. This dictionary is constructed as:

[0227] time_to_space_map = {

[0228] t0: S[t0],

[0229] t1: S[t1], ...

[0231] tn: S[tn]

[0232] }

[0233] Among them, t0, t1, ..., tn is a sequence of time points, and S[t0], S[t1], ..., S[tn] is the corresponding airspace information at these time points. The airspace information can be the geographical location, number or other relevant identification of the insulation node.

[0234] As a further example, instead of manually creating a mapping for each time point, you can use Python to use a loop or other automated method to fill the dictionary. Suppose there is a list times containing time points and a list space_info containing the corresponding space information, you can use the zip function and dictionary comprehension to create a mapping:

[0235] times = [t0, t1, t2, ..., tn]

[0236] space_info = [S[t0], S[t1], S[t2], ..., S[tn]]

[0237] time_to_space_map = {t: s for t, s in zip(times, space_info)}

[0238] Among them, the time_to_space_map dictionary contains the mapping of all time points to spatial domain information.

[0239] When you need to find the corresponding spatial information at a given time point t, you can enter the following command to achieve it: space_info_at_time_t = time_to_space_map.get(t, None)

[0240] Among them, time_to_space_map.get(t, None) will try to get the value of key t from the dictionary. If this key does not exist in the dictionary, it will return None (or you can specify another default value).

[0241] Furthermore, the spatial information is integrated into a single data structure, and each time point t has a corresponding spatial index s. The three-dimensional array is as follows:

[0242] ;

[0243] X'[k,t,s] is treated as a three-dimensional array in this expression, where for each valid (k,t) pair, s is the spatial index obtained by mapping the time point t to the index s corresponding to the spatial information S[t]. If the [k,t,s] combination is invalid (for example, because s is not obtained by a valid mapping), then X'[k,t,s] can be set to 0, NaN coefficients, or other placeholders representing invalid data.

[0244] Preferably, the mapping mechanism of the spatial index (s) obtained by mapping is:

[0245] P1. Data preparation: The system first collects or generates a data set containing timestamps and spatial information (such as insulation node numbers). These data come from a sensor network deployed along the railway line, and each sensor can detect and report the signal of the insulation node passing by.

[0246] P2. Time-spatial mapping table: A "time-spatial mapping" table is constructed based on the collected data. This table maps each time point t at which an insulation node signal is detected to a unique spatial index s. The spatial index can be an insulation node identifier (ID), location coordinates, or other information that can uniquely refer to a specific insulation node.

[0247] P3, mapping function f: This function accepts a time point t as input and returns the corresponding spatial index s. Figure 7 , the steps are:

[0248] P300, first use the binary search method as the initial search step, the binary search method is applied to the sorted list or array, and the element is quickly located by continuously halving the search range. In this case, the specific implementation of the mapping function f is described as the following steps:

[0249] P3000. Determine the search range: Initialize two pointers, one pointing to the starting position of the time-space mapping table and the other pointing to the end position.

[0250] P3001, Binary search: Calculate the time point at the middle position and compare it with the target time point t. If the middle time point is equal to t, return the corresponding spatial index s. If the middle time point is less than t, continue searching in the right half; if it is greater than t, search in the left half.

[0251] P3002. Repeat P3000 to P3001 until an exact match is found or it is determined that no exact match exists.

[0252] Among them, the execution program (MATLAB) of P300 is as follows:

[0253] function idx = binarySearch(arr, target)

[0254] % arr: ordered array

[0255] % target: the target value to be found

[0256] % idx: The index of the target value in the array, if not found, returns empty

[0257] % Initialize array bounds

[0258] low = 1;

[0259] high = numel(arr);

[0260] % Loop through binary searches until the target is found or determined not to exist

[0261] while low<= high

[0262] % Calculate the middle position

[0263] mid = floor((low + high) / 2);

[0264] % Check the middle element

[0265] if arr(mid) == target

[0266] idx = mid; % Find the target and return the index

[0267] return;

[0268] elseif arr(mid) <target

[0269] % If the target value is greater than the middle element, search the right half

[0270] low = mid + 1;

[0271] else

[0272] % If the target value is less than the middle element, search the left half

[0273] high = mid - 1;

[0274] end

[0275] end

[0276] % If the target value is not found at the end of the loop, return empty

[0277] idx = [];

[0278] end

[0279] When using the above function, you need to pass in an ordered array arr and the target value target you want to find. The function returns the index idx of the target value in the array, or an empty array [] if the value is not found. The input array arr is sorted in ascending order. If the array is sorted in descending order, you need to adjust the comparison operator accordingly.

[0280] When the binary search does not find an exact matching time point, the following P301 uses polynomial interpolation to estimate the spatial domain index. Of course, you can also directly ignore the above step P300 and directly execute P301.

[0281] P301, Polynomial interpolation: If there are only two known time-space data points and ,in , use linear interpolation to estimate the spatial index s corresponding to time point t:

[0282] ;

[0283] P302, Lagrange interpolation method: If there are more known data points and you want to perform more accurate interpolation, use a higher-order Lagrange interpolation method, that is, for a given n+1 data points ,in :

[0284] ;

[0285] Then, the corresponding spatial index s is estimated by substituting the time point t into L(t): s = L(t);

[0286] In this example, see Figure 6 , regarding step S2, analyze the power spectral density (PSD): based on the new data structure X'[k,t,s], calculate the spatiotemporal power spectral density through the PSD algorithm to form the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P'[k,t]. Specifically, it includes three sub-steps, namely steps S200~S202. Power spectral density is a measurement value that describes the power distribution of a signal or time series in the frequency domain. In the railway mileage measurement, by analyzing the power spectral density of the insulation node signal, important information about the frequency components of the signal can be obtained, which can be used to judge the railway mileage and determine the location of the insulation node.

[0287] Specifically, step S200, calculate the power spectrum density: for a given discrete data structure X'[k,t,s], first calculate its power spectrum density. Here X'[k,t,s] is regarded as a three-dimensional array, where k represents the frequency index, t represents the time index, and s represents the spatial index (i.e., the position index of the insulation node), and its power spectrum density is estimated as:

[0288] ;

[0289] in,

[0290] P[k] represents the estimated value of the power spectral density at frequency index k. This is a scalar value that describes the power strength of the signal at that particular frequency;

[0291] is the square of the modulus of the frequency domain signal. The modulus square calculation here is to take the modulus (i.e. amplitude) of each complex element in X'[k,t,s] and then square it. This is done to obtain the power value corresponding to each frequency component;

[0292] N is a normalization factor in this formula only, which is used to ensure the correct normalization of the power spectrum density. In discrete signal processing, N represents the number of data points or other relevant parameters to ensure that the calculation results are within a reasonable range.

[0293] It is understandable that the above formula is a standard and classic power spectrum density calculation formula, but it may not be perfectly adapted to the spatial characteristics considered in this embodiment. Because the expression focuses on the average power distribution of the signal within a certain frequency range, rather than the power of a single frequency point. Therefore, if necessary, the technical solution provided in the full text of Embodiment 2 can also be used to obtain a more accurate estimate P[k].

[0294] Specifically, step S201 calculates the spatiotemporal power spectral density: starting from the data structure X'[k,t,s] under the estimated value P[k]. X'[k,t,s] in this step is regarded as a three-dimensional array. The goal of this step is to calculate the power distribution of the signal in time and space, form the spatiotemporal power distribution P[k]', and analyze the energy changes of the signal at different positions and times:

[0295] ;

[0296] Among them, the time-space power distribution P[k]' is the time-frequency power distribution at frequency index k and window t (regarded as a time point). It gives the signal power at a specific time, frequency and spatial position.

[0297] It should be pointed out that the spatiotemporal power distribution P[k]' at this time still contains the spatial domain index s, so P[k]' at this time should be more appropriately regarded as a time-frequency power distribution in the Fourier correlation concept.

[0298] It is understandable that, since the embodiment needs to calculate the total spatiotemporal power distribution, if the spatial information is to be further and more perfectly obtained, it is necessary to consider how to process the spatial information, that is, to obtain an overall spatiotemporal power distribution P'[k, t] that is not dependent on a specific spatial position, so all spatial positions can be averaged to meet the above requirements:

[0299] ;

[0300] in, is the total number of spatial points.

[0301] It can be understood that this expression calculates the average value of the signal power at all spatial locations for each time point t and frequency index k. In this way, P'[k,t] becomes a true overall spatiotemporal power distribution, which no longer depends on a specific spatial location.

[0302] Preferably, the overall spatiotemporal power distribution P'[k, t] can further provide a two-dimensional matrix, where each element represents the signal power at a given time point and frequency index:

[0303] ;

[0304] It can be understood that for the matrix dimension: the dimension of the matrix is ,in is the total number of frequency indices, is the total number of time points. Each row of the matrix represents a specific frequency and each column represents a specific time point.

[0305] It can be understood that for the matrix elements: Each element P′[k,t] in the matrix represents the signal power at a given time point t and frequency index k. This value is obtained by calculating the average value of the signal power at all spatial locations, as described above.

[0306] Specifically, step S202, consider the propagation characteristics of the signal: This step is essentially an optional step. This step considers that the signal of the insulation node may attenuate with the propagation distance. If it is decided to take this into account when calculating the power spectrum density, step S201 can be modified to further improve the average accuracy of step S201. Among them, the attenuation of the wireless signal can be expressed as a function of distance. This technical solution regards it as an exponential decay model, which means that the signal strength decreases exponentially with the increase of distance. In this way, the overall spatiotemporal power distribution P'[k, t] of step S201 is modified:

[0307] Let the attenuation function A(d) represent the attenuation factor of the signal after propagation distance d, where d is the distance from the signal source to the receiving point. The exponential attenuation model is:

[0308] ;

[0309] Among them, α is the attenuation coefficient; e is the base of the natural logarithm (taken as 2.71828).

[0310] Preferably, the attenuation coefficient α is assigned by the solutions provided in the following three documents:

[0311] (1) Xu Leying, Cheng Linfang, Li Yinan, et al. Influence of the installation method of mechanical insulation joint lead connection of track circuit in station on locomotive signal reception [J]. China Railway, 2014 (4): 52-56;

[0312] (2) Zhao Kai. Design of railway track insulation online test system based on DSP[D]. Xi'an University of Science and Technology, 2018;

[0313] (3) Wang Xiaoming, Shang Wenbin, Shang Wenqing. Research on dynamic contact impedance of insulated wheel-rail in high-speed railway station [J]. Railway Communication Signal, 2017 (8): 10-13;

[0314] Preferably, the attenuation coefficient α can also be assigned by a person skilled in the art using conventional signal propagation analysis techniques and at his / her own initiative.

[0315] It is understandable that in order to apply this decay function to the overall spatiotemporal power distribution P'[k, t], it is necessary to know the distance corresponding to each spatial point s . Then use the decay function (decay factor) To correct the power value of each spatial point.

[0316] It should be pointed out that since the overall spatiotemporal power distribution P'[k, t] has been averaged for all spatial points in step S201, it is necessary to combine this averaging process with the attenuation correction. Therefore, further, instead of directly performing spatial averaging in step S201, the weighted power of each spatial point is calculated first and then averaged. In this way, the attenuation function can be directly applied to each spatial point. The corrected overall spatiotemporal power distribution P''[k, t] is calculated as:

[0317] ;

[0318] The power at each spatial point is The corresponding attenuation factor The weights are then averaged. Therefore, the contribution of farther spatial points to the overall power is reduced, reflecting the characteristic that the signal attenuates with distance.

[0319] In this example, see Figure 1 Regarding step S3, determining a preset threshold T: based on the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P'[k,t], calculating the preset threshold T; there are two schemes:

[0320] (Scheme 1) Single threshold scheme:

[0321] A preset threshold T is set based on the distribution characteristics of the overall spatiotemporal power distribution P''[k, t]. The preset threshold T is set based on the statistical characteristics (mean or standard deviation) in the power spectrum density. These include:

[0322] 1) Preset threshold based on mean: Assume that the background level of the power spectrum density is approximated by its mean. The preset threshold T can be set to the mean of the overall spatiotemporal power distribution P''[k, t] plus an offset, which can be fixed or based on a percentage of the mean:

[0323] T = μ + offset;

[0324] Where μ is the mean of the overall spatiotemporal power distribution P''[k, t], which reflects the background level of the power spectral density and is calculated as:

[0325] ;

[0326] in, is the total number of frequency indices, is the total number of time points.

[0327] Next, set a positive offset; the offset is a fixed percentage of the mean (for example, 10%), or it can be a fixed constant:

[0328] offset=0.1μ;

[0329] 2) Preset threshold based on standard deviation: This method takes into account the volatility of the power spectral density. The preset threshold T can be set to the mean of the overall spatiotemporal power distribution P''[k, t] plus several times the standard deviation. This can capture more significant signal changes relative to the background level:

[0330] T = μ + σn;

[0331] where μ is the mean of the overall spatiotemporal power distribution P''[k, t], as described above; σ is the standard deviation of the overall spatiotemporal power distribution P''[k, t], calculated as:

[0332] ;

[0333] In this formula, n is a positive coefficient, for example, it can be 2 or 3, indicating that the preset threshold T is the mean plus 2 or 3 times the standard deviation.

[0334] (Scheme 2) Dynamic threshold scheme:

[0335] The preset threshold T is dynamically adjusted according to the real-time statistical characteristics of the signal, that is, the mean value within the sliding window.

[0336] Assume a sliding window with a width of W, which contains the power spectral density data of the most recent W time points. For a given current overall spatiotemporal power distribution P''[k, t], calculate the mean within the sliding window , and then use this mean to dynamically set the dynamic preset threshold T[t]:

[0337] S3010, calculating the mean value in the sliding window:

[0338] ;

[0339] in, is the total number of frequency indexes k, and τ is the time point operation traversed within the sliding window. The power spectral density data at all time points and frequency indices in the sliding window are summed and divided by the total number of data points. , get the mean value in the sliding window .

[0340] It can be understood that the mean here reflects the average level of the power spectral density within the sliding window.

[0341] S3011, set dynamic preset threshold: based on the calculated mean value μ_w[t] in the sliding window, set the dynamic preset threshold T[t]. There are two setting methods to choose from:

[0342] Based on the mean within the sliding window , set the dynamic preset threshold T[t] as:

[0343] ;

[0344] or,

[0345] ;

[0346] in, is a positive offset, which is assigned by the staff according to the actual situation. The selection of the offset should ensure that the dynamic preset threshold T[t] is sensitive enough to capture effective signal changes, while being robust enough to avoid too many false alarms.

[0347] in, is a positive percentage value, which is also assigned by the staff according to the actual situation. This setting method takes into account the relative change of the power spectrum density in the sliding window, so that the dynamic preset threshold T[t] can be adaptively adjusted according to the real-time characteristics of the signal.

[0348] Preferably, the offset and percentage values According to the scheme given in the fourth embodiment, the overall steps S1 to S3 of this embodiment can obtain the function of adaptive correction.

[0349] In this example, see Figure 1 , about step S4: suppose the train mileage coordinates at time t are , the train mileage coordinates at time t+1 are , the absolute mileage of the current insulation node obtained through the database is , and then get the mileage correction value C:

[0350] ;

[0351] If the absolute value of C is greater than or equal to the preset threshold T, the positioning error is judged to be too large, and an alarm is issued to notify ground personnel or on-board engineers to intervene. If the absolute value of C is less than the preset threshold T, C is used to correct the current absolute mileage.

[0352] Example 2: Please refer to Figure 6 Based on the first embodiment, this embodiment further provides a more accurate method for obtaining the estimated value P[k]:

[0353] In the railway mileage determination method based on Fourier analysis, in order to more accurately capture the characteristics of the isolation node signal, this technical solution proposes a method of independently calculating and combining the power spectral density (PSD) of the spatial point. This method takes into account the independence of the isolation node signal in the spatial domain and calculates the PSD for each spatial point (isolation point) separately, and finally combines these PSDs to form a comprehensive PSD or keeps them separate for subsequent analysis.

[0354] The implementation method of this embodiment is: each spatial point (isolation point) is considered as an independent signal source, then the power spectrum density (PSD) is calculated independently for each spatial point, and finally these PSDs are combined to form a comprehensive PSD or kept separate for subsequent analysis. This implementation method specifically includes three sub-steps in total, namely, steps S2000 to S2002.

[0355] In this embodiment, regarding step S2000, independent initialization: a specific spatial point s is selected, which corresponds to an isolation point. Since each isolation point is regarded as an independent signal source, the PSD will be calculated separately for each isolation point. To achieve this, it is necessary to fix the spatial point s and traverse the frequency index k and time index t to extract the signals of all frequencies and time points corresponding to the selected spatial point s from the multidimensional data structure X'[k,t,s].

[0356] It is understandable that the purpose of step S2000 is to extract the signal of a specific spatial point from the complex three-dimensional data structure so as to perform subsequent power spectrum density calculation. The advantage of this is that the signal characteristics of each insulation point can be analyzed more accurately without being disturbed by the signals of other spatial points.

[0357] In this embodiment, regarding step S2001, for the signal of the extracted specific spatial point s, the power spectrum density at each frequency index k is calculated using the power spectrum density calculation formula:

[0358] ;

[0359] in, is the power spectral density estimate for spatial point s at frequency index k;

[0360] is the total number of time points used to calculate the PSD (i.e. the length of the signal in time or the length of the time window used for averaging);

[0361] is the square of the signal modulus at a given frequency k (frequency domain index k), time t (window t), and spatial point s (spatial domain index s).

[0362] It should be pointed out that the calculation process of this step is similar to the standard power spectrum density calculation, but the key difference is that the calculation here is performed for a single spatial point. The advantage of this is that the frequency characteristics of the signal at each insulation point can be captured more accurately, providing more precise information for subsequent analysis and railway mileage judgment.

[0363] In this embodiment, regarding step S2002, combining the power spectrum densities of all spatial points, the following three solutions are provided for reference only:

[0364] (Scheme 1) Subjective selection scheme: Based on the subjective initiative of the technician in the field, keep the PSD of each spatial point separate and analyze them independently. Applicable to the case where the spectral characteristics of a specific insulation point are of interest. Then independently formulate an estimate P[k].

[0365] (Scheme 2) Average power spectral density scheme: In order to obtain an overall PSD estimate rather than one specific to a certain insulation point, the average of the power spectral density of all spatial points can be calculated:

[0366] ;

[0367] in, is the total number of points in space, is the average power spectral density at each spatial point.

[0368] This scheme is suitable for situations where the overall spectral characteristics are of interest and can provide a more comprehensive estimate of railway mileage.

[0369] (Scheme 3) Weighted Power Spectral Density Scheme: If some spatial points have different importance or uncertainty, the power spectral density can be combined using weighted averaging, where the contribution of each spatial point is adjusted according to its weight. The final output is a weighted estimated value P[k].

[0370] In this embodiment, the MATLAB execution program of the above steps S2000-S2002 is as follows:

[0371] % The size of the three-dimensional array X_prime is [N_k, N_t, N_s],

[0372] % Where N_k is the number of frequency indices, N_t is the number of time indices, and N_s is the number of spatial points (insulation points).

[0373] % Step S2000, independent initialization: initialize variables

[0374] N_k = size(X_prime, 1); % Frequency index number

[0375] N_t = size(X_prime, 2); % Number of time indexes

[0376] N_s = size(X_prime, 3); % Number of spatial point indices (number of insulation points)

[0377] P_s = zeros(N_k, N_s); % Pre-allocate the power spectral density array for each spatial point

[0378] P = zeros(N_k, 1); % Pre-allocate the average power spectral density array

[0379] % Step S2001: Calculate the power spectrum density of a single spatial point

[0380] for s = 1:N_s

[0381] for k = 1:N_k

[0382] % Extract the signal corresponding to all frequencies and time points of the selected spatial point s and calculate the PSD

[0383] signal_s = X_prime(k, :, s);

[0384] P_s(k, s) = mean(abs(signal_s).^2); % Use mean function to calculate time average

[0385] end

[0386] end

[0387] % Step S2002: Combine the power spectral densities of all spatial points (using the average power spectral density scheme)

[0388] for k = 1:N_k

[0389] % Calculate the average value of the power spectral density of all spatial points

[0390] P(k) = mean(P_s(k, :));

[0391] end

[0392] The PSD calculation in the above code uses time averaging instead of the more standard frequency domain calculation. In actual applications, if the more advanced FFT (Fast Fourier Transform) is used to calculate the power spectral density in the frequency domain, step S2001 in the above program is transformed into:

[0393] % Perform FFT on each spatial point

[0394] for s = 1:N_s

[0395] X_freq(:, s) = fft(X_time(:, s)); % Calculate FFT

[0396] end

[0397] % Calculate power spectral density and average

[0398] for k = 1:N_t

[0399] for s = 1:N_s

[0400] % Calculate the power spectral density of each spatial point

[0401] P_s(k, s) = abs(X_freq(k, s))^2 / N_t; % Normalization

[0402] end

[0403] % Calculate the average value of the power spectral density of all spatial points

[0404] P(k) = mean(P_s(k, :));

[0405] end

[0406] % It is recommended to take half of the frequency points after FFT (due to symmetry) and normalize the frequency axis.

[0407] In the above procedure, power spectral density (PSD) is a measurement that describes the power distribution of a signal or time series in the frequency domain. In the railway mileage determination method, by analyzing the PS of the insulation node signal, information about the frequency characteristics of the signal can be obtained, which may be crucial for determining the railway mileage and other related parameters.

[0408] In step S2000 of the above program: The program does not actually perform explicit initialization operations because MATLAB does not require all variables to be declared in advance like other programming languages. However, conceptually, space is allocated for the PSD values ​​to be calculated (i.e., arrays are preallocated), which is done by creating zero matrices P_s and P in the above code, which will store the PSD and average PSD for each spatial point, respectively. Preallocating arrays is an important step because it can help improve the performance of the code, especially when dealing with large data sets. By allocating enough space for the results in advance, you can avoid dynamically resizing the arrays in the loop.

[0409] In step S2001 of the above program: it iterates through each spatial point (isolation point) and each frequency index through two layers of loops. For each spatial point, it extracts the signal corresponding to that point (in the sample code, this is done by directly indexing the three-dimensional array X_prime), and then calculates the power spectral density of the signal. FFT is performed on the signal at each spatial point, and then the power spectral density of the first N_f frequency points is calculated. This is because for real signals, the FFT output is symmetrical, so only the first half of the frequency points need to be considered.

[0410] In step S2002 of the above procedure: the average value of the power spectral density of all spatial points is calculated to obtain an overall PSD estimate. This is done by traversing each frequency index and averaging the PSD values ​​of all spatial points at that frequency.

[0411] Embodiment 3: Based on step S3 of embodiment 1, this embodiment further provides another scheme for assigning a preset threshold value T, namely, a multi-threshold scheme:

[0412] For different time domain and space domain conditions, multiple preset thresholds T can be set to form a threshold set. In practical applications, a suitable preset threshold T is selected according to the current state of the signal. The steps include:

[0413] P1. Define the preset threshold T as a set T (hereinafter referred to as "set T" in this embodiment): Set T contains N preset thresholds Each preset threshold is associated with specific temporal and spatial conditions.

[0414] P2. Extract signal features: Extract features related to threshold selection from the current overall spatiotemporal power distribution P''[k, t]. These features include power level at a specific frequency or time interval, signal-to-noise ratio, dynamic range of the signal, etc.

[0415] P3. Determine the current state: Based on the extracted features, determine the current state of the signal. This can be achieved by comparing with previously defined state conditions.

[0416] For example, if the power level of a specific frequency exceeds a frequency threshold limit (which needs to be formulated according to the actual railway environment), the signal is considered to be in a specific state.

[0417] P4. Select a preset threshold: According to the current state of the signal, select the most appropriate preset threshold from the threshold set T . This is achieved through conventional table lookups, conditional statements, or logical structures such as decision trees.

[0418] P5. Apply preset threshold: Apply the selected preset threshold Applied in subsequent signal processing or decision-making steps.

[0419] In this embodiment, the above steps are abstracted into the logical structure of MATLAB pseudo code as follows:

[0420] %Define the threshold set T = {T1, T2, ..., TN} and the related condition set Conditions = {Condition1, Condition2, ..., ConditionN}

[0421] %For each time point t and frequency point k:

[0422] Extract features features = ExtractFeatures(P''[k, t])

[0423] % For each condition in Conditions:

[0424] If condition matches features:

[0425] Select the threshold T_i, where i is the index of the condition in Conditions

[0426] Break out of the loop

[0427] % If there is no matching condition:

[0428] Select the default threshold T_default

[0429] % Apply the selected threshold T_i or T_default to subsequent processing

[0430] Embodiment 4: Based on step S3011 of embodiment 1, this embodiment further provides a value assignment scheme. The value assignment scheme of this embodiment enables steps S1 to S3 as a whole to obtain a function of adaptively adjusting according to changes in the railway environment.

[0431] In the first embodiment, offset and are positive offsets, which are assigned by staff based on their own initiative according to the actual situation. The selection of the offset should ensure that the dynamic preset threshold T[t] is sensitive enough to capture effective signal changes, while being robust enough to avoid excessive false alarms. is a positive percentage value, which is also assigned by the staff according to the actual situation. This setting method takes into account the relative change of the power spectrum density in the sliding window, so that the dynamic preset threshold T[t] can be adaptively adjusted according to the real-time characteristics of the signal.

[0432] The above subjective assignment strategy limits the preset threshold T of the final output. There are still certain subjective limitations. In this embodiment, it is hoped that the different background noise levels of each insulation node signal or the signal strength required for detection can be considered as the basis for correction, and the offset or percentage value Assign a value. For example, the technical solution of this embodiment is implemented as follows:

[0433] The DS evidence theory algorithm is introduced (based on the distinction in the first embodiment, this embodiment is called step S5), and the different background noise levels of each isolation node signal and the signal strength required for detection are used as references for confidence, and a joint confidence function Bel(A∪B) is formed using the Dempster's combination principle, and finally the percentage value is formed in the form of sigmoid function mapping. The specific steps include:

[0434] In this embodiment, regarding step S500, a definition set is initialized: all the insulation nodes and their numbers along a section of railway are divided into group A and group B in the form of odd-even grouping, each group has a total of n insulation nodes, and the insulation nodes of group A and group B are interspersed along the railway;

[0435] According to the background noise intensity or signal intensity of each insulation node fed back by the ground simulation recognition system, it is set as a set:

[0436] ;

[0437] ;

[0438] in, It is an insulation joint or Respective background noise intensity or signal strength vector or parameter value.

[0439] The above-mentioned group A and group B are regarded as "evidence set A" and "evidence set B" (hereinafter referred to as evidence A and evidence B) in the DS evidence theory algorithm, and then the subsequent steps are performed.

[0440] In this embodiment, regarding step S501, Dempster's combination rule: for a given A and evidence B, Dempster's combination rule combines the two to form a combined trust distribution set Bel(C), where C is a possible hypothesis or state space. Dempster's combination rule also calculates the combined uncertainty distribution Pl(C):

[0441] ;

[0442] ;

[0443] Among them, ∅ represents the empty set, and are a subset of evidence A and evidence B respectively (corresponding to an isolation node and its background noise intensity or signal intensity vector); C is a "possible set";

[0444] Specifically, W is the weight. If the staff believes that a certain insulation joint is special, the weight W value of the insulation joint can be adjusted subjectively; if it is believed that all insulation joint signals have the same status, the weight W can be assigned a value of 0.5.

[0445] In this embodiment, regarding step S502, the joint trust function Bel(A∪B) is calculated:

[0446] ;

[0447] Where Bel(A) and Pl(A) are the trust and uncertainty distribution of the set of evidence A;

[0448] Where Bel(B) and Pl(B) are the trust and uncertainty distribution of the set of evidence B;

[0449] in, and are subsets of evidence A and evidence B respectively. The numerator is the sum of the mass of evidence supporting A∩B, and 1 minus this sum gives Pl(A∩B), which indicates the degree of uncertainty about A∩B.

[0450] Among them, (Pl(A)∩Pl(B)) is the uncertainty distribution of the intersection of evidence A and evidence B:

[0451] ;

[0452] Specifically, among them, Represents any subset of evidence A and a subset of evidence B The intersection of is not empty. That is, it seeks a subset that exists in the objective , the subset W in A is not zero, and there is a subset , the W of this subset in B is also not zero, that is, the two subsets have an intersection.

[0453] It is understandable that Responsible for finding some elements in hypothesis set C that are supported by both evidence A and evidence B. , then when calculating the combined trust distribution Bel(C), these common elements must be taken into account because they contribute to the support of the hypothesis set C.

[0454] Furthermore, after the above concept is concretized, it is determined whether the "influence" of the background noise intensity or signal strength of a certain insulating node itself is "supported" by the "influence" of the background noise intensity or signal strength of another insulating node adjacent to it; if so, it is considered that evidence A and evidence B support each other.

[0455] This means that the isolation node will affect the offset due to its background noise intensity or signal strength. or percentage value If this has an adverse effect, the interval value of the subsequent mapping will be biased, which will lead to an offset or percentage value The assignment is more inclined to the effect of "causing a deviation between the preset threshold value T and the mileage correction value C".

[0456] In this embodiment, regarding step S503, the joint trust function Bel(A∪B) is mapped to an interval value between [0,1] as a percentage value through a sigmoid function. :

[0457] ;

[0458] e is the base of natural logarithms;

[0459] It is understandable that the percentage value of the output The size of will be affected by Bel(A∪B), thereby adjusting the state transition. The output value of the sigmoid function approaches 1 when Bel(A∪B) is large, and approaches 0 when Bel(A∪B) is small. Therefore, the percentage value of this embodiment is The size reflects the size of Bel(A∪B). This mapping preserves the causal relationship, that is, the percentage value The change of is consistent with the change of Bel(A∪B). Through this mapping relationship, this embodiment can introduce a correction for uncertainty when considering the credibility of the background noise intensity or signal strength of each isolation node, so as to adapt to different situations more flexibly. This design helps the system to be more adaptive and robust when facing new data and changes, so as to achieve the purpose of adaptive adjustment.

[0460] It can be understood that for steps S1 to S3 of Example 1, the solution provided in this embodiment is that when the background noise intensity or signal strength of the insulating joint in the railway environment changes, the overall system can make corresponding operation modes similar to weight assignment according to the changes, thereby realizing the adaptive adjustment of the overall method according to the changes in the railway environment.

[0461] For demonstration, combined with the content of step S3 of the first embodiment, if the mean is 100, and the percentage value output by the above step S503 is ,but:

[0462] ;

[0463] That is, the preset threshold T is regarded as 110.

[0464] In summary, for this embodiment, the advantage of the adaptive correction scheme is that it can use the different background noise levels of each insulation node signal or the signal strength required for detection as the basis for correction, and automatically adjust the threshold according to the dynamic range of the signal, making the threshold setting more flexible and adaptable.

[0465] Example 5: Please refer to Figure 8 , this embodiment provides a railway mileage measurement system based on Fourier analysis:

[0466] This system is a railway mileage measurement solution that integrates advanced signal processing technology and geographic positioning technology. It is designed to accurately measure mileage by precisely analyzing the electrical signals of railway insulation joints.

[0467] 1. System composition:

[0468] 1) Insulation node layout: Multiple insulation nodes are strategically arranged along the railway line. These insulation nodes will generate specific electrical signals when a train passes through, and these signals contain the time domain and space domain information of the train passing through.

[0469] 2) Ground simulation identification system: This system is equipped with a ground simulation identification system, which is connected to the insulation joint through wireless communication. The ground simulation identification system is responsible for receiving the electrical signal generated by the insulation joint and transmitting it to the processor for further analysis.

[0470] Preferably, the ground simulation identification system is the ground simulation identification system described in the document (Yang Ji, Xu Qingyang, Duan Hehui, et al. Research on dynamic automatic detection system of railway signals [J]. Railway Technology Innovation, 2021 (6): 98-103).

[0471] 3) Processor and Registers: The ground simulation identification system is connected to a high-performance processor equipped with registers dedicated to signal processing. The processor is responsible for executing program instructions stored in the registers, which implement the railway mileage determination method based on Fourier analysis.

[0472] (II) Signal processing flow:

[0473] 1) Time domain signal analysis: The processor first focuses on the voltage parameter changes generated when the train passes through the insulation joint, that is, the time domain signal. The change pattern of these signals is directly related to the train passing event, providing basic data for mileage determination.

[0474] 2) Spatial signal analysis: In addition to the time domain signal, the processor also receives the geographic coordinate location information of the insulation node, that is, the spatial signal. These coordinate points accurately identify the location of different insulation nodes on the railway line and provide a spatial reference for mileage calculation.

[0475] 3) Fourier analysis application: The processor uses Fourier analysis technology to perform spectrum analysis on the received time domain signal. By decomposing the different frequency components in the signal, the system can extract the features most relevant to the train passing event, thereby improving the accuracy of mileage measurement.

[0476] 4) Mileage calculation and calibration: Combining the time domain signal characteristics and spatial coordinate information, the processor can calculate the precise mileage of the train when it passes through each insulation section. The system can also further improve the accuracy and robustness of mileage determination by fusing data with other positioning systems.

[0477] All the above embodiments only express the implementation methods of the relevant practical applications of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be based on the attached claims.

[0478] For those skilled in the art, it can be further appreciated that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, computer software, or a combination of the two. In order to clearly illustrate the interchangeability of hardware and software, the composition and steps of each example have been generally described in the above description according to function. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of the present invention.

[0479] At the same time, those skilled in the art can understand that all or part of the processes in all the above-mentioned embodiments can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, storage, database or other media provided in this application and used in the embodiments may include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM) or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. As an illustration and not limitation, RAM is available in many forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double-speed data rate SDRAM (SSRSDRAM), enhanced SDRAM (ESDRAM), synchronous link (Synchlink) DRAM (SLDRAM), memory bus (Rambus) direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM).

Claims

1. A method for measuring railway mileage based on Fourier analysis, comprising obtaining the time domain signal and the space domain signal generated by a train passing through an insulating section by a ground simulation recognition system, characterized in that: The determination method further comprises the following steps: S1. Fourier analysis: for the time domain signal, convert it to the frequency domain and obtain the spectrum data X[k]; for the spatial domain signal, associate it with the spectrum data X[k] to obtain the data structure X'[k, t, s]; In S1, the converting, associating and obtaining steps are S100 to S102 respectively: S100, Fourier transform: for the discrete time domain signal x[n], obtain the spectrum data X[k]: ； Wherein, N is the number of samples of the time domain signal; i is an imaginary unit; e is the base of the natural logarithm; n is the discrete time index of the time domain signal; S101, time-frequency analysis: performing a short-time Fourier transform operation to correct the spectrum data X[k] to new spectrum data X[k,t]; S102, associating spatial information: when executing S100 and S101, obtaining the spatial information S[t] of the spectrum data X[k] and the new spectrum data X[k, t], and obtaining the data structure X'[k, t, s] by data integration; in S102, the step of data integration includes: S1020, for each time point, record the new spectrum data X[k, t]; S1021, associating the spatial domain information S[t] corresponding to the time point with the new spectrum data X[k, t] and then creating the data structure X'[k, t, s], wherein the data structure X'[k, t, s] includes the original spectrum data and the corresponding spatial domain information; Among them, k is the frequency domain index, t is the starting position of the window, and s is the spatial domain index; S2. Analyze power spectrum density: perform a PSD algorithm on the data structure X'[k, t, s] to calculate the spatiotemporal power spectrum density to form a spatiotemporal power distribution P[k]' or an overall spatiotemporal power distribution P'[k, t]; The PSD algorithm performs the following steps: S200, calculating the estimated value P[k] of the power spectrum density: for the discrete data structure X'[k, t, s], the estimated value of the power spectrum density is: ； Where P[k] represents the estimated value of the power spectral density at frequency index k, is the square of the modulus of the frequency domain signal; S201, calculating the spatiotemporal power spectral density: for the data structure X'[k,t,s] under the estimated value P[k], calculating the power distribution of the time domain signal and the spatial domain signal in time and space to form the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P'[k,t]; The method of calculating the power distribution of the spatiotemporal power distribution P[k]' is: ； Wherein, the spatiotemporal power distribution P[k]' is the time-frequency power distribution at frequency index k and window t; The implementation method of the power distribution calculated by the overall spatiotemporal power distribution P'[k,t] is: ； in, is the total number of points in space; S3, determining a preset threshold value T: calculating the statistical characteristics of the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P'[k,t] to obtain the preset threshold value T; In S3, the preset threshold T is obtained by the statistical characteristics of the mean or standard deviation of the spatiotemporal power distribution P[k]' or the overall spatiotemporal power distribution P'[k,t], or the mean value in the sliding window. , set to the dynamic preset threshold T[t]; S4. Determination: After calculating the train mileage correction value C, based on the preset threshold value T, determine whether the absolute mileage of the train needs to be corrected.

2. The method for measuring railway mileage based on Fourier analysis according to claim 1, characterized in that: In the S101, the short-time Fourier transform operation includes: S1010. Define windowing and window function w[n]: ； Wherein, x[n] is the discrete time domain signal, n is the discrete time index in the window, and t is the starting position of the window; S1011, executing DFT algorithm: ； Where k is the index in the frequency domain; S1012, iteration: move the sliding window to the next position, repeat S1010 and S1011, until the entire signal is analyzed, and the spectrum data X[k, t] is updated.

3. The method for measuring railway mileage based on Fourier analysis according to claim 1, characterized in that: In S4, the mileage correction value C is obtained by: Assume that the train mileage coordinates at time t are , the train mileage coordinates at time t+1 are , the absolute mileage of the current insulation node is obtained through the ground simulation identification system: , and then obtain the mileage correction value C: ； In S4, the judging method is: If the absolute value of the mileage correction value C is greater than or equal to the preset threshold value T, it is determined that there is an error in the train positioning; if the mileage correction value C is less than the preset threshold value T, the mileage correction value C is used to correct the previous absolute mileage of the current train.

4. A railway mileage measurement system based on Fourier analysis, characterized in that: The system includes insulation joints arranged along the railway and arranged at the back-end control end. A ground simulation identification system, the ground simulation identification system is connected to the insulation node via wireless communication; The ground simulation identification system is connected to a processor, and the processor is connected to a register; The processor is used to analyze the electrical signal of the insulation node transmitted by the ground simulation identification system, wherein the electrical signal includes a time domain signal and a space domain signal; The register stores program instructions, and when the program instructions are executed by the processor, the processor executes the railway mileage measurement method based on Fourier analysis as described in any one of claims 1 to 3.

Images