JP2005189179A - 粉粒体の粒度計測方法 - Google Patents

Abstract

【課題】撮像した画像から粒径のばらつきを正確に求めることができる粉粒体の粒度分布計測方法を提供することにある。
【解決手段】撮像した元画像にフィルタリング処理および2値化処理を施した前処理後の画像に対して、複数の拘束条件からなる評価関数を確率論的に最適化処理を行って画像を推定する。この過程で、注目画素の画素変数をその近傍画素に割り当てられた変数からの演算で近似し、更に近傍画素を代表する近傍領域の範囲を調整することによって、活性画素と非活性画素に対する補正作用に偏りを持たせて推定画像を得るようにする。
【選択図】図１

Description

本発明は、粉粒体の粒度分布を画像処理を用いて計測する粒度分布計測方法に関するものであり、特にコークス製造における石炭の粒度分布計測に好適な粒度分布計測方法に関するものである。

粉粒体の粒度分布を計測することは、粉粒体の製造工程では粉粒体の品質管理上、極めて重要である。例えば、石炭を乾留してコークスを製造する工程は、原料（石炭）荷役貯蔵、石炭処理、乾留、コークス処理および化成品処理の各工程より成り立っている。原料（石炭）荷役貯蔵工程では、船舶等で入荷された石炭は、アンローダにより荷揚げされた後、ベルトコンベアで輸送され、スタッカーにより貯炭場に積み付けられる。次の石炭処理工程は、良質コークスが得られるように各種銘柄の石炭を適量に配合し、それらを最適な粒度に粉砕する工程である。先の貯炭場に貯炭されていた石炭は、払い出し設備すなわちローダーなどで払い出され、さらにコンベアで配合槽まで送炭されて各銘柄別に貯蔵される。配合槽に貯蔵される石炭は、石炭の流動性、不活性物質の多少、硬さ、粒度などの炭質を十分考慮して、各銘柄別に配合される。

通常、原料炭として入荷する石炭は、粉状の粉炭であるが、その粒度は、銘柄によって相当に異なっている。とくに不活性物質の石炭は、コークス化性の観点からできる限り細かくすることが好ましく、細かくすることによりコークス強度の向上が期待できる。しかしながら一方では、必要以上に細粒化すると、コークス炉に装入する段階で嵩密度低下をおこし、生産性低下の原因となるとともに、装入作業時の粉塵飛散やコークス炉窯内でのカーボン付着トラブルをひきおこす等の問題点もある。

そこで、互いに相反する要求に対して、両者を満足させる最適な粒度範囲を見出す必要がある。通常、3mm粒径以下のものが、全体の80%前後を占めるような粒度分布になるように操業管理している。具体的には、配合された石炭の粒度分布を計測し、この計測結果に基づき粉砕機の制御を行っている。

このように、コークスの製造工程では、石炭の粒度分布を計測することが操業管理上重要である。従来からの石炭粒度分布の代表的な計測方法として、サンプリングロボットで石炭を自動で採取し、乾燥後に粒径別に篩い分け、秤量することにより粒度分布を求める方法がある。

この方法を用いることにより、精度良く粒度分布を求めることができるものの、バッチ計測であるという計測頻度の制約（1日に数回程度）から、その結果を用いてリアルタイムに石炭粒度分布管理・粉砕機の回転数制御を行うことは困難であるとともに、石炭のサンプリングをするサンプリングロボット装置自体も高価であるという問題があった。

そこで、リアルタイムで粒度分布管理を行える画像解析による粒度分布計測方法の開発が行われてきた。ＣＣＤカメラ等により撮像された石炭画像から画像処理技術により粒度分布を求める方法が考案された。画像処理による粒度分布計測方法の一例として、特開平５−２３１８２１号公報（特許文献１）では、流下する粒状物質を静止画像として映出し、その映出された粒子の体積や重量を計算機による画像処理で算出し、粒度を測定する方法が開示されている。この発明を実施するためには、粒状物質を連続的に流下させる装置が必要であり、コークス粒投入ホッパー、ベルトコンベア、分離スクリーン、分離板、漏斗、反射板、搬送用ホッパー、および分散棒などの構成要素を必要とし、設備化に多大な費用を要する。

もっと簡単な構成による石炭の粒度分布の測定方法として、特開２００２−２２１４８１号公報（特許文献２）の開示がある。ベルトコンベア等の搬送体に載せられて移動中の石炭を撮像して画像処理により粒度を測定する装置であり、撮像装置と被写体との距離変動をなくすために整流板をつける工夫がなされている。

また、粒度分布測定用の画像処理方法に工夫を加えたものとして、特開２００３−８３８６８号公報（特許文献３）の開示がある。これは、撮像された元画像から、当該元画像に「ぼかし処理」を行った「ぼかし画像」を得て、当該「ぼかし画像」を2値化処理することにより、所定粒径以上の測定対象物の分布を測定すると共に、前記撮像された元画像と「ぼかし画像」の差分により形成された「差分画像」を2値化処理することにより、前記所定粒径未満の測定対象物の粒径分布を測定し、これら2種類の粒径測定分布の測定結果に基づいて、全体の粒径分布を測定する粒度分布測定方法である。

さらに、新しい画像処理方法として、Geman ＆ Geman の提案したマルコフ確率場に基づくベイズ推定法（非特許文献１）が提案されている。これは、最適化処理を用いた確率論的画像推定であり、この最適化で用いる評価関数に、複数の拘束条件を採用することにより、データに基づくなめらかな補間を実現している。さらに、この方法を応用した補間法として、特開平５−１５９０７５号公報（特許文献４）の開示がある。これは、マルコフ確率場に基づくベイズ推定法を利用して、まばらな観測データから連続量である各種物理量をなめらかに補間する手法であり、空間的に細かい変化は除去してなめらかな補間を行い、空間的に大きな変化にはなめらかな補間を行わないことを特徴としている。
特開平５−２３１８２１号公報 特開２００２−２２１４８１号公報 特開２００３−８３８６８号公報 S.Geman & D.Geman; "Stochastic Relaxation, Gibbs Distribution, and The Bayesian Restoration of Images", IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol.PAMI−6, No.6,(Nov.1984) 特開平５−１５９０７５号公報

特許文献２で示される方法は、撮像画像を2値化して粒子の像を求め、周知の画像処理方法により粒子の像の大きさを判別し、その大きさの分布から粒度分布を求める方法である。しかしながら、単なる撮像画像には、撮像される粒子の大きさに応じて輝度斑が発生し、この輝度斑が原因となって、2値化を行った場合に、2値化画像の粒子の粒径と実際の粒子の粒径が対応しなくなり、正確な粒度測定が行えないという問題がある。

また、特許文献３で示される方法は、空間フィルタの走査によって「ぼかし処理」を行っているが、実際の撮像では多様な輝度斑やノイズ成分を含んだ元画像が得られるため、その空間フィルタについてのサイズや画素毎の重み、また走査回数の設計は困難となり一意に決められないという問題がある。そして、広い範囲にわたって粒度分布測定を行なおうとする際には、上記問題がより複雑になる。

さらに、非特許文献１および特許文献４で示される方法は、いずれも全てのデータ値について変化に対する作用が同じであるという問題がある。これは例えば2値化画像の場合での、活性画像と非活性画素について変化に対する作用（連結や除去）が同じというものである。活性画素と非活性画素の面積がほぼ同じ場合では問題はないが、実際上は、非活性画素の背景に少量の活性画素が粒子像候補として残っている場合またはその逆の場合が多い。つまり、活性画素と非活性画素の面積には偏りがある場合が多い。このような元画像に対して、上述のマルコフ確率場を利用した画像補正を施すと、粒子像候補である活性画素を消してしまう方向にのみ処理が進みやすい。実際は重要な情報である活性画素の粒子像候補を特に連結しやすくして粒子像として残し、かつ適切な按配で活性画素の孤立領域をノイズとして除去したいといったように、活性画素と非活性画素に対する作用に偏りがある要望が多い。

本発明は上記事情に鑑みてなされたもので、新規な画像処理により粒径のばらつきを正確に求めることができる粉粒体の粒度分布計測方法を提供することにある。

本発明は、撮像した元画像にフィルタリング処理および2値化処理を施した前処理後の画像に対して、当該前処理後の画像と推定画像との距離に関する拘束条件と境界の長さに関する拘束条件からなる評価関数に確率論的最適化処理を行って画像を推定する方法において、前記確率論的最適化処理は、前記前処理後の画像中の注目画素の画素変数をそれ自身とその近傍画素に割り当てられた変数からの演算で近似し、更に近傍画素を代表する近傍領域の範囲を調整することによって、活性画素と非活性画素に対する補正作用に偏りを持たせて推定画像を得て、当該推定画像から外乱に埋もれた粒子像を抽出して粒度分布を得ることを特徴とする粉粒体の粒度分布測定方法である。

また本発明は、請求項１に記載の粉粒体の粒度分布測定方法において、前記最適化処理の解法として、擬似焼鈍(SA:Simulated Annealing)の手法を用いることを特徴とする粉粒体の粒度分布測定方法である。

また本発明は、請求項１または請求項２のいずれか１項に記載の粉粒体の粒度分布測定方法において、異なる評価関数から推定された異なる粒子画像に対応して異なった粒度範囲の粒度分布を測定し、それらの粒度分布を統合して最終的な粒度分布を得ることを特徴とする粉粒体の粒度分布測定方法である。

また本発明は、請求項１ないし請求項３のいずれか１項に記載の粉粒体の粒度分布測定方法において、前記フィルタリング処理は、シェーディング補正処理またはラプラシアンフィルタの走査であることを特徴とする粉粒体の粒度分布測定方法である。

また本発明は、請求項１ないし請求項４のいずれか１項に記載の粉粒体の粒度分布測定方法において、前記粉粒体は、コークス製造過程における石炭であることを特徴とする粉粒体の粒度分布測定方法である。

さらに本発明は、請求項１ないし請求項５のいずれか１項に記載の粉粒体の粒度分布測定方法において、粒度分布の重量換算に、Rosin−Rammlerの式を用いることを特徴とする粉粒体の粒度分布測定方法である。

本発明は、注目画素の画素変数をそれ自身とその近傍画素に割り当てられた変数からの演算で近似し、更に近傍画素を代表する近傍領域の範囲を調整することによって、活性画素と非活性画素に対する補正作用に偏りを持たせて推定画像を得るようにしたので、粒径のばらつきを正確に表わせる粉粒体の粒度分布計測が可能である。

以下、本発明を実施するための最良の形態について、図面と数式を用いて説明する。図２は、本発明を実施するための装置構成の一例を示す図である。図中、１は石炭粒子、２はベルトコンベア、３はエリアセンサカメラ、４はカメラコントローラ、５は記憶装置、６は演算装置、および７は表示装置をそれぞれ示している。まず、カメラコントローラ４によって制御されたエリアセンサカメラ３によって、ベルトコンベア２にて搬送される散在した石炭粒子１の状態を撮像する。

エリアセンサカメラ３には、2次元CCDカメラなどを用い、粒度分布測定範囲内の微細な粒子が充分確認できる位置に設置する。画像採取時期は、ベルトコンベア２の速度と視野範囲から適合する値を選定し、石炭粒子１が完全に砂塵に埋もれてしまっていない場所を選択するのが理想的である。また、表面に露出している粒子は全体のなかのごく一部である可能性があるので、統計的信頼性を増すためには可能な限り多くの画像を得るようにするとよい。

得られた画像は、直ちに記憶装置５に蓄積され、さらに演算装置６において粒度分布を求める演算処理して、表示装置７にて粒度分布計測結果として表示するという流れを取る。演算装置６における演算処理は、本発明の核となるので以下に詳説する。なお、この演算処理した粒度分布は、記憶装置５に蓄積されるとともに、上位のプロセス計算機（（図示せず）にデータ転送し、石炭の粉砕機等への制御指令に用いることもできる。

図１に、演算装置６における画像処理による粒度分布計測処理の流れを説明する図である。初め（S100）で処理をスタートして、まず記憶装置５に蓄積された画像から計測対象となる画像を、画像採取する（S101）。撮像・画像採取された元画像には、光学系に起因する画像のぼけがあり、また散在する砂塵や背景の微細なテクスチャーあるいは電気系などのノイズに覆われている。

そこで、画像の前処理として、フィルタリング処理（S102）を行なう。低周波域の大きな輝度変化を除去して、粒子像など注目する要素を強調する処理であり、例えばシェーディング補正処理、ラプラシアンフィルタの走査など行なう。シェーディング補正処理は、注目画素の近傍領域にて画素変数の平均をとり、注目画素の値から引き算または割り算を行うものである。また、ラプラシアンフィルタの走査は、一例として図10に示したSavitzy−Golayのラプラシアンフィルタカーネルを画像に対して走査して、コンボリューション（たたみ込み積分）をとる。このとき除去される低周波域の大きな輝度変化とは、たとえば照明の具合や撮像角度など光学的要因によるものである。

次に、2値化処理（S103）を行ない、粒子候補の抽出をする。図４に、画像の前処理によって得られた画像の一例を示す。図４（ａ）は、石炭粒子を撮像した元画像、図４（ｂ）は、フィルタリング処理（S102）により低周波域の大きな輝度変化を除去して、粒子像など注目する要素を強調した画像、および図４（ｃ）は、2値化処理後の画像をそれぞれ示し、粒子境界がよりはっきりしてきている様子が分かる。この前処理の対象となる量は、通常は輝度であるが、測定対象物の性質により他の特定の物理量、例えば特定の色の強さ等を選択することもできる。また、2値化閾値等は経験により決定するが、画像が規格化（通常は8bit表現において中心輝度を128にする）されていれば、容易に2値化閾値を決定することができる。
得られた画像のM×N個の画素をI≡[( i,j)|i=1,2,3,・・・,M, j=1,2,3,・・・,N]とし、上に述べた前処理手段により得られた2値化画像X⁰=[x⁰ _i,j|(i,j)∈I]を、候補画像とする。更に、抽出された粒子像候補を構成する画素は活性画素(x⁰ _i,j = +1)、その他は非活性画素(x⁰ _i,j = −1)と呼ぶこととする。なお、8bit表現においては、活性画素の値を255で白、非活性画素の値は0で黒とするとよい。この候補画像には、まだ砂塵などのノイズ要素が含まれており、正確な粒界の認識は困難である。
そこで、この候補画像X⁰=[x⁰ _i,j|( i,j)∈I]から、マルコフ確率場MRF（D）による画像推定（S104）ステップにて、マルコフ確率場を適用した最適化処理により画像を推定する。ここでの処理は、後に詳述するが、異なるマルコフ確率場を適用して候補画像から推定した画像の例を、図11に示す。図11（１）は、候補画像を、（２）から（５）はマルコフ確率場MRF（D）でのパラメータDを０から３まで変化させた推定画像をそれぞれ示している。この例において、大径粒子はマルコフ確率場MRF(３)を適用して推定した画像から、小径粒子はマルコフ確率場MRF(0)を適用して推定した画像から計数するといった測定法が考えられる。このように粒子の特徴に応じたマルコフ確率場を適用して候補画像から画像を推定し、その分布を計測すれば外乱要素を避けることができる。

次に、補正処理（S105）では、既によく知れている画素連結や孤立点除去など確定論的な手法による更なる補正を補助として行う。こうして得られた画像をもとににして、粒径Rmin(D)からRmax(D)までの粒度分布の測定（S106）を行う。ここでは、パラメータDの値によって、最小粒径Rmin(D)と最大粒径Rmax(D)のペアを事前に与えてあり、この範囲での粒度分布解析を行う。処理の詳細は、後述する。

次のステップ−粒度分布範囲全域が包括されたか（S107）にて、希望とする粒度範囲をカバーした処理が終了したかを判定する。処理が終了していなければ、次のステップ−D=D+1（S108）にてパラメータDを増加させて、マルコフ確率場MRF（D）による画像推定（S104）に戻る。

（S107）での処理が終了していれば、粒度分布統合（S109）にて異なった粒度範囲で得られた粒度分布を統合する処理を行う。Rosin−Rammler式補正（S110）では、２次元画像から得られた粒度分布を、重量換算した粒度分布に変換する。そして、結果表示（S111）を行い、終わり（S112）で一連の処理を終了する。

次に、図１のマルコフ確率場MRF（D）による画像推定（S104）の処理について、詳説してゆく。まず、本発明の基礎となっている「マルコフ確率場の例と画像推定の手法」を、以下説明する。

この手法では適当な大きさの粒子像を得るために、候補画像X⁰=[x⁰ _i,j|( i,j)∈I]から画像Φ=[ ψ_i,j|( i,j)∈I]を推定する問題を、以下の（１）式で示される最適化問題として定式化している。
Φ=argmin_xH(X| X⁰, J(D)) ・・・・・・・・（1）
ここで評価関数（エネルギー関数）H(X| X⁰, J(D))は、（２）式および（３）・（４）式で表わされる。
H(X| X⁰, J(D) )≡ d(X| X⁰) +J(D)σ( X) ・・・・・・・・（2）
d(X| X⁰) ≡ 1/4ΣΣ(x_i,j− x⁰ _i,j ) ² ・・・・・・・・（3）
σ(X) ≡ 1/4ΣΣ｛(x_i,j− x_i+1,j ) ²+(x_i,j− x_i,j+1) ²｝ ・・・・・・・・（4）
なお（３）・（４）式での総和は、全ての画素について総計するものである。また、以下の説明では、4近傍の場合についてのマルコフ確率場による画像推定を行なっているが、8近傍の場合への拡張は容易である。

（３）式のd(X| X⁰) は、「候補画像と推定画像の距離」に関する拘束条件であり、画像X⁰と 画像Xにおいて異なっている画素の総数に対応している。「候補画像と推定される画像の距離」は、候補画像と推定画像が全く同じとき最小値0となる。また、（４）式のσ( X)は、「境界の長さ」に関する拘束条件であり、推定画像において活性画素と非活性画素の境界の長さを示している。「境界の長さ」は、推定画像が複雑にいりくんでいるほど大きな値をとることとなる。（１）式の最適化問題においてJ(D)が極端に大きい場合は、「境界の長さ」を小さくするように全てが活性画素（粒子像）であるか非活性画素である画像が推定結果として得られ、逆にJ(D)が0であれば「候補画像と推定される画像の距離」を0にするように候補画像そのものが得られる。このような拘束条件の性質を踏まえて、測定する粒塊の大きさの範囲Dに応じてJ(D)の値を経験的に求める。J(D)が大きいほど大まかな像の推定に適しており、J(D)が小さいほど微細な像の推定に適している。またJ(D)が大きいほど画像がぼけており、J(D)が小さいほど画像が鮮明であるとも言える。

コークス製造過程における石炭の場合は、J(D)が大きいほど大径の粒子に注目した画像推定に適しており、J(D)が小さいほど小径の粒子に注目した画像推定に適している。測定する粒子の大きさが広い範囲にわたっている場合は、J(D)の値を変えて複数の推定画像を得るとよい。図８は、粒径範囲DによりJ(D)を決定するテーブルを示す図であり、この例では、４つの粒径範囲Dをとっている。候補画像に対する推定画像は評価主体がどれを粒子像と想定するかによって変わるものであるから、注目する粒子の大きさによりマルコフ確率場を変えていかなければならない。この画像補正においては、注目する粒子の大きさにより調整するパラメータは、上に述べたJ(D)だけであり、簡便な画像補正が可能である。
x_i,j もx⁰ _i,j も+1（活性）か−1（非活性）しかとらないことから, （２）式の評価関数は下の（５）式のように、磁性体の物性を説明する2次元イジングモデルのハミルトニアンになっている。
H(X| X⁰, J(D))≡1/2ΣΣ｛(x⁰ _i,j x_i,j+J(D)x_i,j x_i+1,j +J(D)x_i,j x_i,j+1) +（J(D)+1/2）MN・・・（5）
ここで粒子画像の推定は、各スピンに異なる外場のかかったイジングモデルの基底状態を探索する問題に帰着される。基底状態の探索に後述する擬似焼鈍(SA:Simulated Annealing)という手法を用いるため、温度Tという量を新たに導入する。擬似焼鈍(SA:Simulated Annealing)は、熱した金属を徐々に冷却していくと、内部熱エネルギーが最低の状態に落ち着く自然界の最適化プロセスを、人工的に真似たものである。適切な冷却計画による擬似焼鈍(SA:Simulated Annealing)は、局所的な解に捕らえられず大局的な最適解を得られることが知られている。以下の説明においては、（５）式について設定されるパラメータJ(D)を交換エネルギーと呼ぶことにする。
上の最適化問題はギブス分布を用いた計算により、（６）・（７）式のような確率最大化問題に置き換えることができる。
Φ=argmax_xρ(X| X⁰, T,J(D)) ・・・・・・・・（6）
ρ(X| X⁰, T,J(D)) ≡exp[− H(X| X⁰, J(D)) /T ]/ Σ_xexp[− H(X| X⁰, J(D)) /T ]・・・・・・（7）
温度Tの値は正であり、最適化問題（１）式と確率最大化問題（６）式は等価である。

ここで確率最大化問題（６）式をそのまま解くのは計算量が多いので、統計力学の手法である平均場近似を用いて解く。平均場近似とは、「画素間の相互作用は、画素変数の平均値に影響を与えるだけで、すべての画素変数は独立な確率変数である」という仮定をおくものである。この仮定により、最隣接する画素間にしか働かない相互作用が平均場というかたちで全ての画素に作用する。すなわち、ある画素変数が平均場を媒介として非常に離れた画素変数に作用することになる。平均場近似は確率最大化問題の解法としては厳密には正しくはないが、石炭粒子分布画像の補正に関する本課題には充分適用できるものである。

これにより画素( i,j)∈I について活性(n=+1)・非活性(n=−1)の確率は、それぞれ（８）式および（９）式のようになる。
ρ_i,j(n| X⁰, T,J(D))=1/2(1+ξ_i,j n), (i,j)∈I, n=±1, Σ_n=±₁ρ_i,j(n| X⁰, T,J(D))=1・・・（8）
ξ_i,j=tanh[1/T(x⁰ _i,j +J(D)ξ_i+1,j +J(D)ξ_i-1,j+J(D)ξ_i,j+1+J(D)ξ_i,j-1)] ・・・・・（9）
充分に小さい正の温度Tにおいて、平均場近似における確率最大化問題は、（１０）式のように帰着される。
ψ_i,j = argmax _n=±₁ρ_ij(n| X⁰, T,J(D)) ・・・・・・・・・・・・（10）
すなわち、J(D)によって特徴づけられる周辺確率分布ρ_ij(n| X⁰, T,J(D))に従って各画素の活性(n=+1)・非活性(n=−1)を決定して、Φ=[ψ_i,j |( i,j)∈I]を推定画像とすればよいこととなる。
画像X^mよりの画像X^m+1の最適化の実際においては、後に説明する模擬焼鈍での冷却段階mがわかるように（９）式に添字を加えて（１１）式に書きあらためる。
ξ^m _i,j=tanh[1/T(x^m _i,j +J(D)ξ^m _i+1,j +J(D)ξ^m _i-1,j+J(D)ξ^m _i,j+1+J(D)ξ^m _i,j-1)]・・・・・・（11）
またξ^m _i,jの初期値は、ρ_i,j(n|X, T,J(D))が確率1で候補画像 X⁰ =[ x⁰ _i,j|( i,j)∈I]を生起するように、｛ξ⁰ _i,j = x⁰ _i,j |( i,j)∈I｝としておく。
x^m+1 _i,j= argmax _n=±₁ρ_ij(n|X^m, T,J(D)) ・・・・・・・・（12)
ρ_i,j(n|X^m,T,J(D))=1/2(1+ξ^m _i,j n), ( i,j)∈I, n=±1, Σ_n=±₁ρ_i,j(n|X^m, T,J(D))=1 ・・・・・・・・（13)
冷却段階m＋１の画像X^m+1は、（１２）および（１３）式のように表わされ，各画素についてその近傍で周辺確率分布を最大化していき、初期画像X⁰=[ｘ⁰ _i,j|( i,j)∈I]から粒子画像を得ることになる。

すなわち図５で示すように、注目画素の画素変数は、そのもともとの値x^m _i,j と近傍画素に割り与えられた変数（以降は平均場近似変数と呼ぶ）ξ^m _i+1,j,ξ^m _i-1,j,ξ^m _i,j+1,ξ^m _i,j-1からの（１１）・（１２）・（１３）式に従う確率論的な計算により変換される。ここで遠くの画素は、平均場近似変数を通して注目画素に影響することになる。
（６）式および（７）式からの演算に比べて、平均場近似によれば大幅に演算量を減らすことができる。特に石炭粒度分布測定においては、統計的信頼性を増すために可能な限り多くの画像を処理することが好ましく、また上流にある石炭粉砕機の制御性を良くすることからも、このような演算の高速化は効果的である。ここで、2次元空間の各画素について（１２）式に従って活性(n=+1)・非活性(n=−1)交換を試行する順序は、各画素位置と対となる一様乱数を発生させ昇順または降順に並び替えるなどの手法でランダムに決めることができる。
上述のようなマルコフ確率場の演算においては、活性画素と非活性画素に対する連結性や除去性が同じである。ところが実際の粒子分布画像補正問題においては、図4（ｃ）の2値化画像にみるように、非活性画素の背景に少量の活性画素が粒子像候補として残っている場合が多い。つまり活性画素と非活性画素の面積には偏りがある。このような元画像に対して、要素の連続性を得るため上述のマルコフ確率場を利用した画像補正を施すと、粒子像候補である活性画素を消してしまう方向にのみ処理が進みやすい。実際は重要な情報である活性画素の粒子像候補を特に連結しやすくして粒子像として残し、かつ適切な按配で活性画素の孤立領域をノイズとして除去したいといったように、活性画素と非活性画素に対する作用に偏りがある要望が多い。
そこで発明者は、図6に示すような新しいマルコフ確率場を提案するものである。ここで近傍画素には、適当な大きさで定義した近傍領域のなかから注目画素を最も活性画素にしやすい平均場近似変数をもったものが代表として選ばれる。容易に理解できるように、この近傍領域を大きくとると活性画素は連結しやすくかつ除去しにくくなる。そこで図6に示すように、この近傍領域を最適化計算の進行につれて収縮させてゆく操作を行なう。この近傍領域の範囲の操作により、最適化演算の初期は分離領域連結の作用が強く、後期には孤立領域除去の作用が強くなる。すなわち、初期には近接した活性画素領域を連結して大きな領域にし、後期には孤立して残った砂塵等による小さな活性画素領域をノイズとして除去するように処理が進み、最終的に良好に補正された石炭粒子分布画像を得ることができる。このように近傍領域を調整したマルコフ確率場による画像補正では、活性画素と非活性画素に対する作用に偏りを持たせることができる。以降は、非対称なマルコフ確率場と呼ぶことにする。このように、本発明では近傍画素を近傍領域から適切に代表して選ぶようにしているので、広く画像補正問題に対応することができる。
（１０）式の確率最大化問題は、温度T（正値）を一定値として解けるが、ここではk+1段階の温度T_i（T₀> T₁> T₂>・・・> T_k>0）からなる冷却計画の擬似焼鈍(SA:Simulated Annealing)によることとして説明する。このことにより大局的な推定画像から次第に最適な粒子画像に収束させ、局所的最適解に捕らえられる危険性を避けることができる。ここではランダムに選んだ各画素について、図7に示したとおりの計算を行う。

図７は、 図2に示したマルコフ確率場MRF（D）による画像推定（S104）の処理の流れを示す図である。まず、始め（S200）で処理をスタートして、2値化処理後の候補初期温度T₀にて画像を初期画像X⁰として設定する（S201）。次に、冷却計画の温度ステージｍを、m=0（最初の温度ステージ）（S202）とする。

さらに、注目画素に対する近傍領域の設定を行なう（近傍領域の範囲をN(m)とする（S203））。近傍領域の設定の仕方は、前述の図６に示す通りである。すなわち、冷却計画の初期の温度ステージでは、図の左のように近傍領域を広く設定し、温度ステージが進む（図中の矢印の方向）に従って近傍領域を狭くする。このことにより、初期は分離領域連結の作用が強く、後期には孤立領域除去の作用が強くなる。すなわち、初期には近接した活性画素領域を連結して大きな領域にし、後期には孤立して残った砂塵等による小さな活性画素領域をノイズとして除去するように処理が進む。

次に注目画素を選択するために、ランダムにk番目の画素I_k=(i_k,j_k)を取り上げる（S204）。そして、確率論的に画素値を決定して、画素の活性・非活性を決定する。すなわち、0から1までの一様乱数RANDの発生（S205）をさせて、この乱数RANDの値と確率ρ_i,j(+1| X^m,T,J(D))を比べる（S206）。ここで、RAND ＜ ρ_i,j(+1| X^m,T,J(D))であれば、値x^m+1 _i,j=+1(活性画素とする)（S207）とし、そうでなければ、x^m+1 _i,j =−1(非活性画素とする)（S208）とする。

各画素について交換を試行した（S209）の判断にて、全画素についての試行が終わったかどうか判定する。試行が終わってなければ、k=k+1 (次の画素へ)（S210）として、（S204）に戻る。全画素についての試行が終われば、冷却計画が終了したか（S211）の判断にて、冷却計画の終了を判定する。通常、冷却計画実行の前に終了温度ステージを与えておく。

冷却計画が終了していなければ、m=m+1(次の温度ステージへ)（S212）として、温度ステージを一段進めて、近傍領域の範囲設定（S203）に戻る。ここで、ｍの温度ステージが進むにしたがい、図６に示したように近傍領域を狭くしてゆく。そして、冷却計画が終了すれば、一連の処理の終わり（S213）となる。

以上の手順によれば、エネルギー増分が正である遷移もゆらぎとして、温度（正値）の増加関数であるギブス分布に従った確率で許すことになる。擬似焼鈍(SA:Simulated Annealing)の冷却計画については初期条件によらず大局的最適解に収束するものが知られているが、実用的には遅すぎる冷却計画であるため、ヒューリスティックに代わりの早い冷却計画を見つける必要がある。発明者が石炭粒子分布画像の補正のためにシミュレーションした結果によれば、T_m+1= T_m− ΔT（ΔTは正の定数）と温度が線形に降下するものでも充分に画像処理効果が得られた。温度を大きくとるとエネルギー増分が正である遷移を許容しやすくなり元画像の性質を反映しない傾向が現れ、小さくとると確定論的な画像処理に近い傾向が現れる。最大温度、最低温度（通常は0に近い正値とする）、冷却回数で決まる冷却計画は、これらの性質を踏まえて、いくつかのサンプル画像により調整する。

また、活性と非活性の判断においては、0から1での区間で適切な刻みで一様乱数を計算機内部で発生させてρ_i,j(+1| X^m,T,J(D))を越える場合と越えない場合に分けて処理を割り振ればよい。また確率ρ_i,j(+1| X^m,T,J(D)) の計算は（１３）式に基づいて行われるが、画像全領域ではなく図５のように注目する画素の近傍領域に限った計算で求められる。さらに、活性(n=+1)・非活性(n=−1)交換の試行は、説明した擬似焼鈍(SA:Simulated Annealing)の冷却段階毎に必ずしも全画素について行う必要はない。
このような手法によれば、注目する粒子の大きさに応じて交換エネルギーJ(D)の値を調整することにより、図１１（１）候補画像で見られる孤立して活性化された画素や活性化領域の複雑な境界がノイズ成分として除去され、連続性の高い粒子像を選択的に生成・抽出する（図１１（２）〜（５））ことができる。すなわち、大径粒子については大きな交換エネルギーJ(D)で推定した画像（図１１（５））によりその分布を計測し、小径粒子については小さな交換エネルギーJ(D)で推定した画像（図１１（２））によりその分布を計測すればよい。
そして、最終的に補正された粒子像に含まれる活性画素の総和を求めて、粒子の断面積とする。ここで粒子の断面積を求めるには、たとえば図９に示す隣接画素の値を調べるラベリングの手法に従えばよい。図９は、初めの（１）白抜きの活性画素が得られた状態〜最終の（８）同じ粒子としてラベリングされるまでの様子を順次模式的に示している。
まず、（２）では、探索の順序と振り返る領域を示している。X方向では左から右へ、Y方向では上から下へと活性画素を探索し、活性画素が探索されれば、振り返る領域は、斜上・隣・斜下・後の４領域であることを示している。（３）では、ラベリングのルール１「活性画素に出会い、まだその斜上、隣、斜下、後のいずれにもドメインNo．が与えられて
なければ新しいドメインNo．を昇順にて与える」の様子（ドメインNo．１、２、３）が、さら
にラベリングのルール２「活性画素に出会い、その斜上、隣、斜下、後のいずれかに既にドメインNo．が与えられていればその古いほうをドメインNo．として与える」の様子（ドメイ
ンNo．１）を示している。探索が進んで、（４）では（ドメインNo．２）、（５）では（ドメ
インNo．１、２、３）がそれぞれラベリングされていることが分かる。
次に、（６）では、既に生成していたドメインNo．2領域がNo．1領域に連結され、かつ領域
No．2が領域No．1に更新される様子を表わしている。これは、「隣接画素に与えられていた
ドメインNo．は新しいドメインNo．に更新される」というラベリングのルール３によるもので
ある。（７）では、最終的に2つのラベル（ドメインNo．１、３）が生成された様子を示し
ている。最終の（８）では、後の処理を容易にするためドメインNo．を昇順にて振りなおし
ている（ドメインNo．１、２）。断面積は、それぞれの同じドメインNo．を振られた画素の総
和をとればよい。以上は、ラベリングの一手法として示したものであり、これに限られるものでなく他の手法を用いてもよい。
この面積Sについて、例えば断面が円と見なした式R_p=2(S/π)^1/2からR_pを粒塊の径の指標とし粒度分布を求める。R_pが所定範囲の外であるとされた粒子についてはノイズ成分であるとして粒度分布の評価に使用しない。異なる交換エネルギーJ(D)の値により複数の補正画像を得る場合は、計数が重複なく、かつ知りたい粒度分布範囲全域を包括させるよう、あらかじめ図8の表のように計数を行う粒径範囲Dについて、それぞれマルコフ確率場を生成する交換エネルギーJ(D)の値を予め決めておく。このように本手法によれば注目する粒子の異なった形状に応じて適切な画像処理を行い、正確な粒度分布を計測することができる。
異なる粒径範囲Dごとに計測された粒度分布を統合した後、各粒子径に対応する石炭の重さの分布を求める。すなわち粒度分布の重量換算をおこなう。比重ρの粒子径R_pの粒子数がn個観察されたとすると、その粒子径の石炭の重さW_rは、以下の（１４）式のように求めることができる。

W_r =nρπR_p ³／６ ・・・・・・・・（14)
そして石炭の粉砕後粒度をRosin−Rammlerの式にフィッティングさせることにより、測定値のバラツキを除去した粒度分布を求める。すなわち、横軸にlogR_pを縦軸にlog(log100/W_r )をとって計測値をプロットする。次に回帰分析によりx=logR_pとy=log(log100/W_r )の関係を表した1次式を求める。すなわち(x,y)=( logR_p, log(log100/W_r ))のデータを、以下の（１５）式にあてはめる。
ｙ＝ｂｘ+c ・・・・・・・・（15)
次に、（１５）式よりW_r＝36.8%のときのR_pの値を求めこれをaとすると、石炭粒子の直径R_pと粒径R_pの石炭粒子の重量割合W_rとの関係は、以下の（１６）式となる。
W_r=100exp(−R_p/a)^b ・・・・・・・・（16)
このようにしてRosin−Rammlerの式に従った重量分布が求まるが、通常は重量割合の総和が100%にならないために過不足分を按分加減する必要がある。すなわち、重量割合の総和をΣW_rとすると、各重量割合に100/ΣW_rを最終的な重量割合とする。このようにして最終的な粒度分布を求めて、結果は例えば図3のようなヒストグラムにした後に、図２の（7）表示装置にてオペレータに表示する。

なお、粒度分布計測の適用対象はこれまで説明した石炭とは別に、鉱石、石灰石、コークス、プラスチック粉、トナーなどの各種粉粒体にも適用することが可能である。

画像処理による粒度分布計測処理の流れを説明する図である。 本発明を実施するための装置構成の一例を示す図である。 粒度分布ヒストグラムの例を示す図である。 画像の前処理の変遷の一例を示す図である。 マルコフ確率場による画像補正における平均場近似を説明する図である。 非対称なマルコフ確率場の演算を示す図である。 マルコフ確率場による画像推定の処理の流れを示す図である。 粒径範囲DによりJ(D)を決定するテーブルを示す図である。 粒界推定画像から粒子の面積を求める方法を示す図である。 ラプラシアンフィルタカーネルの一例を示す図である。 異なるマルコフ確率場での画像推定の比較を示す図である。

符号の説明

１ 石炭粒子
２ ベルトコンベア
３ エリアセンサカメラ
４ カメラコントローラ
５ 記憶装置
６ 演算装置
７ 表示装置

Claims (6)

1. 撮像した元画像にフィルタリング処理および2値化処理を施した前処理後の画像に対して、当該前処理後の画像と推定画像との距離に関する拘束条件と境界の長さに関する拘束条件からなる評価関数に確率論的最適化処理を行って画像を推定する方法において、
   前記確率論的最適化処理は、前記前処理後の画像中の注目画素の画素変数をそれ自身とその近傍画素に割り当てられた変数からの演算で近似し、
   更に近傍画素を代表する近傍領域の範囲を調整することによって、活性画素と非活性画素に対する補正作用に偏りを持たせて推定画像を得て、
   当該推定画像から外乱に埋もれた粒子像を抽出して粒度分布を得ることを特徴とする粉粒体の粒度分布測定方法。
2. 請求項１に記載の粉粒体の粒度分布測定方法において、
   前記最適化処理の解法として、擬似焼鈍(SA:Simulated Annealing)の手法を用いることを特徴とする粉粒体の粒度分布測定方法。
3. 請求項１または請求項２のいずれか１項に記載の粉粒体の粒度分布測定方法において、
   異なる評価関数から推定された異なる粒子画像に対応して異なった粒度範囲の粒度分布を測定し、それらの粒度分布を統合して最終的な粒度分布を得ることを特徴とする粉粒体の粒度分布測定方法。
4. 請求項１ないし請求項３のいずれか１項に記載の粉粒体の粒度分布測定方法において、
   前記フィルタリング処理は、シェーディング補正処理またはラプラシアンフィルタの走査であることを特徴とする粉粒体の粒度分布測定方法。
5. 請求項１ないし請求項４のいずれか１項に記載の粉粒体の粒度分布測定方法において、
   前記粉粒体は、コークス製造過程における石炭であることを特徴とする粉粒体の粒度分布測定方法。
6. 請求項１ないし請求項５のいずれか１項に記載の粉粒体の粒度分布測定方法において、
   粒度分布の重量換算に、Rosin-Rammlerの式を用いることを特徴とする粉粒体の粒度分布測定方法。

Images