RRUFFのラマンデータで遊ぶ の履歴(No.24)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
RRUFFはアリゾナ大学のDowns教授（今は名誉教授になっている）らが作成している鉱物のラマンスペクトルのオープンデータベースで、ラマンスペクトルだけでなく同試料の組成分析、IRスペクトルや粉末XRDパターンも公開されています。最近サイトがリニューアルされました。このサイトでは鉱物名や化学組成で検索して、出てきた試料のラマンスペクトル等を見ることができますが、手持ちの観察スペクトルデータやピーク位置（ラマンシフト）を使っての検索機能まではありません。ラマンスペクトルのデータベース検索については、別途CrystalSleuthというプログラムが同じサイトで古くから公開されています（現在はOther databasesのSoftwareのところ）。ただしCrystalSleuthはWindows用で、また2008年で更新が止まっていて、付属データベースもそのころの古いままのようです。さらに私はMacユーザーなので利用しづらい問題もあります。そのため検索プログラムを自作することを以前から考えてはいたのですが、この機会に作ってみることにしました。また最終的には自分の持っているスペクトルデータも追加して、合わせて利用できるようにしたいと思ってます。
RRUFFプロジェクトは元々Appleの初代CEOのMichael Scottがお金を出したそうです（RRUFFの名付け親でもあるらしい）。その辺りの話は引用すべき論文の中に詳しいです。Downs教授とはIMA神戸の時にちょっとだけ話した記憶があります。京都大学にしばらく滞在したことがあると言ってました。記憶は曖昧ですが、1990年くらいのGoldschmidtの時にもちょっと話したことがあったかも。
以下は前半の方ではピーク位置のindexファイルを作ってそれを使って検索する方法を、後半ではスペクトル自体を使ったbrute forceな検索方法について書いてます。
この件についてはJpGU2026で自作簡易ラマン分光器と絡めてポスター発表に出そうかなと予定しています。
RRUFFのラマンスペクトル集はRRUFFのサイトから自由にダウンロードすることができます。いくつかのデータセットがありますが、私は最もサイズの大きいexcellent_unoriented.zipをダウンロードしました。解凍してみると2025/12/03の時点で11,413点のファイル（テキストファイル）がありました。excellent_unorientedのファイルを詳しく見ると、測定したデータ(Data_RAW)と処理したデータ(Data_Processed)があるので、実際に検索に使うスペクトルとしてはその約半分となります。このデータセットの相対波数範囲は100くらいから1300 cm-1くらい（ダイヤモンドやグラファイトなどはもう少し高い波数まで）になります。4000 cm-1くらいまでの広い波数範囲が欲しい場合はLR-Raman.zipの方をダウンロードしてください。後で他のファイルも合わせてindexファイルを作ろうかと思ってますが…
ファイルをいくらか読んでみると、フォーマットは大体は統一されているのですが、全く同じではないので、読み込む時に注意が必要なので、以下に注意点をまとめました。
プログラムはpythonで作成し、いつものようにGUIにはTkinterを使うのですが、以下の方針で作ることにしました。
最初にRRUFFのラマンデータファイルを開いてプロットして、平滑化処理やピーク検出を試すテストプログラムを作成しました。それを使っていくらかのファイルを処理して分かったことは、find_peaks()関数によるピーク検出はバックグラウンド（傾き）があるスペクトルではあまりうまくいきません（ピーク自体は弱くてもバックグラウンドが高いと拾われる）。そういう意味でもバックグラウンド処理済みのProcessedを使うことになります。また検索プログラムの自分のスペクトルの前処理でもバックグラウンド除去が必須となります。
Savitzky-Golay法で平滑化した方がより正しいピーク検出が行われるようでした。しかし平滑化は場合によっては元のスペクトルを大きく変化させてしまうことがあります。そのために大きく変化した場合（最大ピーク強度が半分以下になる）はパラメータを変えて再度計算しなおすようにしてます。ピーク検出では検出感度のパラメータを自動で調整して、8個ピークが得られるようにしてます。元々ピーク数が８個もない鉱物もあるため、index.txtの一部は本当のピークを拾ってない場合もありますが、そこはあまり大事にはならないだろうと考えてます。この辺りをもう少しrefineするには自分でピーク検出部分を書く必要がありそうです。
ピーク強度は取り入れてませんが、強度の順番にindex.txtファイルに入れているので、それを定性的な強度情報として検索で利用してます。ラマンの強度は方位に依存するので、検索で強度にあまり重点を置かない方がいいと考えました。
Processedを全部処理するプログラムを作成して、それを実行しました。結構速くて１分以内で終了します。index.txtのファイルサイズは815 kb程度となりました。見ると5779行あります。RAWデータを無視しているので、ファイル数のほぼ半数ですが、対応するRAW自体がないものもあるので、半数よりは少し多くなってます。また、ガラススペクトルに関して５つほどピーク検出がうまくいってなくてその分少なくなってます。
index.txtファイルは先頭８つがピーク位置を高いものから低いものの順番で並べて、最後にそのファイル名を入れてあります。ファイル名は後でオリジナルのファイルにアクセスするためと、ファイル名から鉱物名を拾うためです。なおピーク位置はピーク関数を使ってfitはしてなくてfind_peaks関数で得られたものをそのまま使ってます。ピーク位置は校正の方法などに依存するので実際的には使った装置が違えばある程度のずれが存在すること、鉱物の場合は固溶体が多く、ピーク位置は組成でかなり変わりうるので、同定目的では非常に正確なピーク位置を求める必要性はないと考えるからです（その分、検索時には幅を広げて検索する必要があります）。
続いて検索するプログラムを作成しました。同定したいスペクトルを読んでプロットし、平滑化、必要ならバックグラウンドを引く処理をして8個ピーク検出して、それをlistbox1に入れます。そのリストとindex.txtを比較して、一致度がいいものをlistbox2にリストアップします。それらを実行するにはプログラムのwinodwで、ボタンを上から順番にクリックすることで達成できます。最終的にリストアップされたものを選択して、スペクトルとして表示して、同定したいスペクトルと一致しているかをユーザー自身が判断することで同定を行います。
実測スペクトルによってはピークが１つ（ダイヤモンド）とか２つ（グラファイト）しかない場合があります。そういう場合には無理矢理8個拾ってきた位置で検索すると余分なものまで候補として拾ってきます。また宇宙線由来のピークやレイリー散乱線をピークとして拾ってくる場合もあります。そこでlistbox1の8個のピークで検索で使わないピークをユーザーが選択できるようにしました（複数選択可）。そうやって不要なピークを除外して検索するとより正しい候補が得られます。
今の所、一致度（gof）には非常に単純な指標を使ってます。toleranceで指定した値の範囲内にあるindex.txtにあるピークを一致したとして、強度の高い順に8から1のポイントを与えて、その合計ポイントを計算しているだけです。ラマンの強度は方位にもよるので強度の強弱は必ずしも一定ではありませんが、意外にもこの単純な指標で正しい候補が得られてます。ただし正解ではない候補のスペクトルと比べてみると全然似てないスペクトルの場合が多いので、逆に候補側の強度が高いものからのピークが存在するかどうかも同様に調べて一致度に加えるように直しました。これがあると候補の確度が上がります。ピーク位置のずれによる指標をもう少し複雑な関数にしてもいいのですが、測定データが校正されていない場合や固溶体でピーク位置がズレるなどもあるので、複雑な関数を使っても検索がうまくいく訳ではなさそうです。
バックグラウンド除去には移動平均を利用する方法を使いました。これはS.Bruckner (J.Appl.Cryst.33,977-979,2000)で提案されている方法で、こちらで紹介されていてpythonのサンプルコードも公開されてます。これを実行しておくとピークの検出がより正確になります。
自分で測定したデータ数十個でテストしてみましたが、ほとんどの場合は出てきた候補の中に正解が含まれていて同定できてました。うまくいかないケースも数件ありましたが、それはデータベースにない、データベースにあるスペクトルがよくない、測定スペクトルがよくない（ブロードでピーク数が少ない）などが理由でした。
面白いことに石英の測定スペクトルを使って検索したところ、石英が候補リストに出てくるのは当然なのですが、それ以外でredmonditeも候補リストに出てきました。そこでRRUFFのredmonditeのスペクトルを見てみたところ、石英のスペクトルと完全に一致しました。これはRRUFFに対応を報告済みです（既に修正されている）。さらにもう１件、おかしいものを見つけます。
自分の持っているラマンデータも検索に使えるようにしたいのですが、その前に使うスペクトルの吟味や情報追加などが必要で少しづつ処理してます。
まずはプログラムのGUI windowが出てきます。全然洗練されてませんが、基本的には上からボタンをクリックしていけば処理が進みます。赤くなっているボタンが次の処理を示してます。結果はプロットされるので、あまりうまくない時はパラメータを変えて、再度ボタンをクリックすると改善されるはずです。
まずはLoadボタンをクリックして、データを読み込みます。例は隕石中の鉱物のラマンスペクトルです。読み込むとMatPlotLibを使ってスペクトルを表示します。
次にSmoothボタンをクリックすると平滑化したスペクトルが表示されます。あまりスペクトルの変形がひどい時はwindowの値を小さくして再度ボタンをクリックします。
バックグラウンドがある場合はbackgroundボタンをクリックすると、バックグラウンドを引いたスペクトルが表示されます。
次にピーク検出を行います。パラメータは普通はそのままで問題ないはずですが、近接する２つのピークで片側が検出できてない場合はDistanceを小さくします。find peaksボタンをクリックすると、スペクトルとピークを検出した位置が示されます。この場合はバックグラウンドを引いたスペクトルでなくsmoothingしただけのスペクトルが表示されます。また、下のリストボックスに検出したピーク位置のリストが表示されます。
スペクトルとリストボックスのピークリストを吟味して、本来のピークでないものがあった場合（例えば宇宙線由来ピーク）はリストボックスのその行をクリックすると選択されます。選択されたものは検索で使いません。それでSearch indexボタンをクリックすると検索が始まって、一致のよい候補を一番下のリストボックスに表示します。より多くの候補を見たい場合はthresholdを小さくしてください。測定データが波数校正してない等の場合はtoleranceを少し大きくした方がいいかもしれません。なお、RRUFFのデータは約50から100 cm-1以下の低波数データがないので、私のスペクトルではレイリー散乱やantistokes側を含むものが多いのですが、50 cm-1以下は使わないようにして、プロットでも使わないようにしてます。
候補のリストボックスから見たいものを選んでPlotボタンをクリックするとそのスペクトルと測定されたスペクトルが同時に表示されるので、２つが一致しているかどうかを判断できます。また隣のOpenボタンをクリックするとテキストエディタでデータファイルが開きます。その試料の細かい情報がデータファイルのヘッダ部分にあります。この例の場合は候補にはほとんどenstatiteのみがリストされており、enstatiteのスペクトルをプロットすると測定スペクトルとよく一致していることが分かりますので、enstatiteと同定できました。
最新版では鉱物名で検索できるようにしました。その結果はリストボックスに入るので、先ほどと同じように１つ選んで、Plotボタンをクリックするとそのスペクトルと比較することができます。鉱物名で検索して比較を行う場合には最低観測スペクトルを読み込んでおく必要があります。
この例の場合はenstatiteなので、orthoenstatiteかclinoenstatiteかどちらかの可能性があります。この例の場合はorthoenstatiteと同定しました。
可能性としてはprotoenstatiteもあり得るのですが（特に常圧合成の場合）、残念ながらRRUFFにはprotoenstatiteのデータがありません。protoenstatiteは670 cm-1付近のピークが２本でなく１本になるので、簡易的にはそれで判定できます（実際にはスペクトルを全体的に評価しますが）。今回のは明らかにprotoenstatiteではありません。隕石からprotoenstatiteは見つかったことがないはずなので、もし見つけられたら論文書けます。
参考のためにコードを公開してます（Mac上で作成）。クリックするとファイル自体が開いてしまうので（その場合、セーブしてください）、リンク上で右クリックでダウンロードするを選んでください。また実際に使うためにはpythonライブラリの追加（matplotlib, numpy,scipyなど）、ファイルのディレクトリをPCに合わせて変更するなどが必要になります。一応Windowsでも動くはずです。
上記のindexファイルを使った方法で十分同定できるのですが、ケースによってはあまりうまく同定できない場合もあります。例えば相手がピークの非常に多いスペクトルだと結構一致が生じてしまうためにindexファイルを使った方法では候補の上位に来ますが、スペクトルを目で比べてみると全く合ってません。
そこでスペクトル自体を使う方法としてコサイン類似度を考えてみました。コサイン類似度はスペクトルデータの値をベクトル要素と見て、それが作る多次元ベクトルの方向をcosine関数で評価して、２つが一致しているかどうかを判定する方法です。強度は標準化する必要はありませんが、横軸は同じサイズに合わせる必要があります。
まずは２つのラマンスペクトルを読み込んで、それのコサイン類似度を計算してみました。コサイン類似度は２つのスペクトルが完全一致だと1に、完全に不一致だと0になる量です。上の例のようにolivine同士のスペクトルだと0.8以上になってます（類似度はプロット中に示してます）。一方、下のオリビンとエンスタタイトの場合は0に近い値になりました。色々試してみると同じ鉱物同士のスペクトルだと0.8以上になっていました。同定に十分使えそうです。
１スペクトルの処理時間（ファイル読み込み時間も込み）を調べると10 ms以下なので、RRUFFのファイル（excellent_unorientedのprocessedのみで5800個くらい）を全部処理しても１分以内で終わりそうです（実際はもっと速かった）。
これをベースに検索プログラムを作ってみました。最初同定したいラマンデータを処理しますが、そこはindexファイルを使うプログラムとほぼ同じです。それから１つづつRRUFFデータを読み込んで、横軸を揃えて、コサイン類似度等を計算するだけです。指標のよいものから100個をリストボックスに入れてます。そのリストボックスで選択して、スペクトルを比較することができます。この辺りもindexファイルを使うプログラムとほぼ同じです。
コサイン類似度の計算はscipyのライブラリを使ってます。以下をプログラム先頭に書いておきます。
from scipy.spatial.distance import cosine
そしてコサイン類似度の計算はスペクトル（numpy1次元行列）を下のようにcosine()に与えるだけです。cosine()は距離を与えるので、類似度にする場合は1から距離を引きます。
cs = 1.0 - cosine(spectrum1,spectrum2)
このようにコサイン類似度の計算はライブラリを使うと簡単ですが、面倒なところは比較するスペクトルはその範囲と横軸が共通ではないので、それを合わせてからcosine()に入れる必要があります。そのために片側の横軸に合わせた強度を再計算するために内挿法を使ってます。これにはscipyのinterpolateライブラリを使ってます。計算速度を上げるために最も単純な線形近似を使ってますが、特に問題はないようでした。
コサイン類似度だけだとちょっと正解でないケースが多かったので、ピアソン相関係数も入れた指標を使うように今はなってます。ピアソン相関係数も以下で導入できます。
from scipy.stats import pearsonr
そして以下で計算できます。ccの方にピアソン相関係数が入ってます。
cc,p = pearsonr(spectrum1,spectrum2)
ピアソン相関係数は1から-1の値を取るので、ピアソン相関係数に１を足して、２で割ったものとコサイン類似度を掛けたものの平方根（sqrt）を現在は指標に使ってます。これが最適かどうかはまだチェックしているところですが、コサイン類似度やピアソン相関係数を単独で使うよりはいいようです。
cx = math.sqrt(cs*(cc+1.0)/2.0)
それで検索プログラムを作ってみました（最後の図の中にコサイン類似度、ピアソン相関係数と指標が書いてあります）。チェックした限りは大体は正しい候補を示してくれるようです。検索にかかる時間は10秒以内なのでこのbrute forceプログラムでも十分実用的だと思いました。RRUFFにはガラスのラマンデータもいくらかあるのですが、ガラスはブロードなので比較的何にでも一致してしまうためにガラスが上位に来る場合がありました。現在ガラスはファイル読むところで除外してます。
indexファイルを使った方で正解が出て、こちらで出ない場合もその逆もあるので両方使ってみるのがいいのかもしれません。現状diamondのようなかなり高波数に出るピークもうまく検索できてません。indexファイルの方だと不要なピークを使わない選択ができますが、こちらではそれができないために間違った候補がでる場合がありました。
これをテストしていたらRRUFFの間違いをまた１つ見つけました。apatiteの実測スペクトルで検索するとapatiteが候補で出るのは当然ですが、なぜかmonazite-(La)も出るので、詳しくみるとどうもmonaziteでなくapatiteのようでした。どちらもPO4があるので、それの伸縮振動が強くほぼ同じ位置に出るので間違ったのかも。RRUFFには他の方法でまだ同定していないデータが載っていることがあるので注意が必要です。
これをテストしていると何となくCrystalSleuthの動作と似た感じがあるので、CrystalSleuthも同様な検索方法を使っているのかもしれません。CrystalSleuthについては、どう検索しているかとかの情報が出てないのでよく分からないのですが…
参考のためにコードを公開してます（Mac上で作成）。クリックするとファイル自体が開いてしまうので（その場合、セーブしてください）、リンク上で右クリックでダウンロードするを選んでください。また実際に使うためにはpythonライブラリの追加（matplotlib, numpy,scipyなど）、ファイルのディレクトリをPCに合わせて変更するなどが必要になります。Windowsでも動くはずです（文字フォントやボタンのサイズがあまり合ってませんが）。
