JP6795453B2 - Measuring device and measuring method - Google Patents
Measuring device and measuring method Download PDFInfo
- Publication number
- JP6795453B2 JP6795453B2 JP2017100683A JP2017100683A JP6795453B2 JP 6795453 B2 JP6795453 B2 JP 6795453B2 JP 2017100683 A JP2017100683 A JP 2017100683A JP 2017100683 A JP2017100683 A JP 2017100683A JP 6795453 B2 JP6795453 B2 JP 6795453B2
- Authority
- JP
- Japan
- Prior art keywords
- vivo
- error
- value
- estimation
- weighting
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
Description
本発明は、非侵襲な生体内成分濃度の測定技術に関する。 The present invention relates to a technique for measuring non-invasive in vivo component concentration.
高齢化が進み、成人病に対する対応が大きな課題になりつつある。血糖値などの検査は血液の採取が必要なために患者にとって大きな負担である。そのため、血液を採取しない非侵襲な生体内成分濃度の測定装置が注目されている。 As the population ages, dealing with adult diseases is becoming a major issue. Testing such as blood glucose level is a heavy burden for patients because it requires blood sampling. Therefore, a non-invasive device for measuring the concentration of in vivo components that does not collect blood is drawing attention.
非侵襲な生体内成分濃度の測定装置として、例えば、光音響法を用いた測定装置や、誘電分光法を用いた測定装置が提案されている。光音響法は、皮膚内に電磁波を照射し、測定対象とする血液成分、例えば、グルコース分子に電磁波を吸収させ、グルコース分子からの熱の放射によって局所的に熱膨張を起こし、熱膨張によって生体内から発生した音波を観測する（例えば、非特許文献１参照。）。 As a non-invasive measuring device for the concentration of components in the living body, for example, a measuring device using a photoacoustic method and a measuring device using a dielectric spectroscopy have been proposed. In the photoacoustic method, an electromagnetic wave is irradiated into the skin, the blood component to be measured, for example, a glucose molecule absorbs the electromagnetic wave, and the radiation of heat from the glucose molecule causes local thermal expansion, which is generated by the thermal expansion. Observe the sound waves generated from inside the body (see, for example, Non-Patent Document 1).
図８は、従来の光音響法を用いた生体内成分濃度の測定装置１００の構成例を示す図である。従来の測定装置１００は、光信号発生装置１０１と、光信号出射装置１０２と、音響素子１０３と、音響信号計測装置１０４と、信号処理装置１０５と、濃度演算装置１０６と、データベース１０７とから構成される。 FIG. 8 is a diagram showing a configuration example of the in-vivo component concentration measuring device 100 using a conventional photoacoustic method. The conventional measuring device 100 includes an optical signal generator 101, an optical signal emitting device 102, an acoustic element 103, an acoustic signal measuring device 104, a signal processing device 105, a concentration calculation device 106, and a database 107. Will be done.
従来の測定装置１００では、被検体である生体に光が照射され、照射された光は、生体の内部に光音響信号である音波を発生させる。この音波は電気信号に変換され、信号解析がなされて、背景成分及び対象成分が混合されている生体における、グルコース値等の生体内成分濃度が測定される（例えば、特許文献１参照。）。しかし、グルコースと電磁波の相互作用は小さく、また生体に安全に照射しうる電磁波の強度には制限があり、従来の測定装置１００ではグルコース値についての十分な精度を得ることが困難であった。 In the conventional measuring device 100, the living body as a subject is irradiated with light, and the irradiated light generates a sound wave which is a photoacoustic signal inside the living body. This sound wave is converted into an electric signal, signal analysis is performed, and the concentration of an in-vivo component such as a glucose value in a living body in which the background component and the target component are mixed is measured (see, for example, Patent Document 1). However, the interaction between glucose and electromagnetic waves is small, and the intensity of electromagnetic waves that can safely irradiate the living body is limited, and it has been difficult to obtain sufficient accuracy for the glucose value with the conventional measuring device 100.
そこで、生体内グルコース値等の測定値における精度をより向上させるために、図８に示す従来の光音響法を用いた生体内成分濃度の測定装置２００では、連続的に強度変調した光源が用いられる（例えば、特許文献１参照）。測定装置２００は、発振器２０１と、駆動回路２０３ａ、２０３ｂと、遅延調整器２０２と、第１光源２０４ａと、第２光源２０４ｂと、光合波器２０５と、音波検出器２０６と、波形観測器２０７と、記録器２０８とから構成される。本例では、二つの光源を用いており、第１光源２０４ａは、波長λ１の測定光を発生し、第２光源２０４ｂは、波長λ２の参照光を発生する。 Therefore, in order to further improve the accuracy of the measured values such as the in-vivo glucose value, the continuously intensity-modulated light source is used in the in-vivo component concentration measuring device 200 using the conventional photoacoustic method shown in FIG. (See, for example, Patent Document 1). The measuring device 200 includes an oscillator 201, drive circuits 203a and 203b, a delay adjuster 202, a first light source 204a, a second light source 204b, an optical combiner 205, a sound wave detector 206, and a waveform observer 207. And a recorder 208. In this example, two light sources are used, the first light source 204a generates the measurement light of the wavelength λ1, and the second light source 204b generates the reference light of the wavelength λ2.
発振器２０１は、第１光源２０４ａ及び第２光源２０４ｂから出力される光を強度変調するための変調信号を出力する。遅延調整器２０２は、発振器２０１からの変調信号のうち一方を反転して出力する。駆動回路２０３ａは第１光源２０４ａを駆動させる。駆動回路２０３ｂは、遅延調整器２０２で反転された変調信号を基に第２光源２０４ｂを駆動させる。 The oscillator 201 outputs a modulation signal for intensity-modulating the light output from the first light source 204a and the second light source 204b. The delay regulator 202 inverts and outputs one of the modulated signals from the oscillator 201. The drive circuit 203a drives the first light source 204a. The drive circuit 203b drives the second light source 204b based on the modulated signal inverted by the delay regulator 202.
第１光源２０４ａは、駆動回路２０３ａからの信号により波長λ１の測定光を強度変調して出力する。第２光源２０４ｂは、駆動回路２０３ｂからの信号により波長λ２の参照光を強度変調して出力する。これにより、異なる２波長λ１及びλ２光のそれぞれを、同一周波数で逆位相の信号により電気的に強度変調して出力する。 The first light source 204a intensity-modulates the measurement light having the wavelength λ1 according to the signal from the drive circuit 203a and outputs it. The second light source 204b intensity-modulates the reference light having the wavelength λ2 by the signal from the drive circuit 203b and outputs it. As a result, each of the two different wavelengths λ1 and λ2 light is electrically intensity-modulated by signals having opposite phases at the same frequency and output.
ここで、図１０にグルコース水溶液の吸光度スペクトルを示す。２つの波長λ１１及び波長λ２１は、対象成分の呈する吸収の差が、背景成分の呈する吸収の差よりも大きい波長である。また、波長λ１１は、対象成分が特徴的な吸収を呈する波長に設定する。波長λ１１及び波長λ２１は、対象成分の呈する吸収の差がそれ以外の成分の呈する吸収の差よりも大きい２波長であってもよい。これにより、水や測定対象の生体内成分以外の成分による吸収の影響を少なくして測定装置の測定精度をより向上させることができる。 Here, FIG. 10 shows the absorbance spectrum of the aqueous glucose solution. The two wavelengths λ11 and λ21 are wavelengths in which the difference in absorption exhibited by the target component is larger than the difference in absorption exhibited by the background component. Further, the wavelength λ11 is set to a wavelength at which the target component exhibits characteristic absorption. The wavelength λ11 and the wavelength λ21 may be two wavelengths in which the difference in absorption exhibited by the target component is larger than the difference in absorption exhibited by the other components. As a result, the influence of absorption by components other than water and components in the living body to be measured can be reduced, and the measurement accuracy of the measuring device can be further improved.
また、図１１にグルコース差分吸収スペクトルを示す。例えば、波長１６００ｎｍおよび２１００ｎｍにグルコースの吸収ピークが観測できる。したがって、波長λ１１はこれらの波長に設定する。 In addition, FIG. 11 shows a glucose differential absorption spectrum. For example, glucose absorption peaks can be observed at wavelengths of 1600 nm and 2100 nm. Therefore, the wavelength λ11 is set to these wavelengths.
図９に戻り、従来の生体内成分濃度の測定装置２００では、２つの波長λ１及び波長λ２の各々を電気的に強度変調する変調周波数を、被測定物で発生する音波の検出に関わる共鳴周波数と同一の周波数で変調することにより、音波の測定値における吸収係数に関わる非線形性に配慮して選択された２波長の光に対する音波を測定する。そして、２波長差分の音波から一方の波長の音波を規格化することで、一定に保ちがたい多数のパラメータの影響を排除して、より高精度に被測定物内に発生する音波を検出することができる。 Returning to FIG. 9, in the conventional in-vivo component concentration measuring device 200, the modulation frequency for electrically intensifying each of the two wavelengths λ1 and the wavelength λ2 is the resonance frequency related to the detection of the sound wave generated in the object to be measured. By modulating at the same frequency as, the sound wave for two wavelengths of light selected in consideration of the non-linearity related to the absorption coefficient in the measured value of the sound wave is measured. Then, by standardizing the sound wave of one wavelength from the sound wave of the difference of two wavelengths, the influence of many parameters that cannot be kept constant is eliminated, and the sound wave generated in the object to be measured is detected with higher accuracy. be able to.
２つの波長λ１及び波長λ２の光により被測定物の内部で発生した光音響信号は、それぞれ音響センサにより検出され、音圧に比例した電気信号に変換され、波形観測器２０９によって観測される。２波長に対応する光音響信号の強度の差は、血液中に含まれるグルコースの量に対応した電気信号として測定される。 The photoacoustic signal generated inside the object to be measured by the light of the two wavelengths λ1 and λ2 is detected by the acoustic sensor, converted into an electric signal proportional to the sound pressure, and observed by the waveform observer 209. The difference in the intensity of the photoacoustic signal corresponding to the two wavelengths is measured as an electric signal corresponding to the amount of glucose contained in the blood.
また、図１２に示すように、例えば、従来の測定装置２００における２つの選択波長では、それぞれの波長で水の吸光度スペクトルは温度依存係数が正負の傾向を示す。このように、吸光度変化の影響が大きい場合、温度変化による吸光度変化Δαｗに対して、温度を直接計測して既知の吸光度スペクトルから求める、若しくは、光の吸収量変化から温度変化ΔＴを推定して既知の吸光度スペクトルからΔαｗを求める。この場合、Δα＝Ｍαｇ＋Δαｗと記述できる。したがって、モル濃度Ｍ＝（Δα−Δαｗ）／αｇを演算し求めることができる。 Further, as shown in FIG. 12, for example, at two selected wavelengths in the conventional measuring apparatus 200, the absorbance spectrum of water at each wavelength tends to have a positive or negative temperature dependence coefficient. In this way, when the influence of the absorbance change is large, the temperature is directly measured and obtained from the known absorbance spectrum with respect to the absorbance change Δαw due to the temperature change, or the temperature change ΔT is estimated from the change in the amount of light absorbed. Obtain Δαw from a known absorbance spectrum. In this case, it can be described as Δα = Mαg + Δαw. Therefore, the molar concentration M = (Δα−Δαw) / αg can be calculated and obtained.
このように、従来の測定装置１００、２００では、光音響信号変化と生体内成分濃度との相間を予め測定することによって検量モデルをデータベースとして構築し、計測した光音響信号の変化から生体内成分濃度の検量を行っていた。 As described above, in the conventional measuring devices 100 and 200, a calibration model is constructed as a database by measuring the phase between the change in the photoacoustic signal and the concentration of the in-vivo component in advance, and the in-vivo component is obtained from the measured change in the photoacoustic signal. The concentration was calibrated.
また、図１３に、別の従来の非侵襲な生体内成分濃度の測定装置である、誘電分光法を用いた測定装置３００の構成例を示す。測定装置３００は、同軸プローブ３０１と、高周波信号計測装置３０２と、信号処理装置３０３と、濃度演算装置３０４と、データベース３０５とから構成される。 Further, FIG. 13 shows a configuration example of a measuring device 300 using dielectric spectroscopy, which is another conventional non-invasive measuring device for in vivo component concentration. The measuring device 300 includes a coaxial probe 301, a high frequency signal measuring device 302, a signal processing device 303, a concentration calculation device 304, and a database 305.
従来の誘電分光法による生体内成分濃度の測定装置３００は、同軸プローブを介して電磁波を皮膚内に照射し、測定対象の血液成分、例えば、グルコース分子と水の相互作用に従い、電磁波を吸収させ、電磁波の周波数に対する振幅及び位相を観測する（例えば、非特許文献２参照。）。 The conventional device 300 for measuring the concentration of components in a living body by dielectric spectroscopy irradiates the skin with electromagnetic waves via a coaxial probe, and absorbs the electromagnetic waves according to the interaction between the blood component to be measured, for example, a glucose molecule and water. , Observe the amplitude and phase of the electromagnetic wave with respect to the frequency (see, for example, Non-Patent Document 2).
観測される電磁波の周波数に対する振幅及び位相から、図１４Ａ及び図１４Ｂに示す誘電緩和スペクトル（濃度０〜１．４６２Ｍ）を算定する。一般的には、Ｃｏｌｅ−Ｃｏｌｅ式に基づき緩和カーブの線形結合として表現し、複素誘電率を算定する。 The dielectric relaxation spectrum (concentration 0 to 1.462M) shown in FIGS. 14A and 14B is calculated from the amplitude and phase of the observed electromagnetic wave with respect to the frequency. Generally, it is expressed as a linear combination of relaxation curves based on the Core-Cole equation, and the complex permittivity is calculated.
血液成分の定量では、例えば血液中に含まれるグルコースやコレステロール等の血液成分の量に複素誘電率は相間があるため、その変化に対応した電気信号（振幅、位相）として測定される。したがって、従来の測定装置３００では、複素誘電率変化と生体内成分濃度との相間を予め測定することによって検量モデルをデータベースとして構築し、計測した誘電緩和スペクトルの変化から生体内成分濃度の検量を行っていた。 In the quantification of blood components, for example, since the complex permittivity is between phases in the amount of blood components such as glucose and cholesterol contained in blood, it is measured as an electric signal (amplitude, phase) corresponding to the change. Therefore, in the conventional measuring device 300, a calibration model is constructed as a database by measuring the phase between the change in the complex permittivity and the concentration of the in-vivo component in advance, and the in-vivo component concentration is calibrated from the changed change in the measured dielectric relaxation spectrum. I was going.
このように、従来の光音響法を用いた測定装置１００、２００や、誘電分光法を用いた測定装置３００では、測定対象に応じた検量モデルをデータベースとして構築し、センサのキャリブレーションを行っていた。 In this way, in the measuring devices 100 and 200 using the conventional photoacoustic method and the measuring device 300 using the dielectric spectroscopy, a calibration model corresponding to the measurement target is constructed as a database and the sensor is calibrated. It was.
しかし、実際の生体グルコース等の生体内成分の測定では、様々な要因が影響を与える。生体内グルコースの測定に影響を与える要因としては、例えば、測定環境（温湿度、風）、発汗（電解質）や精神状態、皮膚部位間の毛細血管密度のバラツキ、皮膚の凹凸や堅さ、測定時の接触圧と他パラメータとの相関、皮膚下の骨の有無がプローブ接触面積に与える影響、測定者依存性、皮膚滲出液の有無、機器の皮膚内到達深度等の特性の設定等が挙げられる。 However, in the actual measurement of biological components such as biological glucose, various factors have an influence. Factors that affect the measurement of in vivo glucose include, for example, measurement environment (temperature and humidity, wind), sweating (electrolyte) and mental state, variation in capillary density between skin parts, unevenness and hardness of skin, and measurement. Correlation between contact pressure at time and other parameters, influence of presence or absence of bone under the skin on probe contact area, measurer dependence, presence or absence of skin exudate, setting of characteristics such as depth of reach of the device in the skin, etc. Be done.
したがって、従来におけるデータベースとして予め構築された検量モデルに基づくセンサのキャリブレーションでは、変化する測定環境等の影響が十分に反映されず、正確なキャリブレーションをすることが困難な場合があり、測定される生体内成分濃度の定量精度が十分でないことがあった。 Therefore, the calibration of the sensor based on the calibration model constructed in advance as a conventional database does not sufficiently reflect the influence of the changing measurement environment, etc., and it may be difficult to perform accurate calibration. In some cases, the quantification accuracy of the concentration of components in the body was not sufficient.
本発明は、様々な状況下で生体内成分濃度の定量精度を向上させることができる測定装置及び方法を提供することを目的とする。 An object of the present invention is to provide a measuring device and a method capable of improving the quantification accuracy of a biological component concentration under various circumstances.
上述した課題を解決するために、本発明に係る測定装置において、被験者の生体内成分値に関する情報の時系列信号に基づいて複数の特徴量の情報を抽出する特徴量抽出回路と、前記特徴量抽出回路により抽出された前記特徴量に基づいて、前記特徴量ごとに前記生体内成分値の事前推定値を算出する事前推定回路と、前記事前推定回路により算出された前記事前推定値に含まれるノイズを濾した事後推定値を算出するカルマンフィルタと、前記カルマンフィルタで算出されたそれぞれの前記事後推定値を重み付け平均化して最終推定値を出力する統合処理回路と、を備えることを特徴とする。 In order to solve the above-mentioned problems, in the measuring device according to the present invention, a feature amount extraction circuit that extracts information on a plurality of feature amounts based on a time-series signal of information on the in-vivo component values of a subject, and the feature amount extraction circuit. Based on the feature amount extracted by the extraction circuit, the pre-estimation circuit that calculates the pre-estimated value of the in-vivo component value for each feature amount and the pre-estimated value calculated by the pre-estimation circuit are used. It is characterized by including a Kalman filter that calculates a post-estimated value obtained by filtering out the included noise, and an integrated processing circuit that outputs the final estimated value by weighting and averaging each of the post-estimated values calculated by the Kalman filter. To do.
また、本発明に係る測定装置において、前記生体内成分値に関する前記情報の信頼度を示すエラー係数を測定するエラー測定回路と、前記エラー係数に基づいて重み付け行列の成分である重み付け定数を求める重み付け定数生成回路と、前記カルマンフィルタの推定誤差が最小となるように前記カルマンフィルタの共分散推定誤差行列を更新する適応推定回路と、前記重み付け行列と前記共分散推定誤差行列とに基づいて前記カルマンフィルタのカルマン利得を更新する重み付け演算回路と、をさらに備える。 Further, in the measuring device according to the present invention, an error measuring circuit for measuring an error coefficient indicating the reliability of the information regarding the in-vivo component value and a weighting for obtaining a weighting constant which is a component of the weighting matrix based on the error coefficient. The Kalman of the Kalman filter based on the constant generation circuit, the adaptive estimation circuit that updates the covariance estimation error matrix of the Kalman filter so that the estimation error of the Kalman filter is minimized, and the weighting matrix and the covariance estimation error matrix. a weighting calculation circuit for updating the gain, Ru further comprising a.
また、本発明に係る測定装置において、前記生体内成分値は生体内グルコース値であってもよい。 Further, in the measuring device according to the present invention, the in-vivo component value may be an in-vivo glucose value.
また、本発明に係る測定装置において、前記生体内成分値に関する前記情報の時系列信号は、光音響センサにより検出される光音響信号と、誘電分光センサにより検出される透過信号又は反射信号であってもよい。 Further, in the measuring device according to the present invention, the time-series signal of the information regarding the in-vivo component value is a photoacoustic signal detected by a photoacoustic sensor and a transmission signal or a reflection signal detected by a dielectric spectroscopy sensor. You may.
また、本発明に係る測定装置において、前記特徴量抽出回路により抽出される前記特徴量は、前記光音響信号の位相と振幅に基づいて算定される音速と、光吸収係数と、前記透過信号又は前記反射信号に基づいて算出される複素誘電率と、を含んでもよい。 Further, in the measuring device according to the present invention, the feature amount extracted by the feature amount extraction circuit is a sound velocity calculated based on the phase and amplitude of the photoacoustic signal, a light absorption coefficient, and the transmitted signal or. It may include a complex dielectric constant calculated based on the reflected signal.
また、本発明に係る測定方法は、被験者の生体内成分値に関する情報の時系列信号に基づいて複数の特徴量の情報を抽出する特徴量抽出ステップと、前記特徴量抽出ステップで抽出された前記特徴量に基づいて、前記特徴量ごとに前記生体内成分値の事前推定値を算出する事前推定ステップと、前記事前推定ステップで算出された前記事前推定値に含まれるノイズをカルマンフィルタにより濾して事後推定値を算出するフィルタリングステップと、前記フィルタリングステップで算出されたそれぞれの前記事後推定値を重み付け平均化して最終推定値を出力する統合処理ステップと、を備えることを特徴とする。 Further, the measurement method according to the present invention includes a feature amount extraction step of extracting information on a plurality of feature amounts based on a time-series signal of information on the in-vivo component values of the subject, and the feature amount extraction step extracted by the feature amount extraction step. Based on the feature amount, the pre-estimation step of calculating the pre-estimated value of the in-vivo component value for each feature amount and the noise contained in the pre-estimated value calculated in the pre-estimation step are filtered by a Kalman filter. It is characterized by including a filtering step for calculating the ex-post estimated value and an integrated processing step for weighting and averaging each ex-post estimated value calculated in the filtering step and outputting the final estimated value.
また、本発明に係る測定方法において、前記生体内成分値は生体内グルコース値であってもよい。 Further, in the measuring method according to the present invention, the in-vivo component value may be an in-vivo glucose value.
本発明によれば、被験者の生体内成分値に関するデータから複数の特徴量を抽出し、それぞれの特徴量から生体内成分値の事前推定をし、事前推定値各々についてカルマンフィルタでノイズを濾した事後推定値を算出し、複数の事後推定値を重み付け平均化するため、様々な状況下で生体内成分濃度の定量精度を向上させることができる。 According to the present invention, a plurality of feature amounts are extracted from the data on the in-vivo component values of the subject, the in-vivo component values are pre-estimated from each feature amount, and the noise is filtered by the Kalman filter for each of the pre-estimated values. Since the estimated values are calculated and a plurality of ex-post estimated values are weighted and averaged, the quantification accuracy of the in vivo component concentration can be improved under various situations.
以下、本発明の好適な実施の形態について、図１から図７を参照して詳細に説明する。また、以下の実施の形態では、本発明に係る測定装置を血糖自己測定（ＳＭＢＧ：Ｓｅｌｆ−ｍｏｎｉｔｏｒｉｎｇ ｏｆ ｂｌｏｏｄ ｇｌｕｃｏｓｅ）装置に適用し、生体内グルコース値を測定する場合について説明する。また、各図について共通する部分には、同一の符号が付されている。 Hereinafter, preferred embodiments of the present invention will be described in detail with reference to FIGS. 1 to 7. Further, in the following embodiment, a case where the measuring device according to the present invention is applied to a self-monitoring blood glucose (SMBG: Self-monitoring of blood glucose) device to measure an in vivo glucose level will be described. Further, the same reference numerals are given to the common parts of each figure.
＜実施の形態＞
本発明の実施の形態に係る測定装置は、測定対象である生体内グルコースに関する音響信号の時系列信号と、透過信号又は反射信号の時系列信号のそれぞれに基づいて複数の特徴量の情報を抽出する信号抽出部と、抽出された特徴量に基づいて、特徴量ごとに生体内グルコースの事前推定値を算出する生体内グルコース推定部と、生体内グルコース推定部により算出された各事前推定値に含まれるノイズを濾した事後推定値を算出するカルマンフィルタと、各事後推定値を統合して生体内グルコースの最終推定値を出力する統合処理部とを有する。
<Embodiment>
The measuring device according to the embodiment of the present invention extracts information on a plurality of feature quantities based on each of a time-series signal of an acoustic signal relating to glucose in the living body to be measured and a time-series signal of a transmitted signal or a reflected signal. To the signal extraction unit, the in-vivo glucose estimation unit that calculates the pre-estimated value of in-vivo glucose for each feature amount based on the extracted feature amount, and each pre-estimated value calculated by the in-vivo glucose estimation unit. It has a Kalman filter that calculates a post-estimated value by filtering out the contained noise, and an integrated processing unit that integrates each post-estimated value and outputs a final estimated value of glucose in vivo.
図１は、本発明の実施の形態に係る測定装置１の構成例を示すブロック図である。測定装置１は、光音響センサ２と、誘電分光センサ３と、信号抽出部４と、生体内グルコース推定部５と、カルマンフィルタ６と、統合処理部７と、エラー測定部８と、重み付け定数生成部９と、適応推定部１０と、重み付け演算部１１とを含んで構成される。 FIG. 1 is a block diagram showing a configuration example of the measuring device 1 according to the embodiment of the present invention. The measuring device 1 includes a photoacoustic sensor 2, a dielectric spectroscopy sensor 3, a signal extraction unit 4, an in vivo glucose estimation unit 5, a Kalman filter 6, an integrated processing unit 7, an error measurement unit 8, and a weighting constant generation. A unit 9, an adaptive estimation unit 10, and a weighting calculation unit 11 are included.
測定装置１による生体内グルコース値の測定は、生体内グルコース値の推定値を求めることにより行われる。測定装置１によって測定される生体内グルコース値は、後述する生体内グルコースの事前推定値のノイズが濾された事後推定値３をさらに統合した最終推定値ということになる。また、本実施の形態において、生体内グルコース値とは、グルコースの血中濃度を意味する。 The measurement of the in-vivo glucose level by the measuring device 1 is performed by obtaining an estimated value of the in-vivo glucose level. The in-vivo glucose value measured by the measuring device 1 is a final estimated value in which the noise of the in-vivo glucose pre-estimated value, which will be described later, is filtered out and the post-estimated value 3 is further integrated. Further, in the present embodiment, the in vivo glucose level means the blood concentration of glucose.
光音響センサ２は、生体内グルコースが光を吸収することにより生ずる光音響信号を検出し、光音響の時系列信号列を出力する。誘電分光センサ３は、生体内グルコースを透過した電磁波の透過又は生体内グルコースを反射した電磁波を検出し、透過又は反射信号の時系列信号列を出力する。 The photoacoustic sensor 2 detects a photoacoustic signal generated by absorption of light by glucose in the living body, and outputs a time-series signal sequence of photoacoustic. The dielectric spectroscopy sensor 3 detects the transmission of an electromagnetic wave transmitted through the in-vivo glucose or the electromagnetic wave reflected by the in-vivo glucose, and outputs a time-series signal sequence of the transmitted or reflected signal.
信号抽出部４は、光音響センサ２により検出される光音響信号の位相と振幅の平均から変化量を規定し、特徴量である実効音速、及び実効光吸収係数を算定する。また、信号抽出部４は、誘電分光センサ３により検出される透過信号又は反射信号のスペクトルにおける位相と振幅により、特徴量である実効複素誘電率を算定する。信号スペクトルに関しては、信号抽出部４は、Ｓａｖｉｔｓｋｙ−Ｇｏｌａｙフィルタ処理や微分処理等を用いて、スペクトル変動量を求める。 The signal extraction unit 4 defines the amount of change from the average of the phase and amplitude of the photoacoustic signal detected by the photoacoustic sensor 2, and calculates the effective sound velocity and the effective light absorption coefficient, which are characteristic quantities. Further, the signal extraction unit 4 calculates the effective complex permittivity, which is a feature quantity, from the phase and amplitude in the spectrum of the transmitted signal or the reflected signal detected by the dielectric spectroscopy sensor 3. Regarding the signal spectrum, the signal extraction unit 4 obtains the amount of spectrum fluctuation by using Savitzky-Goray filter processing, differential processing, or the like.
また、信号抽出部４は、実効音速、実効光吸収係数、及び実効複素誘電率の各時系列信号について、後述のダウンサンプリングを行うサンプリング回路、さらに、これらの信号について帯域制限するバンドパスフィルタ回路として機能する。 Further, the signal extraction unit 4 includes a sampling circuit that performs downsampling of each time-series signal of effective sound velocity, effective light absorption coefficient, and effective complex permittivity, which will be described later, and a bandpass filter circuit that limits the band of these signals. Functions as.
生体内グルコース推定部５は、信号抽出部４によりダウンサンプリング及び帯域制限された実効音速、実効光吸収係数、及び実効複素誘電率の信号について、周波数解析により生体内グルコースの事前推定値を算出する。 The in-vivo glucose estimation unit 5 calculates a pre-estimated value of in-vivo glucose by frequency analysis for signals of effective sound velocity, effective light absorption coefficient, and effective permittivity downsampled and band-limited by the signal extraction unit 4. ..
カルマンフィルタ６は、生体内グルコース推定部５により算出された実効音速、実効光吸収係数、及び実効複素誘電率の信号に基づいて推定された生体内グルコースの事前推定値各々についてノイズを濾した生体内グルコースの事後推定値を算出する。 The Kalman filter 6 filters out noise for each of the pre-estimated values of in-vivo glucose estimated based on the signals of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity calculated by the in-vivo glucose estimation unit 5. Calculate the ex-post estimate of glucose.
統合処理部７は、カルマンフィルタ６で算出された生体内グルコースの事後推定値を重み付け平均化する、重み付け定数生成処理部７０１と、重み付け平均化成処理部７０２とを有し、データを統合する。 The integrated processing unit 7 has a weighted constant generation processing unit 701 that weights and averages the ex post facto estimated value of glucose in the living body calculated by the Kalman filter 6, and a weighted averaging processing unit 702, and integrates the data.
エラー測定部８は、実効音速、実効光吸収係数、及び実効複素誘電率の各々の信頼度を示すエラー係数を測定する。 The error measuring unit 8 measures an error coefficient indicating the reliability of each of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity.
重み付け定数生成部９は、エラー測定部８により測定されたエラー係数に基づいて重み付け行列の成分である重み付け定数を求める。 The weighting constant generation unit 9 obtains a weighting constant which is a component of the weighting matrix based on the error coefficient measured by the error measurement unit 8.
適応推定部１０は、カルマンフィルタ６の推定誤差が最小となるようにカルマンフィルタ６の共分散推定誤差行列を更新する。 The adaptive estimation unit 10 updates the covariance estimation error matrix of the Kalman filter 6 so that the estimation error of the Kalman filter 6 is minimized.
重み付け演算部１１は、重み付け行列と共分散推定誤差行列に基づいてカルマン利得を更新する。 The weighting calculation unit 11 updates the Kalman gain based on the weighting matrix and the covariance estimation error matrix.
図２は、本実施の形態に係る測定装置１の動作を説明するフローチャートである。まず、信号抽出部４は、光音響センサ２と誘電分光センサ３各々について、事前のキャリブレーションを行う（ステップＳ１）。具体的には、信号抽出部４は、初期パラメータとして測定値と温度（測定物温度、外気温）を用い、実測グルコース値から検量直線を生成する。 FIG. 2 is a flowchart illustrating the operation of the measuring device 1 according to the present embodiment. First, the signal extraction unit 4 pre-calibrates each of the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3 (step S1). Specifically, the signal extraction unit 4 uses the measured value and the temperature (measured object temperature, outside air temperature) as initial parameters, and generates a calibration straight line from the measured glucose value.
さらに、信号抽出部４は、実測グルコース値に基づいて生成された検量直線から生体内グルコース値を算定し、記憶部（図示しない）に記憶する。このとき、信号抽出部４は、光音響センサ２と誘電分光センサ３による観測値、及び温度（測定物温度、外気温）を用いる。これにより、光音響センサ２と誘電分光センサ３の事前のキャリブレーションが行われる。 Further, the signal extraction unit 4 calculates the in-vivo glucose value from the calibration straight line generated based on the actually measured glucose value, and stores it in a storage unit (not shown). At this time, the signal extraction unit 4 uses the observed values by the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3 and the temperature (measured object temperature, outside air temperature). As a result, the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3 are pre-calibrated.
次に、信号抽出部４は、光音響センサ２により検出された光音響信号の位相と振幅を検出する（ステップＳ２）。また、信号抽出部４は、誘電分光センサ３により検出された透過信号又は反射信号スペクトルの位相と振幅を検出する（ステップＳ３）。 Next, the signal extraction unit 4 detects the phase and amplitude of the photoacoustic signal detected by the photoacoustic sensor 2 (step S2). Further, the signal extraction unit 4 detects the phase and amplitude of the transmitted signal or reflected signal spectrum detected by the dielectric spectroscopy sensor 3 (step S3).
そして、信号抽出部４は、ステップＳ２で検出した光音響信号の位相と振幅の平均から変化量を規定した後に、実効音速を算定し（ステップＳ４）、さらに、実効光吸収係数を算定する（ステップＳ５）。また、信号抽出部４は、ステップＳ３で検出した透過信号又は反射信号スペクトルの位相と振幅とにより、実効複素誘電率を算定する（ステップＳ６）。本実施の形態では、信号抽出部４により算定された実効音速、実効光吸収係数、及び実効複素誘電率に基づいて、最終的な生体内グルコース値（最終推定値）が推定される。 Then, the signal extraction unit 4 determines the amount of change from the average of the phase and amplitude of the photoacoustic signal detected in step S2, then calculates the effective sound velocity (step S4), and further calculates the effective light absorption coefficient (step S4). Step S5). Further, the signal extraction unit 4 calculates the effective complex permittivity based on the phase and amplitude of the transmitted signal or reflected signal spectrum detected in step S3 (step S6). In the present embodiment, the final in-vivo glucose value (final estimated value) is estimated based on the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity calculated by the signal extraction unit 4.
次に、信号抽出部４は、ステップＳ４〜Ｓ６でそれぞれ算定した実効音速、実効光吸収係数、および実効複素誘電率の各時系列信号についてダウンサンプリングを行う（ステップＳ７）。具体的には、信号抽出部４は、例えば、０．３Ｈｚ間隔で、実効音速、実効光吸収係数、および実効複素誘電率の時系列信号のサンプリング周波数を下げる。さらに、実効複素誘電率の時系列信号については、信号抽出部４は、０．３Ｈｚ間隔へのダウンサンプリングを行った後に、１分間ステップサイズで５分間の時間領域について平均化する。 Next, the signal extraction unit 4 downsamples each time-series signal of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity calculated in steps S4 to S6, respectively (step S7). Specifically, the signal extraction unit 4 lowers the sampling frequency of the time-series signal of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity, for example, at intervals of 0.3 Hz. Further, for the time-series signal having an effective complex permittivity, the signal extraction unit 4 performs downsampling to 0.3 Hz intervals and then averages the time domain for 5 minutes with a 1-minute step size.
次に、信号抽出部４は、ステップＳ７でダウンサンプリングした実効音速、実効光吸収係数、および実効複素誘電率の時系列信号の各々を帯域制限する（ステップＳ８）。生体信号である実効音速、実効光吸収係数、および実効複素誘電率は、低周波のみに限られるため、バンドパスフィルタを用いる。このバンドパスフィルタの典型的な通過帯域は、例えば、０．１５〜０．４Ｈｚである。 Next, the signal extraction unit 4 band-limits each of the time-series signals of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity downsampled in step S7 (step S8). Since the effective sound velocity, effective light absorption coefficient, and effective complex permittivity, which are biological signals, are limited to low frequencies, a bandpass filter is used. The typical pass band of this bandpass filter is, for example, 0.15 to 0.4 Hz.
なお、実効音速の時系列信号、実効光吸収係数の時系列信号、実効複素誘電率の時系列信号に対するバンドパスフィルタの通過帯域を異なる値に設定してもよい。例えば、光音響センサ２からの実効音速と実効光吸収係数の各時系列信号の通過帯域の値は０．１５〜０．４Ｈｚ、誘電分光センサ３からの実効複素誘電率の時系列信号の通過帯域の値は０．１〜０．５Ｈｚとしてもよい。 The pass band of the bandpass filter for the time-series signal of the effective sound velocity, the time-series signal of the effective light absorption coefficient, and the time-series signal of the effective complex permittivity may be set to different values. For example, the pass band values of the effective sound velocity and the effective light absorption coefficient of each time-series signal from the photoacoustic sensor 2 are 0.15 to 0.4 Hz, and the pass of the time-series signal of the effective complex permittivity from the dielectric spectroscopy sensor 3. The band value may be 0.1 to 0.5 Hz.
次に、生体内グルコース推定部５は、ステップＳ７〜Ｓ８で信号抽出部４によりダウンサンプリングおよび帯域制限された実効音速の時系列信号、実効光吸収係数の時系列信号、実効複素誘電率の時系列信号から生体内グルコースの事前推定値を各々算出する（ステップＳ９）。より詳細には、生体内グルコース推定部５は、実効音速対生体内グルコース値の関係、実効光吸収係数対生体内グルコース値の関係、および実効複素誘電率対生体内グルコース値の関係を求める。その後に、生体内グルコース推定部５は、温度と特徴量（実効音速、実効光吸収係数、および実効複素誘電率）から、特徴量ごとの生体内グルコースの事前推定値を算出する。 Next, the in-vivo glucose estimation unit 5 performs a time-series signal with an effective sound velocity downsampled and band-limited by the signal extraction unit 4 in steps S7 to S8, a time-series signal with an effective light absorption coefficient, and an effective complex permittivity. Pre-estimated values of in vivo glucose are calculated from the series signals (step S9). More specifically, the in-vivo glucose estimation unit 5 obtains the relationship between the effective sound velocity and the in-vivo glucose value, the relationship between the effective light absorption coefficient and the in-vivo glucose value, and the relationship between the effective complex permittivity and the in-vivo glucose value. After that, the in-vivo glucose estimation unit 5 calculates a pre-estimated value of in-vivo glucose for each feature from the temperature and the feature (effective sound velocity, effective light absorption coefficient, and effective complex permittivity).
ここで、本実施の形態において特徴量として用いる各物理パラメータとグルコース濃度変化の関係を表１に示す。 Here, Table 1 shows the relationship between each physical parameter used as a feature amount in the present embodiment and the change in glucose concentration.
表１によれば、吸光度α、音速ｖ、および誘電率実部εは、それぞれグルコース濃度変化と相関がある。また、これらの物理パラメータは温度変化とも相関を有することが分かる。よって、ステップＳ４〜Ｓ６で信号抽出部４により得られた特徴量であるこれらの物理パラメータは、例えば、行列Ｈを用いて次式で表される。 According to Table 1, the absorbance α, the speed of sound v, and the real dielectric constant ε are each correlated with the change in glucose concentration. It can also be seen that these physical parameters also correlate with temperature changes. Therefore, these physical parameters, which are the feature quantities obtained by the signal extraction unit 4 in steps S4 to S6, are represented by, for example, the following equation using the matrix H.
ここで、ｖ（ｔ）は音速、α（ｔ）は吸光度、ｅ（ｔ）は誘電率、Ｇ（ｔ）は生体内グルコース値、Ｔ（ｔ）は測定部体温、Ｒ（ｔ）は残差、ｔは時間である。さらに、上式（１）〜（３）から逆行列Ｗを求めることにより、生体内グルコース値Ｇ（ｔ）、測定部体温Ｔ（ｔ）、残差Ｒ（ｔ）は次の推定モデル式で表される。 Here, v (t) is the speed of sound, α (t) is the absorbance, e (t) is the permittivity, G (t) is the in vivo glucose value, T (t) is the body temperature of the measuring part, and R (t) is the residual. The difference, t, is time. Furthermore, by obtaining the inverse matrix W from the above equations (1) to (3), the in vivo glucose value G (t), the body temperature of the measuring part T (t), and the residual R (t) are obtained by the following estimation model equations. expressed.
残差Ｒを最も小さくするように行列Ｈの係数を調整してもよい。また、周波数情報を利用して、それぞれの物理パラメータから生体内グルコース値Ｇを求めてもよい。すなわち、ある波長λにおける吸光度α（ｔ、λ）より、式（１）を次式のように表すことができる。 The coefficients of the matrix H may be adjusted to minimize the residual R. Further, the in-vivo glucose value G may be obtained from each physical parameter by using the frequency information. That is, the equation (1) can be expressed as the following equation from the absorbance α (t, λ) at a certain wavelength λ.
本実施の形態では、従来の主成分分析やＰＬＳ分析に代表される多変量解析により、予め求めた検量モデルから生体内グルコース値Ｇ（ｔ）や測定部体温Ｔ（ｔ）を推定する。また、音速ｖ（ｔ）や誘電率ε（ｔ）に対しても、同様な方法により生体内グルコース値Ｇ（ｔ）や測定部体温Ｔ（ｔ）を推定することができる。 In the present embodiment, the in-vivo glucose value G (t) and the body temperature T (t) of the measurement part are estimated from a calibration model obtained in advance by multivariate analysis represented by conventional principal component analysis and PLS analysis. Further, for the sound velocity v (t) and the dielectric constant ε (t), the in-vivo glucose value G (t) and the body temperature T (t) of the measuring part can be estimated by the same method.
生体内グルコース推定部５は、光音響センサ２と誘電分光センサ３による情報である音速ｖ（ｔ）、吸光度α（ｔ）、および誘電率ε（ｔ）から推定した生体内グルコース値Ｇ（ｔ）の時系列データに対して、予め定めた生体内グルコース変動モデルに基づき生体内グルコース値Ｇ（ｔ）の推定を行う。生体内グルコース変動モデルは、グルコース推移曲線の理想モデルとして、人体の生理的なグルコース値変化を数理モデルで仮定する。 The in-vivo glucose estimation unit 5 estimates the in-vivo glucose value G (t) from the sound velocity v (t), the absorbance α (t), and the permittivity ε (t), which are information obtained by the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3. ), The in-vivo glucose value G (t) is estimated based on a predetermined in-vivo glucose fluctuation model. The in-vivo glucose fluctuation model assumes the physiological glucose level change of the human body as an ideal model of the glucose transition curve by a mathematical model.
生体内グルコース値の時間推移は様々なモデル化手法によりモデル化できることが知られている。例えば、非線形自己回帰モデル、確率モデル等のパラメトリックモデルを用いることができる。本実施の形態では、糖代謝メタボリズムをシステム工学的に記述した生体内グルコース変動モデルを用いる（例えば、非特許文献３参照）。このモデルにおいて、インシュリンと独立なグルコースＧ値の推移は次式のように時刻ｔに対する１次微分方程式により表される。 It is known that the time course of glucose level in vivo can be modeled by various modeling methods. For example, a parametric model such as a nonlinear autoregressive model or a stochastic model can be used. In this embodiment, an in vivo glucose fluctuation model in which glucose metabolism metabolism is described in a system engineering manner is used (see, for example, Non-Patent Document 3). In this model, the transition of the glucose G value independent of insulin is expressed by a first-order differential equation with respect to time t as shown in the following equation.
ここで、ｐ１、ｐ２は定数である。ｐ１＝−ｋ１、ｐ２＝Ｐとなり、ｋ１はグルコースＧの消費速度を代表する定数、Ｐは肝臓におけるグルコース生成に関する定数である。これらの定数は、初期的に採血等による生体内グルコース値により校正されて、実用に供される。 Here, p1 and p2 are constants. p1 = −k1 and p2 = P, where k1 is a constant representing the rate of glucose G consumption and P is a constant relating to glucose production in the liver. These constants are initially calibrated by the in vivo glucose level by blood sampling or the like and put into practical use.
このように定められる生体内グルコース変動モデルは、図３に示す、測定装置１による生体内グルコース値の推定を説明する模式図において、破線で示される。また、生体内グルコース値の推定開始前における「白丸」の点で示すキャリブレーション値は、ステップＳ１で事前のキャリブレーションが行われた光音響センサ２および誘電分光センサ３からの実測値に基づいて、生体内グルコース変動モデルのパラメータフィッティングが行われた、実際の生体内グルコース値を示している。 The in-vivo glucose fluctuation model defined in this way is shown by a broken line in the schematic diagram shown in FIG. 3 for explaining the estimation of the in-vivo glucose value by the measuring device 1. Further, the calibration value indicated by the “white circle” before the start of estimation of the in-vivo glucose value is based on the measured value from the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3 that were pre-calibrated in step S1. , The actual in-vivo glucose level is shown after the parameter fitting of the in-vivo glucose fluctuation model.
本実施の形態における測定装置１による生体内グルコース値の測定は、前述のとおり生体内グルコースの推定値を求めることにより行われる。図３に示すように、生体内グルコース値の推定開始後においては、生体内グルコース変動モデルによる推定値と、各センサから得られる特徴量である、実効音速、実効光吸収係数、および実効複素誘電率に基づいて特徴量ごとに算出される生体内グルコースの事前推定値の両方が用いられる。各センサによる生体内グルコースの事前推定値は、推定モデル式（４）より算出される。 The measurement of the in-vivo glucose value by the measuring device 1 in the present embodiment is performed by obtaining the estimated value of the in-vivo glucose as described above. As shown in FIG. 3, after the start of estimation of the in-vivo glucose value, the estimated value by the in-vivo glucose fluctuation model and the feature quantities obtained from each sensor, such as the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity. Both pre-estimated values of in vivo glucose calculated for each feature based on the rate are used. The pre-estimated value of glucose in the living body by each sensor is calculated from the estimation model formula (4).
図２に戻り、カルマンフィルタ６は、生体内グルコース推定部５が実効音速の時系列信号、実効光吸収係数の時系列信号、および実効複素誘電率の時系列信号から推定した生体内グルコースの事前推定値の各々についてノイズを濾して事後推定値を求める（ステップＳ１０）。本実施の形態では、時系列自己回帰（ＡＲ）モデルで適応カルマンフィルタを用いる場合について説明する。 Returning to FIG. 2, the Kalman filter 6 pre-estimates the in-vivo glucose estimated by the in-vivo glucose estimation unit 5 from the effective sound velocity time-series signal, the effective light absorption coefficient time-series signal, and the effective complex permittivity time-series signal. Noise is filtered for each of the values to obtain a post-estimation value (step S10). In this embodiment, a case where an adaptive Kalman filter is used in a time series autoregressive (AR) model will be described.
まず、本実施の形態では、カルマンフィルタ６は、生体におけるグルコース推移システムを次式のように定める。 First, in the present embodiment, the Kalman filter 6 defines the glucose transition system in the living body as follows.
式（１１）において、Ｇvは生体内グルコース推定部５が実効音速の時系列信号から推定した生体内グルコースの事前推定値、Ｇαは、生体内グルコース推定部５が実効光吸収係数の時系列信号から推定した生体内グルコースの事前推定値、およびＧεは生体内グルコース推定部５が実効複素誘電率の時系列信号から推定した生体内グルコース値である。また、上式（１２）において、ｇv、ｇα、ｇεはそれぞれＧv、Ｇα、Ｇεに対応するシステムに入力値である。 In equation (11), G v is a pre-estimated value of in-vivo glucose estimated by the in-vivo glucose estimation unit 5 from the time series signal of the effective sound velocity, and G α is when the in-vivo glucose estimation unit 5 has an effective light absorption coefficient. The in-vivo glucose pre-estimated value estimated from the series signal and G ε are the in-vivo glucose values estimated by the in-vivo glucose estimation unit 5 from the time-series signal of the effective complex dielectric constant. Further, in the above equation (12), g v , g α , and g ε are input values to the system corresponding to G v , G α , and G ε , respectively.
図４は、本実施の形態に係る測定装置１におけるカルマンフィルタのブロック線図である。ベイジアン推定と自動再帰推定に基づくカルマンフィルタはノイズが無相関で、ゼロ平均、白色雑音の際に、最適線形推定器であることが知られている。この条件が真でなくても、最適ではないが最良の線形推定器と考えられる。 FIG. 4 is a block diagram of the Kalman filter in the measuring device 1 according to the present embodiment. The Kalman filter based on Basilian estimation and automatic recursive estimation is known to be an optimal linear estimator for zero mean and white noise with no correlation of noise. Even if this condition is not true, it is considered to be the best linear estimator, though not optimal.
生体内グルコース推定部５が実効音速の時系列信号から推定した生体内グルコースの事前推定値、実効光吸収係数の時系列信号から推定した事前推定値、および、実効複素誘電率の時系列信号から推定した事前推定値は、それぞれカルマンフィルタ６に入力される。 From the in-vivo glucose estimation unit 5 estimated from the time-series signal of the effective sound velocity, the pre-estimated value estimated from the time-series signal of the effective light absorption coefficient, and the time-series signal of the effective complex dielectric constant. The estimated pre-estimated values are input to the Kalman filter 6, respectively.
また、カルマンフィルタ６は、式（１０）について、次のようにサンプル時間Ｔとした離散時間状態方程式を求める。 Further, the Kalman filter 6 obtains the discrete-time state equation with the sample time T as follows for the equation (10).
したがって、生体システムを記述する物理量Ｇkは再帰的に決定されており，次式で表される。 Therefore, the physical quantity G k that describes the biological system is recursively determined and is expressed by the following equation.
ここで，ｇkはシステム入力、Ｇkは雑音無しの物理量である。また、Ａはシステムモデルを示すｎ×ｍ行列（ｍは測定数、ｎは生体システムからの信号数）、Ｈは測定系モデルを示すｍ×ｎ行列である。物理量Ｇk+1は、前時刻における物理量ＡＧkとＣkとｗkの生体システム雑音が含まれる。システム入力ｇkには、測定値ＨＧkとｖkの測定系雑音が含まれる。システムの更新は下記式により行う。 Here, g k is a system input and G k is a physical quantity without noise. Further, A is an n × m matrix indicating a system model (m is the number of measurements, n is the number of signals from the biological system), and H is an m × n matrix indicating the measurement system model. The physical quantity G k + 1 includes the biological system noise of the physical quantities AG k , C k, and w k at the previous time. The system input g k includes measurement system noise of the measured values HG k and v k . The system is updated by the following formula.
ここで、ｄkはカルマンイノベーション、ハットＧkは生体内グルコースの事前推定値Ｇkについての事後推定値である。以下、同様に文字上に付した「∧」をハットと呼ぶ。ここで、Ｓkは従来のｎ×ｍ行列のカルマン利得とは異なる。変更カルマン利得Ｓkは次式より求める。 Here, d k is Kalman innovation, and hat G k is a post-estimation value for the pre-estimated value G k of in vivo glucose. Hereinafter, the “∧” attached to the letters is referred to as a hat. Here, S k is different from the Kalman gain of a conventional n × m matrix. Change Kalman gain S k is determined from the following equation.
ここで、Ｐk -はハットＧkの共分散推定誤差行列であり、ＨTは行列Ｈの転置行列であり、Ｄkは測定系雑音ｖkの共分散行列であるＲkに基づくパラメータ（行列）である。式（１７）によれば、既知の情報を用いて、共分散行列Ｒkを調整することによりパラメータＤkを求め、変更カルマン利得Ｓkを決定することができる。既知の情報は、例えば、外部からのインパルス雑音やセンサ固有の雑音であり、重み付け演算部１１からの定数および定常パラメータからなる。 Here, P k - is the covariance estimation error matrix hat G k, H T is the transpose of the matrix H, D k is based on R k is the covariance matrix of the measurement system noise v k parameters ( Matrix). According to the equation (17), the parameter D k can be obtained by adjusting the covariance matrix R k using the known information, and the modified Kalman gain S k can be determined. The known information is, for example, external impulse noise or sensor-specific noise, and includes constants and steady-state parameters from the weighting calculation unit 11.
パラメータＤkには従来のカルマンフィルタで用いることができない情報が含まれているので、共分散推定誤差行列Ｐk -を更新するのに用いるのは適切ではない。パラメータＤkの推定方法は後述するが、本実施の形態では、共分散推定誤差行列Ｐk -の更新には、共分散行列Ｒkを用いた標準的な算定法を採用する。従来のカルマン利得行列Ｋkは次のように定義される。 Since the parameter D k includes information that can not be used in a conventional Kalman filter covariance estimation error matrix P k - used to update is not appropriate. The method for estimating the parameter D k will be described later, but in the present embodiment, a standard calculation method using the covariance matrix R k is adopted for updating the covariance estimation error matrix P k − . The conventional Kalman gain matrix K k is defined as follows.
ここで、Ｐk +はＡハットＧkの共分散推定誤差行列であり、Ｑkは生体システム雑音ｗkの共分散行列、Ｉは単位行列、ＡTは行列Ａの転置行列である。式（１８）、式（２０）に用いる共分散行列Ｑk、Ｒkの値を適応的に求める手法は複数存在する。本実施の形態では、適応推定部１０において、次式を用い、共分散行列Ｒkの値を推定する。 Here, P k + is the covariance estimation error matrix of the A hat G k , Q k is the covariance matrix of the biological system noise w k , I is the identity matrix, and AT is the transposed matrix of the matrix A. There are a plurality of methods for adaptively obtaining the values of the covariance matrices Q k and R k used in the equations (18) and (20). In the present embodiment, the adaptive estimation unit 10 estimates the value of the covariance matrix R k by using the following equation.
ここで、ｍは推定に用いる窓長さである。式（２１）の方法では、カルマンイノベーションｄkが窓のなかで定常的であることが必要である。この前提は、非定常的な生体システムでは近似となる。また、カルマンフィルタが最適の際に、カルマンイノベーションｄkはゼロ平均、白色雑音となる。もしカルマンイノベーションｄkが非ゼロ平均で有色であるならば、共分散行列Ｒkを調整することで、白色化を試みる例がある。この際に、定常システムでは、最適な共分散行列Ｒkは最終的に一定値へ到達する。本実施の形態では、非定常的な生体システムを前提とするため、ある限られた窓長さｍで、定常状態となることを仮定する。本実施の形態では、同様の近似をカルマン利得等の他の定常パラメータにも適用する。 Here, m is the window length used for estimation. The method of equation (21) requires that the Kalman innovation d k be stationary in the window. This assumption is an approximation for non-stationary biological systems. Further, when the Kalman filter is optimal, the Kalman innovation d k becomes zero average and white noise. If the Kalman innovation d k is colored with a non-zero mean, there are examples of attempting whitening by adjusting the covariance matrix R k . At this time, in the stationary system, the optimum covariance matrix R k finally reaches a constant value. In this embodiment, since a non-stationary biological system is premised, it is assumed that a steady state is achieved with a certain limited window length m. In this embodiment, the same approximation is applied to other stationary parameters such as Kalman gain.
生体システム雑音ｗkの共分散行列Ｑkに関しては、共分散行列Ｒkと同様に推定手法が知られているが、定常状態が仮定できないため、誤差を含むことが考えられる。そこで、本実施の形態では、共分散行列Ｑkを、予め人工的に高い値に設定された成分を有する所定の行列とし、モデルの誤差が大きいものという前提を課すことで、新たな推定値の信頼性が古い推定値よりも高くなるようにする。 For the covariance matrix Q k of the biological system noise w k, but estimation method similar to the covariance matrix R k is known, because the steady state can not be assumed, it is conceivable to include an error. Therefore, in the present embodiment, the covariance matrix Q k is set to a predetermined matrix having a component artificially set to a high value in advance, and a new estimated value is imposed by imposing the premise that the error of the model is large. Make sure that the reliability of is higher than the old estimate.
次に、重み付け定数生成部９が実行する処理ついて説明する。前述したとおり、既知の情報をカルマンフィルタ６に注入する必要がある。本実施の形態では、光音響センサ２および誘電分光センサ３から得られた信号のエラーを測定することにより、重み付け行列Ｗkの成分である重み付け定数を求める。情報の信頼度に応じて、システムにフィードバックする。 Next, the process executed by the weighted constant generation unit 9 will be described. As mentioned above, it is necessary to inject known information into the Kalman filter 6. In the present embodiment, the weighting constant which is a component of the weighting matrix W k is obtained by measuring the error of the signal obtained from the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3. Feedback is given to the system according to the reliability of the information.
このようなフィードバックする定数（エラー係数）を−１から＋１の間で定める。情報の信頼度が低い場合、エラー係数を−１とし、情報の信頼度が高い場合、エラー係数を＋１とする。エラーが未知の場合は、エラー係数を０とする。エラー測定部８によるエラー係数の具体的な決定方法については後述する。 Such a feedback constant (error coefficient) is defined between -1 and +1. If the reliability of the information is low, the error coefficient is set to -1, and if the reliability of the information is high, the error coefficient is set to +1. If the error is unknown, the error coefficient is set to 0. The specific method for determining the error coefficient by the error measuring unit 8 will be described later.
エラー係数はカルマンフィルタ状態と関連する重み付け行列Ｗkの対角成分に配置される。例えば実効音速についてのエラー係数を重み付け行列Ｗkの対角成分Ｗ11として配置し、実効光吸収係数についてのエラー係数を重み付け行列Ｗkの対角成分Ｗ22として配置し、実効複素誘電率についてのエラー係数を重み付け行列Ｗkの対角成分Ｗ33として配置する。 The error coefficients are placed on the diagonal components of the weighting matrix W k associated with the Kalman filter state. For example arranged an error factor for the effective speed of sound as diagonal elements W 11 of the weighting matrix W k, place the error coefficients for effective optical absorption coefficient as a diagonal W 22 of the weighting matrix W k, the effective complex permittivity The error coefficient of is arranged as the diagonal component W 33 of the weighting matrix W k .
本実施の形態では、信号間の重み付け共分散も必要である。２つの信号は同様に誤差が小さい、誤差が大きい場合には似たようなエラー係数となるが、２つの信号が依存していることを意味していない。しかしながら、それらの信号の誤差の大小は、独立の雑音源が関与しているはずである。 In this embodiment, weighting covariance between signals is also required. The two signals have similarly small errors, and when the errors are large, they have similar error coefficients, but this does not mean that the two signals are dependent. However, the magnitude of the error of those signals should be related to the independent noise source.
本実施の形態では、エラー係数間の差を用いて、重み付け行列Ｗkの非対角成分Ｗij，Ｗji（ｉ，ｊは重み付け定数Ｗの次元を示すインデックスであり、ｉ≠ｊ）として配置される重み付け定数を次式により計算する。 In the present embodiment, the difference between the error coefficients is used as the off-diagonal components Wij and W ji of the weighting matrix W k (i and j are indexes indicating the dimension of the weighting constant W, and i ≠ j). The weighting constant to be placed is calculated by the following formula.
なお、時間インデックスを示すｋは式（２２）に記述していない。こうして、重み付け定数生成部９は、エラー測定部８から出力されるエラー係数に基づいて重み付け行列Ｗkを決定する。 Note that k indicating the time index is not described in the equation (22). In this way, the weighting constant generation unit 9 determines the weighting matrix W k based on the error coefficient output from the error measurement unit 8.
次に、式（１７）に用いるパラメータＤkを求めるためには、単純に重み付け行列Ｗkを用いて、式（２３）の計算を行う。 Next, in order to obtain the parameter D k used in the equation (17), the calculation of the equation (23) is performed simply by using the weighting matrix W k .
式（２３）における記号「○」はアダマール積を意味する。αはエラー係数の強さを調整するスケールファクタである（０≦α≦１）。αは収束性を高めるために、例えば、０．５と設定しておく。 The symbol "○" in equation (23) means the Hadamard product. α is a scale factor that adjusts the strength of the error coefficient (0 ≦ α ≦ 1). α is set to 0.5, for example, in order to improve the convergence.
図５は重み付け定数生成部９と適応推定部１０と重み付け演算部１１の動作を説明するフローチャートである。上述したように、ある窓長さｍ（測定数）におけるエラー測定で得られたエラー係数に基づいて、重み付け定数生成部９は、短期的なシステム推定を行い、重み付け行列Ｗkを決定する（ステップＳ１０２）。 FIG. 5 is a flowchart illustrating the operations of the weighting constant generation unit 9, the adaptation estimation unit 10, and the weighting calculation unit 11. As described above, the weighting constant generation unit 9 performs short-term system estimation and determines the weighting matrix W k based on the error coefficient obtained by the error measurement at a certain window length m (measurement number). Step S102).
適応推定部１０は、カルマンフィルタ状態と過去の状態から、理想的な雑音状態（カルマンイノベーションが定常状態）を仮定して、周知の手法により長期的なシステム推定を行い、共分散行列Ｒkを推定し（ステップＳ１０３）、さらに式（１８）〜式（２０）により定常パラメータである共分散推定誤差行列Ｐk -を決定する（ステップＳ１０４）。式（１８）〜式（２０）はハットＧkの推定誤差が最小となるようにカルマン利得Ｋkを再帰的に決定する周知の手法を表しており、この手法により共分散推定誤差行列Ｐk -，Ｐk +を更新することができる。 The adaptive estimation unit 10 assumes an ideal noise state (Kalman innovation is a steady state) from the Kalman filter state and the past state, performs long-term system estimation by a well-known method, and estimates the covariance matrix R k . (step S103), further covariance estimation error matrix P k which is a constant parameter by equation (18) to (20) - determining (step S104). Equations (18) to (20) represent a well-known method for recursively determining the Kalman gain K k so that the estimation error of the hat G k is minimized, and the covariance estimation error matrix P k is calculated by this method. - , P k + can be updated.
そして、重み付け演算部１１は、重み付け定数生成部９が決定した重み付け行列Ｗkを用いて式（２３）によりパラメータＤkを算出し、このパラメータＤkと適応推定部１０が決定した共分散推定誤差行列Ｐk -とを用いて式（１７）により変更カルマン利得Ｓkを算出する（ステップＳ１０５）。以上により、カルマン利得Ｓkを適応的に変更することができる。 Then, the weighting calculation unit 11 calculates the parameter D k by the equation (23) using the weighting matrix W k determined by the weighting constant generation unit 9, and the covariance estimation determined by this parameter D k and the adaptive estimation unit 10. The modified Kalman gain S k is calculated by Eq. (17) using the error matrix P k − (step S105). Thus, it is possible to change the Kalman gain S k adaptively.
本実施の形態は、実効音速の時系列信号から推定した生体内グルコース値の事後推定値と、実効光吸収係数の時系列信号から推定した生体内グルコース値の事後推定値と、実効複素誘電率の時系列信号から推定した生体内グルコース値の事後推定値とを統合することを特徴としている。すなわち、パラメータＤk、変更カルマン利得Ｓkおよび下式（２４）の行列Ａに基づいて、３つの生体内グルコースの事後推定値は式（１３’）により１つの最終推定値に結合される。なお、式（２５）は測定系モデルを表す。 In this embodiment, an ex post facto estimate of the in vivo glucose value estimated from the time series signal of the effective sound velocity, an ex post facto estimate of the in vivo glucose value estimated from the time series signal of the effective light absorption coefficient, and an effective complex dielectric constant It is characterized by integrating the ex-post estimate of the in-vivo glucose value estimated from the time series signal of. That is, based on the parameter D k , the modified Kalman gain Sk, and the matrix A of equation (24) below, the three ex post facto estimates of glucose in vivo are combined into one final estimate by equation (13'). Equation (25) represents a measurement system model.
本実施の形態では、上述のとおり、共分散行列Ｒkについては適応推定部１０で推定し、一方、共分散行列Ｑkについては推定せずに固定値とするため、共分散行列Ｑkをあらかじめ定めておく必要がある。システムを単純化するために、プロセス雑音は独立と仮定する。すなわち、共分散行列Ｑkの対角成分のみが値を有し、非対角成分が０となるようにする。 In the present embodiment, as described above, the covariance matrix R k is estimated by the adaptive estimation unit 10, while the covariance matrix Q k is not estimated and is set to a fixed value. Therefore, the covariance matrix Q k is used. It is necessary to decide in advance. To simplify the system, process noise is assumed to be independent. That is, only the diagonal component of the covariance matrix Q k has a value, and the off-diagonal component is 0.
例えば実効音速に対応する成分を共分散行列Ｑkの対角成分Ｑ11として配置し、実効光吸収係数に対応する成分を共分散行列Ｑkの対角成分Ｑ22として配置し、実効複素誘電率に対応する成分を共分散行列Ｑkの対角成分Ｑ33として配置する。これら対角成分の値は、それぞれの信号（実効音速、実効光吸収係数、実効複素誘電率）で同じ値とする。 For example arranged components corresponding to the effective speed of sound as a diagonal component Q 11 of the covariance matrix Q k, place the component corresponding to the effective optical absorption coefficient as a diagonal component Q 22 of the covariance matrix Q k, the effective complex permittivity The component corresponding to the rate is arranged as the diagonal component Q 33 of the covariance matrix Q k . The values of these diagonal components are the same for each signal (effective sound velocity, effective light absorption coefficient, effective complex permittivity).
図２に戻り、統合処理部７は、カルマンフィルタ６から出力された、実効音速の時系列信号に基づく生体内グルコースの事後推定値と、実効光吸収係数の時系列信号に基づく生体内グルコースの事後推定値と、実効複素誘電率の時系列信号に基づく生体内グルコースの事後推定値とを、カルマンフィルタ６の自乗推定誤差に基づく重みを用いて重み付け平均化処理することにより、これら生体内グルコースの事後推定値を統合する（ステップＳ１１）。 Returning to FIG. 2, the integrated processing unit 7 performs a post-mortem glucose value based on the time-series signal of the effective sound velocity output from the Kalman filter 6 and a post-mortem glucose based on the time-series signal of the effective light absorption coefficient. The estimated value and the ex post facto estimated value of in vivo glucose based on the time series signal of the effective complex dielectric constant are weighted and averaged by using the weight based on the square estimation error of the Kalman filter 6, and the ex post facto of these in vivo glucose are processed. The estimates are integrated (step S11).
図６は本実施の形態に係る統合処理部７の構成を示すブロック図である。統合処理部７は、実効音速の時系列信号に基づく生体内グルコースの事後推定値、実効光吸収係数の時系列信号に基づく生体内グルコースの事後推定値、実効複素誘電率の時系列信号に基づく生体内グルコースの事後推定値の各々について、カルマンフィルタ６の自乗推定誤差から統合処理のための重み付け定数を算出する重み付け定数生成処理部７０１と、この重み付け定数を用いて、複数の呼吸数の推定値を重み付け平均化する重み付け平均化処理部７０２とから構成される。上記のカルマンフィルタ処理（図２ステップＳ１０）では、時刻ごとに呼吸数の推定値が得られると共に変更カルマン利得Ｓkおよび共分散推定誤差行列Ｐk -，Ｐk +が更新される。 FIG. 6 is a block diagram showing the configuration of the integrated processing unit 7 according to the present embodiment. The integrated processing unit 7 is based on an ex post facto estimate of in vivo glucose based on a time series signal of effective sound velocity, an ex post facto estimate of in vivo glucose based on a time series signal of effective light absorption coefficient, and a time series signal of effective complex dielectric constant. For each of the ex post facto estimates of in-vivo glucose, a weighted constant generation processing unit 701 that calculates a weighted constant for integrated processing from the square estimation error of the Kalman filter 6 and an estimated value of a plurality of breathing numbers using this weighted constant. Is composed of a weighted averaging processing unit 702 that weights and averages. In the above Kalman filter processing (step S10 in FIG. 2), the estimated value of the respiratory rate is obtained for each time, and the modified Kalman gain Sk and the covariance estimation error matrix P k − and P k + are updated.
カルマンフィルタ６の自乗推定誤差は、共分散推定誤差行列Ｐk +から得ることができる。本実施の形態では、実効音速の時系列信号に基づく生体内グルコース値の自乗推定誤差をσ１、実効光吸収係数の時系列信号に基づく生体内グルコース値の自乗推定誤差をσ２、実効複素誘電率の時系列信号に基づく生体内グルコース値の自乗推定誤差をσ３とする。 The square estimation error of the Kalman filter 6 can be obtained from the covariance estimation error matrix P k + . In the present embodiment, the squared estimation error of the in vivo glucose value based on the time series signal of the effective sound velocity is σ1, the squared estimation error of the in vivo glucose value based on the time series signal of the effective light absorption coefficient is σ2, and the effective complex dielectric constant. Let σ3 be the square estimation error of the in vivo glucose value based on the time series signal of.
重み付け定数生成処理部７０１は、実効音速、実効光吸収係数、実効複素誘電率のそれぞれの自乗推定誤差σｉから重み付け定数αｉを式（２６）により算出する。 The weighted constant generation processing unit 701 calculates the weighted constant αi from the squared estimation errors σi of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity by the equation (26).
重み付け定数生成処理部７０１は、式（２６）の算出を実効音速、実効光吸収係数、実効複素誘電率の各々について行う。そして、重み付け平均化処理部７０２は、重み付け定数生成処理部７０１が算出した重み付け定数αｉ（ｉ＝１，２，３）を用いて、実効音速の時系列信号に基づく生体内グルコースの事後推定値ハットＧ１と、実効光吸収係数の時系列信号に基づく生体内グルコースの事後推定値ハットＧ２と、実効複素誘電率の時系列信号に基づく生体内グルコースの事後推定値ハットＧ３とを重み付け平均化処理した生体内グルコースの最終推定値ハットＧｆkを式（２７）のように算出する。統合処理部７は、以上のような統合処理を時刻ごとに行う。 The weighting constant generation processing unit 701 calculates the equation (26) for each of the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity. Then, the weighted averaging processing unit 702 uses the weighting constant αi (i = 1, 2, 3) calculated by the weighting constant generation processing unit 701 to post-estimate the glucose in vivo based on the time-series signal of the effective sound velocity. Weight averaging processing of hat G1, post-estimated value of in-vivo glucose based on time-series signal of effective light absorption coefficient hat G2, and post-estimated value of in-vivo glucose based on time-series signal of effective complex dielectric constant hat G3 The final estimated value hat Gf k of in vivo glucose is calculated by the formula (27). The integrated processing unit 7 performs the above integrated processing for each time.
図３に示す生体内グルコース値の推定を説明する模式図において、カルマンフィルタ６から出力された事後推定値は、「×」で示す点に対応している。また、統合処理部７により出力される最終推定値は「黒丸」の点で示される。 In the schematic diagram for explaining the estimation of the in-vivo glucose value shown in FIG. 3, the ex-post estimation value output from the Kalman filter 6 corresponds to the point indicated by “x”. Further, the final estimated value output by the integrated processing unit 7 is indicated by a "black circle".
次に、図５に戻り、エラー測定部８によるエラー係数の決定方法について説明する。エラー測定部８は、光音響センサ２と誘電分光センサ３各々の信号の時間推移を測定して、例えば、標準偏差等を用いてその平均値から大きく外れた場合にエラーとして、また、平均値から変わらなければエラーはないとして検出する。そして、エラー測定部８は、検出したエラーを−１〜＋１に正規化することによりエラー係数を定める（図５ステップＳ１０１）。 Next, returning to FIG. 5, a method of determining the error coefficient by the error measuring unit 8 will be described. The error measuring unit 8 measures the time transition of each signal of the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3, and if the signal greatly deviates from the average value by using, for example, a standard deviation, it is regarded as an error and the average value. If it does not change from, it is detected as no error. Then, the error measuring unit 8 determines the error coefficient by normalizing the detected error to −1 to +1 (FIG. 5, step S101).
具体的には、音速のエラー係数に関しては、エラー測定部８は、音響波の位相変動の振幅を検出する。そして、エラー測定部８は、経験的に０〜±１％をエラー係数±１とエラー係数−１と対応させる（図５ステップＳ１０１）。実効光吸収係数のエラー係数に関しては、エラー測定部８は、経験的に０〜±２％をエラー係数±１とエラー係数−１と対応させる（図５ステップＳ１０１）。また、実効複素誘電率については経験的に０〜±０．５％をエラー係数±１とエラー係数−１と対応させる（図５ステップＳ１０１）。 Specifically, regarding the error coefficient of the speed of sound, the error measuring unit 8 detects the amplitude of the phase fluctuation of the acoustic wave. Then, the error measuring unit 8 empirically associates 0 to ± 1% with the error coefficient ± 1 and the error coefficient -1 (step S101 in FIG. 5). Regarding the error coefficient of the effective light absorption coefficient, the error measuring unit 8 empirically associates 0 to ± 2% with the error coefficient ± 1 and the error coefficient -1 (step S101 in FIG. 5). As for the effective complex permittivity, 0 to ± 0.5% is empirically associated with the error coefficient ± 1 and the error coefficient -1 (step S101 in FIG. 5).
以上説明したように、本実施の形態に係る測定装置１は、生体内グルコース変動モデルの生体内グルコース推定値と、光音響センサ２により検出される音響信号の時系列信号列から抽出される、実効音速と実効光吸収係数各々の時系列信号と、誘電分光センサ３により検出される透過信号又は反射信号の時系列信号列から抽出される実効複素誘電率の時系列信号のそれぞれから推定される生体内グルコースの事前推定値とに基づいて、生体内グルコースの事前推定値に重み付けして、事後推定値を算出する。さらに、それぞれの事後推定値を統合して一つの最終推定値を求める。 As described above, the measuring device 1 according to the present embodiment is extracted from the in-vivo glucose estimated value of the in-vivo glucose fluctuation model and the time-series signal sequence of the acoustic signal detected by the photoacoustic sensor 2. Estimated from each of the time series signals of the effective sound velocity and the effective light absorption coefficient, and the time series signal of the effective complex dielectric constant extracted from the time series signal sequence of the transmitted signal or the reflected signal detected by the dielectric spectroscopy sensor 3. A post-estimation value is calculated by weighting the pre-estimated value of in-vivo glucose based on the pre-estimated value of in-vivo glucose. Furthermore, each ex post facto estimate is integrated to obtain one final estimate.
このように、複数のセンサデータに重み付けして統合することで、生体内グルコース値に含まれる雑音を低減することができ、安定した生体内グルコース値の推定が可能である。 By weighting and integrating a plurality of sensor data in this way, noise contained in the in-vivo glucose value can be reduced, and a stable in-vivo glucose value can be estimated.
また、本実施の形態に係る測定装置１は、エラー測定部８と重み付け定数生成部９と適応推定部１０と重み付け演算部１１とを設けることにより、カルマンフィルタ６のカルマン利得に、光音響信号に関するデータおよび透過信号又は反射信号に関するデータの信頼度を反映させることができる。 Further, the measuring device 1 according to the present embodiment is provided with the error measuring unit 8, the weighting constant generation unit 9, the adaptive estimation unit 10, and the weighting calculation unit 11, so that the Kalman gain of the Kalman filter 6 is related to the photoacoustic signal. The reliability of the data and the data regarding the transmitted or reflected signal can be reflected.
このように、センサデータの信頼度を反映させることで、環境変動に対応して、適応的にカルマン利得を変更しながら複数のセンサデータを統合し、最良のデータを抽出して、被験者の生体内グルコース値を推定する。したがって、測定装置１における生体内グルコース値の推定精度をより高めることができ、生体内グルコース濃度測定の定量精度を向上させることができる。 In this way, by reflecting the reliability of the sensor data, it is possible to integrate multiple sensor data while adaptively changing the Kalman gain in response to environmental changes, extract the best data, and raw the subject. Estimate the glucose level in the body. Therefore, the estimation accuracy of the in-vivo glucose value in the measuring device 1 can be further improved, and the quantification accuracy of the in-vivo glucose concentration measurement can be improved.
本実施の形態で説明した測定装置１は、図７に示すように、バス１０１ａを介して接続されるＣＰＵ１０２ａ、記憶装置１０３ａ、およびＩ／Ｆ１０４ａを備えるコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。ＣＰＵ１０２ａは、記憶装置１０３ａに格納されたプログラムに従って本実施の形態で説明した処理を実行する。 As shown in FIG. 7, the measuring device 1 described in the present embodiment controls a computer including a CPU 102a, a storage device 103a, and an I / F 104a connected via a bus 101a, and their hardware resources. It can be realized by a program. The CPU 102a executes the process described in the present embodiment according to the program stored in the storage device 103a.
以上、本発明の測定装置１における実施形態について説明したが、本発明は説明した実施形態に限定されるものではなく、請求項に記載した発明の範囲において当業者が想定し得る各種の変形を行うことが可能である。 Although the embodiment of the measuring device 1 of the present invention has been described above, the present invention is not limited to the described embodiment, and various modifications that can be assumed by those skilled in the art within the scope of the invention described in the claims are made. It is possible to do.
例えば、本実施の形態では、センサとして、光音響センサ２と誘電分光センサ３とを用いて、３つの異なる特徴量である実効音速、実効光吸収係数、および実効複素誘電率から生体内グルコースの推定値を算出する場合について説明した。しかし、本発明においてセンサはこれらに限られず、磁場センサや加速度センサ、温度センサ、超音波センサ、インピーダンスセンサなど、その他のセンサを組合せて特徴量を抽出してもよい。 For example, in the present embodiment, the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3 are used as sensors, and the effective sound velocity, the effective light absorption coefficient, and the effective complex permittivity, which are three different feature quantities, are used to determine the in vivo glucose. The case of calculating the estimated value has been described. However, in the present invention, the sensor is not limited to these, and a feature amount may be extracted by combining other sensors such as a magnetic field sensor, an acceleration sensor, a temperature sensor, an ultrasonic sensor, and an impedance sensor.
また、本実施の形態では、測定装置１をＳＭＢＧ装置に適用して生体内グルコース値を測定する場合について説明したが、測定装置１は持続式連続血糖値モニタに適用してもよい。測定装置１を持続式連続血糖値モニタに適用することで、光音響センサ２と誘電分光センサ３のキャリブレーションはより正確になるため、生体内グルコース値の測定精度はより高まる。 Further, in the present embodiment, the case where the measuring device 1 is applied to the SMBG device to measure the glucose level in the living body has been described, but the measuring device 1 may be applied to the continuous continuous blood glucose level monitor. By applying the measuring device 1 to the continuous continuous blood glucose level monitor, the calibration of the photoacoustic sensor 2 and the dielectric spectroscopy sensor 3 becomes more accurate, so that the measurement accuracy of the in-vivo glucose level is further improved.
１、１ａ…測定装置、２…光音響センサ、３…誘電分光センサ、４…信号抽出部、５…生体内グルコース推定部、６…カルマンフィルタ、７…統合処理部、７０１…重み付け定数生成処理部、７０２…重み付け平均化処理部、８…エラー測定部、９…重み付け定数生成部、１０…適応推定部、１１…重み付け演算部。 1, 1a ... Measuring device, 2 ... Photoacoustic sensor, 3 ... Dielectric spectroscopy sensor, 4 ... Signal extraction unit, 5 ... In vivo glucose estimation unit, 6 ... Kalman filter, 7 ... Integrated processing unit, 701 ... Weighting constant generation processing unit , 702 ... Weighting averaging processing unit, 8 ... Error measuring unit, 9 ... Weighting constant generation unit, 10 ... Adaptation estimation unit, 11 ... Weighting calculation unit.
Claims (6)
前記特徴量抽出回路により抽出された前記特徴量に基づいて、前記特徴量ごとに前記生体内成分値の事前推定値を算出する事前推定回路と、
前記事前推定回路により算出された前記事前推定値に含まれるノイズを濾した事後推定値を算出するカルマンフィルタと、
前記カルマンフィルタで算出されたそれぞれの前記事後推定値を重み付け平均化して最終推定値を出力する統合処理回路と、
前記生体内成分値に関する前記情報の信頼度を示すエラー係数を測定するエラー測定回路と、
前記エラー係数に基づいて重み付け行列の成分である重み付け定数を求める重み付け定数生成回路と、
前記カルマンフィルタの推定誤差が最小となるように前記カルマンフィルタの共分散推定誤差行列を更新する適応推定回路と、
前記重み付け行列と前記共分散推定誤差行列とに基づいて前記カルマンフィルタのカルマン利得を更新する重み付け演算回路と、
を備えることを特徴とする測定装置。 A feature extraction circuit that extracts information on multiple features based on time-series signals of information on the in-vivo component values of the subject, and a feature extraction circuit.
A pre-estimation circuit that calculates a pre-estimated value of the in-vivo component value for each feature amount based on the feature amount extracted by the feature amount extraction circuit.
A Kalman filter that calculates a post-estimated value by filtering out noise included in the pre-estimated value calculated by the pre-estimated circuit, and a Kalman filter.
An integrated processing circuit that weights and averages each of the post-estimated values calculated by the Kalman filter and outputs the final estimated value.
An error measurement circuit that measures an error coefficient indicating the reliability of the information regarding the in vivo component value, and an error measurement circuit.
A weighting constant generation circuit that obtains a weighting constant that is a component of the weighting matrix based on the error coefficient,
An adaptive estimation circuit that updates the covariance estimation error matrix of the Kalman filter so that the estimation error of the Kalman filter is minimized.
A weighting arithmetic circuit that updates the Kalman gain of the Kalman filter based on the weighting matrix and the covariance estimation error matrix.
A measuring device characterized by comprising.
前記特徴量抽出ステップで抽出された前記特徴量に基づいて、前記特徴量ごとに前記生体内成分値の事前推定値を算出する事前推定ステップと、
前記事前推定ステップで算出された前記事前推定値に含まれるノイズをカルマンフィルタにより濾して事後推定値を算出するフィルタリングステップと、
前記フィルタリングステップで算出されたそれぞれの前記事後推定値を重み付け平均化して最終推定値を出力する統合処理ステップと、
前記生体内成分値に関する前記情報の信頼度を示すエラー係数を測定するエラー測定ステップと、
前記エラー係数に基づいて重み付け行列の成分である重み付け定数を求める重み付け定数生成ステップと、
前記カルマンフィルタの推定誤差が最小となるように前記カルマンフィルタの共分散推定誤差行列を更新する適応推定ステップと、
前記重み付け行列と前記共分散推定誤差行列とに基づいて前記カルマンフィルタのカルマン利得を更新する重み付け演算ステップと
を備えることを特徴とする測定方法。 A feature extraction step that extracts information on multiple features based on time-series signals of information about the in vivo component values of the subject, and a feature extraction step.
A pre-estimation step of calculating a pre-estimated value of the in-vivo component value for each feature amount based on the feature amount extracted in the feature amount extraction step, and a pre-estimation step.
A filtering step in which noise included in the pre-estimated value calculated in the pre-estimated step is filtered by a Kalman filter to calculate a post-estimated value, and
An integrated processing step that weights and averages each of the ex-post estimates calculated in the filtering step and outputs the final estimate.
An error measurement step for measuring an error coefficient indicating the reliability of the information regarding the in vivo component value, and
A weighting constant generation step for obtaining a weighting constant which is a component of the weighting matrix based on the error coefficient, and
An adaptive estimation step that updates the covariance estimation error matrix of the Kalman filter so that the estimation error of the Kalman filter is minimized.
A measurement method comprising: a weighting operation step for updating the Kalman gain of the Kalman filter based on the weighting matrix and the covariance estimation error matrix .
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2017100683A JP6795453B2 (en) | 2017-05-22 | 2017-05-22 | Measuring device and measuring method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2017100683A JP6795453B2 (en) | 2017-05-22 | 2017-05-22 | Measuring device and measuring method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2018192182A JP2018192182A (en) | 2018-12-06 |
| JP6795453B2 true JP6795453B2 (en) | 2020-12-02 |
Family
ID=64571177
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2017100683A Active JP6795453B2 (en) | 2017-05-22 | 2017-05-22 | Measuring device and measuring method |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP6795453B2 (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111141695B (en) * | 2019-12-24 | 2022-11-25 | 中国船舶重工集团公司第七一八研究所 | Non-dispersive infrared multi-component Freon gas detection system |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7020506B2 (en) * | 2003-11-06 | 2006-03-28 | Orsense Ltd. | Method and system for non-invasive determination of blood-related parameters |
| JP4989331B2 (en) * | 2007-06-21 | 2012-08-01 | 三菱電機株式会社 | Wake integration apparatus, program, and wake integration method |
| JP5089281B2 (en) * | 2007-07-26 | 2012-12-05 | 三菱電機株式会社 | State estimation device and state estimation method |
| GB0719969D0 (en) * | 2007-10-12 | 2007-11-21 | Cambridge Entpr Ltd | Substance monitoring and control in human or animal bodies |
| JP2012004181A (en) * | 2010-06-14 | 2012-01-05 | Sharp Corp | Apparatus, method and program for characteristic prediction, and program recording medium |
| JP5709651B2 (en) * | 2011-06-03 | 2015-04-30 | 三菱電機株式会社 | Tracking device |
-
2017
- 2017-05-22 JP JP2017100683A patent/JP6795453B2/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| JP2018192182A (en) | 2018-12-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Blank et al. | Clinical results from a noninvasive blood glucose monitor | |
| US6157041A (en) | Methods and apparatus for tailoring spectroscopic calibration models | |
| US6671542B2 (en) | Non-invasive method of determining skin thickness and characterizing layers of skin tissue in vivo | |
| US6122535A (en) | Method and apparatus for determining the concentration of a component | |
| Maruo et al. | Near-infrared noninvasive blood glucose prediction without using multivariate analyses: introduction of imaginary spectra due to scattering change in the skin | |
| KR20170126310A (en) | Method and apparatus for predicting analyte concentration | |
| JP2019037752A (en) | Measuring apparatus and measuring method | |
| JP5186044B2 (en) | Blood glucose level estimation device | |
| US10966642B2 (en) | Photoplethysmography apparatus | |
| EP1592339A2 (en) | Non-invasive determination of direction and rate of change of an analyte | |
| EP1279025A2 (en) | Methods and apparatus for spectroscopic calibration model transfer | |
| JP2006126219A (en) | Method and apparatus for multi-spectral analysis in noninvasive infrared spectroscopy | |
| WO2004081524A2 (en) | Measurement site dependent data preprocessing method for robust calibration and prediction | |
| Menon et al. | A survey on non-invasive blood glucose monitoring using NIR | |
| JP6795453B2 (en) | Measuring device and measuring method | |
| JP6807300B2 (en) | Component concentration measuring method and component concentration measuring device | |
| Squara et al. | Metrology part 2: procedures for the validation of major measurement quality criteria and measuring instrument properties | |
| US11428672B2 (en) | Component-concentration measuring apparatus and method | |
| US20010041829A1 (en) | Non-invasive method of determining skin thickness and characterizing layers of skin tissue in vivo | |
| JP2010082246A (en) | Method for processing measurement data of biological spectrum | |
| JP6813848B2 (en) | Component concentration measuring device and method | |
| CN113692530A (en) | Portable device and method for non-invasive blood glucose level estimation | |
| KR20170063039A (en) | Apparatus for estimating blood level of an ingredient using specrum analysis | |
| JP6572169B2 (en) | Component concentration measuring apparatus and component concentration measuring method | |
| CN112914564A (en) | Infant blood oxygen saturation monitoring method and intelligent monitoring device |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20190903 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20200831 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20200908 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20201019 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20201110 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20201112 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 6795453 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |