以下本発明の実施の形態を図面を参照して� ��明する。
 [第1の実施形態]
 図1は、本発明の一実施形態にかかる内燃機 関(以下「エンジン」という)及びその制御装� ��の全体構成図であり、例えば4気筒のエンジ ン1の吸気管2の途中にはスロットル弁3が配さ れている。スロットル弁3にはスロットル弁� �度THを検出するスロットル弁開度センサ4が� �結されており、センサ4の検出信号は、電子� ��御ユニット(以下「ECU」という)5に供給され� ��。

 燃料噴射弁6はエンジン1とスロットル弁3� ��の間かつ吸気管2の図示しない吸気弁の少し 上流側に各気筒毎に設けられており、各噴射 弁は図示しない燃料ポンプに接続されている と共にECU5に電気的に接続されて当該ECU5から� ��信号により燃料噴射弁6の開弁時間が制御さ れる。エンジン1の各気筒には、点火プラグ7� ��設けられており、点火プラグ7はECU5に接続� �れている。ECU5は、点火プラグ7に点火信号を 供給する。

 スロットル弁3の下流側には吸気圧PBAを検 出する吸気圧センサ8、及び吸気温TAを検出す る吸気温センサ9が設けられている。エンジ� �1の本体には、エンジン冷却水温TWを検出す� �冷却水温センサ10及び非共振型のノックセン サ11が装着されている。センサ8~11の検出信号 は、ECU5に供給される。ノックセンサ11として は、例えば5kHzから25kHzまでの周波数帯域の振 動を検出可能なものが使用される。

 吸気管2のスロットル弁3の上流側には、吸� �空気流量GAを検出する吸入空気流量センサ13� ��設けられており、その検出信号はECU5に供給 される。
 ECU5には、エンジン1のクランク軸(図示せず) の回転角度を検出するクランク角度位置セン サ11が接続されており、クランク軸の回転角� ��に応じた信号がECU5に供給される。クランク 角度位置センサ11は、エンジン1の特定の気筒 の所定クランク角度位置でパルス(以下「CYL� �ルス」という)を出力する気筒判別センサ、� ��気筒の吸入行程開始時の上死点(TDC)に関し� �定クランク角度前のクランク角度位置で(4気 筒エンジンではクランク角180度毎に)TDCパル� �を出力するTDCセンサ及びTDCパルスより短い� �定クランク角周期(例えば6度周期)で1パルス( 以下「CRKパルス」という)を発生するCRKセン� �から成り、CYLパルス、TDCパルス及びCRKパル� �がECU5に供給される。これらのパルスは、燃� ��噴射時期、点火時期等の各種タイミング制� ��、エンジン回転数(エンジン回転速度)NEの検 出に使用される。

 エンジン1は、吸気弁(図示せず)の弁リフ� ��量及び開角(開弁期間)を連続的に変更する� �1弁作動特性可変機構と、吸気弁を駆動する� ��ムの、クランク軸回転角度を基準とした作� ��位相を連続的に変更する第2弁作動特性可変 機構とを有する弁作動特性可変装置20を備え� ��いる。ECU5は、弁作動特性可変装置20にリフ� ��量制御信号及び作動位相制御信号を供給し� ��吸気弁の作動制御を行う。第1及び第2弁作� �特性可変機構の構成は、それぞれ例えば特� �2008-25418号公報及び特開2000-227013号公報に示� �れている。

 ECU5は、各種センサからの入力信号波形を 整形し、電圧レベルを所定レベルに修正し、 アナログ信号値をデジタル信号値に変換する 等の機能を有する入力回路、中央演算処理ユ ニット(以下「CPU」という)、該CPUで実行され� ��各種演算プログラム及び演算結果等を記憶� ��る記憶回路(メモリ)、燃料噴射弁6及び点火� ��ラグ7に駆動信号を供給する出力回路等から 構成される。

 次に本実施形態におけるノッキング判定� ��法の概要を説明する。図2(a)は、ノックセン サ11の出力信号波形を示し、図2(b)は、同図(a) の期間TSの波形を拡大して示す図である。本� ��施形態では、サンプリング周期20マクロ秒� �検出される50個のデータを対象として、高速 フーリエ変換による周波数成分分析を行う。 その周波数成分分析の結果が図2(c)に示され� �いる。図2(c)の縦軸はエネルギ強度PSであり� �本実施形態では、5kHzから25kHzまでの周波数� �域における1kHz毎の周波数(5,6,7,…,24,25kHz)に� �応するエネルギ強度PSが、クランク角度6度� �に算出される。

 これにより、周波数成分の時系列データ� ��図3(a)に示すように2次元マトリクス(以下「� ��ペクトル時系列マップ」という)として算出 される。スペクトル時系列マップは縦方向が 周波数f[kHz]であり、横方向がクランク角度(� �焼行程が開始する上死点後のクランク角度)C A[deg]である。ここで、縦方向及び横方向のイ ンデクスとしてそれぞれ「j」(j=0~20)及び「i� �(i=0~14)を用い、スペクトル時系列マップの要 素を強度パラメータKMAP(j,i)と表示する。例え ば強度パラメータKMAP(0,0)は、左下端の「34」� ��強度パラメータKMAP(0,14)は右下端の「31」、� ��度パラメータKMAP(20,0)は左上端の「56」、強� ��パラメータKMAP(20,14)は右上端の「30」に相当 する。なお強度パラメータKMAPは相対強度を� �す無次元量である。

 図3(a)のスペクトル時系列マップを二値化 すると図3(b)の二値化スペクトル時系列マッ� �が得られる。図3(b)のマップの要素である二� ��化強度パラメータNKMAP(j,i)は、強度パラメー タKMAP(j,i)が「50」以上であるとき「1」、「50� ��未満であるとき「0」をとる。

 図3(b)に示す二値化スペクトル時系列マッ プは、図30に破線で示すノッキングが発生し� ��ときに得られるマップである。一方図30に� �線で示す着座ノイズの二値化スペクトル時� �列マップは、図4に示すようになり、図3(b)に 示すノッキング発生時のマップと明確に相違 するため、着座ノイズをノッキング発生と誤 判定することがなく、正確なノッキング判定 を行うことができる。

 本実施形態ではさらに、ノッキングが発� ��していないと判定されたときの二値化スペ� ��トル時系列マップに基づいて、ノイズに対� ��する二値化スペクトル時系列マップ(以下単 に「二値化ノイズマップ」という)を学習演� �により生成し、二値化ノイズマップ上のノ� �ズ学習値を、検出データから算出される二� �化スペクトル時系列マップ上の対応する値� �ら差し引くことにより、ノイズの影響を除� �するようにしている。

 図5(a)は、学習演算により得られた二値化 ノイズマップの一例を示す図であり、クラン ク角度6度のタイミングで6~25kHzの帯域でノイ� ��成分が表れている。図5(b)は、図3(b)に示す� �値化スペクトル時系列マップから図5(a)に示� ��ノイズ成分を減算することによりノイズ除� ��処理を行った後の二値化スペクトル時系列� ��ップを示す。このようにノイズ除去処理を� ��うことにより、ノイズの影響を除き、より� ��確な判定が可能となる。

 また本実施形態では、検出データから得� ��れる二値化スペクトル時系列マップを、図6 (a)に示すマスタパターンマップと比較するこ とにより、ノッキングが発生したか否かを判 定するようにしている。マスタパターンマッ プは、ノッキングが発生したときに得られる 典型的な二値化スペクトル時系列マップに相 当するものである。

 上記比較は具体的には以下のように行わ� ��る。二値化スペクトル時系列マップ上の二� ��化強度パラメータNKMAP(j,i)と、マスタパター ンマップ上のマスタパラメータMMAP(j,i)との積 を積算することにより、強度積算値SUMKを算� �し、マスタパラメータMMAP(j,i)そのものの積� �値である基準積算値SUMMを算出し、基準積算� ��SUMMに対する強度積算値SUMKの比率として適� �率PFIT(=SUMK/SUMM)を算出する。そして、適合率P FITが判定閾値SLVLを超えるとノッキングが発� �したと判定する。

 図6(b)は、図5(b)に示す二値化強度パラメ� �タNKMAP(j,i)と、図6(a)に示すマスタパターンマ ップ上のマスタパラメータMMAP(j,i)との積(NKMAP (j,i)×MMAP(j,i))のマップを示す。この例では、� ��6(b)の積マップは、図5(b)に示すもとのマッ� �と全く同一となり、適合率PFITは、「0.863」� �なる。

 なお、後述する実際のノッキング判定処� ��においては、上記強度積算値SUMK及び基準積 算値SUMMを算出する際にエンジン運転状態に� �じた重み付けを行うようにしている。この� �み付けを行うための重み付けパラメータWMAP( j,i)が設定されている重み付けマップの一例� �図7に示す。図7に示す重み付けマップは、周 波数6kHzの成分についてクランク角12度から72� ��までのパラメータ値、周波数10kHzの成分に� �いてクランク角36度から66度までのパラメー� ��値、並びに周波数11kHz、13kHz、15kHz、及び20kH zの成分についてクランク角12度から30度まで� ��パラメータ値に対して重みが付けられるよ� ��に設定されている。

 図7に示すような重み付けマップによる重 み付けは、エンジン運転状態に依存して二値 化スペクトル時系列マップの周波数に対する 特性が変化することを補償するために行われ るものである。エンジン運転状態に応じて重 み付けを行うことにより、エンジン運転状態 の変化に拘わらず正確な判定を行うことが可 能となる。

 本実施形態におけるノックセンサ出力VS� �サンプリング周期は20マイクロ秒とし、50個� ��サンプリング値を用いて周波数成分分析を� ��う。これにより、25kHzまでの周波数成分を1k Hz間隔で算出することが可能となる。なお、� ��ランク軸が6度回転するのに要する時間TD6は 、エンジン回転数NEが1000rpmのとき1ミリ秒で� �り50個のサンプリング値が得られる。また例 えばエンジン回転数NEが2000rpmのときは時間TD6 が500マイクロ秒となり、サンプリングデータ の数は25となるので、周波数成分分析は、直� ��に検出された、クランク角12度分のサンプ� �ング値を用いてクランク角6度毎に行われる� ��

 周波数成分分析(高速フーリエ変換)は、� �示しない処理において実行され、クランク� �度6度毎に、演算の結果得られる周波数毎の� ��度値STFT(k)(k=0~(KN-1))がECU5の記憶回路(RAM)に格 納される。KNは、上死点後6度から90度までの� ��ランク角度範囲で得られる強度値の数であ� ��、本実施形態では315(図3~6に示すマップ上の データ数)である。

 図8は、上述した手法によりノッキング判 定を行う処理のフローチャートであり、この 処理はECU5のCPUでTDCパルスの発生に同期して� �行される。

 ステップS11では、図10に示す二値化デー� �マップ算出処理を実行し、上述した二値化� �ペクトル時系列マップの算出を行う。ステ� �プS12では、図13に示すノイズ除去処理を実行 し、図5(a)に例示するような二値化ノイズマ� �プを用いてノイズ成分を除去する処理を行� �。

 ステップS13では、図14に示す適合率算出� �理を実行し、ノイズ成分が除去された二値� �スペクトル時系列マップと、マスタパター� �マップとを用いて適合率PFITを算出する。

 ステップS14では、エンジン回転数NE及び� �気圧PBAに応じて例えば図9に示すように設定� ��れたSLVLマップを検索し、判定閾値SLVLを算� �する。SLVLマップに設定されている格子点以� ��のエンジン回転数NE及び吸気圧PBAについて� �、補間演算により判定閾値SLVLを算出する。

 ステップS15では、ステップS13で算出され� ��適合率PFITが判定閾値SLVLより大きいか否か� �判別し、その答が肯定(YES)であるときは、ノ ッキングが発生したと判定し、ノッキングフ ラグFKNOCKを「1」に設定する(ステップS16)。

 ステップS15でPFIT≦SLVLであるときは、ノ� �キングは発生していないと判定し、ノッキ� �グフラグFKNOCKを「0」に設定する(ステップS17 )。次いで図16に示すノイズ学習処理を実行し 、二値化ノイズマップ(図5(a)参照)の更新を行 う。

 図10は、図8のステップS11で実行される二値� ��データマップ算出処理のフローチャートで� ��る。
 ステップS21では、クランク角インデクスi,� �波数インデクスj,及びメモリアドレスインデ クスkをいずれも「0」に初期化する。ステッ� ��S22では、周波数インデクスjが周波数データ 数JN(本実施形態では21)から「1」を減算した� �より大きいか否かを判別する。最初はこの� �は否定(NO)であるので、ステップS23に進み、� ��ランク角インデクスiがクランク角データ数 IN(本実施形態では15)から「1」を減算した値� �り大きいか否かを判別する。

 最初はステップS23の答も否定(NO)であるので 、メモリアドレスインデクスkを下記式(1)に� �り算出する(ステップS24)。したがって、メモ リアドレスインデクスkは、クランク角イン� �クスi及び周波数インデクスjの増加に伴って 、0から(JN×IN-1)まで変化する。
 k=i+j×IN                (1)

 ステップS25では、強度パラメータKMAP(j,i)� ��、メモリに格納されている強度値STFT[k]に設 定し、ステップS26では、図11に示す二値化処� ��を実行する。次いでクランク角インデクスi を「1」だけインクリメントする(ステップS27) 。

 図11のステップS31では、エンジン回転数NE 及び吸気圧PBAに応じて図12に示すBLVLマップを 検索し、二値化閾値BLVLを算出する。BLVLマッ� ��は、エンジン回転数NEが増加するほど二値� �閾値BLVLが増加し、かつ吸気圧PBAが増加する� ��ど二値化閾値BLVLが増加するように設定され ている。ラインL1で示される設定値は、第1所 定吸気圧PBA1(例えば53kPa(400mmHg))から第2所定吸 気圧PBA2(例えば80kPa(600mmHg))の範囲で適用され� ��。ラインL2及びL3は、それぞれ第3所定吸気� �PBA3(例えば93kPa(700mmHg))及び第4所定吸気圧PBA4( 例えば107kPa(800mmHg))に対応する。

 ステップS32では、強度パラメータKMAP(j,i)� ��二値化閾値BLVLより大きいか否かを判別し、 その答が肯定(YES)であるときは、二値化強度� ��ラメータNKMAP(j,i)を「1」に設定する(ステッ� ��S33)。一方、ステップS32でKMAP(j,i)≦BLVLであ� �ときは、二値化強度パラメータNKMAP(j,i)を「0 」に設定する(ステップS34)。

 図10に戻り、ステップS23の答が否定(NO)で� ��る間はステップS24~S27を繰り返し実行し、ク ランク角インデクスiが(IN-1)を超えると、ス� �ップS28に進み、クランク角インデクスiを「0 」に戻すとともに、周波数インデクスjを「1� ��だけインクリメントし、ステップS22に戻る� ��

 ステップS22の答が否定(NO)である間は、ス テップS23~S28を繰り返し実行し、周波数イン� �クスjが(JN-1)を超えると、本処理を終了する� ��

 図10の処理により、メモリに格納されて� �る、周波数成分分析の結果として得られた� �度値STFT[k]が、図3(a)に示すスペクトル時系列 マップの形式に変換されるとともに、強度パ ラメータKMAP(j,i)の二値化が行われ、二値化強 度パラメータNKMAP(j,i)が算出される。すなわ� �、図3(b)に示す二値化スペクトル時系列マッ� ��が生成される。

 図13は、図8のステップS12で実行されるノイ� ��除去処理のフローチャートである。
 ステップS41では、クランク角インデクスi及 び周波数インデクスjをともに「0」に初期化� ��る。ステップS42では、周波数インデクスjが 周波数データ数JNから「1」を減算した値より 大きいか否かを判別する。最初はこの答は否 定(NO)であるので、ステップS43に進み、クラ� �ク角インデクスiがクランク角データ数INか� �「1」を減算した値より大きいか否かを判別� ��る。

 最初はステップS43の答も否定(NO)であるので 、下記式(2)により二値化強度パラメータNKMAP( j,i)を補正し、補正二値化強度パラメータJKMAP (j,i)を算出する。式(2)のNNMAP(j,i)は、学習処理 により更新される二値化ノイズマップ上の二 値化ノイズパラメータである。
 JKMAP(j,i)=NKMAP(j,i)-NNMAP(j,i)     (2)

 ステップS45では、補正二値化強度パラメ� ��タJKMAP(j,i)が負の値であるか否かを判別し、 その答が否定(NO)であるときは直ちにステッ� �S47に進む。JKMAP(j,i)<0であるときは、JKMAP(j, i)を「0」に設定し(ステップS46)、ステップS47� ��進む。

 ステップS47ではクランク角インデクスiを 「1」だけインクリメントし、ステップS43に� �る。ステップS43の答が否定(NO)である間はス� ��ップS44~S47を繰り返し実行し、クランク角イ ンデクスiが(IN-1)を超えると、ステップS48に� �み、クランク角インデクスiを「0」に戻すと ともに、周波数インデクスjを「1」だけイン� ��リメントし、ステップS42に戻る。ステップS 42の答が否定(NO)である間は、ステップS43~S48� �繰り返し実行し、周波数インデクスjが(JN-1)� ��超えると、本処理を終了する。

 図13の処理により、ノイズが除去された� �正二値化強度パラメータJKMAP(j,i)が得られる( 図5(b)参照)。

 図14は、図8のステップS13で実行される適合� ��算出処理のフローチャートである。
 ステップS51では、エンジン回転数NE及び吸� �圧PBAに応じてマスタパターンマップを選択� �、ステップS52では、エンジン回転数NE及び吸 気圧PBAに応じて重み付けマップを選択する。 重み付けマップは、エンジン運転状態に依存 して二値化スペクトル時系列マップの周波数 に対する特性が変化することを補償するため に設けられている。エンジン回転数NEまたは� ��気圧PBA(エンジン負荷)が変化すると、燃焼� �内の温度が変化し、二値化スペクトル時系� �マップが変化する。したがって、エンジン� �転数NE及び吸気圧PBAに応じてマスタパターン マップ及び重み付けマップを選択することに より、エンジン運転状態の変化に拘わらず正 確な判定を行うことが可能となる。

 本実施形態では、図15に示すようにエン� �ン回転数NE及び吸気圧PBAにより定義される9� �のエンジン運転領域に対応して、9個のマス� ��パターンマップ及び9個の重み付けマップが 予め設定されており、ステップS51では9個の� �スタパターンマップうちの1つが選択され、� ��テップS52では、9個の重み付けマップのうち の1つが選択される。図15において低回転領域 は例えばエンジン回転数NEが2000rpm以下の領域 とされ、中回転領域は2000rpmから4000rpmまでの� ��域とされ、高回転領域は4000rpmを超える領域 とされる。また低負荷領域は例えば吸気圧PBA が67kPa(500mmHg)以下の領域とされ、中負荷領域� ��67kPaから93kPa(700mmHg)までの領域とされ、高負 荷領域は93kPaを超える領域とされる。

 ステップS53では、クランク角インデクスi 及び周波数インデクスjをともに「0」に初期� ��するとともに、強度積算値SUMK及び基準積算 値SUMMを「0」に初期化する。強度積算値SUMK及 び基準積算値SUMMは、後述するステップS57で� �新され、ステップS60で適合率PFITの算出に適� ��される。

 ステップS54では、周波数インデクスjが周 波数データ数JNから「1」を減算した値より大 きいか否かを判別する。最初はこの答は否定 (NO)であるので、ステップS55に進み、クラン� �角インデクスiがクランク角データ数INから� �1」を減算した値より大きいか否かを判別す� ��。

 最初はステップS55の答も否定(NO)であるので 、ステップS56に進み、下記式(3)及び(4)により 、重み付けマスタパラメータMMW及び重み付け 積パラメータKMWを算出する。下記式のWMAP(j,i) は、重み付けマップに設定されている重み付 けパラメータである。重み付け積パラメータ KMWは、マスタパラメータMMAP(j,i)と補正二値化 強度パラメータJKMAP(j,i)との積を、重み付け� �ラメータWMAP(j,i)によって重み付けしたもの� �ある。
 MMW=MMAP(j,i)×WMAP(j,i)            (3)
 KMW=MMAP(j,i)×JKMAP(j,i)×WMAP(j,i)    (4)

 ステップS57では、下記式(5)及び(6)により、� ��み付けマスタパラメータMMW及び重み付け積� ��ラメータKMWを積算し、基準積算値SUMM及び強 度積算値SUMKを算出する。
 SUMM=SUMM+MMW                 (5)
 SUMK=SUMK+KMW                 (6)

 ステップS58では、クランク角インデクスiを 「1」だけインクリメントし、ステップS55に� �る。ステップS55の答が否定(NO)である間はス� ��ップS56~S58を繰り返し実行し、クランク角イ ンデクスiが(IN-1)を超えると、ステップS59に� �み、クランク角インデクスiを「0」に戻すと ともに、周波数インデクスjを「1」だけイン� ��リメントし、ステップS54に戻る。ステップS 54の答が否定(NO)である間は、ステップS55~S59� �繰り返し実行し、周波数インデクスjが(JN-1)� ��超えると、ステップS60に進み、下記式(7)に� ��り適合率PFITを算出する。
 PFIT=SUMK/SUMM                (7)

 図16は、図8のステップS18で実行されるノイ� ��学習処理のフローチャートである。
 ステップS71では、クランク角インデクスi、 周波数インデクスj、加算学習パラメータLK、 及び減算学習パラメータLMをいずれも「0」に 初期化する。ステップS72では、周波数インデ クスjが周波数データ数JNから「1」を減算し� �値より大きいか否かを判別する。最初はこ� �答は否定(NO)であるので、ステップS73に進み� ��クランク角インデクスiがクランク角データ 数INから「1」を減算した値より大きいか否か を判別する。

 最初はステップS73の答も否定(NO)であるの で、ステップS74に進み、二値化強度パラメー タNKMAP(j,i)が二値化ノイズパラメータNNMAP(j,i)� ��等しいか否かを判別する。この答が肯定(YES )であるときは、直ちにステップS80に進む。

 ステップS74の答が否定(NO)であって、二値 化強度パラメータNKMAP(j,i)が二値化ノイズパ� �メータNNMAP(j,i)と異なるときは、二値化強度� ��ラメータNKMAP(j,i)が二値化ノイズパラメータ NNMAP(j,i)より大きいか否かを判別する(ステッ� ��S75)。この答が肯定(YES)であるときは、加算� ��習パラメータLKを「1」に設定し、減算学習� ��ラメータLMを「0」に設定する(ステップS76)� �一方、NKMAP(j,i)<NNMAP(j,i)であるときは、加� �学習パラメータLKを「0」に設定し、減算学� �パラメータLMを「1」に設定する(ステップS77) 。

 ステップS78では、下記式(8)及び(9)により、� ��算学習パラメータLK及び減算学習パラメー� �LMを修正する。式(8)、(9)のDSNOISEは、例えば0. 1に設定されるノイズ学習係数である。
 LK=DSNOISE×LK                 (8)
 LM=DSNOISE×LM                 (9)

 ステップS79では、下記式(10)に加算学習パラ メータLK及び減算学習パラメータLMを適用し� �ノイズパラメータNMAP(j,i)を更新する。次に� �明する図17の処理で、ノイズパラメータNMAP(j ,i)を二値化することにより、二値化ノイズパ ラメータNNMAP(j,i)が算出される。
 NMAP(j,i)=NMAP(j,i)+LK-LM        (10)

 ステップS80ではクランク角インデクスiを 「1」だけインクリメントし、ステップS73に� �る。ステップS73の答が否定(NO)である間はス� ��ップS74~S80を繰り返し実行し、クランク角イ ンデクスiが(IN-1)を超えると、ステップS81に� �み、クランク角インデクスiを「0」に戻すと ともに、周波数インデクスjを「1」だけイン� ��リメントし、ステップS72に戻る。ステップS 72の答が否定(NO)である間は、ステップS73~S81� �繰り返し実行し、周波数インデクスjが(JN-1)� ��超えると、ステップS82に進み、図17に示す� �イズマップ更新処理を実行する。

 図17のステップS91では、クランク角イン� �クスi及び周波数インデクスjをともに「0」� �初期化する。ステップS92では、周波数イン� �クスjが周波数データ数JNから「1」を減算し� ��値より大きいか否かを判別する。最初はこ� ��答は否定(NO)であるので、ステップS93に進み 、クランク角インデクスiがクランク角デー� �数INから「1」を減算した値より大きいか否� �を判別する。

 最初はステップS93の答も否定(NO)であるの で、ステップS94に進み、ノイズパラメータNMA P(j,i)がノイズ二値化閾値NLVL(例えば0.8)より大 きいか否かを判別する。この答が肯定(YES)で� ��るときは、二値化ノイズパラメータNNMAP(j,i) を「1」に設定する(ステップS95)。一方、NMAP(j ,i)≦NLVLであるときは、二値化ノイズパラメ� �タNNMAP(j,i)を「0」に設定する(ステップS96)。

 ステップS97ではクランク角インデクスiを 「1」だけインクリメントし、ステップS93に� �る。ステップS93の答が否定(NO)である間はス� ��ップS94~S97を繰り返し実行し、クランク角イ ンデクスiが(IN-1)を超えると、ステップS98に� �み、クランク角インデクスiを「0」に戻すと ともに、周波数インデクスjを「1」だけイン� ��リメントし、ステップS92に戻る。ステップS 92の答が否定(NO)である間は、ステップS93~S98� �繰り返し実行し、周波数インデクスjが(JN-1)� ��超えると、本処理を終了する。

 図16のステップS75でNKMAP(j,i)>NNMAP(j,i)で� �るときは、ステップS78でLK=0.1,LM=0となり、式 (10)により、ノイズパラメータNMAP(j,i)が「0.1� �だけインクリメントされる。一方、ステッ� �S75でNKMAP(j,i)<NNMAP(j,i)であるときは、ステ� �プS78でLK=0,LM=0.1となり、式(10)により、ノイ� �パラメータNMAP(j,i)が「0.1」だけデクリメン� �される。そして、図17の処理でノイズパラメ ータNMAP(j,i)がノイズ二値化閾値NLVLより大き� �とき、二値化ノイズパラメータNNMAP(j,i)が「1 」に設定される一方、ノイズパラメータNMAP(j ,i)がノイズ二値化閾値NLVL以下であるときは� �二値化ノイズパラメータNNMAP(j,i)が「0」に設 定される。

 図16及び図17の処理により、ノッキングが 発生していないと判定されたときの二値化強 度パラメータNKMAP(j,i)に応じて二値化ノイズ� �ップが更新され、例えば吸気弁の着座ノイ� �のように定常的に発生するノイズが二値化� �イズマップに反映される。その結果、ノイ� �の影響を除いて高精度の判定を行うことが� �能となる。

 なお、ノイズ除去処理(図8,ステップS12)を 行わずに、二値化強度パラメータNKMAPをその� ��ま用いて適合率PFITの算出を行うようにして よい。ノイズ除去処理を行わない場合には、 ノイズの影響を受ける可能性が高くなるが、 図30を参照して説明したようにノイズと区別� ��てノッキングの判定を行うことができるか� ��である。

 以上詳述したように本実施形態では、ク� ��ンク角6度間隔でノックセンサ11の出力信号� ��高速フーリエ変換演算(周波数成分分析)が� �われ、その結果得られる、5kHz~25kHzの周波数� ��分の強度の時系列データであるスペクトル� ��系列マップが生成される。すなわち、スペ� ��トル時系列マップの要素が二次元配列デー� ��である強度パラメータKMAP(j,i)としてメモリ� ��格納される。そして、強度パラメータKMAP(j, i)を二値化することにより、二値化強度パラ� ��ータNKMAP(j,i)が算出され、該二値化強度パラ メータNKMAP(j,i)に基づいてノッキングが発生� �たか否かが判定される。二値化強度パラメ� �タNKMAP(j,i)には、エンジンの回転に伴う周波� ��成分分布の変化が反映されるので、二値化� ��度パラメータNKMAP(j,i)と、ノッキング発生時 に特有の変化パターンに対応するマスタパタ ーンマップ上のマスタパラメータMMAP(j,i)とを 比較することにより、ノッキングの発生を正 確に判定することができる。また周波数成分 分析の結果得られる強度パラメータKMAP(j,i)を 二値化することにより、データ量が減少する とともに時系列データの変化パターンが単純 化されるので、メモリ容量を低減するととも に演算速度を高めることができる。

 また二値化強度パラメータNKMAP(j,i)(JKMAP(j, i))が、マスタパラメータMMAP(j,i)に近い変化パ ターンを示すときに、ノッキングが発生して いる可能性が高いので、二値化強度パラメー タNKMAP(j,i)(JKMAP(j,i))とマスタパラメータMMAP(j,i )との類似性(相関性)を示すパラメータを算出 することにより、正確な判定を行うことがで きる。

 本実施形態ではこの類似性(相関性)を示� �パラメータとして適合率PFITを用い、適合率P FITが判定閾値SLVLを超えたときにノッキング� �発生したと判定される。適合率PFITを用いる� ��とにより、二値化強度パラメータNKMAP(j,i)(JK MAP(j,i))とマスタパラメータMMAP(j,i)との類似性 (相関性)を比較的簡単な演算で的確に評価し� ��正確な判定を行うことができる。

 また二値化強度パラメータNKMAP(j,i)に基づ いて、ノイズ成分の時系列データNNMAP(j,i)が� �出され、二値化強度パラメータNKMAP(j,i)がノ� ��ズ成分の時系列データであるノイズパラメ� ��タNNMAP(j,i)により補正され、補正二値化強度 パラメータJKMAP(j,i)に基づいてノッキング判� �が行われる。したがって、上述した着座ノ� �ズのように定常的に表れるノイズ成分を除� �してより正確な判定を行うことが可能とな� �。

 また二値化強度パラメータNKMAP(j,i)(JNKMAP(j ,i))及びマスタパラメータMMAP(j,i)に対して、� �波数に応じて設定された重み付けパラメー� �WMAP(j,i)を乗算して、適合率PFITが算出される� ��ノッキングが発生したときに大きくなる周� ��数成分は予め判明しているので、その周波� ��の近傍の周波数に対応するパラメータに大� ��な重みを付けることにより、判定精度を高� ��ることができる。

 本実施形態では、ECU5がノッキング判定手 段、周波数成分分析手段、データ格納手段、 二値化手段、ノイズ成分算出手段、ノイズ補 正手段、及び重み付け手段を構成する。具体 的には、図8のステップS12~S18がノッキング判� ��手段に相当し、ステップS11(図10の処理)がデ ータ格納手段及び二値化手段に相当する。ま たステップS18がノイズ成分算出手段に相当し 、ステップS12がノイズ補正手段に相当し、図 14のステップS52及びS56が重み付け手段に相当� ��る。

 [第2の実施形態]
 本実施形態は、第1の実施形態におけるノッ キング判定処理(図8)を、図18に示す処理に変� ��したものである。以下に説明する点以外は� ��第1の実施形態と同一である。

 図18に示す処理では、先ず二値化する前� �強度パラメータKMAPについて、ノイズ除去処� ��を実行して、補正強度パラメータJKMAPIを算� ��し(ステップS12a)、補正強度パラメータJKMAPI� ��二値化して、補正二値化強度パラメータJKMA Pを算出する(ステップS11b)。

 図19(a)は、本実施形態におけるスペクト� �時系列マップの一例を示し、図19(b)は、本実 施形態におけるノイズマップの一例を示す図 であり、二値化されていないノイズ学習値が 設定されている。ノイズ除去処理は、図19(a)� ��スペクトル時系列マップの各マップ値から� ��19(b)のノイズマップ値を減算することによ� �行われ、図20(a)に示す補正スペクトル時系列 マップが得られる。そして、図20(a)の補正ス� ��クトル時系列マップを二値化することによ� ��、図20(b)に示す補正二値化強度パラメータ� �ップが得られる。

 図21は、図18のステップS11aで実行される� �ータマップ算出処理のフローチャートであ� �。この処理は、図10に示す二値化データマッ プ算出処理のステップS26を削除したものであ る。すなわち、強度値STFT(k)をスペクトル時� �列マップの形式に変換する処理のみ行われ� �。

 図22は、図18のステップS12aで実行される� �イズ除去処理のフローチャートである。こ� �処理は、図13に示す処理のステップS44~S46を� �れぞれステップS44a~S46aに変更したものであ� �。

 ステップS44aでは、下記式(2a)により、補正� �度パラメータJKMAPIを算出する。式(2a)のNLMAP(j ,i)は、学習処理により更新されるノイズマッ プ(図19(b))上のノイズパラメータである。
 JKMAPI(j,i)=KMAP(j,i)-NLMAP(j,i)     (2a)

 ステップS45aでは、補正強度パラメータJKM API(j,i)が負の値であるか否かを判別し、その� ��が否定(NO)であるときは直ちにステップS47に 進む。JKMAPI(j,i)<0であるときは、JKMAPI(j,i)を 「0」に設定し(ステップS46a)、ステップS47に� �む。

 図22の処理により図20(a)に示す補正スペク トル時系列マップが得られる。

 図23は、図18のステップS11bで実行される� �値化処理のフローチャートである。この処� �は図11に示す二値化処理のステップS32~S34を� �れぞれステップS32a~34aに変更し、さらにすべ てのマップ値について二値化を行うためのス テップS41~S43,S47及びS48(図22に示すステップと� ��じ処理を行うステップ)を追加したものであ る。

 ステップS32aでは、補正強度パラメータJKM API(j,i)が二値化閾値BLVLより大きいか否かを判 別し、その答が肯定(YES)であるときは、補正� ��値化強度パラメータJKMAP(j,i)を「1」に設定� �る(ステップS33a)。一方、ステップS32aでJKMAPI( j,i)≦BLVLであるときは、補正二値化強度パラ� ��ータJKMAP(j,i)を「0」に設定する(ステップS34a )。

 図23の処理により、図20(a)に示すマップか ら図20(b)に示す補正二値化強度パラメータマ� ��プが得られる。

 図24は、図18のステップS18aで実行される� �イズ学習処理のフローチャートである。こ� �処理は、図16に示す処理のステップS74~S78,及� ��S82を削除し、ステップS79をステップS79aに変 更したものである。

 ステップS79aでは、ノイズパラメータNLMAP(j,i )の更新を行う。具体的には、ノッキングが� �出されないときに得られた強度パラメータKM AP(j,i)を下記式(31)に適用し、ノイズパラメー� ��NLMAP(j,i)を更新する。右辺のノイズパラメー タNLMAP(j,i)は、更新前のマップ値であり、DSNOI SEaは例えば「0.5」に設定されるなまし係数で ある。
 NLMAP(j,i)=NLMAP(j,i)×DSNOSEa
        +KMAP(j,i)×(1-DSNOSEa)   (31)

 なお、式(31)のなまし係数DSNOISEaは、以下に� ��明するようにエンジン運転状態に応じて0か ら1の範囲で変更するようにしてもよい。
 1)エンジン回転数NEが急激に上昇したとき(� �えばエンジン回転数の今回値NE(m)と前回値NE( m-1)との差(NE(m)-NE(m-1))が200rpm以上であるとき)� ��、なまし係数DSNOISEaをより小さい値に変更� �る。これにより、強度パラメータの今回値KM AP(j,i)の重みを増加させ、ノイズパラメータNL MAP(j,i)の学習速度が高められる。上記「m」は 、図18の処理の実行周期で離散化した離散化� ��刻である。

 2)吸気弁の開弁時期CAIが急激に変化した� �き(例えば開弁時期の今回値CAI(m)と前回値CAI( m-1)との差の絶対値が5度以上であるとき)は、 なまし係数DSNOISEaをより小さい値に変更する� ��

 3)吸気弁のリフト量LFTが急激に変化した� �き(例えばリフト量の今回値LFT(m)と前回値LFT( m-1)との差の絶対値が0.5mm以上であるとき)は� �なまし係数DSNOISEaをより小さい値に変更する 。

 本実施形態によれば、ノイズマップのマ� ��プ値(NLMAP)が、二値化前の強度パラメータKMA Pを用いて更新されるので、第1の実施形態と� ��べて、より精度の高いノイズパラメータNLMA P(j,i)が得られる。

 本実施形態では、図18のステップS11aがデ� ��タ格納手段に相当し、ステップS12aがノイズ 補正手段に相当し、ステップS18aがノイズ成� �算出手段に相当し、ステップS11bが二値化手� ��に相当し、ステップS13~S17がノッキング判定 手段に相当する。

 [第3の実施形態]
 本実施形態は、5kHzから25kHzまでの21個の周� �数を、3個の周波数からなる7つの周波数グル ープG1~G7に分割(グループ化)し、1つサンプル� ��イミングに対応する強度パラメータKMAPの代 表値であるグループ強度パラメータGKMAPを、� ��波数グループ内の最大値に設定するように� ��たものである。このようにグループ化を最� ��に行うことにより、ノイズ除去処理、二値� ��処理、及び適合率算出処理の対象となるデ� ��タ数が1/3に減少し、演算負荷を大幅に低減� ��ることができる。なお、以下に説明する点� ��外は、第2の実施形態と同一である。

 図25は、グループ化を説明するための図� �ある。図25(a)に示すように、周波数5,6,7kHzが� ��波数グループG1に対応し、周波数8,9,10kHzが� �波数グループG2に対応し、周波数11,12,13kHzが� ��波数グループG3に対応し、周波数14,15,16kHzが 周波数グループG4に対応し、周波数17,18,19kHz� �周波数グループG5に対応し、周波数20,21,22kHz� ��周波数グループG6に対応し、周波数23,24,25kHz が周波数グループG7に対応する。各サンプル� ��イミング(クランク角度)におけるグループ� �の最大値がハッチングを付して示されてい� �。この最大値を抽出することにより、図25(b) に示すグループ化したスペクトル時系列マッ プが得られる。

 本実施形態では、図25(b)に示すスペクト� �時系列マップを用いて、ノッキング判定が� �われる。図25(b)に示すスペクトル時系列マッ プの設定値を、「グループ強度パラメータGKM AP(n,i)」という。「n」は、周波数グループを� ��すグループインデクスであり、本実施形態� ��は周波数グループG1からG7に対応して「0」� �ら「6」の値をとる。

 図26(a)は、図25(b)と同一のマップを示し、 図26(b)は、第2の実施形態で示した学習処理に より得られるノイズマップである。図26(b)に� ��すノイズマップの設定値を、「グループノ� ��ズパラメータGNLMAP(n,i)」という。

 グループ強度パラメータGKMAP(n,i)からグル ープノイズパラメータGNLMAP(n,i)を減算するこ� ��により、グループ補正強度パラメータGJKMAPI (n,i)が算出される。図26(c)は、グループ補正� �度パラメータGJKMAPI(n,i)が設定されたスペク� �ル時系列マップを示す。

 さらにグループ補正強度パラメータGJKMAPI (n,j)を二値化することにより、二値化グルー� ��補正強度パラメータGJKMAP(n,j)が算出される� �図26(d)は、二値化グループ補正強度パラメー タGJKMAP(n,i)が設定されたスペクトル時系列マ� ��プを示す。またマスタパターンマップは例� ��ば図26(e)で与えられる。

 二値化グループ補正強度パラメータGJKMAP� ��び図26(e)に示すマスタパターンマップに基� �いて算出される適合率PFITに応じたノッキン� ��判定は、上述した実施形態と同様にして行� ��ことができる。

 図27は、本実施形態におけるノッキング� �定処理のフローチャートである。図18に示す 処理にグループ化処理を行うステップS100を� �加し、ステップS12a,S11b,S13~S17,S18aをそれぞれ� ��テップS101~S108に変更したものである。

 ステップS101~S103及びS108の処理が対応する 図18のステップの処理と異なる点は、処理対� ��となるマップの行数が21行ではなく7行であ� ��点のみであるので、詳細な説明を省略する� ��

 図28は、図27のステップS100で実行されるグ� �ープ化処理のフローチャートである。
 ステップS111では、クランク角インデクスi,� ��波数インデクスj,及びグループインデクスn� ��すべて「0」に初期化する。ステップS112で� �、周波数インデクスjが周波数データ数JNか� �「1」を減算した値より大きいか否かを判別� ��る。最初はこの答は否定(NO)であるので、ス テップS113に進み、クランク角インデクスiが� ��ランク角データ数INから「1」を減算した値� ��り大きいか否かを判別する。

 最初はステップS113の答も否定(NO)である� �で、ステップS114に進み、パラメータa,b,cを� �れぞれ周波数方向に連続する強度パラメー� �KMAP(j,i),KMAP(j+1,i),KMAP(j+2,i)に設定する。ステ� �プS115では、グループ強度パラメータGKMAP(n,i) をパラメータa,b,cの最大値に設定する。

 ステップS116ではクランク角インデクスi� �「1」だけインクリメントし、ステップS113に 戻る。ステップS113の答が否定(NO)である間は� ��テップS114~S116を繰り返し実行し、クランク� ��インデクスiが(IN-1)を超えると、ステップS11 7に進み、クランク角インデクスiを「0」に戻 すとともに、周波数インデクスjを「3」だけ� ��ンクリメントし、ステップS112に戻る。ステ ップS112の答が否定(NO)である間は、ステップS 113~S117を繰り返し実行し、周波数インデクスj が(JN-1)を超えると、本処理を終了する。

 図28の処理により、グループ強度パラメ� �タGKMAP(n,i)が設定されたスペクトル時系列マ� ��プ(図25(b))が得られる。

 本実施形態によれば、グループ強度パラ� ��ータGKMAPを対象として、ノイズ除去処理、� �値化処理、及び適合率算出処理が行われる� �で、処理対象となるデータ数が減少し、演� �負荷を低減することができる。

 本実施形態では、図28の処理が周波数グ� �ープ時系列データ生成手段に相当し、図27の ステップS108がグループノイズ成分算出手段� �相当し、ステップS101がグループノイズ補正� ��段に相当する。

 なお本発明は上述した実施形態に限るも� ��ではなく、種々の変形が可能である。例え� ��、ノックセンサ出力のサンプリング周期や� ��波数成分分析を行うクランク角度間隔は上� ��したもの(20マイクロ秒、6度)に限るもので� �なく、本発明の目的が達成される範囲内に� �いて変更可能である。また、二値化スペク� �ル時系列マップ(上述した実施形態では21行× 15列のマトリクスで構成)も同様に変更可能で ある。

 また上述した第1の実施形態では、補正二 値化強度パラメータJKMAP(j,i)及びマスタパラ� �ータMMAP(j,i)に重み付けパラメータWMAP(j,i)を� �算して、適合率PFITを算出したが、重み付け� ��ラメータWMAP(j,i)を乗算せずに、すなわち重� ��付けを行わずに算出するようにしてもよい� ��

 また上述した実施形態では、判定閾値SLVL 、二値化閾値BLVL、マスタパターンマップ、� �び重み付けマップをエンジン回転数NE及び吸 気圧PBAに応じて、算出または選択するように したが、予め設定された値または1つのマッ� �に固定しておいてもよい。

 また適合率PFITは、図29に示す処理により� ��出するようにしてもよい。図29に示す処理� �、図14のステップS53,S56,S57,S60をステップS53a,S 56a,S57a,S57b,及び60aに変更したものである。

 ステップS53aでは、クランク角インデクス i及び周波数インデクスjを初期化するととも� ��、強度積算値SUMKa及び基準積算値SUMMaを「0� �に初期化する。

 ステップS56aでは、補正二値化強度パラメー タJKMAP(j,i)が対応するマスタパラメータMMAP(j,i )と等しいか否かを判別し、その答が肯定(YES) であるときは、下記式(11)により、強度積算� �SUMKaを更新する(ステップS57a)。
 SUMKa=SUMKa+WMAP(j,i)           (11)

 ステップS56aでJKMAP(j,i)が対応するマスタパ� �メータMMAP(j,i)と等しくないときは、直ちに� �テップS57bに進む。
 ステップS57bでは、下記式(12)により、基準� �算値SUMMaを更新する。
 SUMMa=SUMMa+WMAP(j,i)           (12)

 ステップS60aでは、下記式(13)に強度積算値SU MKa及び基準積算値SUMMaを適用し、適合率PFITを 算出する。
 PFIT=SUMKa/SUMMa              (13)

 図29の処理により算出される適合率PFITも� ��補正二値化強度パラメータJKMAP(j,i)とマスタ パラメータMMAP(j,i)の類似性(相関性)が高いほ� ��、大きな値となり、類似性(相関性)を示す� �ラメータとして使用することができる。

 また上述した実施形態では、補正二値化強� ��パラメータJKMAP(j,i)とマスタパラメータMMAP(j ,i)との類似性(相関性)を示すパラメータとし� ��、適合度PFITを用いたが、例えば下記式(21)� �算出される補正二値化強度パラメータJKMAP(j, i)とマスタパラメータMMAP(j,i)の差の絶対値の� ��計である差分積算値SUMDIFを用いることもで� ��る。

 差分積算値SUMDIFは、その値が小さいほど� ��似性(相関性)が高いことを示すので、差分� �算値SUMDIFが差分判定閾値SDLVLより小さいとき 、ノッキングが発生したと判定する。

 また上述した第3の実施形態では、グループ 強度パラメータGKMAPを
グループノイズパラメータGNLMAPにより補正す ることによりグループ補正強度パラメータGJK MAPIを算出し、グループ補正強度パラメータGJ KMAPIを二値化することにより、二値化グルー� ��補正強度パラメータGJKMAPを算出するように� ��たが、第1の実施形態と同様に、グループ強 度パラメータGKMAPを二値化して二値化グルー� ��強度パラメータGNKMAPを算出し、二値化グル� ��プ強度パラメータGNKMAPを二値化グループノ� ��ズパラメータNNMAPで補正することにより二� �化グループ補正強度パラメータGJKMAPを算出� �るようにしてもよい。

 また上述した実施形態では、ノッキング� ��定を常時行うようにしているが、ノッキン� ��の発生し難い気筒の有無を判定し、そのよ� ��な気筒(以下「気筒CYLXK」という)が有る場合 には、気筒CYLXKの爆発行程におけるノッキン� ��判定処理を行わないようにしてもよい。こ� ��によって、ノッキング判定処理の演算負荷� ��低減することができる。

 気筒CYLXKの選択は例えば以下のようにし� �行う。エンジン回転数NEが比較的低い運転状 態で全気筒についてノッキング判定を実行し 、そのときに点火時期の遅角補正値IGKNOCKを� �筒毎に算出する。そして遅角補正値IGKNOCKの� ��小値が所定閾値より小さい気筒を、気筒CYLX Kとして選択する。そして、エンジン回転数NE が所定回転数より高い高回転領域において、 気筒CYLXKに対応するノッキング判定処理を停� ��する。点火時期の遅角補正値IGKNOCKは、ノッ キングが検出されたとき、第1所定量だけ増� �させ、ノッキングが検出されないときは点� �時期毎に漸減することにより算出される(こ� ��算出手法は、例えば特開2004-353473号公報に� �されている)。

 また本発明は、クランク軸を鉛直方向と� ��た船外機などのような船舶推進機用エンジ� ��などのノッキング検出にも適用が可能であ� ��。