A Method for Identifying the s-Domain Parameters of a High-Order Complex Pulse

Abstract

The invention discloses a high-order complex pulse s A domain parameter identification method is used to identify high-order complex pulse signals output by the front-end simulation system of nuclear radiation detectors. s The parameter identification of the domain characteristic function provides guarantee for the subsequent waveform shaping and accurate measurement of the energy spectrum. This method pre-identifies the structure of the characteristic function based on the spectrum information of the pulse signal, and obtains the initial value of the parameter of the characteristic function, which narrows the parameter search range and overcomes the adverse consequences of the parameter optimization process caused by the uncertainty of the function structure and the unknown range of parameter values ​​(such as non-convergence or non-global optimality), ensuring that the parameter search process of the characteristic function converges quickly and the result is the global optimal; using this parameterized form, it can not only realize signal filtering, but also reduce the dimension of the subsequent shaping algorithm.

Description

Technical Field

[0001] The invention relates to an S-domain parameter identification method for high-order complex pulses. Background Art

[0002] Due to the diversity or complexity of the front-end analog system of nuclear radiation detectors, the pulse signal it outputs is not just a simple typical form (such as exponential signal, step signal, etc.); especially when there are multiple conditioning links in the circuit, the pulse signal shows high-order complexity; in addition, affected by the test conditions and environment, as well as the change of the parameters of the circuit components over time, the output pulse waveform will also change. In the measurement of nuclear energy spectrum, it is necessary to extract the amplitude after shaping these high-order complex pulses; if they are fitted or approximated with simple typical signals (such as exponential signal, step signal, etc.), the pulses after shaping will be distorted and even invalid; in addition, even if the waveform shaping algorithm is accurate in a certain measurement, due to the change of system parameters over time, the waveform shaping algorithm used may cause the pulses after shaping to be inconsistent with the early ones in the next measurement. Therefore, it is necessary to accurately identify the parameters of the high-order complex pulses before shaping, and use the identified parameters to construct the waveform shaping algorithm, so as to achieve accurate shaping and finally obtain a high-precision nuclear energy spectrum.

[0003] The purpose of this method is to accurately identify the parameters of the pulse signal through the form of characteristic functions, so as to provide guarantee for the subsequent waveform shaping and accurate measurement of the energy spectrum. It is assumed that the structure of the characteristic function is determined according to the spectrum information of the pulse signal and the initial value of the parameter is obtained. This initial value of the parameter has the characteristic of being quite close to the true value, so as to narrow the parameter search range and overcome the complexity, drift and non-convergence of the parameter optimization process caused by the uncertainty of the function structure and the unknown range of parameter values, and even the parameters finally searched are not globally optimal. Summary of the invention

[0004] The object of the present invention is to disclose a method for identifying the s-domain parameters of high-order complex pulses, which is used to identify the parameters of the s-domain characteristic function of the high-order complex pulse signal output by the front-end analog system of a nuclear radiation detector, so as to provide guarantee for subsequent waveform shaping and accurate energy spectrum measurement. According to the frequency spectrum information of the pulse signal, this method pre-judges the framework of the characteristic function (i.e., the types and quantities of typical links included), and obtains the initial parameter values of the typical links (integral, first-order differential, second-order differential, inertia link, oscillation link, etc.) contained in the characteristic function. Such initial parameter values are quite close to the true values, narrowing the parameter search range and overcoming the adverse consequences of the parameter optimization process caused by the uncertainty of the function architecture and the unknown parameter value range: namely, the optimization process is complex, drifting and non-convergent, and the finally searched parameters are not globally optimal; ensuring that the parameter search process of the characteristic function converges quickly and the result is globally optimal; adopting this parameterized form can not only realize signal filtering but also reduce the dimension of the subsequent shaping algorithm.

[0005] The identification of the s-domain parameters of the high-order complex pulse signal output by the nuclear radiation detection channel of the present invention is realized through the following steps ① to ⑤.

[0006] Step ①: Use the s-domain expression as the characteristic function of the pulse signal, which includes v integral links, N inertia links, M first-order differential links, P second-order differential links, and Q oscillation links.

[0007] Step ②: Obtain the frequency characteristic curve of the pulse signal and its asymptote, correct the slope of the asymptote, and obtain the intersection points of the asymptote and the intersection points with the coordinate axes, which are carried out according to the following steps A to B:

[0008] A. Obtain the frequency characteristic curve of the pulse signal, and obtain the asymptote of the frequency characteristic curve. The asymptotes are numbered from left to right along the abscissa.

[0009] B. Correct the slope of the asymptote, and obtain the intersection coordinates of the asymptote and the intersection coordinates with the coordinate axes.

[0010] Step ③: According to the intersection coordinates and the corrected asymptote slope in Step ②, obtain the initial parameter values of the characteristic function: (1) Solving the parameters of the integral link; (2) Solving the parameters of the first-order differential link; (3) Solving the parameters of the inertia link; (4) Solving the parameters of the second-order differential link; (5) Solving the parameters of the oscillation link.

[0011] Step ④: Substitute the parameter values obtained in Step ③ into the characteristic function, and obtain the proportional coefficient, zeros and poles, and use them as the initial values of the subsequent recursive algorithm.

[0012] In step ⑤, iterative recursion is used to search and update each parameter of the characteristic function; when the search process ends, if the error meets the requirements, the final parameters of the characteristic function are obtained.

[0013] The parameters of the characteristic function are identified through steps ① to ⑤.

[0014] Taking the trapezoidal shaping of the pulse waveform as an example, an algorithm for trapezoidal shaping of the pulse signal output by the front-end detection system of the nuclear instrument using the discrete form of the identified characteristic function is given.

[0015] The beneficial effects of the present invention are:

[0016] Due to the diversity or complexity of the front-end analog system of the nuclear radiation detector, the pulse signal it outputs is not just a simple typical form (such as exponential signal and step signal, etc.); especially in the case where there are multiple conditioning links in the circuit, the pulse signal shows high-order complexity; in addition, affected by the test conditions and environment, and the change of the parameters of its own circuit components over time, the pulse waveform it outputs will also change. In nuclear energy spectrum measurement, it is necessary to shape these high-order complex pulses and then extract the amplitude; if simple typical signals (such as exponential signal and step signal, etc.) are used for fitting or approximation, the shaped pulse will be distorted, and even the shaping will be invalid; in addition, even if the waveform shaping algorithm is accurate and error-free during a certain measurement, due to the change of system parameters over time, the waveform shaping algorithm used may cause the shaped pulse to be inconsistent with the early stage during the next measurement. Therefore, it is necessary to accurately identify the parameters of the high-order complex pulses before shaping, and use the identified parameters to construct the waveform shaping algorithm, so as to achieve accurate shaping and finally obtain a high-precision nuclear energy spectrum.

[0017] Aiming at the high-order complexity of the pulse signal output by the detector system, through the form of the characteristic function, the accurate identification of the pulse signal parameters is realized, providing guarantee for the subsequent waveform shaping and accurate energy spectrum measurement. The advantages of this method are: (1) This method pre-judges the framework of the characteristic function (that is, the types and quantities of the typical links included) according to the spectrum information of the pulse signal, and obtains the initial values of the parameters of the typical links (integral, first-order differential, second-order differential, inertia link, oscillation link, etc.) contained in the characteristic function. These initial parameter values are quite close to the true values, narrowing the parameter search range and overcoming the adverse consequences of the parameter optimization process caused by the uncertainty of the function framework and the unknown parameter value range: that is, the optimization process is complex, drifting and non-convergent, and the finally searched parameters are not globally optimal; (2) Ensure that the parameter search process of the characteristic function converges quickly and the result is globally optimal; (3) Using this parameterized form can not only realize signal filtering, but also reduce the dimension of the subsequent shaping algorithm. Description of the Drawings

[0018] Figure 1 Schematic diagram of the frequency characteristics curve of the characteristic function and its asymptote Figure 2 Flow chart of this recognition method Specific implementation mode

[0019] The following will describe in detail the embodiments of the present invention with reference to the accompanying drawings. These embodiments are implemented on the premise of the technical solution of the present invention, and detailed implementation manners and processes are given. However, the protection scope of the present invention is not limited to the following embodiments.

[0020] The present invention relates to a method for identifying the s-domain parameters of high-order complex pulses, which is realized through the following steps ① to ⑤, and is used to identify the parameters of the s-domain characteristic function of the pulse signal output by the front-end analog system of a nuclear radiation detector.

[0021] Step ① Use the following s-domain expression as the characteristic function of the pulse signal

[0022]

[0023] Each link it contains is as follows

[0024] v integral links N inertia links

[0025] M first-order differential links P second-order differential links

[0026] Q oscillation links

[0027] Step ② Obtain the frequency characteristics curve of the pulse signal and its asymptote, correct the slope of the asymptote, and obtain the intersection points of the asymptote and the intersection points with the coordinate axes. Perform according to the following steps A to B

[0028] A. Obtain the frequency characteristics curve of the pulse signal and obtain the asymptote of the frequency characteristics curve

[0029] Assume that the measured pulse signal is h(t), and obtain its frequency domain expression as follows

[0030]

[0031] Let y = 20lgH(jω), x = lgω; draw the x-y curve with x as the abscissa and y as the ordinate, that is, the amplitude-frequency characteristic curve of H(jω); and draw the asymptote of the straight line segment on the characteristic curve with a slope of 20n (n is an integer). The asymptotes are all straight lines. The slopes of the asymptotes are as follows

[0032]

[0033] where n is an integer, K 0i and K i are the slope of the i-th straight line segment and the slope of the asymptote corresponding to this straight line segment respectively; the straight line segments are numbered from left to right along the abscissa, that is, the i-th one is on the left of i + 1, and the slopes of the asymptotes corresponding to the straight line segments are marked as K1, K2,..., K M+N+Q+P+1 .

[0034] An example is given to illustrate. Suppose Figure 1 (a) is the frequency characteristic curve of h(t), and the asymptotes corresponding to 5 straight line segments on its frequency characteristic curve are drawn as Figure 1 (b) shows; the intersection point of the first asymptote and the ordinate is A1, and its slope K1 is -40. In fact, it can be deduced that it corresponds to the integral link K / s 2 ; the slope K2 of the second asymptote is 0, and the intersection point with the first asymptote is A2. It can be deduced that it corresponds to the second-order differential link τ1s 2 +γ1s + 1; the slope K3 of the third asymptote is -40, and the intersection point with the second asymptote is A3. It can be deduced that it corresponds to the oscillation link The slope K4 of the fourth asymptote is -20, and the intersection point with the third asymptote is A4. It can be deduced that it corresponds to the first-order differential link b1s + 1; the slope K5 of the fifth asymptote is -40, and the intersection point with the fourth asymptote is A5. It can be deduced that it corresponds to the inertia link From this, the s-domain expression H(s) of h(t) can be obtained as

[0035] B. Correct the slopes of the asymptotes, and obtain the intersection coordinates of the asymptotes and the intersection coordinates with the coordinate axes:

[0036] The corrected slopes of the asymptotes are K1′, K2′,..., K′ M+N+Q+P+1 :

[0037]

[0038] That is: K i ′ = K i - K i-1 , 1 < i ≤ M + N + Q + P + 1 (5)

[0039] The intersection point of the first asymptote and the ordinate is A1, and the intersection point of the i-th asymptote and the (i - 1)-th asymptote is A i , where 2 ≤ i ≤ M + N + P + Q + 1; that is, these intersection points are A1, A2,..., A M+N+Q+P+1 , and obtain their corresponding coordinates as:

[0040] For subsequent convenience, let M + N + Q + P + 1 = SUM.

[0041] Step ③: According to the intersection coordinates and the corrected asymptote slope in Step ②, find the initial values of the parameters of the characteristic function G(s), and proceed as follows in (1) to (5):

[0042] (1) Solve the parameters K and v of the integral link as follows:

[0043] v = -K1 / 20 (6)

[0044] From we get:

[0045]

[0046]

[0047] Let K1″ = K1′, and the parameter set of the integral link is expressed as

[0048] (2) Solve the parameters M and b of the first-order differential link n as follows in steps a - c:

[0049] a. Let m = 0, i = 2, j = 1;

[0050] b. If K i ′ = 20, then m = m + 1, j = j + 1, and modify the parameter values as follows:

[0051]

[0052] c. i = i + 1; if i > SUM, stop; otherwise return to b.

[0053] After solving, the parameter set of each first-order differential link is expressed as:

[0054]

[0055] That is: (K j ″, X(j), Y(j), b j-1 s + 1), 2 ≤ j ≤ M + 1 (11)

[0056] (3) Solve the parameters N and a of the inertia link n as follows in steps a - c:

[0057] a. Let n = 0, i = 2, j = M + 1;

[0058] b. If K i ′ = -20, then n = n + 1, j = j + 1, and modify the parameter values as follows:

[0059]

[0060] c. i = i + 1; If i > SUM, stop; otherwise return to b.

[0061] The parameter set of each inertia link after solution is expressed as:

[0062]

[0063] That is:

[0064] (4) Parameters τ p , γ p of the second-order differential link are solved according to the following steps a-c:

[0065] a. Let p = 0, i = 2, j = M + N + 1;

[0066] b. If K i ′ = 40, then p = p + 1, j = j + 1, and modify the parameter values as follows:

[0067]

[0068] c. i = i + 1; If i > SUM, stop; otherwise return to b.

[0069] The parameter set of each second-order differential link after solution is expressed as:

[0070]

[0071] That is:

[0072] (K j ″, X(j), Y(j), τ j-M-N-1 s 2 + γ j-M-N-1 s + 1), M + N + 2 ≤ j ≤ M + N + P + 1 (17)

[0073] (5) Parameters Q, T q , λ q of the oscillation link are solved according to the following steps a-c:

[0074] a. Let q = 0, i = 2, j = M + N + P + 1;

[0075] b. If K i ′ = -40, then q = q + 1, j = j + 1, and modify the parameter values as follows:

[0076]

[0077] c. i = i + 1; If i > SUM, stop; Otherwise, return to b.

[0078] The parameter set of each oscillation link after solution is expressed as:

[0079]

[0080] That is:

[0081]

[0082] In step ④, substitute the parameter values obtained in step ③ into the following formula G(s), and calculate k, zero point δ j and pole η i .

[0083]

[0084] G(s) has a total of N + 2Q + v poles: η i , i = 1,..., N + 2Q + v; Among them, v poles are all 0. For convenience, assume the first v poles are 0, that is, η i = 0 (i = 1, 2,..., v).

[0085] Calculate k, δ j , η i After that, for the convenience of subsequent formula derivation, take them as the initial values, and let:

[0086]

[0087] In step ⑤, use iterative recursion to search and update each parameter of the characteristic transfer function G(s); Let the number of searches be l, and assume the initial value l = 0; k(l), δ j (l), η i (l) represents the latest values of each parameter corresponding to G(s) after the l-th search is completed, and G l (s) represents the expression of G(s) after the l-th search is completed; Perform search and update according to the following steps 5S1 to 5S2:

[0088] In step 5S1, the update algorithms corresponding to the parameters k(l), δ j (l), η i (l) are shown in formulas (22) to (24):

[0089] k(l + 1) = k(l) + △J k (l)μ k (l) (22)

[0090]

[0091]

[0092] μ in Formulas (22) to (24) k (l), △J k (l), and as shown in Formulas (25) to (28):

[0093] μ k (l) = k(l)μ0; μ δj (l) = δ j (l)μ0; μ ηi (l) = η i (l)μ0 (25)

[0094]

[0095]

[0096]

[0097] g in Formulas (26) to (28) l (t), g l_Δk (t), and as shown in Formulas (29) to (32):

[0098]

[0099]

[0100]

[0101]

[0102] G in Formulas (29) to (32) l (s), G l_Δk (s), and as shown in Formulas (33) to (36):

[0103]

[0104]

[0105]

[0106]

[0107] In the above formulas: μ0 = 0 to 0.1; Δk = (0 to 0.1)k(0); Δη i = (0 to 0.1)ηi (0),

[0108] Δδ i =(0 to 0.1)δ i (0). It should be noted that: we assume that there are no other repeated roots in the actual system; if there are repeated roots, the corresponding time-domain components are processed in the same way as the time-domain components corresponding to the repeated roots when s = 0, as shown in formulas (29) to (32).

[0109] Step 5S2 If the error e(l) ≤ E0 or the search times l ≥ l max , then the search process for the parameters of G(s) ends; when the search process ends, if e(l) ≤ E0, then the parameters k(l), δ j (l), η i (l) of G(s) are obtained. Otherwise, set l = l + 1 and return to Step 5S1 to continue; E0, l max are the end conditions of the search process set as needed; e(l) is calculated according to the following formula (37).

[0110]

[0111] The parameters of the characteristic function G(s) are identified through Steps ① to ⑤, and G l (s) at the end of the search process is the final result of G(s), that is

[0112]

[0113] The above is the characteristic function in the form of poles and zeros, and it can also be obtained through the following conversion formula to the form of the characteristic function composed of typical links (integral, first-order differential, second-order differential, inertia link, oscillation link, etc.) in Step ①.

[0114]

[0115] The method for identifying the characteristic function of the pulse signal output by the detection system as described above first analyzes the trend of the frequency characteristic curve of the pulse signal, and determines the framework of the characteristic function through the asymptote of the frequency characteristic curve, that is, the types and quantities of the typical links included; then, according to the intersection points of the asymptotes and the slope of the corrected asymptote, the initial values of the parameters of the typical links (integral, first-order differential, second-order differential, inertia link, oscillation link, etc.) included in the characteristic function are obtained; then, the characteristic function is converted into a function form characterized by zeros and poles; finally, through iterative recursion, each parameter such as the proportionality coefficient, zeros, and poles of the characteristic function is searched and updated until the optimal parameters are finally searched. By such a method, the framework of the characteristic function of the signal (that is, the types and quantities of the typical links included) is pre-determined, and the initial values of the parameters of the typical links (integral, first-order differential, second-order differential, inertia link, oscillation link, etc.) included in the characteristic function are obtained. Such initial parameter values are quite close to the true values, narrowing the parameter search range and overcoming the adverse consequences of the parameter optimization process caused by the uncertainty of the function framework and the unknown parameter value range: that is, the optimization process is complex, drifting, and non-convergent, and the parameters finally searched are not globally optimal. Ensure the rapid convergence of the characteristic function parameter search process and the global optimality of the results.

[0116] If the z-domain expression G(z) is obtained by performing the following transformation on the searched G(s), it is convenient for the shaping of the nuclear pulse

[0117]

[0118] The algorithm provides great convenience. Taking trapezoidal shaping as an example, assuming that the pulse signal output by the front-end detection system of the nuclear instrument is X(z), and the set trapezoid is Y(z), then the trapezoidal shaping algorithm of X(z) is as follows:

[0119]

[0120] The pulse rise time corresponding to the trapezoidal pulse Y(z) is t a , and the sum of the rise time and the flat top time is t b , and the width of the pulse is t c , and its z-domain expression is as follows:

[0121]

[0122] In formula (42), A0 is the flat top amplitude of the trapezoidal pulse, n a =t a / T s , n b =t b / T s , n c =t c / Ts , T S is the sampling period.

[0123] In the above embodiments of the present invention, a method for identifying the pulse signal characteristic function parameters output by the detection system has been described in detail. However, it should be noted that the above is only one embodiment of the present invention. When other similar characteristic function identification methods involve using the method proposed in this article, the present invention is still valid. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for identifying S-domain parameters of high-order complex pulses, characterized in that: The parameter identification of the s-domain characteristic function of the pulse signal is achieved through the following steps ①~⑤: Step ① Use the following s-domain expression as the characteristic function of the pulse signal The various links it includes are as follows: v points points: N inertia links: M first-order differential links: P second-order differential links: Q oscillation links: Step ② obtain the frequency characteristic curve and asymptote of the pulse signal, correct the slope of the asymptote, and obtain the intersection of the asymptote and the intersection with the coordinate axis, according to the following steps A to B: A. Obtain the frequency characteristic curve of the pulse signal and the asymptote of the frequency characteristic curve: Assume that the measured pulse signal is h(t), and obtain its frequency domain expression as follows: Let y = 20lg|H(jω)|, x = lgω; with x as the abscissa and y as the ordinate, draw the xy curve, that is, the amplitude-frequency characteristic curve of H(jω); and draw the asymptotes of the straight line segment on the characteristic curve with a slope of 20n. The asymptotes are all straight lines; the slopes of the asymptotes are as follows: Where n is an integer, K 0i and K i are the slope of the i-th straight line segment and the slope of the asymptote corresponding to the straight line segment; the straight line segments are numbered from left to right along the horizontal axis, that is, the i-th straight line is to the left of i+1, and the slopes of the asymptotes corresponding to the straight line segments are marked as K1, K2, ..., K M+N+Q+P+1 ; B. Correct the slope of the asymptote and find the coordinates of the intersection of the asymptote and the intersection with the coordinate axis: The slopes of the corrected asymptotes are K1′, K2′, ..., K′ M+N+Q+P+1 : That is: K i ′=K i -K i-1 , 1<i≤M+N+Q+P+1(5) The intersection of the first asymptote and the ordinate is A1, and the intersection of the i-th asymptote and the i-1-th asymptote is A i , where 2≤i≤M+N+P+Q+1; that is, these intersection points are A1, A2, ..., A M+N+Q+P+1 , and find its corresponding coordinates: For the convenience of the following, let M+N+Q+P+1=SUM; Step ③: According to the intersection coordinates in step ② and the corrected asymptote slope, calculate the initial value of the parameter of the characteristic function G(s), as follows (1) to (5): (1) The integration phase parameters K, v are solved as follows: v=-K1 / 20 (6) Depend on have to: Let K1″=K1′, and the parameter set of the integral link is expressed as (2) First-order differential link parameter M,b n To solve, follow the following ac steps: a. Let m = 0, i = 2, j = 1; b. If K i ′=20, then m=m+1, j=j+1, and the parameter values ​​are modified as follows: c. i=i+1; if i>SUM, stop; otherwise return to b; After solving, the parameter set of each first-order differential link is expressed as: That is: (K j ″,X(j),Y(j),b j-1 s+1),2≤j≤M+1(11)(3) Inertia link parameter N,a n To solve, follow the following ac steps: a. Let n = 0, i = 2, j = M + 1; b. If K i '=-20, then n=n+1, j=j+1, and the parameter values ​​are modified as follows: c. i=i+1; if i>SUM, stop; otherwise return to b; After solving, the parameter set of each inertia link is expressed as: Right now: (4) Second-order differential link parameter τ p , γ p To solve, follow the following ac steps: a. Let p = 0, i = 2, j = M + N + 1; b. If K′ i =40, then p=p+1, j=j+1, and the parameter values ​​are modified as follows: c. i=i+1; if i>SUM, stop; otherwise return to b; After solving, the parameter set of each second-order differential link is expressed as: Right now: (K′′ j ,X(j),Y(j),τ j-M-N-1 s 2 +γ j-M-N-1 s+1),M+N+2≤j≤M+N+P+1 (17) (5) Oscillation parameters Q, T q ,λ q To solve, follow the following ac steps: a. Let q = 0, i = 2, j = M + N + P + 1; b. If K′ i =-40, then q=q+1, j=j+1, and the parameter values ​​are modified as follows: c. i=i+1; if i>SUM, stop; otherwise return to b; After solving, the parameter set of each oscillation link is expressed as: Right now: Step ④: Substitute the parameter values ​​obtained in step ③ into the following equation G(s), and calculate k and zero point δ according to the following equations: j and the pole η i ; G(s) has a total of N+2Q+v poles: η i , i=1,...,N+2Q+v; where v poles are all 0; for convenience, the first v poles are set to 0, that is, η i =0(i=1, 2, ..., v); Find k,δ j ,η i After that, for the convenience of subsequent formula derivation, it is used as the initial value and set: Step ⑤ Use iterative recursion to search and update the parameters of the characteristic recursive function G(s); set the number of searches to l, and set the initial value l = 0; k(l),δ j (l),η i (l) represents the latest value of each parameter of G(s) after the lth search is completed, G l (s) represents the expression of G(s) after the lth search is completed; the search is updated according to the following steps 5S1-5S2: Step 5S1 parameter k(l),δ j (l),η i (l) The corresponding update algorithm is shown in formulas (22) to (24): <h2 style=";text-align:left;direction:ltr">k(l+1) = k(l)+ΔJ<h2 style=";text-align:left;direction:ltr"> k <h2 style=";text-align:left;direction:ltr"> (l)μ<h2 style=";text-align:left;direction:ltr"> k <h2 style=";text-align:left;direction:ltr"> (l) (22) μ in formulas (22) to (24) k (l), ΔJ k (l), and As shown in formulas (25) to (28): In formulas (26) to (28), g l (t), g l_ Δ k (t), g l_ηi (t) and As shown in formulas (29) to (32): In formulas (29) to (32), G l (s), G l_Δk (s), and As shown in formulas (33) to (36): In the above formula: μ0 = 0 ~ 0.1; Δk = (0 ~ 0.1) k (0); Δη i =(0~0.1)η i (0), Δδ i =(0~0.1)δ i (0); It should be noted that we assume that there are no other repeated roots in the actual system; if there are repeated roots, the corresponding time domain components are processed in the same way as the time domain components corresponding to the repeated roots when s = 0, as shown in formulas (29) to (32); Step 5S2 If error e(l)≤E0 or search times l≥l max , then the search process for G(s) parameters ends; At the end of the search process, if e(l)≤E0, the parameters k(l),δ of G(s) are obtained. j (l),η i (l); otherwise, set l = l + 1 and return to step 5S1 to continue; E0, l max is the end condition of the search process set as needed; e(l) is calculated according to the following formula (37); Through steps ① to ⑤, the parameters of the characteristic function G(s) are identified. At the end of the search process, G l (s) is the final result of G(s), as shown in formula (38).

Images