A method for evaluating the running stability of straight-moving vehicles at intersections

Abstract

The present invention discloses a method for evaluating the running stability of through vehicles at an intersection, which is used to evaluate the running safety of through vehicles affected by left-turning vehicles in a permitted phase at an intersection. The method comprises the following parts: parameter initialization, obtaining the time series of the headway of through traffic, phase space reconstruction, Lyapunov stability analysis, outlier removal and outlier rate calculation. The only data required by the evaluation method is the headway sequence of the through traffic. The algorithm framework is automatically implemented by a written program, and the running safety of through traffic at four intersections in Changchun is analyzed. The contribution of the invention is to provide a quantitative and easy-to-implement method to evaluate the running safety of through traffic at an intersection. In practice, it can be applied to the effect evaluation before and after the improvement of signal control schemes or channelization schemes.

Description

Technical Field

[0001] The invention belongs to the field of traffic safety and specifically relates to a method for evaluating the running stability of straight-moving vehicles at an intersection. Background Art

[0002] At a signalized intersection in the left-turn permission phase, according to laws and regulations, left-turning vehicles must give way to straight-moving vehicles. This means that the passage of straight-moving vehicles should not be affected by left-turning vehicles in the permission phase. However, this assumption is too ideal and cannot reflect the actual traffic flow operation status. In fact, some drivers in my country do not have much awareness of the right of way. They prefer to obtain priority in the permission phase through vehicle game and make decisions based on actual traffic conditions. This behavior is called non-strict priority phenomenon.

[0003] Some scholars have paid attention to the phenomenon of non-strict priority. Wang Wei found that conflicting traffic at unsignalized intersections took turns to pass through the intersection during peak hours. He abandoned the gap acceptance analysis method and proposed a fleet-based single-lane intersection capacity and delay model. After that, Meng Yongping and Li Aizeng did a lot of work to calculate the capacity of two-lane intersections. They assumed that left-turning traffic was equivalent to straight traffic with a large critical gap. In addition, Kaysi and Abbany developed a binary probit-based behavioral model to predict the probability of drivers engaging in non-strict priority behaviors, such as competing to obtain priority at an unsignalized intersection.

[0004] For intersections with a left turn permission phase, previous studies mostly discussed left-turn behavior under non-strict priority. Qu Zhaowei's research found that compared with vehicles in the left turn protection phase, vehicles in the left turn permission phase always have greater acceleration to make themselves reach potential conflict points earlier than oncoming straight vehicles. Because drivers pay more attention to the accelerator pedal, the time for braking reaction increases. In order to alleviate the competition for priority, Bai Qiaowen tried to optimize the design form of the left turn guide line. He set three key points on the left turn guide line. His method is to constrain the path of left-turn vehicles in the permission phase through the left turn guide line and force left-turn vehicles to arrive at potential conflict points later than oncoming straight vehicles. These research results show that under non-strict priority conditions, the operating safety of straight vehicles in the left turn permission phase is seriously affected. However, few studies have conducted quantitative analysis on the operating safety of straight vehicles under the influence of left turn permission vehicles. Therefore, we make up for this deficiency through this invention research. Summary of the invention

[0005] The purpose of the present invention is to overcome the deficiencies in the above-mentioned prior art and to provide a method for evaluating the running stability of straight-moving vehicles at an intersection. The evaluation method is used to evaluate the running stability and safety of straight-moving vehicles affected by left-turning vehicles in the permitted phase at the intersection.

[0006] To achieve the above object, the technical solution of the present invention is: a method for evaluating the running stability of straight-moving vehicles at an intersection, characterized in that it comprises the following steps:

[0007] S1: parameter initialization;

[0008] S2: Collect the headway time parameters of the through traffic flow and obtain the headway time sequence;

[0009] S3: Reconstruct the phase space according to the acquired headway sequence;

[0010] S4: Perform Lyapunov exponent stability analysis;

[0011] S5: Remove outliers and calculate the outlier rate.

[0012] Furthermore, the parameter initialization in the above step S1 specifically includes the following steps:

[0013] S11: Define the outlier rate γ, which is expressed as is the number of outliers N removed from the original headway sequence r The ratio of the number of samples N0 of the original headway sequence;

[0014] S12: Initialize the parameters, set the number of iterations ρ = 0, the number of outliers N r =0.

[0015] Furthermore, the specific step S21 of obtaining the headway sequence in step S2 is: collecting the headway parameters of the straight traffic flow, using Record the original headway time series collected during the field survey.

[0016] Furthermore, the specific steps of the phase space reconstruction in step S3 above include:

[0017] S31: The through traffic flow is regarded as a dynamic system, and all headway data form a one-dimensional time series;

[0018] S32: Determine the 4 key parameters of phase space reconstruction, namely the hysteresis time Embedding window τ w , embedding dimension m, average period p;

[0019] S33: The CC method is used to calculate the hysteresis time and the embedded window value, and the hysteresis time and the embedded window value are estimated simultaneously through the correlation integral. The specific calculation process is as follows:

[0020] Assume that the headway sequence is h = {h i |i=1,2,...N}, N is the sample size, a new phase space H={H i |H i =[h i ,h i+t ,h i+(m-1)t ]} will be determined by the hysteresis time The embedding dimension m is constructed with two variables, and the correlation integral formula of the embedded time series is as follows:

[0021]

[0022] Where: M is the number of embedded points in the m-dimensional space, calculated by M = N-(m-1)t, d ij Yes i -h j The supernormal state, θ(rd ij ) is determined by the following equation, which gives a pair of formulas for embedding points in the time sequence h,

[0023]

[0024] The following formula determines The test statistic formula ΔS(m,t):

[0025] ΔS(m,t)=max{S(m,r,t)-minS(m,r,t)}

[0026] The following formula gives the calculation formula of S(m,r,t) when N→∞:

[0027]

[0028] Optimal delay time is the first minimum value of the function ΔS(m,t)~t,

[0029] The following formula is given to determine τ w The test statistic formula S cor (t):

[0030]

[0031] in

[0032] Embedding window τ w is the function S c o r(t)~t minimum point,

[0033] Based on the above, m can be calculated by the following formula:

[0034] τ w =(m-1)·t

[0035] To calculate the average period P, we must first obtain the amplitude of the headway sequence A1, A2, ...An through Fourier transform. Then, the average period P can be calculated by the following formula:

[0036]

[0037] Where: f i Indicates frequency, A i Indicates amplitude.

[0038] Furthermore, the specific steps of the Lyapunov exponent stability analysis in the above step S4 include:

[0039] S41: Determine the stability of headway sequences by calculating the Lyapunov exponent, which quantifies the stability of headway sequences with hysteresis time. and the divergence of adjacent initial trajectories in the phase space of embedding dimension m, which is expressed as:

[0040]

[0041] S42: The Lyapunov exponent is calculated by the Wolf algorithm. The Wolf algorithm gives the maximum Lyapunov exponent when reaching the end of the trajectory. The basic function is given by the following formula:

[0042]

[0043] Where: L m represents the maximum Lyapunov exponent, represents the system evolution time, is the total number of iterations, ξ(τ0) is the Euclidean distance between the reference point and the adjacent point, ξ'(τ κ ) is the Euclidean distance after ξ(τ0) evolves with a step size of κ, which is expressed by the maximum Lyapunov exponent L m The value of L determines the stability of the time series. m When L is less than 0, the system will converge to a fixed point and can be identified as a stable system. On the contrary, when L m When it is greater than 0, it indicates that the system is unstable. In addition, when the Lyapunov exponent is 0, the system will be in a critical state.

[0044] S43: Determine the termination condition of the loop for calculating the Lyapunov exponent. Each time when an outlier is deleted to obtain a new time interval sequence, the maximum Lyapunov exponent is continuously calculated. Once the maximum value of the Lyapunov exponent is negative, the loop is terminated.

[0045] Furthermore, the specific steps of deleting outliers in step S5 and calculating outlier rates include:

[0046] S51: From the time series Get its maximum value h (ρ) max , minimum value h (ρ) min 、average value h (ρ) mean ;

[0047] S52: Compare the maximum value h (ρ) max , minimum value h (ρ) min With the average h (ρ) mean The largest difference is considered as an outlier, and the outlier sample is deleted from the time interval sequence, and the number of outlier samples N deleted in iteration ρ is output r (ρ) , on this basis, a new time interval sequence is established and the Lyapunov stability analysis is performed again;

[0048] S53: Calculate the outlier rate γ,

[0049] Compared with the prior art, the technical solution of the present invention has the following technical effects:

[0050] 1. Among the existing technologies for evaluating the operational safety of intersections under the permitted phase, most focus on left-turning vehicles under the permitted phase. The present invention focuses on the impact on straight-moving vehicles and can quantitatively analyze the operational safety of intersections through the headway data of straight-moving vehicles.

[0051] 2. The data acquisition method of the present invention is simple, economical and efficient, and can well solve the defects of the current difficulty in acquiring traffic data and data accuracy. By simply collecting the headway data of straight-moving vehicles and obtaining a one-dimensional headway sequence, the intersection operation stability evaluation can be performed, avoiding the complexity of acquiring traffic data, such as vehicle trajectory, vehicle speed, vehicle conflict probability and other data. At the same time, the headway data acquisition accuracy is high, which can be accurate to 0.01s, ensuring the basic data accuracy of this study.

[0052] 3. The present invention verifies the proposed outlier rate algorithm framework through an example, and the verification shows that the present invention provides a quantitative, operational, and effective intersection safety evaluation indicator. In particular, when evaluating the effect of the signal control scheme or channelization scheme before and after improvement, the outlier rate can provide data support for further evaluation. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 is a flow chart of the method of the present invention;

[0054] Figure 2 It is the technical roadmap of the present invention;

[0055] Figure 3 The survey cross-sections of the two research sites, (a) is the survey cross-section of the intersection of Jiefang Road and Tongzhi Street, and (b) is the survey cross-section of the intersection of Haoyue Road and Heping Street;

[0056] Figure 4 This is the phase reconstruction image of the intersection of Haoyue Road and Heping Street;

[0057] Figure 5 Phase reconstruction diagram for the intersection of Jiefang Road and Tongzhi Street;

[0058] Figure 6 The maximum Lyapunov index changes with the number of outliers deleted, where (a) is the maximum Lyapunov index of the intersection of Haoyue Road and Heping Street, and (b) is the maximum Lyapunov index of the intersection of Jiefang Road and Tongzhi Street. DETAILED DESCRIPTION

[0059] See attached Figure 1-2 The present invention provides a method for evaluating the running stability of straight-moving vehicles at an intersection, which is characterized by comprising the following steps:

[0060] S1: parameter initialization;

[0061] S2: Collect the headway time parameters of the through traffic flow and obtain the headway time sequence;

[0062] S3: Reconstruct the phase space according to the acquired headway sequence;

[0063] S4: Perform Lyapunov exponent stability analysis;

[0064] S5: Remove outliers and calculate the outlier rate.

[0065] Furthermore, the parameter initialization in the above step S1 specifically includes the following steps:

[0066] S11: Define the outlier rate γ, which is expressed as The outlier rate is introduced to evaluate the operational stability of the through-traffic flow in the left-turn permission phase. It is the number of outliers N removed from the original headway sequence. r The ratio of the number of samples N0 of the original headway sequence;

[0067] S12: Initialize the parameters, set the number of iterations ρ = 0, the number of outliers N r =0.

[0068] Furthermore, the specific step S21 of obtaining the headway sequence in step S2 is: collecting the headway parameters of the straight traffic flow, using Record the original headway time series collected during the field survey.

[0069] Furthermore, the specific steps of the phase space reconstruction in step S3 above include:

[0070] S31: The through traffic flow is regarded as a dynamic system, and all headway data form a one-dimensional time series;

[0071] S32: Determine the 4 key parameters of phase space reconstruction, namely the hysteresis time Embedding window τ w , embedding dimension m, average period p;

[0072] S33: The CC method is used to calculate the hysteresis time and the embedded window value, and the hysteresis time and the embedded window value are estimated simultaneously through the correlation integral. The specific calculation process is as follows:

[0073] Assume that the headway sequence is h = {h i |i=1,2,...N}, N is the sample size, a new phase space H={H i |H i =[h i ,h i+t ,h i+(m-1)t ]} will be determined by the hysteresis time The embedding dimension m is constructed with two variables, and the correlation integral formula of the embedded time series is as follows:

[0074]

[0075] Where: M is the number of embedded points in the m-dimensional space, calculated by M = N-(m-1)t, d ij Yes i -h j The supernormal state, θ(rd ij ) is determined by the following equation, which gives a pair of formulas for embedding points in the time sequence h,

[0076]

[0077] The following formula determines The test statistic formula ΔS(m,t):

[0078] ΔS(m,t)=max{S(m,r,t)-minS(m,r,t)}

[0079] The following formula gives the calculation formula of S(m,r,t) when N→∞:

[0080]

[0081] Optimal delay time is the first minimum value of the function ΔS(m,t)~t,

[0082] The following formula is given to determine τ w The test statistic formula S cor (t):

[0083]

[0084] in

[0085] Embedding window τ w is the function S c o r (t)~t minimum point,

[0086] Based on the above, m can be calculated by the following formula:

[0087] τ w =(m-1)·t

[0088] To calculate the average period P, we must first obtain the amplitude of the headway sequence A1, A2, ...An through Fourier transform. Then, the average period P can be calculated by the following formula:

[0089]

[0090] Where: f i Indicates frequency, A i Indicates amplitude.

[0091] Furthermore, the specific steps of the Lyapunov exponent stability analysis in the above step S4 include:

[0092] S41: Determine the stability of headway sequences by calculating the Lyapunov exponent, which quantifies the stability of headway sequences with hysteresis time. and the divergence of adjacent initial trajectories in the phase space of embedding dimension m, which is expressed as:

[0093]

[0094] S42: The Lyapunov exponent is calculated by the Wolf algorithm. The Wolf algorithm gives the maximum Lyapunov exponent when reaching the end of the trajectory. The basic function is given by the following formula:

[0095]

[0096] Where: L m represents the maximum Lyapunov exponent, represents the system evolution time, is the total number of iterations, ξ(τ0) is the Euclidean distance between the reference point and the adjacent point, ξ'(τ κ ) is the Euclidean distance after ξ(τ0) evolves with a step size of κ, which is expressed by the maximum Lyapunov exponent L m The value of L determines the stability of the time series. m When L is less than 0, the system will converge to a fixed point and can be identified as a stable system. m When it is greater than 0, the system will not stabilize at a fixed point, nor will it have a periodic solution, which indicates that the system is unstable but chaotic. In addition, when the Lyapunov exponent is 0, the system will be in a critical state.

[0097] S43: Determine the termination condition of the loop for calculating the Lyapunov exponent. Each time when an outlier is deleted to obtain a new time interval sequence, the maximum Lyapunov exponent is continuously calculated. Once the maximum value of the Lyapunov exponent is negative, the loop is terminated.

[0098] Furthermore, the specific steps of deleting outliers in step S5 and calculating outlier rates include:

[0099] S51: From the time series Get its maximum value h (ρ) max , minimum value h (ρ) min 、average value h (ρ) mean ;

[0100] S52: Compare the maximum value h (ρ) max , minimum value h (ρ) min With the average h (ρ) mean The largest difference is considered as an outlier, and the outlier sample is deleted from the time interval sequence, and the number of outlier samples N deleted in iteration ρ is output r (ρ), on this basis, a new time interval sequence is established and the Lyapunov stability analysis is performed again;

[0101] S53: Calculate the outlier rate γ,

[0102] The present invention is further described below in conjunction with specific embodiments and drawings.

[0103] Taking an intersection in a certain city as an example, a method for evaluating the running stability of straight-moving vehicles at an intersection according to the present invention comprises the following steps:

[0104] 1. First, complete the basic data acquisition work. Define the key variables, obtain the headway data of the intersection traffic flow, and initialize the parameters in the iterative process.

[0105] First, a field survey was conducted at four signalized intersections in Changchun City between 5pm and 6pm. The basic information of these locations is listed in Table 1. Among them, two intersections with similar geometric features, Haoyue Avenue-Heping Avenue and Jiefang Avenue-Tongzhi Street, were selected to capture significant differences, such as Figure 3 The proposed calculation model is further verified using two intersections, Jianshe Street-Beian Road and Tongzhi Street-Tongguang Road.

[0106] Table 1 Basic information of the study sites

[0107]

[0108] 2. The algorithm framework was automatically implemented through MATLAB programming. First, the phase space of the two intersections of Haoyue Avenue-Heping Avenue and Jiefang Avenue-Tongzhi Street was reconstructed.

[0109] Figure 4 and Figure 5 This is the phase space reconstruction process of the original data (outliers not deleted) at the intersection of Haoyue Avenue-Heping Street and Jiefang Avenue-Tongzhi Street. The first minimum value of ΔS(m,t) is when the t value is 3. When the t value is 141, S cor (t) reaches its minimum value. From this, we can get τ w =141, m=48. Similarly, for the Jiefang Avenue-Tongzhi Street intersection with a protection phase, we can get τ w =94, m=24.5.

[0110] 3. Lyapunov exponent stability analysis, outlier removal and outlier rate calculation were performed on the two intersections of Haoyue Avenue-Heping Avenue and Jiefang Avenue-Tongzhi Street, and further verification was carried out on the two intersections of Jianshe Street-Bei'an Road and Tongzhi Street-Tongguang Road.

[0111] Then, we calculated the maximum Lyapunov exponent of the original time sequence in the phase space reconstruction. The calculation shows that the L of the two intersections m The values ​​are all positive, with a permitted phase of 0.5403 and a protected phase of 0.5477. They all fail the Lyapunov stability analysis. At this point, the outlier removal procedure is activated. Figure 6 The change of the maximum Lyapunov exponent when we delete outliers is given. It can be found that when 39 samples are deleted in the permitted phase and 10 samples are deleted in the protected phase, the maximum Lyapunov exponent is negative. Finally, according to the outlier rate calculation formula, the outlier rate of the intersection of Haoyue Avenue and Heping Street in the permitted phase is 13.98%, and the outlier rate of the intersection of Jiefang Avenue and Tongzhi Street in the protected phase is 5.38%.

[0112] Although the geometric features of the two intersections are similar, the safety of through traffic under the two control modes is significantly different. The outlier rate of through traffic with the protection phase is much smaller than that of through traffic under the permission phase. Through traffic under the protection phase has higher safety. It should be noted that some sample values ​​under the protection phase still need to be deleted. These outliers mainly come from vehicles leaving during the green light-on phase. This is also consistent with the description in the Highway Capacity Manual, that is, due to the driver's reaction time at the beginning of the green light, the first few vehicles in the queue will leave with a large headway. The situation is completely different under the permission phase. More vehicle samples that do not leave during the green light-on phase are deleted because the through traffic is always forced to give way to the left turn behavior under the permission phase. This may seriously increase the headway value. In addition, some minimum values ​​are deleted as outliers, which come from the last few through vehicles in the queue, which like to move forward with a small headway to avoid being interrupted by the left-turning vehicles in the opposite direction.

[0113] In addition to the above two intersections, two intersections, Jianshe Street-Beian Road and Tongzhi Street-Tongguang Road, were selected for further analysis. The results are shown in Table 3. The outlier rate of Jianshe Street-Beian Road intersection in the permitted phase is much higher than that of other intersections. Through observation, it is found that there is a large proportion of left-turning vehicles at this intersection. In addition, the green light phase duration of this intersection is very short. It is difficult for straight traffic to achieve a stable operating state. The outliers deleted at the Tongzhi Street-Tongguang Road intersection also mainly come from vehicles leaving at the beginning of the green light. However, compared with the Jiefang Road-Tongzhi Street intersection, the outlier rate is slightly larger. This is because its signal cycle is shorter and more vehicles are released during the green light on stage.

[0114] Table 2 Results of motion stability at four intersections

[0115]

[0116] In a stable headway sequence, each element should not be far from the mean. A large dispersion of the headway sequence is a necessary condition for traffic instability. Therefore, we further performed descriptive statistics on the original data of the four intersections. Table 3 gives the mean and standard deviation of the headway of through traffic at the four intersections. It can be seen that the headway of through vehicles in the permitted phase is more widely distributed than that of through vehicles in the protected phase. In fact, through traffic in the permitted phase is more likely to operate in an unstable state. The results reflected by the descriptive statistics are consistent with the results of the proposed method.

[0117] Table 3 Descriptive statistics of the headway distribution of through traffic

[0118]

[0119] It should be pointed out that the description of the above embodiments is only used to help understand the method of the present application and its core idea. For ordinary technicians in this technical field, several improvements and modifications can be made to the present application without departing from the principles of the present application. These improvements and modifications are also within the scope of protection of the claims of the present application.

Claims

1. A method for evaluating the running stability of straight-moving vehicles at an intersection, characterized in that: The steps include: S1: parameter initialization; S2: Collect the headway time parameters of the through traffic flow and obtain the headway time sequence; S3: Reconstruct the phase space according to the acquired headway sequence; S4: Perform Lyapunov exponent stability analysis; S5: Delete outliers and calculate the outlier rate; The specific steps of the phase space reconstruction in step S3 above include: S31: The through traffic flow is regarded as a dynamic system, and all headway data form a one-dimensional time series; S32: Determine the 4 key parameters of phase space reconstruction, namely the hysteresis time Embedding window τ w , embedding dimension m, average period p; S33: The CC method is used to calculate the hysteresis time and the embedded window value, and the hysteresis time and the embedded window value are estimated simultaneously through the correlation integral. The specific calculation process is as follows: Assume that the headway sequence is h = {h i |i=1,2,...N}, N is the sample size, a new phase space H={H i |H i =[h i ,h i+t ,h i+(m-1)t ]} will be determined by the hysteresis time The embedding dimension m is constructed with two variables, and the correlation integral formula of the embedded time series is as follows: Where: M is the number of embedded points in the m-dimensional space, calculated by M = N-(m-1)t, d ij Yes i -h j The supernormal state, θ(rd ij ) is determined by the following equation, which gives a pair of formulas for embedding points in the time sequence h, The following formula determines The test statistic formula ΔS(m,t): ΔS(m,t)=max{S(m,r,t)-minS(m,r,t)} The following formula gives the calculation formula of S(m,r,t) when N→∞: Optimal delay time is the first minimum value of the function ΔS(m,t)~t, The following formula is given to determine τ w The test statistic formula S cor (t): in Embedding window τ w is the function S cor (t)~t minimum point, Based on the above, m can be calculated by the following formula: t w =(m-1)·t To calculate the average period P, we first need to obtain the amplitude of the headway sequence A1, A2, ...An through Fourier transform; then, the average period P can be calculated by the following formula: Where: f i Indicates frequency, A i represents amplitude; The specific steps of the Lyapunov exponent stability analysis in the above step S4 include: S41: Determine the stability of headway sequences by calculating the Lyapunov exponent, which quantifies the stability of headway sequences with hysteresis time. and the divergence of adjacent initial trajectories in the phase space of embedding dimension m, which is expressed as: S42: The Lyapunov exponent is calculated by the Wolf algorithm. The Wolf algorithm gives the maximum Lyapunov exponent when reaching the end of the trajectory. The basic function is given by the following formula: Where: L m represents the maximum Lyapunov exponent, represents the system evolution time, is the total number of iterations, ξ(τ0) is the Euclidean distance between the reference point and the adjacent point, ξ'(τ κ ) is the Euclidean distance after ξ(τ0) evolves with a step size of κ, which is expressed by the maximum Lyapunov exponent L m The value of L determines the stability of the time series. m When L is less than 0, the system will converge to a fixed point and can be identified as a stable system. On the contrary, when L m When it is greater than 0, it indicates that the system is unstable. In addition, when the Lyapunov exponent is 0, the system will be in a critical state; S43: Determine the termination condition of the loop for calculating the Lyapunov exponent. Each time when an outlier is deleted to obtain a new time interval sequence, the maximum Lyapunov exponent is continuously calculated. Once the maximum value of the Lyapunov exponent is negative, the loop is terminated.

2. A method for evaluating the running stability of straight-moving vehicles at an intersection according to claim 1, characterized in that: The specific steps of initializing the parameters in step S1 above include: S11: Define the outlier rate γ, which is expressed as is the number of outliers N removed from the original headway sequence r The ratio of the number of samples N0 of the original headway sequence; S12: Initialize the parameters, set the number of iterations ρ = 0, the number of outliers N r =0.

3. The method for evaluating the running stability of straight-moving vehicles at an intersection according to claim 1, characterized in that: The specific step S21 of obtaining the headway sequence in step S2 is: collecting the headway parameters of the straight traffic flow, using Record the original headway time series collected during the field survey.

4. The method for evaluating the running stability of straight-moving vehicles at an intersection according to claim 1, characterized in that: The above step S5 deletes outliers and the specific steps of calculating the outlier rate include: S51: From the time series Get its maximum value h (ρ) max , minimum value h (ρ) min 、average value h (ρ) mean ; S52: Compare the maximum value h (ρ) max , minimum value h (ρ) min With the average h (ρ) mean The largest difference is considered as an outlier, and the outlier sample is deleted from the time interval sequence, and the number of outlier samples N deleted in iteration ρ is output r (ρ) , on this basis, a new time interval sequence is established and the Lyapunov stability analysis is performed again; S53: Calculate the outlier rate γ,

Images