JP2004109105A - Flaw type classification boundary setting method in surface defect detection, and defect detection method - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a flaw type classification boundary setting method in surface defect detection allowing even an operator without having a special technique to simply realize complicated defect classification, and enabling evaluation capable of remarkably reducing work and of rapidly responding to change of an operation state. <P>SOLUTION: Feature data are inputted for every flaw type (steps S11-S14); parameters (the average vector, a covariance matrix and a coupling coefficient) of a mixture density distribution for every flaw type are found (step S15); an identification boundary for every flaw is determined by using Bayes' theorem based on the found mixture density distribution for every flaw type (step S16); and the likelihood of the determination is calculated (step S17). <P>COPYRIGHT: (C)2004,JPO

Description

【００２３】
と表される。ここで混合密度分布を特徴付ける各Ｎ次元正規分布の平均ベクトルμ_ｋ，ｍと共分散行列σ_ｋ，ｍ、および結合係数ζ_ｋ，ｍは、ＥＭアルゴリズムによる繰り返し計算によって求める。
【００２４】
混合密度分布モデルとＥＭアルゴリズムについては例えば「クラスター分析入門」（森北出版株式会社）（非特許文献１）に詳しいが、概略を述べると次のようになる。表面欠陥検査装置にて観測されたＮ次元特徴量パターンの不完全データ系列ｘ＝（ｘ_１，ｘ_２，ｘ_３，…，ｘ_ＮＸ）に加えて、仮想的な完全データ系列ｙ＝（ｙ_１，ｙ_２，ｙ_３，…ｙ_ＮＸ）を考える。そして完全データ系列が観測されるのではなく部分的な不完全データ系列ｘ＝（ｘ_１，ｘ_２，ｘ_３，…，ｘ_ＮＸ）が観測されているのだと考える。また完全データ系列ｙ＝（ｙ_１，ｙ_２，ｙ_３，…ｙ_ＮＸ）はパラメータ条件つき確率密度関数ｆ（ｙ｜Φ）に従うものとすると、条件ｘおよびΦ’のもとで、ｆ（ｙ｜Φ）の対数をとったもの、すなわち対数尤度関数の期待値は
Ｑ（Φ｜Φ’）＝Ｅ（ｌｏｇｆ（ｙ｜Φ）｜ｘ，Φ’）
となる。
【００２５】
このときＥＭアルゴリズムは次のようになる。
▲１▼まず各Ｎ次元正規分布とその結合状態についてのパラメータの初期推定値をΦ^（０）とし、ｔ＝０とする。
▲２▼条件つき期待値（Ｅｘｐｅｃｔａｔｉｏｎ）Ｑ（Φ｜Φ^（ｔ））を計算する。
▲３▼上のＱ（Φ｜Φ^（ｔ））を最大化（Ｍａｘｉｍｉｚａｔｉｏｎ）するパラメータΦ^〜を求める。すなわち
Φ^〜＝ａｒｇ　ｍａｘ_ΦＱ（Φ｜Φ^（ｔ））
▲４▼ｔ＝ｔ＋１、Φ^（ｔ）＝Φ^〜とし解が収束するまで▲２▼▲３▼を繰り返す。
【００２６】
具体的には完全データ系列とは、表面欠陥検査装置にて観測されたＮ次元特徴量パターン系列（不完全データ系列）のそれぞれについて、結合するＮ次元正規分布のいずれから生じたかについての情報が加わったものとする。
【００２７】
ここでＮＸ個のＮ次元特徴量パターンのデータ系列ｘ_１，ｘ_２，ｘ_３，…，ｘ_ＮＸが、与えられたとき、結合するＮ次元正規分布をＮＧ個とすると混合分布決定の繰り返し計算は図２に示すようなものとなる。図２は、図１におけるステップＳ１５を更に詳細に示したものである。ここで結合するＮ次元正規分布のうちｚ番目の正規分布の平均ベクトルをμ_ｚ ^（ｔ）、共分散行列をσ_ｚ ^（ｔ）結合係数をζ_ｚ ^（ｔ）とする。
【００２８】
まず、ステップＳ１５１において、パラメータ（結合係数、平均ベクトル、共分散行列）の初期値を入力する（前回値を入力する場合は、後に説明する）。続いて、これらのパラメータを繰り返し計算で収束させることにより求める。
【００２９】
すなわち、図２にて流れを示した以下の（１）→（２）→（３）→（４）→（５）の計算を収束するまで繰り返すものである。
（１）　特徴量パターンｘ_ｊがｚ番目（ｚ＝１，２，…，ＮＧ）のＮ次元正規分布から生じた確率ω_ｚｊ ^（ｔ）を計算する（ステップＳ１５２）。
【００３０】
【数２】

【００３１】
総和Σは結合する全てのＮ次元正規分布についておこなう。ここでＮ次元正規分布の式は（１）式に示したとおりである。
（２）　Ｑ（Φ｜Φ^（ｔ））を最大とするようにｚ番目のＮ次元正規分布についての結合係数ζ_ｚを更新する（ステップＳ１５３）。
ζ_ｚ ^{（ｔ＋１）}＝｛Σω_ｚｊ ^（ｔ）｝／ＮＸ
総和ΣはＮＸ個の全てのデータ列についておこなう。
（３）　Ｑ（Φ｜Φ^（ｔ））を最大とするようにｚ番目のＮ次元正規分布についての平均ベクトルμ_ｚを更新する（ステップＳ１５４）。
μ_ｚ ^{（ｔ＋１）}＝｛Σω_ｚｊ ^（ｔ）ｘ_ｊ｝／（ＮＸ・ζ_ｚ ^{（ｔ＋１）}）
総和ΣはＮＸ個の全てのデータ系列についておこなう。
（４）　Ｑ（Φ｜Φ^（ｔ））を最大とするようにｚ番目のＮ次元正規分布についての共分散行列σ_ｚを更新する（ステップＳ１５５）。
【００３２】
【数３】
σ_ｚ ^{（ｔ＋１）}
＝｛Σω_ｚｊ ^（ｔ）（ｘ_ｊ−μ_ｚ ^{（ｔ＋１）}）（ｘ_ｊ−μ_ｚ ^{（ｔ＋１）}）^ｔ｝／（ＮＸ・ζ_ｚ ^{（ｔ＋１）}）
【００３３】
総和ΣはＮＸ個の全てのデータ系列についておこなう。
（５）　必要に応じて、Ｎ次元正規分布の消滅処理（ステップＳ１５６）、Ｎ次元正規分布の分裂処理（ステップＳ１５７）を行う。これらの処理は、必ずしも行わなくてもよいので、後に詳しく説明する。
（６）　混合密度分布Ｐ（ｘ｜Φ^{（ｔ＋１）}）を計算する（ステップＳ１５８）。
【００３４】
【数４】
Ｐ（ｘ｜Φ^{（ｔ＋１）}）＝Σζ_ｚφ（ｘ：μ_ｚ ^{（ｔ＋１）}，σ_ｚ ^{（ｔ＋１）}）
【００３５】
総和ΣはＮＧ個の全てのＮ次元正規分布についておこなう。
そして、計算が収束したかどうかを判断する（ステップＳ１５９）。この収束は、例えばパラメータ（結合係数、平均ベクトル、共分散行列）のうちの一つの前回との差が、所定値以下となったことにより判断することができる。
【００３６】
以上の計算において、結合するＮ次元正規分布の数を多くとると複雑な識別境界に対応できるが、収束性が悪くなるため事前に経験的に適当な値を決める。また結合するＮ次元正規分布の数は各疵種で必ずしも同じである必要はないので、要求される識別境界の複雑さと収束性から個別に選んでもよい。
【００３７】
また最初の学習で設定するＥＭアルゴリズムの初期値は、まずデータを１つの多次元正規分布に従うものとし、平均ベクトルと共分散行列を求め、共分散行列の初期値はその時の共分散行列そのものとし、平均ベクトルの初期値は、その多次元正規分布に従う乱数を計算機内部で発生させて結合する正規分布の数だけ決めてもよい。また、結合係数の初期値は、それぞれの多次元正規分布について全て同じ１／ＮＧ（ＮＧは結合する正規分布の数）としてもよい。ＥＭアルゴリズムの収束計算では、初期値の選び方によって大局的最適解ではなく、局所的最適解に捕らえられてしまう可能性がある。そのため、最初の学習においては、出力画面に表示された識別境界と確からしさの分布を見て、疵種識別境界が正しく行われなかったと判断される場合は、別の初期値を与えて計算をし直す必要がある。
【００３８】
この方法として、手動で初期値を入力し直す方法や、乱数を発生し直す方法が考えられるが、最初の学習でのＥＭアルゴリズムの収束計算については、いくつかの初期値について試行してみてその中で最適な解を選んでも大きな手間にはならない。
【００３９】
ＥＭアルゴリズムによる混合密度分布の推定においてＮ次元正規分布の共分散行列の行列式が０に近くなれば、（１）式を演算する際にオーバーフローが発生する。共分散行列の行列式が０に近くなるのはＮ次元正規分布のある変数間にて相関が極端に強いことを示し、その最たる場合は分布が点状になることを意味する。このようないびつなＮ次元正規分布の生成も、ある局所的最適解への収束と考えることができる。
【００４０】
これを回避するには別の初期値から出発して適切な解に収束させることも考えられるが、そのようないびつなＮ次元正規分布は消去してしまうことも考えられる（図８（ａ））。ここで複数のＮ次元正規分布のうち、いびつな形になったものを消去した場合、結合するＮ次元正規分布の個数が減少する可能性がある。
【００４１】
またＮ次元特徴量空間に大きく広がりすぎたＮ次元正規分布を逆に分裂させることも考えられる（図８（ｂ））。このようにすれば、最初おおまかにクラスタリングした学習データを細かにクラスタリングし直すという動作を加えることができ、最終的に識別境界設定の精度が向上する。ここで複数のＮ次元正規分布のうち、大きく広がったものを分裂させた場合、結合するＮ次元正規分布の個数が増加する可能性がある。
【００４２】
このようにＮ次元正規分布が消滅と分裂を繰り返すことにより、初期に設定する数にあまり気をつかわなくとも、適切な個数となる。どのとき、どのように消滅又は分裂させるかどうかはＮ次元正規分布の共分散行列の固有値と正規直交固有ベクトルから決定することができる。
【００４３】
図２の「Ｎ次元正規分布の消滅処理」（Ｓ１５６）について図９のフローチャートに従って説明すると、またはＮ次元正規分布の各軸方向の広がりが小さいとき（Ｓ１５６３）、結合するＮ次元正規分布の数を１つ減らし（Ｓ１５６４）、その結合係数を残りの全ての正規分布に分配し（Ｓ１５６５）、該当する正規分布を消滅させる（Ｓ１５６６）。
【００４４】
それぞれのＮ次元正規分布は、その共分散行列σ_ｋ，ｍのＮ個の正規直交固有ベクトル方向に広がる超楕円体のようになる。そして共分散行列σ_ｋ，ｍの固有値と正規直交固有ベクトルの組を（λ_１，ａ_１）、（λ_２，ａ_２）、（λ_３，ａ_３）、・・・、（λ_Ｎ，ａ_Ｎ）とすれば、各軸方向の標準偏差に相当するばらつきは
ｂ_１　＝　λ_１ ^１／２ａ_１
ｂ_２　＝　λ_２ ^１／２ａ_２
ｂ_３　＝　λ_３ ^１／２ａ_３
・・・・・・・・
ｂ_Ｎ　＝　λ_Ｎ ^１／２ａ_Ｎ
であるから、これらいずれかの絶対値が予め決められた閾値（Ｄ_ｍｉｎ）より小さいときに消滅させればよい。
【００４５】
α番目のＮ次元正規分布が消滅したとすると
ζ_１ ^（ｔ）→ζ_１ ^（ｔ）＝ζ_１ ^（ｔ）＋１／（ＮＣ−１）・ζ_α ^（ｔ）
ζ_２ ^（ｔ）→ζ_２ ^（ｔ）＝ζ_２ ^（ｔ）＋１／（ＮＣ−１）・ζ_α ^（ｔ）
ζ_３ ^（ｔ）→ζ_３ ^（ｔ）＝ζ_３ ^（ｔ）＋１／（ＮＣ−１）・ζ_α ^（ｔ）
・・・・・・・・
ζ_α ^（ｔ）→ζ_α ^（ｔ）＝０
・・・・・・・・
ζ_Ｎ ^（ｔ）→ζ_Ｎ ^（ｔ）＝ζ_Ｎ ^（ｔ）＋１／（ＮＣ−１）・ζ_α ^（ｔ）
と結合係数を分配して、以降はα番目のＮ次元正規分布を除いた（ＮＣ−１）個のＮ次元正規分布にてＥＭアルゴリズムによる混合密度分布の推定を行う。
【００４６】
図２の「Ｎ次元正規分布の分裂処理」（Ｓ１５７）について図１０のフローチャートに従って説明すると、Ｎ次元正規分布の各軸方向の広がりが大きいとき（Ｓ１５７３）、結合するＮ次元正規分布の数を１つ増やし（Ｓ１５７４）、その該当するＮ次元正規分布の結合係数の分裂（Ｓ１５７５）、共分散行列の分裂（Ｓ１５７６）、平均ベクトルの分裂（Ｓ１５７７）を後述する規則にて行う。
【００４７】
それぞれのＮ次元正規分布は、その共分散行列σ_ｋ，ｍのＮ個の正規直交固有ベクトル方向に広がる超楕円体のようになる。そして共分散行列σ_ｋ，ｍの固有値と正規直交固有ベクトルの組を（λ_１，ａ_１）、（λ_２，ａ_２）、（λ_３，ａ_３）、・・・、（λ_Ｎ，ａ_Ｎ）とすれば、各軸方向の標準偏差に相当するばらつきは
ｂ_１　＝　λ_１ ^１／２ａ_１
ｂ_２　＝　λ_２ ^１／２ａ_２
ｂ_３　＝　λ_３ ^１／２ａ_３
・・・・・・・・
ｂ_Ｎ　＝　λ_Ｎ ^１／２ａ_Ｎ
であるから、これらいずれかの絶対値が予め決められた閾値（Ｄ_ｍａｘ）より大きいときに分裂させればよい。
【００４８】
β番目のＮ次元正規分布が分裂したとすると、その後の結合係数ζ_β１ ^（ｔ），ζ_β２ ^（ｔ）は元の結合係数ζ_β ^（ｔ）を
ζ_β ^（ｔ）→ζ_β１ ^（ｔ）＝１／２ζ_β ^（ｔ），ζ_β２（ｔ）＝１／２ζ_β ^（ｔ）
のように等分したものとし、
また共分散行列σ_β１ ^（ｔ），σ_β２ ^（ｔ）は
σ_β ^（ｔ）→σ_β１ ^（ｔ）＝ｈσ_β ^（ｔ），σ_β２ ^（ｔ）＝ｈσ_β ^（ｔ）
のように元の共分散行列σ_β ^（ｔ）に定数（ｈ＜１）を乗じたものとし、
また平均ベクトルμ_β１ ^（ｔ），μ_β２ ^（ｔ）は、
μ_β ^（ｔ）→μ_β１ ^（ｔ）＝μ_β ^（ｔ）＋Δμ_β，μ_β２ ^（ｔ）＝μ_β ^（ｔ）−Δμ_β
のように元の平均ベクトルμ_β ^（ｔ）から分裂させる方向ｂ_βの延長線上（Δμ_βとｂ_βは平行）に等しい間隔だけ離れたものとする。
すなわち対応する正規直交固有ベクトルをｂ_β＝（ｂ_βｘ，ｂ_βｙ）とすると
【００４９】
【数５】

【００５０】
であり、
Δμ_β＝Δｋ（ｃｏｓθ_β，ｓｉｎθ_β）
となる。
【００５１】
以上のようにして消滅と分裂を繰り返し、結合する正規分布の個数を変えてＥＭアルゴリズムによる混合密度分布の推定を行う。ここで許容するＮ次元正規分布の広がりの範囲（Ｄ_ｍｉｎ，Ｄ_ｍａｘ）をあまりに狭くとりすぎると、ＥＭアルゴリズムによる混合密度分布の推定の収束性が悪くなることがある。あらかじめ採取データや模擬データなどにて演算をおこない、収束性と結果の精度を確認しておくことが好ましい。
【００５２】
操業状態の変化等のため疵種識別境界設定を修正しなければなくなり、またそのための新たな学習データが蓄積されたときは、前回の混合密度分布を特徴づける平均ベクトル、共分散行列、結合係数などのパラメータを初期値としてＥＭアルゴリズムによって更新し、疵種識別境界設定を修正するようにすることが好ましい。疵種分類の多くの場合、操業状態の変化等による疵種分類境界の修正は大きなものではないので、局所的最適解に捕らえられてしまう可能性は少ない。
【００５３】
実際には各疵種について並行して繰り返し計算を行う。そして、疵種α_ｋが特徴量パターンｕを生起する確率をＰ（ｕ｜α_ｋ）＝Ｐ（ｕ｜Φ_ｋ ^{（ｔ＋１）}）として求め、次の計算をおこなう。
【００５４】
ＥＭアルゴリズムの繰り返し計算の収束は、２次元の場合は画面に表示された識別境界と確からしさの分布を見て確認するが、３次元以上の場合は図３の（２）に代表されるような２次元断面の観察と、それぞれのカテゴリでの対数尤度の変化により確認する。
【００５５】
対数尤度は確率密度関数がＰ（ｘ｜Φ_ｋ ^{（ｔ＋１）}）のとき
Ｌ＝ΣｌｏｇＰ（ｘ_ｉ｜Φ_ｋ ^{（ｔ＋１）}）
（総和は全ての観測データｘ_ｉについておこなう）
と表され、収束が進むにつれて増加し、完了すると変化しなくなる。また対数尤度の代わりに平均ベクトルと共分散行列の各成分の変化を観察してもよい。また大きな次元の場合、カルーネン・レーベ展開や主成分分析などで次元を小さくしてから上の識別境界設定をおこなうとよい。
【００５６】
欠陥種類および有害度について分類するためにベイズの定理に基づいた決定をおこなう（ステップＳ１６）。分類すべき欠陥種類または有害度を広く疵種と呼ぶことにすると特徴量パターンｕが与えられたとき、その疵種をα_ｋとする確率をＰ（α_ｋ｜ｕ）、疵種α_ｋが特徴量パターンｕを生起する確率をＰ（ｕ｜α_ｋ）、疵種α_ｋの生起する事前確率をＰ（α_ｋ）、特徴量パターンｕの生起確率をＰ（ｕ）とすれば、ベイズの定理により
Ｐ（α_ｋ｜ｕ）＝Ｐ（ｕ｜α_ｋ）Ｐ（α_ｋ）／Ｐ（ｕ）
【００５７】
ここでＰ（ｕ）は疵種によらず、Ｐ（ｕ｜α_ｋ）は上で述べた混合分布に対するＥＭアルゴリズムによって求める。また事前確率Ｐ（α_ｋ）はデータ採取がランダムサンプリングと仮定できるのならＰ（α_ｋ）＝Ｎ_ｋ／ΣＮ_ｌ（Ｎ_ｌは系列ｌについて得られたデータ数であり、総和は全ての系列についてとる）とするのが自然であるが、データ採取についてそのような仮定が成り立たない場合はＰ（α_ｋ）＝１／ｃ（ｃは系列数）のように意図的に全ての疵種について同じにしてもよいし、ある特定の誤識別に伴う損失が著しく大きいなど特定の疵種をより正しく識別したい場合はそのときのＰ（α_ｋ）を他の事前確率に対して大きく見積もってもよい。すなわちＰ（α_ｋ）は調整係数として機能する。
【００５８】
Ｐ（ｕ）は疵種によらないため、
Ｐ’（α_ｋ｜ｕ）＝Ｐ（ｕ｜α_ｋ）Ｐ（α_ｋ）
を求め、最もＰ’（α_ｋ｜ｕ）が高いものを疵種とするような識別境界の設定を行い図３（１）のようにテーブルに書き込む。また、このテーブルは計算毎、ＣＲＴ等に描画しオペレータに表示するようにしてもよい。
【００５９】
つぎに特徴量空間全域について、判断の確からしさとして情報のエントロピーを求め、テーブルに書き込む（ステップＳ１７）。情報のエントロピーは事象の予測しやすさを表現する量であり、事象ｌの生成確率をＰ_ｌ（ただしΣＰ_ｌ＝１，総和は起こる事象全てについてとる）とすると
Ｈ＝−ΣＰ_ｌｌｏｇＰ_ｌ　（総和は起こる全ての事象についてとる）
という式で表される。全ての事象が同じ確率で起こる最も予測しにくいとき情報のエントロピーは最大となり、ある事象だけが起こり（すなわちその生成確率＝１）他の事象が全く起こらない（すなわちその生成確率＝０）という、最も予測しやすいとき、情報のエントロピーは最低値０になる。
【００６０】
この量は今回の場合は判断の確からしさを表す、判別境界の近くなど判断に迷う部分では情報のエントロピーは高くなり、確からしさが低く、判別境界の内側など判断に自信のある場所では情報のエントロピーが低くなり確からしさが高くなる。今回は多次元の特徴量パターンｕ_ｎが与えられたとき、上でベイズの定理より計算したＰ’（α_ｋ｜ｕ_ｎ）を
Ｐ”（α_１｜ｕ_ｎ）＝Ｐ’（α_１｜ｕ_ｎ）／ΣＰ’（α_ｌ｜ｕ_ｎ）
Ｐ”（α_２｜ｕ_ｎ）＝Ｐ’（α_２｜ｕ_ｎ）／ΣＰ’（α_ｌ｜ｕ_ｎ）
Ｐ”（α_３｜ｕ_ｎ）＝Ｐ’（α_３｜ｕ_ｎ）／ΣＰ’（α_ｌ｜ｕ_ｎ）
…
Ｐ”（α_ｋ｜ｕ_ｎ）＝Ｐ’（α_ｋ｜ｕ_ｎ）／ΣＰ’（α_ｌ｜ｕ_ｎ）
…
のようにΣＰ”（α_ｌ｜ｕ_ｎ）＝１（総和は全ての疵種についてとる）となるよう規格化してから上式で情報のエントロピーを計算し、図３（２）のように確からしさをテーブルに書き込めばよい。このエントロピーを表示するようにすれば、判断の確からしさを判断することができる。
【００６１】
ステップＳ１８においては、このエントロピーの量を見たり、実際の識別境界を画面で見たりして、識別境界が適当なものであるかを判定する。そして、判別境界が不適当なものであると判断した場合は、ステップＳ２０に移行して、ＥＭアルゴリズムの初期値を変えてステップＳ１５に戻り、再度計算を行う。適当なものであると判断されれば、ステップＳ１９に移行し、平均ベクトル、共分散行列、結合係数等のパラメータ、及び識別境界と確からしさをテーブルに保存して終了する。
【００６２】
ここで確からしさは、Ｎ次元特徴量空間の各点毎の確からしさである。例えば、Ｎ次元特徴量空間と特徴量パターンｕが欠陥カテゴリＡに分類されたとき、その点ｕが欠陥カテゴリＡであることの確からしさを表すものである。具体的には前に述べたエントロピーのような量がこれに相当する。
【００６３】
以上の説明においては、Ｎ次元特徴量空間全域に亘って識別境界を設定する方向で説明した。しかしながら学習データとして入力されるＮ次元の特徴量パターンは有限の範囲において観測される。すなわちＮ次元特徴量空間は有界領域である。Ｎ次元特徴量空間の外部では学習データが観測されないことからＮ次元特徴量空間の縁の近傍ではクラスタリングのためのデータが不足してしまう。そのため図６のように推定される混合密度分布の信頼性が低くなる。
【００６４】
そこで図５に示すように、混合密度分布の推定による手法をＮ次元特徴量空間の縁近傍を除いた内部に限定して適用したほうがよい場合が多い。その際、Ｎ次元特徴量空間の縁の近傍においては、従来からよく知られたｉｆ　ｔｈｅｎ　ｅｌｓｅ文によるツリー型の論理などを使用して識別境界を設定するか、混合密度分布の推定により得られた結果を外挿して識別境界を設定する。だだしこれらの場合でも通常、学習データ系列は特徴量空間全域に渡って得られた学習データを全て採用する。Ｎ次元特徴量空間の縁近傍の混合密度分布によらず識別境界が設定される部分については確からしさの指標である情報のエントロピーを０（最も判断が確からしい領域）としておく。このようにすれば識別境界の設定の柔軟性および安定性の両者で効果的である。
【００６５】
オンラインでの欠陥識別時においては画像処理装置にて観測された特徴量パターンが識別境界テーブル上で対応する点をもとめて疵種を決定する。複数の表面欠陥検出手段がある場合は同様に確からしさテーブル上で対応する点をもとめて比較し、それが最も高いものを信頼できる最終的な疵種分類結果とする。
【００６６】
以上に述べたような、特徴量空間を占める疵種についての新たなデータ列がオペレータの目視判断による分類と共に学習データ保存装置に蓄積され、また操業状態の変化等のため疵種識別境界設定を修正しなければならないとき、学習演算装置内の前回の混合密度分布を特徴づける平均ベクトル、共分散行列、結合係数などのパラメータを初期値としてＥＭアルゴリズムによって更新し、疵種識別境界設定を修正する。これが、図２のステップＳ１５１において、パラメータの前回値を入力することに相当する。
【００６７】
以下、本発明を，鉄鋼分野の熱間圧延プロセスの光学式表面欠陥検出装置に適用した例について述べる。
【００６８】
図７は本実施例における熱間圧延プロセスの概要である．圧延ロールは７段になっており，入り側最初のロールをＮｏ．１，最終ロールをＮｏ．７ロールとする。入り側より１３００℃程度の温度のスラブが挿入され，次第に圧延されて，Ｎｏ．７ロール出口で１ｍｍ程度にまで圧延される．ここでの鋼板移動速度は２０ｍ／ｓｅｃ程度，鋼板の幅は１０００ｍｍである。
【００６９】
撮像のための光源としては，例えばメタルハライド２５０Ｗランプの光をバンドルファイバーで線状に形成し，シリンドリカルレンズで集光して鋼板表面を照らす。受光はラインセンサカメラを用いる。カメラの画素数は１０２４であり，１０００ｍｍの板幅を見るので，横方向の画素分解能は１ｍｍである。カメラコントローラ１は、ロールに設けられたロータリエンコーダより信号を受け取って、鋼板が１ｍｍ進むごとにカメラにトリガ信号を送る。カメラはコイルの全長に亘ってトリガ信号が入った時のみ逐次、１０２４ラインの画像信号を出力するので，鋼板の運転速度が変化しても，画像メモリ内では常に１画素が実際の鋼板での縦横１ｍｍの像に相当している。得られた情報は画像メモリ２に蓄積する。
【００７０】
図７に従って処理の流れを説明すると、まず画像データの欠陥候補から得られる多次元の特徴量を画像処理装置３にて求め、オペレータの目視判断による疵種分類と共に学習データ入力装置４に入力し、学習データ保存装置５に蓄積しておく。ここで多次元の特徴量は例えば欠陥の幅、欠陥の長さ、濃度積算値の３次元である。
【００７１】
全ての疵種についてデータが所定の数以上集まったなら、オペレータが設定盤１１の学習開始ボタンを押すと、学習演算装置６が疵種分類の計算をおこなう。まずおのおのの疵種について混合密度分布をＥＭアルゴリズムによって求め確率論的にクラスタリングする。経験から予想されるＮ次元正規分布の個数を初期に設定して、混合密度分布の推定演算を始める。演算のなかで結合させるＮ次元正規分布がＮ次元空間に充分に広がっていない場合は消滅させ、大きく広がっている場合は分裂させる。以上のようにして消滅と分裂を繰り返し、結合する正規分布の個数を変えながらＥＭアルゴリズムによる混合密度分布を推定する。そして多次元の特徴量空間での各点についてベイズの定理に基づいて疵種を判断し識別境界を設定すると同時に確からしさ（情報のエントロピー）も計算する。
【００７２】
計算した結果は、それぞれ識別境界テーブル７と確からしさテーブル８に書き込む。識別境界と識別の確からしさは、計算の繰り返しごとに結果表示装置１０に表示する。解の収束は結果表示装置１０でのｘ、ｙ、ｚの各軸で切った複数の断面の観察と、それぞれのカテゴリについての対数尤度の変化によって判断し、オペレータがこれらを見て分類できたと判断したなら停止ボタンを押して計算を終了する。対数尤度の変化量が所定値以下となったら、分類できたと判断して自動的に計算を終わらせるようにしてもよい。
【００７３】
解が収束したならば最終的な識別境界設定としては３次元特徴量空間の内側の領域だけに上記の手法の結果を適用し、縁の近傍では従来からよく知られたｉｆ　ｔｈｅｎ　ｅｌｓｅ文によるツリー型の論理などを使用して識別境界を設定するか、上記の手法により得られた結果を外挿して識別境界を設定する。このようにすれば識別境界の設定の柔軟性および安定性の両者で効果的である。
【００７４】
表面欠陥検出手段が複数あり、それらの結果を統合して最終結果を求めたい際は、全ての測定手段について上記の演算を行い、それぞれの判断の確からしさを比較して最終的な判断結果とする。
【００７５】
オンラインでの疵種判別を高速化するために、上の学習計算による最終的な判断結果をテーブルに記憶しておくのが理想的であるが、特徴量の次元数が大きかったり特徴量空間が広かったりして計算機のメモリが足りない場合は平均ベクトル、共分散行列、結合係数などのパラメータだけを記憶しておいてもよい。
【００７６】
学習計算の後にオンラインで欠陥識別装置９にて光学式表面欠陥検出装置で検出した欠陥の多次元特徴量について疵種を判別する。最終的な判断結果をテーブルに記憶しておいた場合は画像処理装置９にて検出した欠陥の特徴量を得て、その特徴量がＮ次元特徴量空間のどの位置を占めるかを求め、その位置に対応する疵種を判別境界テーブル７から、確からしさを情報の確からしさテーブル８から読みこみ、結果表示装置１０にて帳票等に記録すれば、この系での識別は一瞬で終わる。
【００７７】
別の表面欠陥検出手段によっても同様の手法で識別境界が得られているならば、多次元の特徴量空間での各点について全ての測定手段から得られた識別境界テーブルと確からしさテーブルを参照して対応する疵種を求め、それらの判断の確からしさを比較し、最も確からしいものを最終的な判断結果とする。
【００７８】
平均ベクトル、共分散行列、結合係数といったパラメータだけを記憶しておいた場合は、検出した欠陥の特徴量について、それらのパラメータから各表面欠陥検出手段での疵種と確からしさを計算した後、すべての表面欠陥検出手段にて確からしさを比較し、最も確からしさの高い判断結果を最終的な判断結果として出力する。
【００７９】
以上に述べた疵種識別境界設定方法において、同様の特徴量空間を占める複数の疵種についての新たなデータ列がオペレータの目視判断による分類と共に学習データ保存装置５に蓄積され、また操業状態の変化などにより疵種識別境界設定を修正する必要が生じたときは、学習演算装置６内のメモリに保存しておいた前回の混合密度分布を特徴づける平均ベクトル、共分散行列、結合係数などのパラメータを初期値としてＥＭアルゴリズムによって更新し、疵種識別境界設定を修正する。
【００８０】
【発明の効果】
以上説明したように、本発明によれば、特殊な技術をもたないオペレータでも簡便に複雑な欠陥分類が実現され、労力が大幅に削減されるとともに操業状態の変化に迅速に対応した評価が可能となるような表面欠陥検出における疵種分類境界設定方法、及びこの方法を利用した表面欠陥検出方法を提供することができる。
【図面の簡単な説明】
【図１】本発明の実施の形態の１例における表面欠陥検出における、疵種分類境界設定方法の学習計算の流れを示すフローチャートである。
【図２】図１に示したフローチャートにおけるＥＭアルゴリズムの詳細を示す図である。
【図３】識別境界と確からしさの表示例示す図である。
【図４】２次元のデータを１つの正規分布と５つの正規分布でクラスタリングした結果を示す図である。
【図５】混合密度分布の推定による手法をＮ次元特徴量空間の縁近傍を除いた内部に限定して適用した例を示す図である。
【図６】Ｎ次元特徴量空間の縁の近傍で、クラスタリングのためのデータの不足に起因して混合密度分布の信頼性が低くなることを示す図である。
【図７】本発明を鉄鋼分野の熱間圧延プロセスに光学式表面欠陥検出に適用した例における装置構成を示す図である。
【図８】Ｎ次元正規分布の消滅と分裂の様子を示す図である。
【図９】図２に示したフローチャートにおけるＮ次元正規分布の消滅処理の詳細を示す図である。
【図１０】図２に示したフローチャートにおけるＮ次元正規分布の分裂処理の詳細を示す図である。
【符号の説明】
１：カメラコントローラ
２：画像メモリ
３：画像処理装置
４：学習データ入力装置
５：学習データ保存装置
６：学習演算装置
７：識別境界テーブル
８：確からしさテーブル
９：欠陥識別装置
１０：結果表示装置
１１：設定盤[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention provides a method of setting a flaw type classification boundary for classifying defects into categories based on multidimensional feature amounts obtained in surface defect detection based on these, and a surface defect detection method using this method It is about.
[0002]
[Prior art]
In an apparatus for detecting defects generated on the surface of an industrial product by an optical method (for example, Patent Document 1), it is extremely important to classify the flaw type and the degree of harm in quality assurance to a customer. Therefore, conventionally, various characteristic amounts such as a defect length, a defect width, and a density integrated value are extracted from an image signal obtained by imaging, and a tree diagram using an if @ then @ else statement is generated based on these characteristic amounts. The defect has been classified according to the logic or table. Further, a method of applying a neural network to such a classification has been studied (Patent Document 2).
[0003]
[Patent Document 1]
JP 2001-24089 A
[Patent Document 2]
JP-A-6-94643
[Non-patent document 1]
"Introduction to Cluster Analysis" (Morikita Publishing Co., Ltd.)
[0004]
[Problems to be solved by the invention]
However, in a method of classifying defects by a logic or a table of a tree diagram based on an if @ then @ else statement, it is difficult to classify based on complex identification boundaries, and an operator who has a special technique for analyzing and coding collected data. It takes great effort. In addition, since such a determination logic often needs to be readjusted due to a change in the operation state, it is essential to simplify and speed up this processing. In addition, when there are a plurality of surface defect detection means and the judgment results are different, there is confusion as to which one to trust.
[0005]
Further, in the method of applying a neural network to such classification, although complicated classification is possible in principle, a large amount of data and time are required for learning, and it cannot be said that it is simple and quick.
[0006]
The present invention has been made in view of such circumstances, and by adopting a novel method, a complex defect classification can be easily realized even by an operator who does not have a special technique, and labor is greatly reduced. An object of the present invention is to provide a method for setting a flaw type classification boundary in surface defect detection, which enables an evaluation to quickly respond to a change in an operation state, and a surface defect detection method using the method.
[0007]
[Means for Solving the Problems]
A first means for solving the above problem is to collect flaw data for learning, and assume a mixed density distribution of an N-dimensional normal distribution for each of a plurality of flaw type categories occupying an N-dimensional feature amount space. Then, parameters such as an average vector, a covariance matrix, and a coupling coefficient that characterize this are determined by convergence by an EM algorithm so as to be fitted to the learning flaw data, and the determined N-dimensional normal for each flaw category is determined. A flaw type classification boundary setting method for surface defect detection, comprising a step of determining a flaw type classification boundary in surface defect detection using Bayes' theorem based on the mixture density distribution of the distribution. ).
[0008]
A second means for solving the above-mentioned problem is the first means, wherein a new data sequence is obtained for a plurality of flaw type categories occupying the same feature amount space, and the flaw type classification boundary setting is performed. When parameters need to be corrected, parameters updated by EM algorithm convergence are determined using parameters such as the mean vector, covariance matrix, and coupling coefficient that characterize the previous mixture density distribution stored in memory as initial values. And correcting the flaw type classification boundary in the surface defect detection using Bayes' theorem based on the determined mixed density distribution of the N-dimensional normal distribution for each flaw category. (Claim 2).
[0009]
A third means for solving the above-mentioned problem is the first means or the second means, wherein a parameter is determined by convergence by an EM algorithm and, at the same time, information is obtained from a probability at each point in a feature space. The method according to claim 3, further comprising a step of calculating entropy and outputting a certainty of the judgment.
[0010]
A fourth means for solving the above-mentioned problem is that, for an area where the flaw type classification boundary is not stable, the flaw in the surface defect detection according to any one of the first means to the third means is provided. A flaw type classification boundary is set by the type classification boundary setting method, and in a region where the remaining flaw type classification boundary is stable, the obtained flaw type classification boundary is extrapolated to obtain a flaw type classification boundary. This is a flaw type classification boundary setting method (claim 4).
[0011]
A fifth means for solving the above-mentioned problem is that, in a region where the flaw type classification boundary is clear and stable, the flaw type classification boundary is determined by any conventionally used method, and the other areas are determined. A flaw type classification boundary setting method in surface defect detection, wherein the flaw type classification boundary is set by the flaw type classification boundary setting method in any one of the first means to the third means. (Claim 5).
[0012]
A sixth means for solving the above-mentioned problem is that, in an area near an edge of an N-dimensional feature amount space, a flaw type classification boundary is determined by an arbitrary method conventionally used, and only for other areas, A flaw type classification boundary setting method for surface defect detection, wherein the flaw type classification boundary is set by the flaw type classification boundary setting method for surface defect detection according to any one of the first means to the third means. ).
[0013]
A seventh means for solving the above-mentioned problem is the first means to the sixth means, including a step of eliminating when an N-dimensional normal distribution to be combined is not sufficiently spread in an N-dimensional space. (Claim 7).
[0014]
Eighth means for solving the above-mentioned problem is the first means to the seventh means, which includes a step of dividing when an N-dimensional normal distribution to be combined is widely spread in an N-dimensional space. (Claim 8).
[0015]
A ninth means for solving the above problem is to set a flaw type classification boundary by a flaw type classification boundary setting method in surface defect detection according to any one of the first means to the eighth means. A surface defect detection method (claim 9) including a step of identifying a type of the obtained defect according to a type classification boundary.
[0016]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings. In the present embodiment, the surface defect inspection apparatus generally includes an imaging system, an image memory, and a processing system. An imaging system consisting of a light source and a camera controller stores the state of the surface to be inspected in an image memory, and calculates various features such as the length, width, and integrated density of defects detected by the image processing device. I do. In the following description, for example, u = (u₁, U₂, U₃, ..., u_N), A symbol defined somewhere as a set of its elements indicates a vector.
[0017]
As shown in the processing flow of FIG. 1, the flow of learning calculation of the flaw type classification boundary setting method in the surface defect detection according to the present embodiment is input of feature amount data for each flaw type (steps S11 to S14), and for each flaw type. (Step S15), determination of an identification boundary based on Bayes' theorem (step S16), and calculation of certainty (step S17). .
[0018]
First, an N-dimensional feature amount pattern u = (u) of a detected learning defect via the learning data input unit from the image processing apparatus.₁, U₂, U₃, ..., u_N), Eg the length u of the defect₁, Defect width u₂, Concentration integrated value u₃… Is a flaw type α visually observed by the operator_kIs input to the learning data storage device together with the determination result (step S11). It is determined whether the data obtained in this manner is sufficient or sufficient for performing the calculation for setting the flaw type classification boundary (step S12). This can be realized by a method of setting a predetermined number in advance and determining whether or not the obtained learning defect data has reached this number.
[0019]
If the collected learning data is not sufficient, the process returns to step S11 to repeat the collection of learning data. If a sufficient amount of learning data is obtained, the process proceeds to step S13. The identification boundary setting calculation can be started. Then, in step S14, it is determined whether or not there is a discriminating boundary setting calculation start command. The identification boundary setting calculation start command may be instructed by a human using a push button or setting, or may be automatically issued at the same time that the identification boundary setting calculation can be started.
When an identification boundary setting calculation start command is issued, the process proceeds to step S15, and the parameters of the mixture density distribution are calculated for each flaw obtained as learning data.
[0020]
First, α_kTo generate an N-dimensional feature amount pattern u P (u | α_k) Is obtained by stochastic clustering. The feature amount pattern obtained by the surface defect meter often has a complicated distribution, and as shown in FIG. 4A, P (u | α_k) Cannot be obtained and clustered, so that a plurality of N-dimensional normal distributions are approximated as a mixture density distribution as shown in FIG. That is, the average vector μ_{k, m}And the covariance matrix σ_{k, m}Is the N-dimensional normal distribution characterized by φ (u: μ_{k, m}, Σ_{k, m}), Weighted by the determined number and combined, and expressed as in the following equation.
[0021]
Here, k indicates the category of the kth flaw, and m indicates the parameter of the mth N-dimensional normal distribution.
P (u | α_k) = Σζ_{k, m}φ (u: μ_{k, m}, Σ_{k, m})
The range of Σ is the range of the number of mixed density distributions to be combined (Σ for m). Here, the number of N-dimensional normal distributions to be combined is appropriately determined by looking at the distribution of defect information. The N-dimensional normal distribution is
[0022]
(Equation 1)

[0023]
It is expressed as Here, the average vector μ of each N-dimensional normal distribution characterizing the mixture density distribution_{k, m}And the covariance matrix σ_{k, m}, And coupling coefficient ζ_{k, m}Is obtained by repeated calculation using the EM algorithm.
[0024]
The mixture density distribution model and the EM algorithm are described in detail in, for example, "Introduction to Cluster Analysis" (Morikita Shuppan Co., Ltd.) (Non-Patent Document 1). Incomplete data sequence x = (x of the N-dimensional feature pattern observed by the surface defect inspection apparatus₁, X₂, X₃, ..., x_NX), And a virtual complete data sequence y = (y₁, Y₂, Y₃, ... y_NX)think of. Then, a partial incomplete data sequence x = (x₁, X₂, X₃, ..., x_NX) Is observed. Also, complete data series y = (y₁, Y₂, Y₃, ... y_NX) Is assumed to follow the probability density function f (y | Φ) with parameters, and the logarithm of f (y | Φ) under the conditions x and Φ ′, that is, the expected value of the log likelihood function Is
Q (Φ | Φ ′) = E (logf (y | Φ) | x, Φ ′)
Becomes
[0025]
At this time, the EM algorithm is as follows.
{Circle around (1)} First, the initial estimated values of the parameters for each of the N-dimensional normal distributions and their combined state are represented by Φ⁽⁰⁾And t = 0.
(2) Conditional expected value (Expection) Q (Φ | Φ^(T)) Is calculated.
(3) Q (Φ | Φ^(T)) To maximize ()^~Ask for. Ie
Φ^~= Arg @ max_ΦQ (Φ | Φ^(T))
(4) t = t + 1, Φ^(T)= Φ^~And (2) and (3) are repeated until the solution converges.
[0026]
More specifically, the complete data sequence is information about which of the N-dimensional normal distributions to combine for each of the N-dimensional feature pattern sequence (incomplete data sequence) observed by the surface defect inspection apparatus. It shall be added.
[0027]
Here, a data series x of NX N-dimensional feature amount patterns₁, X₂, X₃, ..., x_NXIs given, if the number of N-dimensional normal distributions to be combined is NG, the repetitive calculation of the mixture distribution determination is as shown in FIG. FIG. 2 shows step S15 in FIG. 1 in more detail. Here, the mean vector of the z-th normal distribution among the N-dimensional normal distributions to be combined is μ_z ^(T), The covariance matrix is σ_z ^(T)Coupling coefficient_z ^(T)And
[0028]
First, in step S151, an initial value of a parameter (coupling coefficient, mean vector, covariance matrix) is input (when a previous value is input, it will be described later). Subsequently, these parameters are obtained by converging by repeated calculation.
[0029]
That is, the following calculations (1) → (2) → (3) → (4) → (5) shown in FIG. 2 are repeated until convergence.
(1) Feature amount pattern x_jIs the probability ω generated from the z-th (z = 1, 2,..., NG) N-dimensional normal distribution_zj ^(T)Is calculated (step S152).
[0030]
(Equation 2)

[0031]
The sum Σ is performed for all N-dimensional normal distributions to be combined. Here, the equation of the N-dimensional normal distribution is as shown in equation (1).
(2) Q (Φ | Φ^(T)) To maximize the coupling coefficient ζ for the z-th N-dimensional normal distribution._zIs updated (step S153).
ζ_z ^{(T + 1)}= ｛Σω_zj ^(T)｝ / NX
The sum Σ is performed for all NX data strings.
(3) Q (Φ | Φ^(T)) To maximize the mean vector μ for the z-th N-dimensional normal distribution._zIs updated (step S154).
μ_z ^{(T + 1)}= ｛Σω_zj ^(T)x_j｝ / (NX ・ ζ_z ^{(T + 1)})
The sum Σ is performed for all NX data series.
(4) Q (Φ | Φ^(T)) To maximize the covariance matrix σ for the z-th N-dimensional normal distribution_zIs updated (step S155).
[0032]
(Equation 3)
σ_z ^{(T + 1)}
= ｛Σω_zj ^(T)(X_j−μ_z ^{(T + 1)}) (X_j−μ_z ^{(T + 1)})^t｝ / (NX ・ ζ_z ^{(T + 1)})
[0033]
The sum Σ is performed for all NX data series.
(5) As necessary, an N-dimensional normal distribution annihilation process (step S156) and an N-dimensional normal distribution splitting process (step S157) are performed. These processes are not necessarily performed, and will be described later in detail.
(6) mixture density distribution P (x | Φ^{(T + 1)}) Is calculated (step S158).
[0034]
(Equation 4)
P (x | Φ^{(T + 1)}) = Σζ_zφ (x: μ_z ^{(T + 1)}, Σ_z ^{(T + 1)})
[0035]
The sum Σ is performed for all NG N-dimensional normal distributions.
Then, it is determined whether or not the calculation has converged (step S159). This convergence can be determined, for example, when a difference between one of the parameters (coupling coefficient, mean vector, covariance matrix) from the previous time is equal to or smaller than a predetermined value.
[0036]
In the above calculation, if the number of N-dimensional normal distributions to be combined is large, it is possible to cope with a complicated discrimination boundary. Since the number of N-dimensional normal distributions to be combined does not necessarily need to be the same for each type of flaw, it may be individually selected from the required complexity and convergence of the discrimination boundary.
[0037]
The initial value of the EM algorithm set in the first learning is based on the assumption that the data follows one multidimensional normal distribution, the mean vector and the covariance matrix are obtained, and the initial value of the covariance matrix is the covariance matrix itself at that time. , The initial value of the average vector may be determined by the number of normal distributions that generate and combine random numbers according to the multidimensional normal distribution inside the computer. Also, the initial value of the coupling coefficient may be the same 1 / NG (NG is the number of normal distributions to be coupled) for each multidimensional normal distribution. In the convergence calculation of the EM algorithm, depending on how to select an initial value, there is a possibility that a local optimal solution is captured instead of a global optimal solution. Therefore, in the first learning, the distribution of the identification boundaries and the certainty displayed on the output screen is checked, and if it is determined that the defect type identification boundaries have not been correctly performed, another initial value is given and the calculation is performed. I need to do it again.
[0038]
As this method, a method of manually re-entering the initial value or a method of generating a random number again can be considered, but for the convergence calculation of the EM algorithm in the first learning, try several initial values, and Choosing the best solution inside does not require much effort.
[0039]
If the determinant of the covariance matrix of the N-dimensional normal distribution is close to 0 in the estimation of the mixture density distribution by the EM algorithm, an overflow occurs when calculating the expression (1). The fact that the determinant of the covariance matrix is close to 0 indicates that the correlation is extremely strong between certain variables of the N-dimensional normal distribution, and in the worst case, the distribution is point-like. Generation of such an irregular N-dimensional normal distribution can also be considered as convergence to a certain local optimal solution.
[0040]
To avoid this, it is conceivable to start from another initial value and converge to an appropriate solution, but it is also conceivable that such a distorted N-dimensional normal distribution is eliminated (FIG. 8A). ). Here, when a distorted one of the plurality of N-dimensional normal distributions is deleted, the number of N-dimensional normal distributions to be combined may decrease.
[0041]
It is also conceivable to split the N-dimensional normal distribution, which has spread too much in the N-dimensional feature space, in reverse (FIG. 8B). By doing so, an operation of re-clustering the learning data that has been roughly clustered at first can be added, and finally the accuracy of setting the identification boundary can be improved. Here, when a large one of the plurality of N-dimensional normal distributions is divided, the number of N-dimensional normal distributions to be combined may increase.
[0042]
By repeating the disappearance and division of the N-dimensional normal distribution in this manner, an appropriate number can be obtained without paying much attention to the initially set number. When and how to eliminate or split can be determined from the eigenvalues of the covariance matrix of the N-dimensional normal distribution and the orthonormal eigenvectors.
[0043]
The "N-dimensional normal distribution annihilation process" (S156) of FIG. 2 will be described with reference to the flowchart of FIG. 9, or if the spread of the N-dimensional normal distribution in each axial direction is small (S1563), the number of N-dimensional normal distributions to be combined will be described. Is reduced by one (S1564), the coupling coefficient is distributed to all the remaining normal distributions (S1565), and the corresponding normal distribution is eliminated (S1566).
[0044]
Each N-dimensional normal distribution has its covariance matrix σ_{k, m}The shape becomes a hyperellipsoid spreading in the direction of N orthonormal eigenvectors. And the covariance matrix σ_{k, m}The set of eigenvalues and orthonormal eigenvectors of (λ₁, A₁), (Λ₂, A₂), (Λ₃, A₃), ..., (λ_N, A_N), The variation corresponding to the standard deviation in each axis direction is
b₁= Λ₁ ^1/2a₁
b₂= Λ₂ ^1/2a₂
b₃= Λ₃ ^1/2a₃
・ ・ ・ ・ ・ ・ ・ ・
b_N= Λ_N ^1/2a_N
Therefore, any one of these absolute values is determined by a predetermined threshold (D_minIf it is smaller, it should be extinguished.
[0045]
If the α-th N-dimensional normal distribution disappears
ζ₁ ^(T)→ ζ₁ ^(T)= Ζ₁ ^(T)+ 1 / (NC-1) ζ_α ^(T)
ζ₂ ^(T)→ ζ₂ ^(T)= Ζ₂ ^(T)+ 1 / (NC-1) ζ_α ^(T)
ζ₃ ^(T)→ ζ₃ ^(T)= Ζ₃ ^(T)+ 1 / (NC-1) ζ_α ^(T)
・ ・ ・ ・ ・ ・ ・ ・
ζ_α ^(T)→ ζ_α ^(T)= 0
・ ・ ・ ・ ・ ・ ・ ・
ζ_N ^(T)→ ζ_N ^(T)= Ζ_N ^(T)+ 1 / (NC-1) ζ_α ^(T)
Then, the mixture density distribution is estimated by the EM algorithm using (NC-1) N-dimensional normal distributions excluding the α-th N-dimensional normal distribution.
[0046]
The "split processing of N-dimensional normal distribution" (S157) of FIG. 2 will be described with reference to the flowchart of FIG. 10. When the spread of each axis of the N-dimensional normal distribution is large (S1573), the number of N-dimensional normal distributions to be combined is determined. The number is increased by one (S1574), and the splitting of the coupling coefficient of the corresponding N-dimensional normal distribution (S1575), the splitting of the covariance matrix (S1576), and the splitting of the average vector (S1577) are performed according to the rules described later.
[0047]
Each N-dimensional normal distribution has its covariance matrix σ_{k, m}The shape becomes a hyperellipsoid spreading in the direction of N orthonormal eigenvectors. And the covariance matrix σ_{k, m}The set of eigenvalues and orthonormal eigenvectors of (λ₁, A₁), (Λ₂, A₂), (Λ₃, A₃), ..., (λ_N, A_N), The variation corresponding to the standard deviation in each axis direction is
b₁= Λ₁ ^1/2a₁
b₂= Λ₂ ^1/2a₂
b₃= Λ₃ ^1/2a₃
・ ・ ・ ・ ・ ・ ・ ・
b_N= Λ_N ^1/2a_N
Therefore, any one of these absolute values is determined by a predetermined threshold (D_maxIf it is bigger, it should be split.
[0048]
If the β-th N-dimensional normal distribution is split, then the coupling coefficient ζ_β1 ^(T), Ζ_β2 ^(T)Is the original coupling coefficient ζ_β ^(T)To
ζ_β ^(T)→ ζ_β1 ^(T)= 1 / 2ζ_β ^(T), Ζ_β2(T) = 1 / 2ζ_β ^(T)
And shall be equally divided into
Also, the covariance matrix σ_β1 ^(T), Σ_β2 ^(T) Is
σ_β ^(T)→ σ_β1 ^(T)= Hσ_β ^(T), Σ_β2 ^(T)= Hσ_β ^(T)
The original covariance matrix σ_β ^(T)Multiplied by a constant (h <1),
The mean vector μ_β1 ^(T), Μ_β2 ^(T)Is
μ_β ^(T)→ μ_β1 ^(T)= Μ_β ^(T)+ Δμ_β, Μ_β2 ^(T)= Μ_β ^(T)−Δμ_β
The original mean vector μ_β ^(T)Direction b to split from_β(Δμ_βAnd b_βAre parallel).
That is, the corresponding orthonormal eigenvector is represented by b_β= (B_βx, B_βy)
[0049]
(Equation 5)

[0050]
And
Δμ_β= Δk (cos θ_β, Sin θ_β)
Becomes
[0051]
The annihilation and division are repeated as described above, and the mixed density distribution is estimated by the EM algorithm by changing the number of normal distributions to be combined. Here, the allowable range of the spread of the N-dimensional normal distribution (D_min, D_max) Is too narrow, the convergence of the estimation of the mixture density distribution by the EM algorithm may be poor. It is preferable to perform an operation on collected data or simulation data in advance to confirm the convergence and the accuracy of the result.
[0052]
When the flaw type identification boundary setting must be corrected due to changes in the operating state, etc., and new learning data for that purpose is accumulated, the average vector, covariance matrix, coupling coefficient characterizing the previous mixture density distribution It is preferable to update the parameters such as the initial values by the EM algorithm to correct the flaw type identification boundary setting. In many cases of flaw type classification, the correction of the flaw type classification boundary due to a change in the operation state or the like is not large, so that the possibility of being caught by the local optimum solution is small.
[0053]
In practice, the calculation is repeatedly performed in parallel for each flaw type. And the flaw type α_kIs P (u | α)._k) = P (u | Φ_k ^{(T + 1)}) And perform the following calculation.
[0054]
The convergence of the repeated calculation by the EM algorithm is confirmed by checking the distribution of the identification boundaries and the certainty displayed on the screen in the case of two dimensions, but as represented by (2) in FIG. 3 in the case of three or more dimensions. It is confirmed by observing a proper two-dimensional cross section and changing the log likelihood in each category.
[0055]
The log likelihood is P (x | Φ_k ^{(T + 1)})When
L = ΣlogP (x_i│Φ_k ^{(T + 1)})
(Total is all observation data x_iAbout)
And increases as the convergence proceeds, and does not change when the convergence is completed. Alternatively, the change of each component of the mean vector and the covariance matrix may be observed instead of the log likelihood. In the case of a large dimension, it is preferable to set the upper discrimination boundary after reducing the dimension by Karhunen-Loeve expansion or principal component analysis.
[0056]
A decision is made based on Bayes' theorem in order to classify the defect types and harmfulness (step S16). If the defect type or harmfulness to be classified is broadly referred to as a flaw type, when a feature amount pattern u is given, the flaw type is α_kIs defined as P (α_k| U), flaw type α_kIs P (u | α)._k), Flaw type α_kIs defined as P (α_k), If the occurrence probability of the feature amount pattern u is P (u), Bayes' theorem gives
P (α_k| U) = P (u | α)_k) P (α_k) / P (u)
[0057]
Here, P (u) is P (u | α_k) Is determined by the EM algorithm for the mixture distribution described above. In addition, the prior probability P (α_k) Is P (α) if data collection can be assumed to be random sampling._k) = N_k/ ΣN_l(N_lIs the number of data obtained for the series l, and the sum is taken for all the series.) It is natural to use P (α) when such an assumption does not hold for data collection._k) = 1 / c (c is the number of series), and may be intentionally the same for all flaw types, or it is desired to more accurately discriminate a specific flaw type such as a loss caused by a specific misclassification being extremely large. In that case, P (α_k) May be overestimated with respect to other prior probabilities. That is, P (α_k) Functions as an adjustment factor.
[0058]
Since P (u) does not depend on the type of flaw,
P '(α_k| U) = P (u | α)_k) P (α_k)
And the most P '(α_kAn identification boundary is set such that the one with the higher | u) is regarded as a flaw type, and is written in a table as shown in FIG. This table may be drawn on a CRT or the like for each calculation and displayed to the operator.
[0059]
Next, the entropy of the information is obtained for the entire feature amount space as the certainty of the judgment and is written in the table (step S17). The entropy of information is a quantity expressing the predictability of an event._l(However, ΣP_l= 1, the sum is taken for all events that occur)
H = -ΣP_llogP_l(Summation is for all events that occur)
It is expressed by the formula. When all events occur with the same probability and are the least predictable, the entropy of information is maximized, and only one event occurs (ie, its probability of occurrence = 1) and no other event occurs (ie, its probability of occurrence = 0). When it is easiest to predict, the entropy of the information has a minimum value of zero.
[0060]
This amount represents the certainty of the judgment in this case.The entropy of the information is high in the part where the judgment is lost such as near the discrimination boundary, and the information is low in the place where the certainty is low, such as inside the judgment boundary. Entropy is lower and certainty is higher. This time, the multi-dimensional feature pattern u_nIs given, P ′ (α calculated above from Bayes' theorem_k| U_n)
P "(α₁| U_n) = P ′ (α₁| U_n) / ΣP ’(α_l| U_n)
P "(α₂| U_n) = P ′ (α₂| U_n) / ΣP ’(α_l| U_n)
P "(α₃| U_n) = P ′ (α₃| U_n) / ΣP ’(α_l| U_n)
…
P "(α_k| U_n) = P ′ (α_k| U_n) / ΣP ’(α_l| U_n)
…
ΣP ”(α_l| U_n) = 1 (sum is taken for all flaw types), then calculate the entropy of the information by the above equation, and write the certainty in the table as shown in FIG. 3 (2). If this entropy is displayed, it is possible to determine the certainty of the determination.
[0061]
In step S18, it is determined whether the discrimination boundary is appropriate by looking at the amount of entropy or looking at the actual discrimination boundary on the screen. If it is determined that the discrimination boundary is inappropriate, the process proceeds to step S20, the initial value of the EM algorithm is changed, and the process returns to step S15 to perform the calculation again. If it is determined to be appropriate, the process proceeds to step S19, where the parameters such as the average vector, the covariance matrix, the coupling coefficient, and the identification boundary and certainty are stored in a table, and the processing is terminated.
[0062]
Here, the certainty is certainty for each point in the N-dimensional feature amount space. For example, when the N-dimensional feature amount space and the feature amount pattern u are classified into the defect category A, it indicates the certainty that the point u is the defect category A. Specifically, an amount such as the entropy described above corresponds to this.
[0063]
In the above description, the direction in which the identification boundary is set over the entire N-dimensional feature amount space has been described. However, an N-dimensional feature amount pattern input as learning data is observed in a finite range. That is, the N-dimensional feature amount space is a bounded area. Since no learning data is observed outside the N-dimensional feature amount space, data for clustering runs short near the edges of the N-dimensional feature amount space. Therefore, the reliability of the mixture density distribution estimated as shown in FIG. 6 decreases.
[0064]
Therefore, as shown in FIG. 5, it is often better to apply the method based on the estimation of the mixture density distribution to the inside of the N-dimensional feature amount space except for the vicinity of the edge. At this time, in the vicinity of the edge of the N-dimensional feature amount space, an identification boundary is set by using a tree-type logic based on a conventionally well-known if @ then @ else statement, or obtained by estimating a mixture density distribution. The identification result is extrapolated to set an identification boundary. However, even in these cases, the learning data sequence generally employs all the learning data obtained over the entire feature amount space. Regarding the portion where the identification boundary is set irrespective of the mixture density distribution near the edge of the N-dimensional feature amount space, the entropy of information as an index of certainty is set to 0 (region where the judgment is most certain). This is effective in both flexibility and stability of setting the identification boundary.
[0065]
At the time of online defect identification, a flaw type is determined based on a point on the identification boundary table corresponding to a feature amount pattern observed by the image processing apparatus. If there are a plurality of surface defect detection means, the corresponding points are similarly obtained on the likelihood table and compared, and the highest one is determined as the reliable final defect type classification result.
[0066]
As described above, a new data sequence for the flaw type occupying the feature amount space is stored in the learning data storage device together with the classification based on the visual judgment of the operator, and the flaw type identification boundary setting is performed due to a change in the operation state. When correction is required, parameters such as the average vector, covariance matrix, and coupling coefficient characterizing the previous mixture density distribution in the learning operation device are updated as initial values by the EM algorithm to correct the flaw type identification boundary setting. . This corresponds to inputting the previous value of the parameter in step S151 in FIG.
[0067]
Hereinafter, an example in which the present invention is applied to an optical surface defect detection device in a hot rolling process in the steel field will be described.
[0068]
FIG. 7 is an outline of the hot rolling process in this embodiment. The rolling roll has 7 steps, and the first roll on the entrance side is No. 1, the last roll is No. 7 rolls. A slab having a temperature of about 1300 ° C. was inserted from the entry side, and was gradually rolled. Rolled to about 1 mm at the exit of 7 rolls. Here, the moving speed of the steel sheet is about 20 m / sec, and the width of the steel sheet is 1000 mm.
[0069]
As a light source for imaging, for example, the light of a 250 W metal halide lamp is formed into a linear shape using a bundle fiber and condensed by a cylindrical lens to illuminate the surface of the steel plate. Light is received using a line sensor camera. Since the number of pixels of the camera is 1024 and a board width of 1000 mm is observed, the pixel resolution in the horizontal direction is 1 mm. The camera controller 1 receives a signal from a rotary encoder provided on the roll, and sends a trigger signal to the camera every time the steel plate advances by 1 mm. Since the camera sequentially outputs 1024 lines of image signals only when a trigger signal is input over the entire length of the coil, one pixel is always stored in the image memory in the image memory even if the operating speed of the steel plate changes. This corresponds to an image of 1 mm in length and width. The obtained information is stored in the image memory 2.
[0070]
The flow of the process will be described with reference to FIG. 7. First, a multidimensional feature amount obtained from a defect candidate of image data is obtained by the image processing device 3 and input to the learning data input device 4 together with flaw type classification by visual judgment of the operator. Are stored in the learning data storage device 5. Here, the multidimensional feature amount is, for example, three dimensions of a defect width, a defect length, and a density integrated value.
[0071]
When data of a predetermined number or more has been collected for all flaw types, when the operator presses the learning start button on the setting panel 11, the learning operation device 6 calculates the flaw type classification. First, for each flaw type, a mixture density distribution is obtained by the EM algorithm and clustered stochastically. The number of N-dimensional normal distributions expected from experience is initially set, and the estimation calculation of the mixture density distribution is started. If the N-dimensional normal distribution to be combined in the operation is not sufficiently spread in the N-dimensional space, it is deleted. As described above, the extinction and division are repeated, and the mixed density distribution is estimated by the EM algorithm while changing the number of normal distributions to be combined. Then, for each point in the multi-dimensional feature quantity space, a flaw type is determined based on Bayes' theorem, and an identification boundary is set, and at the same time, certainty (entropy of information) is calculated.
[0072]
The calculated results are written in the identification boundary table 7 and the likelihood table 8, respectively. The identification boundary and the certainty of the identification are displayed on the result display device 10 every time the calculation is repeated. The convergence of the solution can be determined by observing a plurality of cross sections cut along each of the x, y, and z axes on the result display device 10 and changing the log likelihood for each category. If it is determined that the calculation has been completed, the stop button is pressed to end the calculation. If the amount of change in the log likelihood is equal to or less than a predetermined value, it may be determined that the classification has been made, and the calculation may be automatically terminated.
[0073]
If the solution converges, the result of the above method is applied only to the region inside the three-dimensional feature space as the final discrimination boundary setting, and a tree by the well-known if @ then @ else statement is used near the edge. An identification boundary is set using type logic, or the result obtained by the above method is extrapolated to set an identification boundary. This is effective in both flexibility and stability of setting the identification boundary.
[0074]
If there are multiple surface defect detection means and you want to integrate the results to obtain the final result, perform the above calculation for all the measurement means, compare the certainty of each judgment, and compare it with the final judgment result. I do.
[0075]
In order to speed up online flaw type discrimination, it is ideal to store the final judgment result obtained by the above learning calculation in a table, but the dimension number of the feature amount is large or the feature amount space is large. If the computer is too large to have enough memory, only parameters such as the mean vector, the covariance matrix, and the coupling coefficient may be stored.
[0076]
After the learning calculation, the defect type is determined for the multidimensional feature amount of the defect detected by the optical surface defect detection device by the defect identification device 9 online in the defect identification device 9. When the final judgment result is stored in the table, the feature amount of the defect detected by the image processing device 9 is obtained, and the position of the feature amount in the N-dimensional feature amount space is determined. If the flaw type corresponding to the position is read from the discrimination boundary table 7 from the certainty table 8 and recorded in a form or the like on the result display device 10, the identification in this system is completed in an instant.
[0077]
If the identification boundary is obtained by another surface defect detection means in the same manner, refer to the identification boundary table and the probability table obtained from all the measurement means for each point in the multidimensional feature space. Then, the corresponding flaw types are determined, the likelihoods of those determinations are compared, and the most probable one is determined as the final determination result.
[0078]
If only the parameters such as the average vector, the covariance matrix, and the coupling coefficient are stored, for the feature amount of the detected defect, after calculating the flaw type and the likelihood at each surface defect detection unit from those parameters, Probabilities are compared by all the surface defect detection means, and a judgment result with the highest probability is output as a final judgment result.
[0079]
In the flaw type identification boundary setting method described above, new data strings for a plurality of flaw types occupying the same feature amount space are stored in the learning data storage device 5 together with the classification by visual judgment of the operator, and the operation state When it becomes necessary to modify the flaw type identification boundary setting due to a change or the like, the average vector, the covariance matrix, the coupling coefficient, etc., which characterize the previous mixture density distribution stored in the memory of the learning operation device 6 are used. The parameters are updated as initial values by the EM algorithm, and the flaw type identification boundary setting is corrected.
[0080]
【The invention's effect】
As described above, according to the present invention, even an operator who does not have a special technology can easily realize a complicated defect classification, greatly reduce labor, and can quickly evaluate an operation in response to a change in operation state. It is possible to provide a flaw type classification boundary setting method in surface defect detection that can be performed, and a surface defect detection method using this method.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a flow of learning calculation of a flaw type classification boundary setting method in surface defect detection according to an example of an embodiment of the present invention.
FIG. 2 is a diagram showing details of an EM algorithm in the flowchart shown in FIG. 1;
FIG. 3 is a diagram showing a display example of identification boundaries and certainty.
FIG. 4 is a diagram showing a result of clustering two-dimensional data with one normal distribution and five normal distributions.
FIG. 5 is a diagram showing an example in which a technique based on estimation of a mixture density distribution is applied only to the inside of the N-dimensional feature amount space excluding the vicinity of an edge.
FIG. 6 is a diagram showing that near the edge of the N-dimensional feature amount space, the reliability of the mixture density distribution is reduced due to lack of data for clustering.
FIG. 7 is a diagram showing an apparatus configuration in an example in which the present invention is applied to optical surface defect detection in a hot rolling process in the steel field.
FIG. 8 is a diagram showing the state of disappearance and division of an N-dimensional normal distribution.
FIG. 9 is a diagram showing details of the disappearance processing of the N-dimensional normal distribution in the flowchart shown in FIG. 2;
FIG. 10 is a diagram showing details of a splitting process of an N-dimensional normal distribution in the flowchart shown in FIG. 2;
[Explanation of symbols]
1: Camera controller
2: Image memory
3: Image processing device
4: Learning data input device
5: Learning data storage device
6: Learning operation device
7: Identification boundary table
8: Probability table
9: Defect identification device
10: Result display device
11: Setting panel

Claims (9)

The flaw data for learning is collected, and for each of a plurality of flaw type categories occupying the N-dimensional feature amount space, an N-dimensional normal distribution mixture density distribution is assumed for each category, and an average vector, a covariance matrix, The parameters such as the coupling coefficient are determined by convergence by the EM algorithm so as to be fitted to the learning flaw data, and based on the determined N-dimensional normal distribution mixture density distribution for each flaw category, Bayes' theorem is calculated. A method for setting a flaw type classification boundary in surface defect detection, comprising using the step of determining a flaw type classification boundary in surface defect detection. 2. The flaw type classification boundary setting method for surface defect detection according to claim 1, wherein a new data string is obtained for a plurality of flaw type categories occupying the same feature amount space, and the flaw type classification boundary setting is corrected. When it is necessary to determine parameters updated by convergence by the EM algorithm with parameters such as an average vector, a covariance matrix, and a coupling coefficient that characterize the previous mixture density distribution stored in memory as initial values, A surface defect including a step of correcting a defect type classification boundary in surface defect detection using Bayes' theorem based on the determined mixed density distribution of the N-dimensional normal distribution for each category of each defect. A method for setting a flaw type classification boundary in detection. 3. A flaw type classification boundary setting method in surface defect detection according to claim 1 or 2, wherein parameters are determined by convergence by an EM algorithm, and at the same time, entropy of information is calculated from a probability at each point in a feature amount space. And outputting a likelihood of the judgment to determine the flaw type classification boundary in surface defect detection. For a region where the flaw type classification boundary is not stable, a flaw type classification boundary is set by the flaw type classification boundary setting method in surface defect detection according to any one of claims 1 to 3, and the remaining flaw type classification boundary is set. A flaw type classification boundary setting method characterized in that a flaw type classification boundary is obtained by extrapolating the obtained flaw type classification boundary in a region where the flaw type classification boundary is stable. In a region where the flaw type classification boundary is clear and stable, the flaw type classification boundary is determined by an arbitrary method conventionally used, and any one of claims 1 to 3 is determined only for other regions. A flaw type classification boundary setting method for surface defect detection, wherein the flaw type classification boundary is set by the flaw type classification boundary setting method for surface defect detection according to claim 1. In a region near an edge of the N-dimensional feature amount space, a flaw type classification boundary is determined by an arbitrary method conventionally used, and only the other region is determined according to any one of claims 1 to 3. A flaw type classification boundary setting method in surface defect detection, wherein a flaw type classification boundary is set by the flaw type classification boundary setting method in surface defect detection described above. 7. The flaw type classification boundary setting method for surface defect detection according to any one of claims 1 to 6, wherein the N-dimensional normal distribution to be combined is eliminated when the N-dimensional normal distribution is not sufficiently spread in the N-dimensional space. A method for setting a flaw type classification boundary in surface defect detection, the method comprising: 8. A flaw type classification boundary setting method in surface defect detection according to any one of claims 1 to 7, wherein the N-dimensional normal distribution to be combined is split when the N-dimensional normal distribution is greatly spread in an N-dimensional space. And a flaw type classification boundary setting method in surface defect detection. A flaw type classification boundary is set by the flaw type classification boundary setting method in surface defect detection according to any one of claims 1 to 8, and the type of the obtained defect is determined according to the flaw type classification boundary. A method for detecting a surface defect, comprising the step of identifying a surface defect.

Images