A joint estimation method of heart rate and respiratory rate based on MIMO radar

Abstract

The invention discloses a method for joint estimation of heartbeat and respiratory frequency based on MIMO radar, which belongs to the field of radar technology, especially the field of contactless vital signs monitoring technology. By adopting the ML estimation method, the estimation of heartbeat and respiratory frequency is realized. Compared with the traditional method based on time-frequency analysis, the performance of the ML estimation method is more stable, the error is smaller, and the dependence on data duration is smaller; and the present invention provides a more theoretically based estimation performance evaluation method for this field by calculating the CRB of heartbeat and respiratory frequency estimation; based on the derived CRB, the present invention provides a system parameter optimization design method, which can assist technicians in finding the optimal radar parameters, thereby improving the performance of heartbeat and respiratory frequency estimation.

Description

Technical Field

[0001] The invention belongs to the technical field of contactless vital sign monitoring, and particularly relates to a method for calculating a Cramer-Rao bound (CRB) of parameter estimation performance and optimizing system parameters in joint estimation of heartbeat and respiratory frequency. Background Art

[0002] Heart rate and respiratory rate are important parameters in human vital signs and can reflect a person's health status to a large extent. Radar technology can be used to measure heart rate and respiratory rate without contact, which has been widely studied because it greatly improves convenience, safety and privacy. At the same time, multiple-input multiple-output (MIMO) radar has been applied in this field in recent years because it has better parameter estimation performance than traditional radar.

[0003] As the lower bound of unbiased estimation, CRB can be used to evaluate the performance of parameter estimation, and can also be used as a basis for the optimization design of radar system parameters, thereby improving the performance of radar systems. Using CRB for system optimization, the results are more theoretically based and very efficient. In reference 1 (X. Shang, J. Liu, J. Li. Multiple object localization and vital sign monitoring using ir-uwb mimo radar [J]. IEEE Transactions on Aerospace and Electronic Systems, 2020, 56 (6): 4437-4450), the authors used MIMO radar for heart rate and respiratory rate estimation, while in reference 2 (Q. He, RS Blum, AM Haimovich. Noncoherent mimoradar for location and velocity estimation: More Antennas Means Better Performance [J]. IEEE Transactions on Signal Processing, 2010, 58 (7): 3661-3680), the authors gave CRB for position and velocity estimation, and the results are of great guiding significance for the design of radar systems. Therefore, in summary, in the research of heart rate and respiratory rate estimation based on MIMO radar, it is feasible and necessary to use CRB as a tool to conduct performance evaluation and system parameter optimization method research. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention proposes a method for optimizing the parameters of the MIMO radar system in the joint estimation of heartbeat and respiratory frequency to solve the problem that the existing optimization methods lack theoretical guidance. The evaluation index of the estimation performance is obtained by calculating the CRB of the parameter estimation, and this is used as the basis for system optimization. At the same time, a maximum likelihood (ML) method for simultaneous estimation of heartbeat and respiratory frequency is also given.

[0005] In order to solve the above problems and implement the proposed method, the technical solution of the present invention is: a heartbeat and breathing frequency joint estimation method based on MIMO radar, the method comprising:

[0006] Step 1: Construct the received signal vector;

[0007] For the M-transmit and N-receive MIMO radar system, the transmitting and receiving antennas are arranged around the human body and ensure that there is no obstruction; the radar signal sampling values ​​received by the N receivers are arranged into a column vector r according to the number, which can be expressed as:

[0008]

[0009] in

[0010]

[0011] is the column vector of the signal received by the nth receiving end, n=1,…,N, r n [k] is the k-th discrete time received signal value collected by the n-th receiving end, is the transposition operator;

[0012] Step 2: Construct the received signal matrix S;

[0013] S=Diag{S 1 …S n S N} (3)

[0014] Where Diag{·} represents the diagonalization operator;

[0015]

[0016] is the received signal matrix of the nth receiver, and

[0017]

[0018] Among them, d 0mn represents the sum of the distances between the transmitting and receiving antennas and the human target on the transmitting and receiving path mn, E m represents the transmission energy of each transmitting antenna, k is the index of the discrete signal at different sampling times, k = 1, ..., K, T s represents the sampling period, fc Indicates the carrier frequency;

[0019] The reflection coefficient of the signal sent by antenna m and received by antenna n is:

[0020]

[0021]

[0022] Among them, mn is the reflection coefficient on path mn, is the low-pass equivalent transmission signal emitted by the mth transmitting antenna at kT s The sampling value at the time, and the waveform meets the normalization condition τ mn [k] represents the delay on path mn, expressed as:

[0023]

[0024]

[0025] Among them, μ r,mn [k] represents the chest motion model during heartbeat, μ h,mn [k] represents the chest cavity motion model during breathing, c represents the speed of light, a h,mn and a r,mn They represent the amplitude of chest movement caused by heartbeat and breathing on path mn, respectively, h and f r They are heart rate and respiratory rate respectively;

[0026] Step 3: Construct the noise covariance matrix Q;

[0027] Q is expressed as Q = E[ww H ],in(·) H is the conjugate transpose operator, and w is the column vector composed of N receiver noises:

[0028]

[0029]

[0030] Among them, E[·] represents the expectation operation, w n represents the noise at the nth receiving end, w n [k] represents the noise at the kth sampling moment of the nth receiving end;

[0031] Step 4: Calculate the ML estimated value of heart rate and respiratory rate;

[0032] The heart rate and breathing rate that make the log-likelihood function reach the maximum value are calculated by the following formula, as f h and f r ML estimate of

[0033]

[0034] At the same time Substitute into the matrix S to calculate the ML estimate of the reflection coefficient

[0035]

[0036] Step 5: Calculate the matrix U;

[0037]

[0038] It is a 2×KMN dimensional matrix, where

[0039]

[0040] Step 6: Sub-matrices D, G and L are obtained by re-dividing and combining the elements in J(θ). J(θ) represents the Fisher information matrix of the intermediate variable. The division and combination method is:

[0041]

[0042] in j represents the imaginary unit, J τ is a KN×KMN dimensional matrix;

[0043] J τ =[J τ [1] … J τ [K]] (17)

[0044] in

[0045] J τ [k] = S τ [k]diag(ζ) (18)

[0046] is a KN×MN dimensional matrix; diag(ζ) is a diagonal matrix consisting of reflection coefficients; S τ [k] is the result of differentiating the delay of each path at time k, expressed as

[0047]

[0048] in

[0049]

[0050] is a K×M dimensional matrix, where the K-dimensional column vector e k Only the kth element is 1, and the rest are all 0. represents the Kronecker product, and

[0051]

[0052] in

[0053]

[0054] Step 7: Calculate the CRB matrix;

[0055] Using the matrix inversion lemma, the CRB matrix is ​​expressed as

[0056]

[0057] Using the first two elements on the diagonal of the CRB matrix, the heart rate estimate is calculated and respiratory rate estimation Since the ML estimate is asymptotically unbiased, the relationship between the mean square error and CRB is:

[0058]

[0059] Among them, [CRB] 1,1 and [CRB] 2,2 Respectively represent the first and second element values ​​on the diagonal of the CRB matrix;

[0060] Step 8: Optimize system parameters;

[0061] After determining the parameters and constraints that need to be optimized, the optimization problem is solved with the goal of minimizing CRB, thereby achieving the purpose of improving system performance;

[0062] All system parameters to be optimized are represented by vector X. When other parameters are determined, and is a function of the parameters to be estimated, then the optimization problem can be expressed as:

[0063]

[0064] Among them, α 1 and α 2 is a given weight coefficient used to adjust the weight of the objective function's heart rate and respiratory rate estimation performance; g i (X) represents the inequality constraint in the optimization problem, h q(X) represents the equality constraint in the optimization problem, I and Q represent the number of constraints, and the CRB of the heart rate and respiratory rate estimation is minimized under the constraints, and the optimization result of the system parameter X is calculated.

[0065] The present invention realizes the estimation of heart rate and respiratory rate by adopting ML estimation method. Compared with the traditional method based on time-frequency analysis, the performance of ML estimation method is more stable, the error is smaller, and the dependence on data duration is smaller. In addition, by calculating the CRB of heart rate and respiratory rate estimation, the present invention provides a more theoretically based estimation performance evaluation method for this field. Based on the derived CRB, the present invention provides a system parameter optimization design method, which can assist technicians in finding the optimal radar parameters, thereby improving the performance of heart rate and respiratory rate estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 Heart rate and respiratory rate estimation in MIMO radar scenario

[0067] Figure 2 Relationship between ML estimation results of heart rate and respiratory rate and RCRB when L=20,M=3,N=3

[0068] Figure 3 M=3,N=2, the variation of RCRB with frequency difference Δf DETAILED DESCRIPTION

[0069] The specific implementation of the present invention is described below in conjunction with the accompanying drawings:

[0070] Step 1 to Step 4:

[0071] In a two-dimensional Cartesian coordinate system, consider an M-transmitting and N-receiving MIMO radar system, where the transmitting antenna m (m = 1, ..., M) and the receiving antenna n (n = 1, ..., N) are located at and The human target is located at (x B ,y B ). The low-pass equivalent transmission signal emitted by the mth transmitting antenna is kT s The sampling value at time is And the waveform satisfies the normalization condition Where E m is the transmission energy of each transmitting antenna, T s is the sampling period, k (k = 1, ..., K) is the index of the discrete signal at different sampling times. Definition of discrete signal Then the time kT s The discrete signal received by the nth receiving antenna can be expressed as

[0072]

[0073] where w n [k] is a complex Gaussian white noise with independent and identical distribution for n and k, with mean 0 and variance d 0mn =d 0m +d 0n represents the nominal distance on the mn path, which includes the distance between the human chest surface and the transmitting antenna m and the distance from the receiving antenna n ζ mn is the complex reflection coefficient on the mn path, which is determined as unknown and estimated jointly in the present invention; c represents the carrier frequency; τ mn [k] represents the delay on path mn, and its value consists of two parts: one is the nominal distance d that does not change with time. 0mn , and the second is the chest movement μ that changes with time mn [k], and the delay information is expressed as

[0074]

[0075] where μ h,mn [k] and μ r,mn [k] represents the distance change between the chest surface and the radar antenna caused by heartbeat and breathing activities on the mn path, expressed as

[0076]

[0077] where a h,mn and a r,mn They represent the amplitude of chest movement caused by heartbeat and breathing on path mn respectively; f h and f r They are heart rate and respiratory rate, which are also the parameters to be estimated in this paper. Combined with the above model, the scenario of using MIMO radar to estimate heart rate and respiratory rate is as follows: Figure 1 shown.

[0078] The received signal vector of the nth receiver receiving k=1,…,K is expressed as

[0079] r n =S n ζ n +w n (29)

[0080] in

[0081]

[0082] represents the transposition symbol, ζ n is the reflection coefficient vector

[0083]

[0084] w n is the noise vector of the nth receiver

[0085]

[0086] And the signal matrix S n Expressed as

[0087]

[0088] where s mn [k] is the signal sent by antenna m and received by antenna n after removing the reflection coefficient

[0089]

[0090] The signals received by all receivers of the MIMO radar system form an NK-dimensional vector

[0091]

[0092] Where S is a KN×MN-dimensional block diagonal matrix containing all received signals

[0093] S=Diag{S 1 … S N} (36)

[0094] ζ is a vector containing all path reflection coefficients

[0095]

[0096] In this paper, it is considered to be unknown, and the vector w is the sum of all the clutter and noise at the receiving end.

[0097]

[0098] All unknown parameters are represented as a vector Θ = [f h f r ζ 11 … MN ], when the noise is a complex Gaussian random variable with zero mean, that is, w~CN(0,Q), the likelihood function of the received signal vector r is

[0099]

[0100] in(·) H represents the conjugate transpose, the covariance matrix Q = E[ww H ], the unknown parameters are contained in S and ζ. Then the log-likelihood function of the received signal is expressed as

[0101]

[0102] Among them C 1 =-r H Q -1 r-NKln(π)-ln(det(Q)) is a constant that is independent of the unknown parameters. First, solve the ML estimation expression of the reflection coefficient, take the derivative of the log-likelihood function with respect to the reflection coefficient and set it equal to 0

[0103]

[0104] The ML estimate of the reflection coefficient vector can be obtained for

[0105]

[0106] Substituting it into equation (40), we can get the expression of the log-likelihood function:

[0107]

[0108] Then the ML estimation expressions of heart rate and respiratory rate can be expressed as

[0109]

[0110] Assume that the noise vector w is independent in time and between different receivers, then the covariance matrix Q = σ w 2 I, then the maximum likelihood estimate of the parameter to be estimated is expressed as

[0111]

[0112] This equation has no closed-form solution, so and A numerical solution is required.

[0113] Step 5 to Step 8:

[0114] To facilitate calculation, the reflection coefficient in the unknown parameter is divided into real part and imaginary part

[0115] ζ mn =ζ mnR +jζ mnI (46)

[0116] Then all unknown parameters are represented by vectors as

[0117]

[0118] Where f is the parameter vector to be estimated

[0119]

[0120] ζ R and I are vectors consisting of the real and imaginary parts of the reflection coefficient

[0121]

[0122]

[0123] Then the CRB of unknown parameter estimation can be expressed as

[0124] CRB=J(θ) -1 (51)

[0125] The Fisher Information Matrix (FIM) J(θ) is expressed as

[0126]

[0127] Where E(·) represents the mathematical expectation, and the gradient operator

[0128]

[0129] Since the relationship between the likelihood function and the estimated parameters is relatively complex, a new parameter vector is introduced

[0130]

[0131] in

[0132]

[0133] in Contains the parameter to be estimated f h and f r Expressed as

[0134]

[0135] According to the chain rule, we can get FIM as

[0136]

[0137] We first calculate It can be expressed as

[0138]

[0139] Where 0 represents an all-zero matrix, I represents the identity matrix, and

[0140]

[0141] is a 2×KMN matrix, where

[0142]

[0143] Next, we derive the expression of J(θ):

[0144]

[0145] It is a (K+2)MN×(K+2)MN dimensional matrix. According to equation (40), S and ζ are functions of θ. The element in the i-th row and q-th column of the matrix J(θ) can be calculated by the following equation:

[0146]

[0147] According to the form of the intermediate vector, we can use the block matrix to represent J(θ) as

[0148]

[0149] in

[0150]

[0151]

[0152]

[0153] First calculate J τ , is a KN×KMN dimensional matrix. From the form of the delay parameter vector τ, we can know that the derivatives for time k=1 to K can be performed separately. Matrix J τ Can be divided into K block matrices

[0154] J τ =[J τ [1] … J τ [K]] (67)

[0155] in

[0156] J τ [k] = S τ [k]diag(ζ) (68)

[0157] is a KN×MN dimensional matrix, and diag(ζ) is a diagonal matrix consisting of reflection coefficients. τ [k] is the result of differentiating the delay of each path at time k, expressed as

[0158]

[0159] in

[0160]

[0161] is a K×M dimensional matrix, where the K-dimensional column vector e k Only the kth element is 1, and the rest are all 0. represents the Kronecker product, and

[0162]

[0163] in

[0164]

[0165] Substituting equation (58) and equation (63) into equation (57), we can obtain FIM:

[0166]

[0167] Finally, the complete CRB can be expressed as the inverse of FIM. According to the matrix inversion lemma, the CRB of heart rate and respiratory rate are the first two values ​​of the diagonal of the complete CRB matrix. Therefore, we only need to calculate a 2×2 block matrix in the upper left corner of the complete CRB matrix, which can be expressed as

[0168]

[0169] We are concerned with heart and breathing rates, i.e., f h and f r The estimation performance of , so we only need the first two elements on the diagonal of the matrix. Since ML estimation is asymptotically unbiased, the relationship between the MSE and CRB of the two is

[0170]

[0171] So far, we have obtained the MSE lower bound expression for heart rate and respiratory rate estimation, which is used as a performance evaluation indicator and can be applied to solving optimization problems.

[0172] Next, we will conduct simulation verification and illustrate the effectiveness of the method with the attached figure. In the simulation, the radars are evenly arranged on a circle with the target as the center and a radius of 2m. The transmitting end uses a frequency-extended sinusoidal pulse signal. The transmitting signal of the mth antenna is

[0173]

[0174] Where L is the number of transmitted pulses; Δf is the frequency increment between transmitted signals, and its value will affect the orthogonality between transmitted signals. Here, it is taken as 1 / T d ; The values ​​of the remaining parameters are: f s =240Hz, T p =0.1s, Td =0.05s, the duty cycle is defined as D=T d / T p ; The vital sign parameters are set to f h =1.52Hz, f r =0.29Hz, a h =0.8mm, a r =5mm. To simplify the analysis, the heart rate and breathing amplitude parameters of different paths are taken as the same value, i.e., a h,mn =a h ,a r,mn =a r (m=1…M,n=1…N), it has been verified that more complex value selection does not affect the analysis of other parameters; SNR is defined as

[0175]

[0176] Where T = LT p is the total observation time, E 1 =…=E M =1 / MLT d , the change of SNR in numerical simulation is achieved by adjusting the noise variance to achieve.

[0177] like Figure 2 The figure shows the heartbeat frequency f when the pulse number L = 20, the number of transmitting and receiving antennas M = 3, N = 3 h and respiratory rate f r The curves of ML estimation RMSE and CRB square root (RCRB) with SNR change show that the red RMSE dashed line has obvious threshold phenomenon. The threshold is defined as the SNR position where the slope of the RMSE curve changes dramatically. When the SNR is higher than the threshold, the estimation performance and stability will be significantly improved. It can be found that the threshold of heart rate estimation is between -10 and 0dB, and the threshold of respiratory rate estimation is between -20 and -10dB. When the SNR is lower than the threshold, the estimation error is large and cannot be used as an effective estimation. When the SNR is higher than the threshold, the estimated RMSE converges to RCRB, and a stable estimation result can be obtained. At the same time, it can be found that the estimation accuracy of respiratory rate is higher than that of heart rate estimation, which is reflected in the lower threshold and smaller CRB of the former. The reason is that on the one hand, the heart rate itself is greater than the respiratory rate, and on the other hand, the impact of the heartbeat on the chest surface is much smaller than that of breathing, which makes the estimation more difficult.

[0178] Figure 3 The figure shows the variation of the heart rate and respiratory rate estimation CRB with the frequency difference Δf. The value of Δf is 0.2 / T d Increased to 5 / T d , and M=3, N=2, f s=1.2KHz, L=1, SNR=20dB, the rest of the parameters remain the same Figure 2 The simulation is the same. Figure 3 From the results, we can find that with the increase of Δf, the orthogonality of the transmitted signals continues to improve, and while the estimation performance improves overall, there are certain fluctuations. When the frequency difference is high enough, this fluctuation will gradually converge to a stable value, and the curve will begin to change dramatically and then tend to be flat, indicating that Δf does not need to be too large in the actual radar system design. Using CRB as a tool to optimize the frequency difference between the transmitted signals according to different situations can achieve the purpose of saving system resources.

Claims

1. A method for jointly estimating heart rate and respiratory rate based on MIMO radar, the method comprising: Step 1: Construct the received signal vector; For the M-transmit and N-receive MIMO radar system, the transmitting and receiving antennas are arranged around the human body and ensure that there is no obstruction; the radar signal sampling values ​​received by the N receivers are arranged into a column vector r according to the number, which can be expressed as: in is the column vector of the signal received by the nth receiving end, n=1,…,N, r n [k] is the k-th discrete time received signal value collected by the n-th receiving end, is the transposition operator; Step 2: Construct the received signal matrix S; S=Diag{S1 … S n … S N } (3) Where Diag{·} represents the diagonalization operator; is the received signal matrix of the nth receiver, and Among them, d 0mn represents the sum of the distances between the transmitting and receiving antennas and the human target on the transmitting and receiving path mn, E m represents the transmission energy of each transmitting antenna, k is the index of the discrete signal at different sampling times, k = 1, ..., K, T s represents the sampling period, f c Indicates the carrier frequency; The reflection coefficient of the signal sent by antenna m and received by antenna n is: Among them, mn is the reflection coefficient on path mn, is the low-pass equivalent transmission signal emitted by the mth transmitting antenna at kT s The sampling value at the time, and the waveform meets the normalization condition τ mn [k] represents the delay on path mn, expressed as: Among them, μ r,mn [k] represents the chest motion model during heartbeat, μ h,mn [k] represents the chest cavity motion model during breathing, c represents the speed of light, a h,mn and a r,mn They represent the amplitude of chest movement caused by heartbeat and breathing on path mn, respectively, h and f r They are heart rate and respiratory rate respectively; Step 3: Construct the noise covariance matrix Q; Q is expressed as Q = E[ww H ],in(·) H is the conjugate transpose operator, and w is the column vector composed of N receiver noises: Among them, E[·] represents the expectation operation, w n represents the noise at the nth receiving end, w n [k] represents the noise at the kth sampling moment of the nth receiving end; Step 4: Calculate the ML estimated value of heart rate and respiratory rate; The heart rate and breathing rate that make the log-likelihood function reach the maximum value are calculated by the following formula, as f h and f r ML estimate of At the same time Substitute into the matrix S to calculate the ML estimate of the reflection coefficient Step 5: Calculate the matrix U; U=[U f [1] … In f [k] … In f [K]] + (14) It is a 2×KMN dimensional matrix, where Step 6: Sub-matrices D, G and L are obtained by re-dividing and combining the elements in J(θ). J(θ) represents the Fisher information matrix of the intermediate variable. The division and combination method is: in j represents the imaginary unit, J τ is a KN×KMN dimensional matrix; I τ =[J τ [1] … J τ [K]] (17) in J τ [k]=S τ [k]diag(ζ) (18) is a KN×MN dimensional matrix; diag(ζ) is a diagonal matrix consisting of reflection coefficients; S τ [k] is the result of differentiating the delay of each path at time k, expressed as in is a K×M dimensional matrix, where the K-dimensional column vector e k Only the kth element is 1, and the rest are all 0. represents the Kronecker product, and in Step 7: Calculate the CRB matrix; Using the matrix inversion lemma, the CRB matrix is ​​expressed as Using the first two elements on the diagonal of the CRB matrix, the heart rate estimate is calculated and respiratory rate estimation Since the ML estimate is asymptotically unbiased, the relationship between the mean square error and CRB is: Among them, [CRB] 1,1 and [CRB] 2,2 Respectively represent the first and second element values ​​on the diagonal of the CRB matrix; Step 8: Optimize system parameters; After determining the parameters and constraints that need to be optimized, the optimization problem is solved with the goal of minimizing CRB, thereby achieving the purpose of improving system performance; All system parameters to be optimized are represented by vector X. When other parameters are determined, and is a function of the parameters to be estimated, then the optimization problem can be expressed as: Among them, α1 and α2 are given weight coefficients, which are used to adjust the weight of the heart rate and respiratory rate estimation performance of the objective function; g i (X) represents the inequality constraint in the optimization problem, h q (X) represents the equality constraint in the optimization problem, I and Q represent the number of constraints, and the CRB of the heart rate and respiratory rate estimation is minimized under the constraints, and the optimization result of the system parameter X is calculated.

Images