OpticStudio で人間の皮膚および光学心拍数センサーをモデル化する方法

Also available in English, 简体中文

光電式容積脈波記録法（PPG）は、皮膚表面の生理学的計測を行う、低コスト、非侵襲の光学技術です。この方法が最も幅広く適用されている用途の 1 つが、市販のスマートウォッチやスポーツブレスレットなどに搭載された心拍数のウェアラブルセンサーです。日々の活動において心拍を快適かつ継続的に測定できます。この記事では、OpticStudio で生理学的計測の対象として人間の皮膚をモデル化する方法と、PPG による心拍数センサーの時間依存シミュレーションを ZOS-API を使用して実行する方法を紹介します。 

Authored By Csilla Timar-Fulep

添付ファイルのダウンロード

Introduction

PPG デバイスは、赤外または可視光領域の発光ダイオード（LED）と光センサーから構成されます。組織内の血液の体積変化を検出する、シンプルな光学的手段を提供します。血液は、周辺の他の組織よりも強く光を吸収および散乱します。したがって、血液の脈動は、ディテクタに逆位相の信号変化として検出されます。この記事では、OpticStudio で人間の皮膚組織モデルを実装する方法と、PPG デバイスで測定される時間依存信号を ZOS-API アプリケーションによってシミュレートする方法を紹介します。 

基本設計

PPGセンサーは、反射または透過モードのいずれかで使用できます。光の侵入深さは波長によって異なります。このため緑と黄色のLEDは表在血管の血流測定に最適で、通常は反射モードで使用されます。これに対して、赤外および近赤外の波長は深部組織の血流測定に適しており、透過モードでの使用も可能です。この記事では、反射モードのPPGデバイスのユースケースを取り上げます。 

目標は、関連文献で発表されているデータに基づき、現実的な皮膚のモデルを開発することです。したがって、皮膚と血液の両方について光学パラメータが広く入手できる波長の適用を目指しました。そのような波長は、市販デバイスで最も頻繁に使用されている値にも近くなります。検討の結果、モデル化する波長として 575 nm、光源として QSMF-C160 LED（Avago Technologies 社製）を選択しました。このLEDのモデルは Radiant Source Models ライブラリから直接ダウンロードが可能で、Radiant Source Model ファイルから光線を生成することで光源（ファイル）を作成できます。これは、OpticStudio Premium Edition で使用できる機能で、ナレッジベース記事『How to generate a ray set from an RSMX Source Model』で詳細に解説しています。

人間の皮膚のモデル

人体組織からなる媒質内の光伝搬をシミュレートするために、表皮、真皮、皮下脂肪を考慮した、層状の皮膚モデルを実装しました。このユースケースの第一目標は、PPG ベースの心拍数センサーのシミュレーションです。PPG センサーで重要になるのが、血液の脈動によって生じる変化の測定であることから、脈動が観測される層の正確なモデル化に重点を置きました。この方針に従い、真皮を構成する層を、含まれる血液量の値に応じて別々にモデル化しました。具体的には、乳頭層、乳頭下層、上部毛細血管網層、網状層、下部毛細血管網層です。一方、表皮に血液は含まれないことから、モデルをシンプルに保つために、角質層、顆粒層、有棘細胞層、基底細胞層のすべてを 1 つの厚い表皮層と見なすことにしました。最後に、発表されているほとんどの皮膚モデルと同様に、皮下脂肪も 1 つの層で表現しました。 

OpticStudio では、以上の層のすべてを矩形体積としてモデル化します。層の厚みは文献データに基づいて設定し、断面の幅は側面からの光がディテクタ方向に漏れないサイズとしました。後ろの層は、前の層を基準オブジェクトとして配置し、Z 位置の値には、前の層のZ長さの列を参照するピックアップソルブを適用します。この方法により、各層を互いに隙間なく密着して配置できます。 

組織層のカスタマイズ

このケーススタディは、文献データのみを使用しており、研究の過程で新たな測定は一切行っていません。今回のモデルパラメータは公に発表されているデータに基づいてはいるものの、人の皮膚の光学パラメータは母集団間で大きく異なる可能性があることに注意が必要です。つまり、特定の被験者集団には、異なるパラメータが必要になるかもしれないということです。したがって、用途によっては、より正確なデータを入手できる可能性があり、その場合はモデルを遠慮なくカスタマイズしてください。 

皮膚内の個々の血管をすべて記述するには、複雑な空間配置を持つオブジェクトを数百個追加する必要があり、モデルの普遍性が損なわれることから、文献でもそのような手法によるモデル化は一般的ではありません。この研究でもその手法は採用しませんでした。代わりに、血液と周辺組織の光学パラメータの重み付け平均を計算することで、各種皮膚層内の血液量を考慮しました。 

OpticStudio では、以下の生データに基づき、材質ソルブの［モデル］（Model）を使用して皮膚層の材質をモデル化しました。 

組織内のバルク散乱

生体組織などの混濁媒質内の微粒子による光散乱は、Henyey-Greenstein 分布関数で正確に記述できます。Henyey-Greenstein モデルの自由パラメータは、異方性パラメータ g の 1 つだけです。このパラメータの領域は [-1, 1] の範囲にあり、g = -1 は後方散乱、g = 0 は等方性散乱、g = 1 は前方散乱に対応します。散乱光の角度分布は、次式で求められます。 

OpticStudio のノンシーケンシャルモードでは、DLL（Henyey-Greenstein-bulk.DLL）による Henyey-Greenstein バルク散乱モデルを使用できます（DLLはOpticStudio のインストールファイルに含まれています）。このモデルについては、DLL の総合的な分析とともにナレッジベース記事『Henyey-Greenstein分布を使用したバルク散乱をモデル化する方法』に詳述されています。これに加えて、OpticStudio による表面散乱およびバルク散乱モデルの詳細な説明については、ナレッジベース記事『OpticStudioで利用可能な散乱モデル』をお読みください。 

本稿に添付された多層の皮膚モデルでは、文献に発表されている現実的な値に基づいて、各層の散乱パラメータを設定しています。Henyey-Greenstein散乱 DLL の入力パラメータは平均光路（Mean Path）、透過率（Transmission）、異方性パラメータgです。一方、文献では、一般的に、異方性パラメータとともに散乱係数 µs と吸収係数 µa が掲載されています。このため、モデルの入力パラメータを次式によって計算しました。 

前述の屈折率と同様に、各種皮膚層の散乱係数、吸収係数、異方性の各パラメータを、血液とそれ以外の組織に対する重み付け平均として計算しました。計算は、照射する波長 575 nm に対応する、次の生データに基づいています。 

下図は、多層皮膚モデル、および組織内の光伝搬の様子を示す 3D レイアウトプロットです。個々の散乱イベントを表示するために、プロットではセグメントごとに光線の色を変えています。 

プロットに加えて、数値データによる結果も得られるように、設計には 3 個のディテクタ（矩形）オブジェクトも配置しました。ディテクタは、皮膚表面から薄い空気の隙間によって隔てられています。2つのディテクタは断面の幅が皮膚層と同じで、一方は光源、もう一方は皮膚モデルの方向を向いています。これによって前者は入射光、後者は後方散乱光のすべてを参考値として捕捉します。3 番目のディテクタは小さく（2 mm x 2 mm）、これも皮膚の方向を向いています。PPG デバイスの一般的な光ディテクタを表します。 

以上紹介した設計は、この記事の添付ファイル（skinModel.zar）に含まれており、測定/シミュレーションの時間依存性を考慮する必要がない場合は、そのまま皮膚モデルとして使用できます。一方、心拍数センサーのように時間に依存する効果をモデル化する方法については、次のセクションで説明します。 

心拍センサーのシミュレーション

心拍数モニターをシミュレートするために、組織内で脈動する血流を模した ZOS-API を使用しました。皮膚層内の血液量を調整することで、心周期内のさまざまな位相をモデル化しました。その後、検出された後方散乱光を、時間刻みの関数として検討しました。層内の血液量の時間変化は、各層において同時に比例関係で変化するものと仮定し、倍率として考慮しました。今回の例では、.NET を介してOpticStudio に接続された Python API を使用してモデルパラメータを変化させ、微調整された設定によって光線追跡を実行し、最後に結果を分析およびプロットしました。スタンドアロンアプリケーションの Python スクリプトは、記事の添付セクションからダウンロードできます（PPGsimulation.py）。 

API からの組織パラメータの変更

文献によれば、心臓が収縮期に血管に血液を送り出す際、皮膚層の血液量は 2 倍になることがあるといいます。この脈動特性を記述するために経験則による関数を使用しました。関数を 10 ステップ/心周期で 5 心周期にわたってプロットしたのが下図です。 

これに基づいて、まず各層の血液量を計算し、つづいて Henyey-Greenstein 散乱分布の屈折率、平均自由光路、透過率、異方性パラメータ g を血液量に応じて更新しました。対応する Python コードは次のとおりです。 

layer = TheNCE.GetObjectAt(layerNum) 
solver = layer.MaterialCell.CreateSolveType(ZOSAPI.Editors.SolveType.MaterialModel)
solver._S_MaterialModel.IndexNd = n
layer.MaterialCell.SetSolveData(solver)
volPhysData = layer.VolumePhysicsData
volPhysData.ModelSettings._S_DLLDefinedScattering.MeanPath = meanPath
volPhysData.ModelSettings._S_DLLDefinedScattering.SetParameterValue(0, transmission) 
volPhysData.ModelSettings._S_DLLDefinedScattering.SetParameterValue(1, g)

最後に、時間刻みごとに光線追跡を実行して結果を得ます。小さなPPGディテクタで測定した、到達光線の総数と吸収されたパワーです。光線追跡を高速化するために、簡易光線分割を適用しました。この手法の詳細については、ナレッジベース記事『What is Simple Splitting?』をお読みください。光線追跡とデータ抽出は、次のPythonコマンドによって実行できます。 

NSCRayTrace = TheSystem.Tools.OpenNSCRayTrace()
NSCRayTrace.ClearDetectors(0)
NSCRayTrace.SplitNSCRays = True
NSCRayTrace.ScatterNSCRays = True
NSCRayTrace.UsePolarization = True
NSCRayTrace.IgnoreErrors = True
NSCRayTrace.SaveRays = False
NSCRayTrace.Run()
NSCRayTrace.WaitForCompletion()
NSCRayTrace.Close()

hits = -3 # pixel =-3 for total hits
power = 0 # pixel = 0 for total power
Data = 0 
hits_bool_return, total_hits = TheNCE.GetDetectorData(detectorNum, hits, Data, 0)
power_bool_return, total_power = TheNCE.GetDetectorData(detectorNum, power, Data, 0)

結果

血液の散乱および吸収係数は、その他の組織に比べてはるかに大きいため、血液量が多ければ、後方散乱光とディテクタで測定されるパワーが減少し、血液量が少ない場合はその逆になります。今回のシミュレーションによれば、1 W の光源と各時間刻みで 105 本の解析光線を使用した場合、血液量のモデル化された脈動は、ディテクタ信号に約 10～15 %の変動となって表れました。この様子を次のプロットに示します。 

これらの結果は、文献に掲載された実験結果とよく一致しています。市販の心拍数センサーでは通常、このデータに後処理としてわかりやすい信号処理アルゴリズムを適用します。第一段階としてノイズの多い信号に対してスムージングを行い、つづいて事前に定義した一定のしきい値を超えるピークの、1分間あたりの数を数えることで心拍数を計算します。 

References

1. Maeda, et al. Monte Carlo simulation of spectral reflectance using a multilayered skin tissue model. Optical Review (2010)
2. Sinichkin, et al. In vivo fluorescence spectroscopy of the human skin: experiments and models. Journal of Biomedical Optics (1998)
3. Meglinski, Matcher. Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions. Physiological Measurement (2002)
4. Doronin, et al. Assessment of the calibration curve for transmittance pulse-oximetry. Laser Physics (2011)
5. Meglinski. Monte Carlo simulation of reflection spectra of random multilayer media strongly scattering and absorbing light. Quantum Electronics (2001)

KA-01986