特許 7589950 非線形モデルによって支配される少なくとも１つの移動対象の動きを表す状態を包囲する制約付きゾノトープを推定する方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-11-18

(45)【発行日】2024-11-26

(54)【発明の名称】非線形モデルによって支配される少なくとも１つの移動対象の動きを表す状態を包囲する制約付きゾノトープを推定する方法

(51)【国際特許分類】

   G01C 21/28 20060101AFI20241119BHJP

   G01S 19/51 20100101ALI20241119BHJP

【ＦＩ】

G01C21/28

G01S19/51

【請求項の数】 12

(21)【出願番号】P 2021566225

(86)(22)【出願日】2020-05-07

(65)【公表番号】

(43)【公表日】2022-07-13

(86)【国際出願番号】 EP2020062748

(87)【国際公開番号】W WO2020225379

(87)【国際公開日】2020-11-12

【審査請求日】2023-04-20

(31)【優先権主張番号】19305582.9

(32)【優先日】2019-05-07

(33)【優先権主張国・地域又は機関】EP

(73)【特許権者】

【識別番号】521485885

【氏名又は名称】アンスティテュ、シュペリエル、ドゥ、ラエロノティック、エ、ドゥ、レスパス

【氏名又は名称原語表記】ＩＮＳＴＩＴＵＴ ＳＵＰＥＲＩＥＵＲ ＤＥ Ｌ’ＡＥＲＯＮＡＵＴＩＱＵＥ ＥＴ ＤＥ Ｌ’ＥＳＰＡＣＥ

(73)【特許権者】

【識別番号】594016872

【氏名又は名称】サントル、ナショナール、ド、ラ、ルシェルシュ、シアンティフィク、（セーエヌエルエス）

(74)【代理人】

【識別番号】100107342

【弁理士】

【氏名又は名称】横田 修孝

(74)【代理人】

【識別番号】100155631

【弁理士】

【氏名又は名称】榎 保孝

(74)【代理人】

【識別番号】100137497

【弁理士】

【氏名又は名称】大森 未知子

(74)【代理人】

【識別番号】100207907

【弁理士】

【氏名又は名称】赤羽 桃子

(74)【代理人】

【識別番号】100217294

【弁理士】

【氏名又は名称】内山 尚和

(72)【発明者】

【氏名】ジャン、ボルティング

(72)【発明者】

【氏名】ソエイブ、フェルガニ

【審査官】佐々木 佳祐

(56)【参考文献】

【文献】特開２０１０－０５８７６３（ＪＰ，Ａ）

【文献】特開２００７－０９１１８１（ＪＰ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｃ ２１／００－２５／００

Ｇ０１Ｓ １９／５１

(57)【特許請求の範囲】

【請求項1】

少なくとも１つの移動対象の運動状態を包囲する制約付きゾノトープを推定する方法であって、
ａ）第１の瞬間（ｔ_ｋ－１）において前記運動状態を包囲する第１の制約付きゾノトープ（Ｚ_ｘ（ｔ_ｋ－１））の要素における非線形遷移モデルを線形化することであって、それにより、線形遷移モデルを生成する、線形化することと、
ｂ）前記第１の制約付きゾノトープ（Ｚ_ｘ（ｔ_ｋ－１））の異なる要素における遷移線形化誤差を計算することであって、前記第１の制約付きゾノトープの所与の要素における遷移線形化誤差は、前記非線形遷移モデルを前記所与の要素に適用した結果生成される出力と、前記線形遷移モデルを前記所与の要素に適用した結果生成される出力との間の差である、計算することと、
ｃ）前記遷移線形化誤差を全て含む第２の制約付きゾノトープを計算することと、
ｄ）前記第１の制約付きゾノトープを伝播させること（１０２）であって、それにより、前記第１の瞬間（ｔ_ｋ－１）後の第２の瞬間（ｔ_ｋ）において前記運動状態を包囲する事前制約付きゾノトープ

【数1】

を取得し、前記第１の制約付きゾノトープの前記伝播は、
－前記線形遷移モデルを前記第１の制約付きゾノトープに適用して、第３の制約付きゾノトープを生成し、
－前記第２の制約付きゾノトープと前記第３の制約付きゾノトープを合算することを含む、伝播させること（１０２）と、
ｅ）前記事前制約付きゾノトープ

【数2】

を更新すること（１０４）であって、それにより、前記第２の瞬間（ｔ_ｋ）において前記運動状態を包囲する事後制約付きゾノトープ（Ｚ’_ｘ（ｔ_ｋ））を取得し、前記事前制約付きゾノトープを更新することは、観測モデルを前記事前制約付きゾノトープ及び少なくとも１つのセンサによって取得された前記移動対象の観測データに適用することを含む、当該更新すること（１０４）と、
を含む方法。

【請求項2】

前記非線形遷移モデルが線形化される前記要素は、前記第１の制約付きゾノトープを含む最小区間包囲の中心である、請求項１に記載の方法。

【請求項3】

前記第１の制約付きゾノトープ又は前記事前制約付きゾノトープの要素における非線形観測モデルを線形化することであって、それにより、線形観測モデルを生成し、ステップｅ）において使用される前記観測モデルは、前記生成された線形観測モデルである、線形化することを更に含む、請求項１又は２に記載の方法。

【請求項4】

前記非線形観測モデルが線形化される前記要素は、前記事前制約付きゾノトープを含む最小区間包囲の中心である、請求項３に記載の方法。

【請求項5】

－前記第１の制約付きゾノトープ又は前記事前制約付きゾノトープの異なる要素における観測線形化誤差を計算することであって、前記第１の制約付きゾノトープ又は前記事前制約付きゾノトープの所与の要素における観測線形化誤差は、前記非線形観測モデルを前記所与の要素に適用した結果生成される出力と前記線形観測モデルを前記所与の要素に適用した結果生成される出力との間の差である、計算することと、
－計算された前記観測線形化誤差を全て含む第４の制約付きゾノトープ（Ｚ（［ε_ｈ（ｔ_ｋ）］））を計算することと、
－前記第４の制約付きゾノトープ（Ｚ（［ε _ｈ （ｔ _ｋ ）］））を膨張させることであって、それにより、前記第２の瞬間と前記第３の瞬間との間の前記運動状態の可能な全ての変動を考慮に入れるとともに、可能な場合には前記センサによって誘導される測定ノイズも考慮に入れて、状態空間において帯状領域を形成する膨張制約付きゾノトープ（Ｚ（［ε _’ ］））を取得し、前記移動対象の前記観測データは、前記第２の瞬間と異なる第３の瞬間において取得されている、膨張させることと、
－前記事前制約付きゾノトープ

【数3】

と前記膨張制約付きゾノトープ（Ｚ（［ε _’ ］））との交点として前記事後制約付きゾノトープ（Ｚ’ _ｘ （ｔ _ｋ ））を計算することと、
を更に含む、請求項３又は４に記載の方法。

【請求項6】

－前記事後制約付きゾノトープ（Ｚ’_ｘ（ｔ_ｋ））の複雑度を特定することと、
－前記複雑度が閾値よりも大きい場合、簡約化プロセスを前記事後制約付きゾノトープ（Ｚ’_ｘ（ｔ_ｋ））に適用することであって、それにより、前記簡約化前の前記事後ゾノトープ（Ｚ’_ｘ（ｔ_ｋ））よりも低い複雑性を有し、前記簡約化前の前記事後ゾノトープ（Ｚ’_ｘ（ｔ_ｋ））を含む事後ゾノトープ（Ｚ_ｘ（ｔ_ｋ））を取得する、適用することと、
を更に含む、請求項１～５の何れか１項に記載の方法。

【請求項7】

第１の制約付きゾノトープとして前記事後制約付きゾノトープ（Ｚ’_ｘ（ｔ_ｋ））を使用し、ステップａ）～ｅ）の続く実行において第１の瞬間として前記第２の瞬間を使用することを更に含む、請求項１～６の何れか１項に記載の方法。

【請求項8】

－前記第１の瞬間と前記観測データが前記センサによって取得された瞬間との間の最大時間区間を計算することと、
－前記最大時間区間が所定の持続時間よりも大きい場合、前記線形遷移モデルを予め決定された制約付きゾノトープに適用することであって、それにより、前記事前制約付きゾノトープを生成し、その他の場合、前記予め決定された制約付きゾノトープを前記事前制約付きゾノトープとして使用する、適用することと、
を更に含む、請求項１～７の何れか１項に記載の方法。

【請求項9】

前記観測データは、ＧＮＳＳ受信機によって取得される以下のデータ：疑似レンジ、取得された搬送波位相、前記移動対象と別の移動対象と間のレンジの少なくとも１つを含む、請求項１～８の何れか１項に記載の方法。

【請求項10】

前記運動データは、別の移動対象に対する前記移動対象の運動データを含む、請求項１～９の何れか１項に記載の方法。

【請求項11】

前記第１の制約付きゾノトープは、区間分析を介したベクトル集合インバータアルゴリズムを使用して計算される、請求項１～１０の何れか１項に記載の方法。

【請求項12】

コンピュータによって実行されるコンピュータプログラムであって、前記コンピュータに請求項１～１１の何れか１項に記載の方法を実行させる命令を含むコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、集合メンバシップ状態推定技法の技術分野に関する。

【0002】

本発明は、非線形モデルによって支配される少なくとも１つの移動対象の動きを表す状態を含む制約付きゾノトープを推定する方法を扱う。

【背景技術】

【0003】

無人車両は、ここ数年にわたり自律レベルの大幅な上昇を見せ、無人航空機システム（ＵＡＳ）を領空域に統合し、自律車を一般交通に統合する現在進行中の取り組みに繋がっている。もう大分以前から、協働運動戦略が、かなりの経済性を有する幾つかの利点を解き放つことが分かっていた。地上車両の自動隊列走行は、幹線道路の交通容量を上げ、交通渋滞を回避し、個々の+車両の消費を低減することができる。緊密フォーメーション飛行は、有意な航空機レンジ強化を可能にすることが実証されている。最も顕著な例には、緊密航空機フォーメーション及び自動空中給油での伴流渦効果の利用によるレンジ改善がある。これらの協働運動戦略を安全且つ効率的に実行するためには、車両の相対位置をセンチメートルレベル精度とは言わないまでもデシメートルレベルで推定する必要がある。

【0004】

緊密フォーメーション飛行では、全地球的航法衛星システム（ＧＮＳＳ）に基づくスタンドアロン位置決めは必要精度を提供しないため、様々な差動位置決め技法がこれまで提案されてきた。ここ１０年間の文献で、慣性測定ユニット、位相情報を利用した若しくは利用しない差動ＧＮＳＳ、マシンビジョン、又は両方の組合せを利用した多くのセンサ融合手法を見つけることができる。最近、安価なデシメートル精度超広帯域（ＵＷＢ）測距センサが、空中リアルタイムキネマティック（ＲＴＫ）解決策のロバスト性の改善を示すことに成功している。

【0005】

これらの既存の手法の間で共通するのは、アンセンテッドカルマンフィルタ（ＵＫＦ）又は拡張カルマンフィルタ（ＥＫＦ）等の確率的フィルタリング技法への依拠である。周知のように、多くの実用的な非線形推定問題は、ＵＫＦ及びＥＫＦの両方の根本的仮定、特に、プロセスモデル及び観測モデルの不確実性の分布が既知であるという主要仮定に違反する。現実では、通常、不確実性の両ソースの確率論的性質の近似しか利用できない。この固有で周知の問題にも拘わらず、ＥＫＦ及びＵＫＦは両方とも、実際に上手く機能し、正確な平均状態推定を提供することが実証されている。

【0006】

しかしながら、近くでの複数のエージェント（ａｇｅｎｔ）の安全な協働運動に関して言えば、フィルタ状態不確実性の一貫した計算は、正確な点推定よりも更に一層重要である。関わる近似（ＥＫＦによって実行される線形化ステップ及びＵＫＦのアンセンテッド変換）に起因して、フィルタ一貫性の保証は、ＥＫＦ又はＵＫＦを利用する場合、非線形推定問題についての未解決問題である。しかしながら、混合航空交通及び共通道路交通等の異種環境への自律車両の更なる安全な統合の鍵となる要素である車両状態の保証付き包囲を計算することができるアルゴリズムを考える。

【0007】

集合メンバシップ状態推定技法は、これに関して確率的枠組みへの興味深い代替を提供する。集合メンバシップ推定器は、観測及びプロセスモデルの不確実性が既知のコンパクトな集合、すなわち区間に含まれると仮定して、集合値プロセスモデルに依拠する。これらの集合内の誤差分布についての仮定は必要ない。厳密な集合メンバシップ推定問題は普通、非線形システムでは計算的に扱いにくい。研究は、関わる集合演算の実施可能である
が緊密な外近似を見つけることに向けられている。ここ１０年にわたり、楕円体又は区間等の集合形状を利用した多くのアルゴリズムが提案されてきた。全ての集合メンバシップ推定器が直面する共通の問題は、保守性－複雑な集合をより単純で計算的に扱いやすいもので表すことによって導入される過剰近似の程度－と計算労力との間のトレードオフを見つけることからなる。

【0008】

集合メンバシップ手法のメリットにも拘わらず、エージェント相対位置特定への適用は殆ど注目されてこなかった。区間分析アルゴリズムを介した集合反転（ＳＩＶＩＡ：Ｓｅｔ Ｉｎｖｅｒｓｉｏｎ Ｖｉａ Ｉｎｔｅｒｖａｌ Ａｎａｌｙｓｉｓ ａｌｇｏｒｉｔｈｍ）に基づく地上車両の協働集合メンバシップ位置特定アルゴリズムが提案されている。ＳＩＶＩＡ及び関連する非線形集合反転アルゴリズムは、状態推定アルゴリズムに適用される場合、緊密な包囲を提供し、近似の程度は、利用可能な計算リソースによってのみ制限される。しかしながら、これらのアルゴリズムの複雑性はしかしながら、通常、状態数に伴ってどんどん増大し、２Ｄ地上車両位置特定用途を超えての有用性は限られている。

【0009】

楕円体は、一定の複雑性により良好な計算性能を提供し、シミュレーションではエージェント相対位置特定への適用に成功してきている。近年、ゾノトープ（ｚｏｎｏｔｏｐｅ）が、保証付き状態推定に関してますます注目されてきている。ゾノトープは、楕円体よりも大きな柔軟性を提供し、複雑性は設計パラメータである。ゾノトープに固有の一つの問題は対称性であり、これは、実際の状態集合が強い非対称性を示す場合は常に、追加の保守性を導入する。第２の問題は、２つのゾノトープの交叉が常にゾノトープになるわけではなく、したがって、楕円体の場合のように、潜在的に凸最適化問題を解くことを含む中間集合の外近似を各フィルタ再帰で計算しなければならないことである。

【発明の概要】

【発明が解決しようとする課題】

【0010】

本発明は、これらの欠点の両方を解決又は少なくとも軽減することを目的とする。

【課題を解決するための手段】

【0011】

本発明の第１の態様は請求項１に記載の方法である。

【0012】

請求項１に記載の方法で使用される事前（ア・プリオリ）ゾノトープは制約付きゾノトープである。このクラスのゾノトープでは、非対称集合形状をよりよく近似することができる。追加の利点として、このクラスのゾノトープは、交叉下で閉であり、且つミンコフスキー和演算下で閉である。

【0013】

本発明の第１の態様による方法は、以下の任意選択的な特徴を含み得る。

【0014】

非線形遷移モデルが線形化される要素は、第１の制約付きゾノトープを含む最小区間包囲の中心であり得る。

【0015】

本方法は、第１の制約付きゾノトープ又は事前制約付きゾノトープの要素における非線形観測モデルを線形化することであって、それにより、線形観測モデルを生成し、ステップｅ）において使用される観測モデルは、生成された線形観測モデルである、線形化することを更に含み得る。

【0016】

非線形観測モデルが線形化される要素は、事前ゾノトープを含む最小区間包囲の中心であり得る。

【0017】

本方法は、
・第１の制約付きゾノトープ又は事前ゾノトープの異なる要素における観測線形化誤差を計算することであって、第１の制約付きゾノトープ又は事前ゾノトープの所与の要素における観測線形化誤差は、非線形観測モデルを所与の要素に適用した結果生成される出力と線形観測モデルを所与の要素に適用した結果生成される出力との間の差である、計算すること、
・計算された観測線形化誤差を全て含む第４の制約付きゾノトープを計算することであって、事後（ア・ポステリオリ）制約付きゾノトープは第４の制約付きゾノトープに依存する、計算すること
を更に含み得る。

【0018】

移動対象の観測データは、第２の瞬間と異なる第３の瞬間において取得され得、本方法は、
・第４の制約付きゾノトープを膨張させることであって、それにより、第２の瞬間と第３の瞬間との間の運動状態の可能な全ての変動を考慮に入れるとともに、可能な場合にはセンサによって誘導される測定ノイズも考慮に入れて、状態空間において帯状領域を形成する膨張制約付きゾノトープを取得する、膨張させることと、
・事前制約付きゾノトープと膨張制約付きゾノトープとの交点として事後制約付きゾノトープを計算することと
を更に含み得る。

【0019】

本方法は、
・事後制約付きゾノトープの複雑度を特定することと、
・複雑度が閾値よりも大きい場合、簡約化プロセスを事後制約付きゾノトープに適用することであって、それにより、簡約化前の事後ゾノトープよりも低い複雑性を有し、簡約化前の事後ゾノトープを含む事後ゾノトープを取得する、適用することと、
を更に含み得る。

【0020】

本方法は、第１の制約付きゾノトープとして事後制約付きゾノトープを使用し、請求項１において定義されるステップａ）～ｅ）の続く実行において第１の瞬間として第２の瞬間を使用することを更に含み得る。

【0021】

本方法は、
・第１の瞬間と観測データがセンサによって取得された瞬間との間の最大時間区間を計算することと、
・最大時間区間が所定の持続時間よりも大きい場合、線形遷移モデルを予め気低された制約付きゾノトープに適用することであて、それにより、事前制約付きのゾノトープを生成し、その他の場合、予め決定された制約付きゾノトープを事前制約付きゾノトープとして使用する、適用することと、
を更に含み得る。

【0022】

観測データは、ＧＮＳＳ受信機によって取得される以下のデータ：疑似レンジ、取得された搬送波位相、移動対象と別の移動対象と間のレンジの少なくとも１つを含み得る。

【0023】

運動データは、別の移動対象に対する移動対象の運動データを含み得る。

【0024】

第１の制約付きゾノトープは、区間分析を介したベクトル集合インバータアルゴリズムを使用して計算し得る。

【0025】

本発明の第２の態様は、コンピュータによって実行されるコンピュータプログラムであ
って、コンピュータに本発明の第１の態様の方法を実行させる命令を含むコンピュータプログラムである。

【図面の簡単な説明】

【0026】

【図1】移動対象の運動状態を包囲する集合を推定するデバイスを概略的に表す。

【図2】システムの状態ｘを包囲するＣゾノトープを推定する一例の方法のステップを示す。

【図3】図２に示す方法の伝播ステップの一連の動作を示す。

【図4】対象デバイス及びビーコンの位置ｐ及びｘ軸への距離ｄのグラフ表現である。

【図5】図２に示す方法の伝播、更新、及び簡約化ステップのスナップショットを示す。

【図6】対象デバイスの真の状態軌道のグラフ表現である。

【図7】対象デバイスのＣＺＥＳＭＦ包囲のグラフ表現である。

【図8】対象デバイスの位置区間包囲のグラフ表現である。

【図9】対象デバイスの速度区間包囲のグラフ表現である。

【図10】別の対象デバイスの位置のグラフ表現である。

【図11】ＶＳＩＶＩＡによって取得されたサブペイビングの区間包囲を示す。

【図12】事前Ｃゾノトープを図１０の対象デバイスの取得された事後Ｃゾノトープ及び簡約化事後Ｃゾノトープと共に表す。

【図13】観測タイミングノイズが増大する場合の事前及び事後Ｃゾノトープを表す。

【図14】図１０の対象デバイスの完全状態包囲軌道を表す。

【図15】観測タイミングノイズが増大する場合の図１０の対象デバイスの完全状態包囲軌道を表す。

【図16】２つの別個の手法によって取得されたレンジ線形化誤差区間測定の２つの半径のグラフ表現である。

【図17】シミュレーション中の２つの他の対象デバイスの位置を表す。

【図18】ＧＮＳＳ更新又はＧＮＳＳ＋ＵＷＢ更新を使用して速度ノルム予測及びＴ３Ｄ予測によって取得された集合ボリュームを表す。

【図19】ＣＺＥＳＭＦフィルタ及びＩＥＳＭＦフィルタを使用して取得された相対位置集合ボリュームを表す。

【発明を実施するための形態】

【0027】

以下の説明は以下のような構造を有する：表記（セクション１）を定義し、集合メンバシップ状態推定の問題を想起（セクション２）した後、制約付きゾノトープ（Ｃゾノトープ）に基づく計算効率的な状態集合推定法について説明する（セクション３）。次に、方法の３つの実際の適用を例として与える：
・単純な１Ｄ例（セクション４）、
・わずかにより込み入った２Ｄ位置特定問題（セクション５）、
・微分ＧＮＳＳ及びピアツーピア測距観測を使用した３Ｄピアツーピア位置特定ベンチマーク問題（セクション６）。

【0028】

１．表記（Notations)
この複雑性に対応するために、連続時間推定理論（ｃｏｎｔｉｎｕｏｕｓ－ｔｉｍｅ ｅｓｔｉｍａｔｉｏｎ ｔｈｅｏｒｙ）から馴染みのある表記ｘ（ｔ_ｌ）を以下使用する。フィルタ実行の離散性は、時間を索引付けることによって捕捉される。具体的には、ｔ_ｋは、事前（ア・プリオリ）及び事後（ア・ポステリオリ）状態集合が処理ユニットによって計算される時間を示す。ベクトルｘ（ｔ）では、時間ｔはこのベクトルの有効時間を示す。

【0029】

それに対応して、集合Ｚ（ｔ_ｋ）は、時間ｔ_ｋにおいて有効であり、すなわち、ｔ＝ｔ_ｋの場合、ｘ（ｔ）∈Ｚ（ｔ_ｋ）である。有効時間が区間［ｔ_ｋ］である場合、集合はこの区間の全ての時間で有効であり、すなわち、ｔ∈［ｔ_ｋ］である場合、ｘ（ｔ）∈Ｚ（ｔ_ｋ）である。

【0030】

それに加えて、行列Ｊ_ｎｘｍは、全要素が１であるｎ行及びｍ列の行列を示す。サイズｎの恒等行列はＩ_ｎで示され、サイズｎのゼロ正方行列は０_ｎで示される。

【0031】

２つのベクトルａ、ｂのドット積は＜ａ，ｂ＞として示される。

【0032】

区間行列は

【数1】

によって示され、式中、

【数2】

の要素はそれぞれ、［Ｍ］の区間要素の実数下界及び実数上界である。ｍ＝１且つｎ＞１である場合、区間ベクトル、すなわち［ｐ］を指す。スカラー区間は

【数3】

によって示され、式中、下界及び上界は

【数4】

である。間隔［Ｍ］の場合、Ｍは区間の実現値を示し、すなわち、Ｍ∈［Ｍ］である。

【0033】

付加誤差が仮定される場合、区間行列を［Ｍ］＝ｍｉｄ（Ｍ）＋［－ｒａｄ（Ｍ），ｒａｄ（Ｍ）］として中心－レンジ表記で書くことがより好都合であることができ、式中、ｍｉｄ（Ｍ）は［Ｍ］の中心点であり、区間半径行列は

【数5】

によって与えられる。［ｃ］＝ａ＋［ｂ］等、項及び区間項が同じ式に現れる場合は常に、項（ａ）は退化区間（［ａ，ａ］）を表し、その式中の全ての演算は区間演算であると仮定される。

【0034】

Ｃゾノトープは非線形変換下で非閉であるため、非線形系方程式の一次テイラー近似が以下使用される：セクション３．２及び３．３参照。これらの近似が、結果として生成される状態包囲の一貫性を損なわないことを保証するために、他の集合推定手法のように、線形包含の概念が使用される。例えば、Ｅ．Ｓｃｈｏｌｔｅらによる「Ａ ｎｏｎｌｉｎｅａｒ ｓｅｔ－ｍｅｍｂｅｒｓｈｉｐ ｆｉｌｔｅｒ ｆｏｒ ｏｎｌｉｎｅ ａｐｐｌｉｃａｔｉｏｎｓ」参照。非線形関数の線形包含は、特定の入力集合にわたり非線形関数の解を含むことが保証される集合値線形関数である。一例として、ｘ∈Ｘである、すなわち、ｘが何らかのコンパクト集合Ｘに含まれる関数

【数6】

を考える。線形関数

【数7】

は、全てのｘ∈Ｘについて

【数8】

である場合、式１の線形包含であり、すなわち、線形集合値関数の集合値解

【数9】

は常に非線形解を含む。

【0035】

２．問題文（Problem statement)
移動対象（ターゲット）の運動は、運動データを含む運動状態によって記述することができる。運動データは、例えば、対象の位置、対象の速度、対象の加速度、他のデータ、又はそれらの組合せを含み得る。運動状態は、移動対象が移動する場合、時間に伴って進化する。

【0036】

以下、簡潔にするために、移動対象１を「システム」と呼ぶことがあり、運動状態を「状態」と呼ぶことがある。

【0037】

図１を参照すると、移動対象の運動状態を包囲する集合を推定するデバイス１は、処理ユニット２及び観測ユニット４を備える。

【0038】

処理ユニット２は、後述する推定方法を実施するコード命令を含むコンピュータプログラムを実行するように構成される。このコンピュータプログラムは、実際にフィルタのように作用するため、以降「ＣＺＥＳＭＦフィルタ」（制約付きゾノトープ拡張集合メンバシップフィルタ：Ｃｏｎｓｔｒａｉｎｅｄ Ｚｏｎｏｔｏｐｅ Ｅｘｔｅｎｄｅｄ Ｓｅｔ－ｍｅｍｂｅｒｓｈｉｐ Ｆｉｌｔｅｒ）又は「フィルタ」と呼ばれる。

【0039】

フィルタは、所定の遷移モデル（以下、「伝播モデル」又は「ダイナミクスモデル」と呼ばれることもある）及び所定の観測モデルを使用するように構成される。

【0040】

処理ユニット２は、異なる瞬間での移動対象の運動状態を包囲する異なる集合を推定することができるように、上記異なる処理瞬間で動作するように構成される。

【0041】

処理ユニット２は任意のタイプであることができる：１つ又は複数のプロセッサ、１つ又は複数のマイクロプロセッサ、ＦＰＧＡ、ＡＳＩＣ等。

【0042】

観測ユニット４は、移動対象を観測するように構成された少なくとも１つのセンサを含む。したがって、観測ユニット４は、移動対象に関連する観測データを出力するように構成される。観測ユニットはまた、異なる瞬間に異なる観測を取得することができるように、上記異なる観測瞬間で動作するようにも構成される。

【0043】

集合メンバシップ状態推定の目的は、システムの状態ｘ（ｔ）が時間ｔに存在することが保証される領域

【数10】

を計算することである。

【0044】

２．１．状態ダイナミクス（State dynamics)
離散状態推定では、多くの場合、システム状態ダイナミクスが、前状態

【数11】

、システム入力

【数12】

、伝播モデル誤差

【数13】

、及び一定サンプリング時間Ｔ_ｓ＝ｔ_ｋ－ｔ_ｋ－１を有するタイプ

【数14】

の静的モデルによって支配されると仮定される。

【0045】

実際には、不規則な観測時間により不規則なサンプリング時間、すなわちフィルタでの不規則な伝播ホライゾン（ｐｒｏｐａｇａｔｉｏｎ ｈｏｒｉｚｏｎ）が生じる。サンプリング時間に応じて、異なる伝播モデルは緊密性が異なる状態包囲を生成し得る。ナビゲーションからの実際例は慣性推測航法位置特定であり、この場合、観測された加速度の二重積分が非常に緊密な短期位置包囲を生成するが、長期ではセンサバイアスが位置包囲の二次成長を生じさせる。逆に、車両ダイナミクスの単純な速度ノルム有界モデル（例えばセクション５参照）は、時間の経過に伴う位置包囲の一次成長を生じさせ、したがって、場合によっては、伝播ホライゾンが長いほど緊密な包囲を提供する。

【0046】

別の例はセクション６で扱われ、その例では、時間微分ＧＮＳＳ位相観測を使用した密な相対位置伝播が、ＧＮＳＳ受信機のサンプリング時間においてのみ利用可能であり、一方、ＧＮＳＳ観測間では、粗い速度ノルム有界伝播モデルしか利用できない。

【0047】

本発明では、非線形写像

【数15】

である非常に一般的な時変状態ダイナミクスモデルを仮定することにより、複数の伝播モデルに適合する。

【数16】

【0048】

引数ｔ_ｋ及びｔ_ｋ－１は、モデルが、ｔ_ｋ－１からｔ_ｋまでの時間ホライゾンにわたりシステム状態の伝播を実行すること及びｆ（）の時変性を示す。伝播モデルが時変であることが可能なことにより、フィルタは、各伝播ステップ１０２（図２参照）において利用可能な一揃いのモデルから保守性が最小のものを利用することができる。

【0049】

２．２．観測（Observations）
実際には、観測ユニット４によって行われる観測は多くの場合、完全には同期されず、また完全には周期的ではない。更なる複雑性として、例えば、ハードリアルタイム通信プロトコルを欠くことに起因して、観測時間は大まかでしか分からないことがある。

【0050】

観測時間の不確実性は、かすかであるが、既存の集合メンバシップ状態推定方式では対処されない深刻な問題を生じさせる。伝播状態包囲は、選択された伝播時間ｔ_ｋにおいてのみ状態を含むことが保証される。これを所与の状態包囲の有効時間

【数17】

と呼ぶ。観測は区間内の何らかの未知の時間に対応するため、この観測を用いて事前包囲を更新すると、ｔ_ｋが観測時間区間内に入る場合であっても、任意の時間で有効であると保証されない状態包囲が生じる。これは単に、

【数18】

及び観測に対応する状態空間制約が同じ時点を考慮しないことに起因する。したがって、観測時間の不確実性を無視すると、わずか１回の更新後に無効な状態包囲が生じる。この問題への改善策を以下に提示する。

【0051】

観測は、非線形写像

【数19】

によって支配されると仮定される。式中、

【数20】

はセンサ観測のベクトルであり、ｗ（ｔ）∈Ｗは対応する測定誤差であり、

【数21】

は観測された観測タイムスタンプであり、ｗ_ｔ（ｔ）∈［ｗ_ｔ］は観測時間誤差であり、ｔは真の観測時間である。観測タイムスタンプ誤差ｗ_ｔ（ｔ）は、タイミングジッタのみならず、例えばデジタルデータリンクにわたるパッケージ符号化／復号化遅延に起因した通信遅延も網羅する。全ての時間はウォールクロック時間である。

【0052】

なお、何らかの既知の集合Ｗ（ｔ）に含まれること以外の測定モデル誤差についての仮定は行われない。したがって、観測タイムスタンプ誤差は、ハードウェアノイズ、バイス、及びモデリング誤差を網羅する。観測モデルｈ（．．．）及び観測ベクトルのサイズｍ（ｔ）が時変であることができることにも留意する。

【0053】

ｙ（ｔ）はモデルにより提供される観測を示すが、センサから受信される実際も測定値はｚ（ｔ）で示される。

【0054】

３．Ｃゾノトープに基づく拡張集合メンバシップフィルタ（Extended set-membership filter based on C-zonotopes)
集合状態推定で使用される最も普及している計算効率的な集合表現は、楕円体変換単位超球及びゾノトープ変換単位超立方体等の基本幾何学的実体の線形変換である。超球及び超立方体の両方の固有の対称性は、非対称又は一般的に不規則形状の集合を密に近似する
楕円体及びゾノトープの能力を制限する。状況を更に悪化させるように、両方ともフィルタ更新の交叉演算下で非閉である。その結果、各交叉の結果の外近似を計算する必要がある。制約付きゾノトープはこれらの問題の両方を軽減する。制約付きゾノトープは、制約付き単位超立方体の線形変換であり、したがって、非対称であることができる。交叉下で閉であり、すなわち、２つのＣゾノトープの交叉は別のＣゾノトープであるが、複雑性がより高い。

【0055】

Ｃゾノトープ

【数22】

は、生成行列Ｇ、単位超立方体の中心ｃ及び制約パラメータＡ、ｂによって定義され、集合：

【数23】

を記述する。

【0056】

各制約は、単位超立方体と交叉する平面を表す。結果生成される交叉は次に、生成行列Ｇを介して射影され、最終的に中心ベクトルｃによって平行移動する。

【0057】

Ｃゾノトープは、凸多面体の代替のパラメータ化である。それに加えて、制約なしの各Ｃゾノトープは通常のゾノトープであるため、Ｃゾノトープは通常のゾノトープの超集合である。簡潔にするために、既存の文献によれば、凸多面体のＣゾノトープ表現はＧ－ｒｅｐとして示され、その半空間表現はＨ－ｒｅｐとして示される。

【0058】

ＣゾノトープＺの区間包囲は［Ｚ］で示され、区間［ａ］のＧ－ｒｅｐはＺ（［ａ］）で示される。

【0059】

図２を参照すると、システムの状態ｘを包囲するＣゾノトープを推定する方法は、時間ｔ_ｋ－１において状態を包囲する所定のＣゾノトープＺ_ｘ（ｔ_ｋ－１）に基づいて、フィルタによって実行される以下の主なステップを含む：
・所定のＣゾノトープが処理されて、時間ｔ_ｋにおいて状態を包囲する事前Ｃゾノトープ

【数24】

を取得する伝播ステップ１０２；
・事前制約付きゾノトープ

【数25】

が、観測ユニット４によって取得された観測データと結合されて、時間ｔ_ｋにおいて状態を包囲する事後ＣゾノトープＺ’_ｘ（ｔ_ｋ）を取得する更新ステップ１０４。

【0060】

時間ｔ_ｋにおいて状態ｘを包囲し、ｔ_ｋを含めたｔ_ｋまでの全ての観測を考慮に入れた事後Ｃゾノトープは、

【数26】

として示される。

【0061】

一方、更新を実行する前、時間ｔ_ｋにおいて状態ｘを包囲する事前Ｃゾノトープは、

【数27】

によって示される。

【0062】

伝播ステップ１０２及び更新ステップ１０４が線形遷移モデル及び線形観測モデルをそれぞれ使用し得、上記線形モデルが特定の点で非線形モデルを線形化することによって取得されることを以下に詳述する。

【0063】

さらに、方法は、フィルタが事後ＣゾノトープＺ’_ｘ（ｔ_ｋ）を、複雑性はより低いがＣゾノトープＺ’_ｘ（ｔ_ｋ）を含む別のＣゾノトープＺ_ｘ（ｔ_ｋ）に変換する簡約化ステップ１０６を更に含む。

【0064】

伝播ステップ１０２、更新ステップ１０４、及び簡約化ステップ１０６は再帰的に繰り返される。より正確には、方法は複数の再帰を含み、各再帰は伝播ステップ１０２、更新ステップ１０４、及び簡約化ステップ１０６を実行する。

【0065】

方法の最初の再帰には、初期Ｃゾノトープ

【数28】

が特定される初期化ステップ１００が先行する。

【0066】

ステップ１００、１０２、１０４、及び１０６についてこれより詳細に説明する。

【0067】

３．１．初期化（Initialization）
初期状態は、第１の事前状態包囲である何らかのＣゾノトープ

【数29】

に含まれることが分かっている。

【0068】

３．２．伝播（Propagation）
ｋ＞１の場合に実行される伝播ステップ１０２において、フィルタは事前状態集合

【数30】

を計算する。先の事後状態集合Ｚ_ｘ（ｔ_ｋ－１）を所与として、時間ｔ_ｋ＞ｔ_ｋ－１において状態を含むことが保証される。

【0069】

不確実な観測の特殊性として、事前Ｃゾノトープを計算する時間の選び方の問題が生じる。将来、フィルタ再帰が開始されるとき、フィルタ実行時間を補償することさえあり得る。しかしながら、再帰（伝播及び更新）の実行時間に応じて、ｔ_ｋを観測時間の近くに設定し、最後の伝播を実行して実行時間を補償するほうが賢明であり得る。

【0070】

幾つかの選択肢は、観測が公称的に周期的である場合、公称周期的観測時間又は観測時間区間の中心である。最初の選択肢の有利な効果として、フィルタ結果は周期的になり、観測時間不確実性をフィルタ内部に隠す。

【0071】

Ｃゾノトープは非線形写像下で非閉であるため、本発明は、拡張カルマンフィルタ（Ｅ
ＫＦ）と同様に、作用点

【数31】

についての状態ダイナミクスモデル式（５）の一次テイラー近似を実行することで開始する。先の事後状態の区間包囲及び入力包囲の中心が作用点として選ばれる：

【数32】

【0072】

Ｃゾノトープの中心はＣゾノトープ内にある必要さえないため、この選択はより小さな線形化誤差に繋がるが、線形化誤差を最小化はしない。

【0073】

線形化は、馴染みのある一次テイラー展開

【数33】

に繋がる。ここで、線形化誤差はε_ｆ（ｔ_ｋ）∈［ε_ｆ（ｔ_ｋ）］であり、

【数34】

である。

【0074】

対応する線形包含は

【数35】

によって与えられ、時間ｔ_ｋ－１から時間ｔ_ｋへの伝播が次に、

【数36】

として実行される。

【0075】

式１８における行列の乗算及び合算演算は、Ｊ．Ｋ．Ｓｃｏｔｔらによる「Ｃｏｎｓｔｒａｉｎｅｄ ｚｏｎｏｔｏｐｅｓ：Ａ ｎｅｗ ｔｏｏｌ ｆｏｒ ｓｅｔ－ｂａｓｅｄ ｅｓｔｉｍａｔｉｏｎ ａｎｄ ｆａｕｌｔ ｄｅｔｅｃｔｉｏｎ」において定義されるＣゾノトープへの集合演算である。線形化誤差及び非線形状態ダイナミクスモデルのモデル不確実性はＺ_ｑ（ｔ_ｋ－１）によって囲まれ、Ｚ_ｕ（ｔ_ｋ－１）は任意の不確実入力のＣゾノトープ包囲である。伝播不確実性項ｑ（ｔ_ｋ）は実際にＥＫＦのプロセスノイズに対応し、ＣゾノトープＺ_ｑ（ｔ_ｋ）はＥＫＦのプロセスノイズ共分散行列Ｑに対応する。

【0076】

［ε_ｆ（ｔ_ｋ）］は、ＣゾノトープＺ_ｘ（ｔ_ｋ－１）の異なる要素における遷移線形化誤差の集合として見なすことができる。Ｚ_ｘ（ｔ_ｋ－１）の所与の要素における遷移線形化誤差は、非線形遷移モデルを所与の要素に適用した結果生成される出力と線形遷移モデルをその所与の要素に適用した結果生成される出力との間の差である。

【0077】

それに加えて、先に述べた式Ｚ_ｑ（ｔ_ｋ－１）＋Ｚ（［ε_ｆ（ｔ_ｋ－１）］）＋Ｚ（［ｃ_ｆ（ｔ_ｋ－１］）は、全ての線形化誤差を含むＣゾノトープと見なし得る。さらに、式

【数37】

は、線形遷移モデルをＺ_ｘ（ｔ_ｋ－１）に適用した結果生成されるＣゾノトープと見なすことができる。事後ゾノトープは両式の和である。

【0078】

３．２．１．条件付き伝播（Conditional propagation）
再帰フィルタリング方式では、伝播ステップ１０２の目的は、システム状態についての観測が行われる時間－又はここではその時間の近くの時間－において有効である状態包囲を取得することである。

【0079】

本方法では、伝播を実行する必要がなく、フィルタ実行時間を短縮する２つの状況がある。

【0080】

第１の状況は、一連の同期ベクトル観測にわたり更新１０４が実行される場合、生じる。この場合、第１のベクトル要素を用いた更新１０４から生じる事後Ｃゾノトープは、観測ベクトルの残りの要素を用いた続く順次更新の事前Ｃゾノトープとして再使用することができる。

【0081】

それに加えて、観測が近い時間に到着する場合、伝播を実行する代わりに、観測区間を更に膨張させて（セクション３．３．１参照）、観測の有効時間を先の事後Ｃゾノトープの有効時間にわたり拡張することを許容することができる（保守性追加に関して）。この状況に適合するために、閾値

【数38】

がフィルタによって使用される。各伝播前、フィルタは、伝播を実行する必要があるか否かを判断する。

【0082】

以下のアルゴリズム１は、条件付き伝播を実際にいかに実施することができるかを示す：（図３も参照）。

【数39】

【0083】

３．３．更新（Update）
更新を実行するために、非線形観測式はその線形包含を形成するように線形化される。区間算術（ＩＡ）によって残り（すなわち、線形化誤差）を評価することにより、事前状態包囲にわたり有効であることが保証される。次に、観測ベクトルの各要素が、状態空間における帯状領域を定義する。事前Ｃゾノトープは、これら全ての帯状領域と一度に又は順次交叉することができる。楕円体及びゾノトープ等の普及している集合形式は、帯状領域との交叉下で非閉であり、したがって、各交叉結果の外近似を必要とし、順次交叉を準最適にする。他方、順次更新を用いる場合、各交叉後の観測式の再線形化は、緊密線形化誤差区間に起因してはるかに緊密な更新に繋がり得る。これらの２つの対立する効果は、関わる観測の非線形性の程度に依存し、したがって、用途依存である。

【0084】

Ｃゾノトープの独自の特徴として、帯状領域との交叉下で閉である。したがって、観測ベクトル全体を用いた一度での更新又は各スカラー要素との順次更新は同等の動作である。その結果、Ｃゾノトープを使用して、次の両方の効果から恩恵を受けることができる：すなわち、交叉最適性を失わない順次更新の単純性及び再線形化の緊密化効果。

【0085】

したがって、一般性を失うことなく、スカラー観測は以下のように考慮される：同期観測のベクトルは、フィルタ再帰を実行する前、同一のタイムスタンプを有する一連のスカラー観測に変換される。したがって、以下、一例としてより単純なスカラー観測式：

【数40】

が採用される。式中、ｔ_ｙ，ｋは観測ｋの真の観測時間であり、

【数41】

はタイミング誤差を受ける見掛けの観測時間である。

【0086】

区間値誤差を仮定すると、時間ｔ_ｙ，ｋにおいて有効な非線形観測写像の一次近似は、形態：

【数42】

を有する。式中、

【数43】

であり、式中、ε_ｈ（ｔ_ｋ）＝ε_ｈ，ｘ（ｔ_ｋ）＋ε_ｈ，ｕ（ｔ_ｋ）∈［ε_ｈ（ｔ_ｋ）］は事前状態

【数44】

及び入力Ｚ_ｕにわたる線形化誤差である。なお、ｔ_ｋは事前Ｃゾノトープの有効時間である。すると、対応する線形包含は

【数45】

によって与えられる。

【0087】

実際に、線形観測式及び事前Ｃゾノトープ内にある全ての状態に対するこの線形式の誤差を含む集合がここで計算される。この誤差は複数の成分を有する：
・Ｄ（ｔ_ｋ）Ｚ_ｕ（ｔ_ｋ）：システム入力、例えばオートパイロット又は外部摂動によって生成される制御入力の影響。ＣゾノトープＺ_ｕは可能な全ての入力又は摂動を含む。
・Ｚ［ε_ｈ（ｔ_ｋ）］：Ｃゾノトープに変換された線形化誤差区間。このゾノトープは、事前ゾノトープの異なる要素における全ての観測誤差を含む。事前ゾノトープの所与の要素における観測線形化誤差は、非線形観測モデルを所与の要素に適用した結果生成される出力と線形観測モデルをその所与の要素に適用した結果生成される出力との間の差である。
・Ｚ（［ｃ_ｈ（ｔ_ｋ）］）：線形化された観測式及び元の非線形観測式の測定ノイズ区間の全ての定数項を表すＣゾノトープ。

【0088】

更新を実行するために、２つの集合が交叉される。第１の集合は事前Ｃゾノトープによって与えられる（式１８）。第２の集合（式３５参照）は観測によって定義される。本発明では、観測ｚ（ｔ_ｙ，ｋ）が式２６の右辺に含まれることが分かっているため、式２６中のｙ（ｔ_ｙ，ｋ）をｚ（ｔ_ｙ，ｋ）で置換して、式３５を形成することができる。スカ
ラー観測では、式３５は状態空間における帯状領域を表す。図３参照。明らかに、第２の集合は、閉ではない（帯状領域が無限に延在する）ため、Ｃゾノトープではない。しかしながら、事前Ｃゾノトープとのその交叉はＣゾノトープ（事後Ｃゾノトープ）である。

【0089】

なお、線形化されるシステム方程式は近似であるが、対応する線形包含は、対応する非線形解を含むことが保証される。

【0090】

３．３．１．観測区間膨張（Observation interval inflation)
真の観測時間が有界不確実性を有することが分かっており、すなわち、区間［ｔ_ｙ，ｋ］内に入ることが分かっていることを想起する。換言すれば、観測はこの時間区間内のどこかの場所で有効である。観測及び事前Ｃゾノトープの有効時間は一致すると保証されないため、更新は実行されないことがある。

【0091】

この状況を救済する単純な手法として、観測誤差区間

【数46】

を膨張させ、観測の有効時間を拡張する。線形包含式２６から始まり、事前Ｃゾノトープの有効時間と観測の有効時間との間の時間は区間

【数47】

によって境される。式中、状態増分区間

【数48】

は、

【数49】

である場合、

【数50】

を満たす。拡張線形包含

【数51】

を計算することができる。

【0092】

膨張ステップの主な目的は、以下の問題の解決することである：観測の有効時間が、事前Ｃゾノトープの有効時間と一致することが保証されない。これを解決する基本概念は、観測誤差区間を膨張させて、その有効時間が事前Ｃゾノトープの有効時間を含むことを保証することである。これは、観測時間と

【数52】

の有効時間（式２７）との間の状態変動（式２８）の上界を計算することによって得られる。次に、既知の項は全て１つの膨張観測誤差区間に一緒にまとめられて、式３２に繋がり、これは、更新ステップ１０４においてフィルタによって使用することができる（式３６）。

【0093】

既知の項は全て１つの膨張観測誤差区間にひとまとめにされて、

【数53】

に簡易化され、式中、

【数54】

である。

【0094】

なお、線形化誤差ε_ｈ（ｔ_ｋ）はここで、膨張した事前状態集合

【数55】

にわたり計算する必要がある。実施のために、ε_ｈ（ｔ_ｋ）は区間包囲

【数56】

にわたり評価することができ、それにより区間算術を使用することが可能である。Ｃゾノトープ合算演算に関わるＣゾノトープ複雑性の増大に起因して、これは代替の区間

【数57】

を選ぶことよりも概して高速である。

【0095】

スカラー観測ｚ（ｔ_ｙ，ｋ）を用いると、

【数58】

及び

【数59】

を有し、式中、［ε］＝ｚ－［ε’］であり、Ｚ（［ε］）＝｛ｒａｄ（ε），ｍｉｄ（［ε］），０，０｝は区間［ε］のＧ－ｒｅｐである。

【数60】

であり、式中、

【数61】

である。伝播及び更新タイミングの概要については図３参照のこと。

【0096】

３．４．簡約化（Reduction)
式３６～４１から見て取ることができるように、各更新は状態Ｃゾノトープの複雑性－すなわち、生成元の数ｎ_ｇ及び制約の数ｎ_ｃ－を上げる。計算複雑性及びメモリフットプリントを制御するために、フィルタは、Ｚ_ｘ（ｔ_ｋ）の複雑性が、

【数62】

によって定義される選ばれた特定の複雑性閾値を超える場合、更新後に簡約化を実行する。

【0097】

以下のアルゴリズム２は、簡約化ステップ１０６を実際にいかに実施することができるかを示す。

【数63】

【0098】

以下のアルゴリズム３、４は、一実施形態によりフィルタによって実行される異なるステップを示す。

【数64】

【数65】

【0099】

４．１Ｄ適用（1D application）
この最初の例では、運動がデカルト系のｘ軸に制限される対象デバイス（又はエージェント）を考える。図４参照。ビーコンへのレンジ観測のみに基づいて対象デバイスの位置ｐ及び速度ｖを推定する。距離ｄだけｘ軸からずれて位置するビーコンは、観測式を非線形にし、したがって、ＣＺＥＳＭＦ等の非線形推定器を必要とする。この最小例では観測時間に不確かさはない。

【0100】

状態のダイナミクスｘ＝（ｐν）^Ｔは、

【数66】

によって支配され、式中、Ｔ_ｓ＝ｔ_ｋ－ｔ_ｋ－１であり、α（ｔ）は加速度である。

【0101】

レンジ観測は、

【数67】

によって支配され、式中、ｗ_ｋ∈［ｗ］＝［－０．２２４，０．１６８］ｍである。

【0102】

４．１．初期化（Initialization）
まず、対象の位置及び速度が

【数68】

に制限されると仮定する。

【0103】

４．２．伝播（Propagation）
対象加速度は

【数69】

及びＴ_ｓ＝ｔ_ｋ－ｔ_ｋ－１によって境されるため、ｘ（ｔ_ｋ）の包含は

【数70】

によって与えられ、式中、Ｚ_ｑ（ｔ_ｋ－１）＝ＱＺ（［α］）である。

【0104】

包含式４６は既に線形であるため、事前Ｃゾノトープは

【数71】

として直接得られる。

【0105】

４．３．更新（Update）
観測時間はこの例では不確実性を有さないため、ｔ_ｋ＝ｔ_ｙ，ｋであり、すなわち、状態包囲を観測の有効時間まで伝播させる。式４３から、線形包含

【数72】

が形成され、式中、

【数73】

であり、

【数74】

である。

【0106】

線形化誤差の包囲を計算する際、区間算術によってラグランジェの剰余を評価する一般的な手法は多くの場合、従属性に起因して過度に保守的な包囲に繋がる。ε_ｈを明示的に評価することは、より緊密な包囲を提供することができ、単調性を利用する場合、更に一層緊密になる。線形化誤差の明示的な式は、

【数75】

である。

【0107】

単調性をチェックするために、導関数

【数76】

を考える。

【0108】

事前Ｃゾノトープ

【数77】

の区間包囲にわたるε_ｈの緊密包囲を計算するには、その頂点及び

【数78】

によって与えられる可能な変曲点を評価するだけでよい。

【0109】

次に、更新は式３６に従って実行される。

【0110】

４．４．シミュレーション（Simulations）
伝播ステップ１０２、更新ステップ１０４、及び簡約化ステップ１０６のスナップショットを図５に表示する。図６及び図７はそれぞれ真の状態軌道及びＣＺＥＳＭＦ包囲を示す。図８及び図９はそれぞれ、対応する位置及び速度区間包囲を示す。各更新後、状態Ｃゾノトープは

【数79】

に簡約化されている。

【0111】

５．２Ｄ位置特定用途（2D localization application）
この例では、２つの静的測距ビーコンを有する、フィールドで移動する対象デバイスを考える（図１０参照）。これは、指定ゾーン内の花々の刈り取りを回避し、且つフィールドの周縁から出ないようにする必要がある自動芝刈り機又は近代倉庫内の自動リフトトラックであることができる。

【0112】

対象の位置ｘ＝（ｘ_１ｘ_２）^Ｔは

【数80】

によって支配され、式中、

【数81】

であり、対象の最高速度は

【数82】

である。この速度ノルム制約は、内部ダイナミクスのよりよい知識が入手可能ではない場合、移動対象デバイスの保守的であるが容易に取得可能なダイナミックモデルである。

【0113】

対象デバイスは定期的に、不確実な位置

【数83】

に配置された両ビーコンへの測距観測を受信し、式中、

【数84】

である。

【0114】

各ビーコンｉ_ｋへの測距観測は、

【数85】

によって支配される。

【0115】

これらの２つの観測は公称的に周波数１Ｈｚにおいて同期する。ウォールクロック時間との測距センサの不完全な同期は、有界タイミングジッタｗ_ｔ（ｔ）∈［－１０^－２１０^－２］ｓを生じさせる。すなわち、一対の測距観測は概ね毎秒到着するが、厳密に毎秒到着するわけではない。

【0116】

周期性状態包囲を得るために、フィルタ再帰が観測時間区間の上界、すなわち、（ｔ_１，ｔ_２，ｔ_３，ｔ_４，．．．）＝（０．０１ｓ，０．０１ｓ，１．０１ｓ，１．０１ｓ，．．．）においてトリガーされる。

【0117】

５．１．初期化（Update）
対象デバイスは最初、フィールド内部にある。これにより、フィールドを網羅する区間に対応するＣゾノトープを計算することができる。しかしながら、この極めて保守的な初期包囲は極めて保守的なフィルタ解に繋がる。したがって、ＶＳＩＶＩＡを用いて初期状態のより緊密な包囲が計算される。このワンタイムステップは、二次元ではそれほど計算的に高価ではない（０．８ｓ及びε＝０．１ｍで２３３個の処理済みボックス）。より高次の状態での推定問題のために、ＶＳＩＶＩＡを制約伝播と結合して、十分に緊密な初期状態包囲を許容可能な時間で計算することが好ましいことがある。

【0118】

ＶＳＩＶＩＡによって得られるサブペイビング（ｓｕｂｐａｖｉｎｇ）の区間包囲に対応するＣゾノトープ（図１１では「解」（ｓｏｌｕｔｉｏｎ）及び「未定義」（ｕｎｄｅｆｉｎｅｄ）と示される）が、フィルタの初期化に使用される。

【0119】

５．２．伝播（Propagation）
状態において状態ダイナミクスは既に線形である。したがって、伝播は

【数86】

によって直接実行することができる。

【0120】

式５８によって定義される状態ダイナミクス不確実球（ｓｔａｔｅ ｄｙｎａｍｉｃｓ
ｕｎｃｅｒｔａｉｎｔｙ ｂａｌｌ）は、Ｚ_ｓ（ｔｋ－１）を得るためのＣゾノトープによる外近似である（以下の付録Ａ参照）。その結果生成された事前Ｃゾノトープの一例を図１２に表示する。

【数87】

がｔ_９において真の位置をいかに含むかに留意する。２つの測距観測が約１Ｈｚで到着し、あらゆる観測がフィルタ再帰インデックスｋを増大させるため、ｔ_９が、ｔ＝４ｓの近くで実行された最初の伝播更新シーケンスに対応することにも留意する。

【0121】

５．３．更新
これは、対応する測距ビーコンｉ_ｋの公称系において各更新を実行することを簡易化する。ビーコン系における事前Ｃゾノトープは、

【数88】

によって与えられるため、各測距観測は

【数89】

によって支配される。

【0122】

更に簡易化するために、

【数90】

によって小さなビーコン不確実性を観測誤差ｗ_ｒ，ｋに吸収させ、

【数91】

を生成させ、式中、

【数92】

である。

【0123】

式６７から、線形包含

【数93】

を形成し、式中、

【数94】

であり、且つ

【数95】

である。

【0124】

５．３．１．観測区間膨張（Observation interval inflation）
タイミングジッタを補償する前、状態増分区間［Δｘ（ｔ_ｋ）］

【数96】

が速度ノルム界

【数97】

から得られ、膨張観測包含が以下のように計算される。

【数98】

【0125】

式３６～４１を用いてビーコン系における事後Ｃゾノトープ

【数99】

を計算した後、事後Ｃゾノトープはナビゲーション系に平行移動する。

【数100】

【0126】

図１２は、ｔ_９における取得された事後Ｃゾノトープ及び簡約化された事後Ｃゾノトープを示す。図１３は、タイミングノイズ値を増大させて得られた結果を示す。なお、観測区間膨張の増大に応じて事後Ｃゾノトープはより大きくなる。

【0127】

完全状態包囲軌道を図１４に表示し、図１５では、観測タイミングノイズを増大した場合を示す。各更新後、状態Ｃゾノトープは

【数101】

に簡約化されている。Δｔを０．１ｓに設定することにより、伝播が測距観測のあらゆる対で１回だけ実行されることが保証され、したがって、フィルタ実行時間を短縮する（条件付き伝播）。

【0128】

６．３Ｄピアツーピア位置特定への適用（Application to 3D peer-to-peer localization）
この例では、関心のある状態は、密なフォーメーション飛行での無人機等の２つのエージェント間の相対位置である。多くの航空宇宙用途では、空中給油、エネルギー効率的フ
ォーメーション飛行、及び小型衛星群等の相対位置の保証が求められる。この例では、生のＧＮＳＳコード及び搬送波位相観測並びにＵＷＢピアツーピア測距観測が得られる。両方の種類のセンサとも安価で非常に低質量－数十ｇ－であり、且つ民生部品（ＣＯＴＳ）として容易に入手可能である。ＣＺＥＳＭＦを実用的に高い関連性があるこの位置特定問題にいかに適用することができるかを示し、既存の集団メンバシップ状態推定方式と比較して、現実的なシミューレションベンチマークにおいてその保守性を評価する。

【0129】

以下、明らかに読みやすさを改善することが可能な場合は常に、時間インデックスを省略する。

【0130】

６．１．初期化（Initialization）
ＧＮＳＳコード二重位相差観測の高度の線形性に起因して、フィルタは単に、初期相対位置についての非常に保守的な推測を用いて初期化することができる。

【0131】

６．２．伝播（Propagation）
相対位置ｘが保証された集合を２台の車両ＡとＢとの間で伝播させる。
ｘ＝ｐ_Ｂ－ｐ_Ａ

【0132】

中程度の保守性の集団値ダイナミック車両モデルは、大方、確立された集団理論的システム識別法がないことに起因して、利用可能なことは希である。したがって、車両に依存しない伝播モデルがかなり望ましい。集団メンバシップ慣性航法は、更なる研究の価値がある一選択肢であるが、本明細書ではこれ以上掘り下げない。代わりに、利用可能な場合、時間差ＧＮＳＳ搬送波二重位相差観測及びＧＮＳＳ観測がない場合には大まかな速度ノルム制約モデルに頼る。スタンドアロン時間差搬送波位相観測は元々、精密な飛行軌道再構築のために受信、導入されていた。それを足場として、時間差二重位相差（Ｔ３Ｄ：Ｔｉｍｅ－Ｄｉｆｆｅｒｅｎｃｅｄ Ｄｏｕｂｌｅ ｐｈａｓｅ Ｄｉｆｆｅｒｅｎｃｅ）観測がまず、緊密な集合値相対軌道追跡のために最近提案された。ミリメートルレベルの搬送波ノイズにより、ＧＮＳＳ観測間の位置増分の非常に狭い包囲を計算することができる。相対位置推定に適用される場合、伝播状態集合は保守性が低く、更新ステップにおける線形化誤差集合を小さくし、最終的に保守性のより低い事後状態集合を生じさせる。

【0133】

６．２．１．Ｔ３Ｄ観測（T3D observations）
搬送波位相観測の単一差の形成は、大気誤差及び衛星クロック誤差の大部分をなくす。

【0134】

２つの受信機Ａ及びＢ並びに衛星Ｓの単一差搬送波位相観測は、

【数102】

によって与えられ、式中、ε_ΔΦは、ハードウェアノイズ及びマルチパスに起因した残余誤差を捕捉し、ΔＴは差動受信機クロック誤差であり、ｃは光速であり、λは搬送波波長であり、ΔＮは差動搬送波周期曖昧性である。

【0135】

ｄ_ｓ＝ｐ_Ａ－ｐ_Ｓ及びｐ_Ｂ＝ｐ_Ａ＋ｘを使用して、よりコンパクトに

【数103】

と書く。

【0136】

差動受信機クロック誤差は、基準衛星Ｒに関して二重差を形成することによってなくなる。

【数104】

【0137】

したがって、二重差搬送波位相観測は、受信機Ａ及びＢの相対位置並びに絶対受信機位置と衛星軌道位置との間のベクトルの非線形関数である。

【0138】

受信機が２つの後続エポックにわたり位相ロックを維持すると仮定すると、曖昧性は、２つの後続二重差観測の差を形成することによってなくなる定数であり、Ｔ３Ｄ観測

【数105】

を形成する。

【0139】

以下、位置増分Δｘ＝ｘ（ｔ_ｋ）－ｘ（ｔ_ｋ－１）を包囲するＣゾノトープをいかに計算することができるかを示す。

【0140】

直接線形化し（式７８を使用して）、Ｅ．Ｓｃｈｏｌｔｅらによる「Ａ ｎｏｎｌｉｎｅａｒ ｓｅｔ－ｍｅｍｂｅｒｓｈｉｐ ｆｉｌｔｅｒ ｆｏｒ ｏｎｌｉｎｅ ａｐｐｌｉｃａｔｉｏｎｓ」に提案されているように区間算術を使用してラプラス剰余の区間包囲を計算することは、重い従属性に起因して非常に大きな過剰近似をもたらす。

【0141】

この問題は、単一差分式７６を代わりに線形化することによって克服される。

【0142】

次に、ＩＡを使用してラグランジェ剰余が評価され、その結果生成された区間はハードウェア誤差区間ε_ΔΦに追加される。

【0143】

必要なヤコビアン及びヘシアンの計算を簡易化するために、式７６を

【数106】

と書くことが好都合であり、式中、拡張状態ベクトル

【数107】

及び

【数108】

である。

【0144】

次に、式８１のヤコビアン及びヘシアンは、

【数109】

として好都合に形成することができる。

【0145】

したがって、線形化及びパラメトリック不確実性（衛星の不確実な軌道位置）に起因した観測モデル誤差の区間境界は、

【数110】

を評価することによって取得することができる。

【0146】

線形単一差観測の形成に進む。

【数111】

【0147】

分離が小さい場合、［Ｒ_ΔΦ］は一般に非常に狭く、ミリメートルのオーダであるため、約ｘ＝０と線形化し、ＲＴＫ文献から馴染みのあるよりコンパクトな形態

【数112】

を有することができ、式中、

【数113】

である。次に、従属性がない線形二重差観測モデルを形成することができ、

【数114】

式中、

【数115】

はここで、観測及び線形化誤差に一緒にまとめられる。

【0148】

観測時間ｔ_ｋ－１及びｔ_ｋをそれぞれ有する２つの連続搬送波位相観測の差を採用する
と、曖昧性定数項がなくなる。Δｘ（ｔ_ｋ）＝ｘ（ｔ_ｋ）－ｘ（ｔ_ｋ－１）を導入することで、Ｔ３Ｄ観測は、

【数116】

によって与えられ、式中、

【数117】

である。

【0149】

ここで、幾何行列

【数118】

が、受信機と衛星との間の距離が大きいことに起因して２つの連続時間ステップ間で略変わらず、非常に狭い増分項

【数119】

をもたらすことを利用する。２つの時間ステップ間の増分変位が線形の観測を得る。この観測を利用して、以下のように増分変位の非常に正確な包囲を計算する。項を一緒にまとめて、式９６を線形包含

【数120】

として書く。

【0150】

Ｔ３Ｄ観測

【数121】

は、半空間表現が

【数122】

によって与えられるポリトープＰにΔｘを制約する。

【0151】

次に、対応するＣゾノトープＺΔｘ（ｔ_ｋ）が、Ｊ．Ｋ．Ｓｃｏｔｔらによる「Ｃｏｎｓｔｒａｉｎｅｄ ｚｏｎｏｔｏｐｅｓ：Ａ ｎｅｗ ｔｏｏｌ ｆｏｒ ｓｅｔ－ｂａｓｅｄ ｅｓｔｉｍａｔｉｏｎ ａｎｄ ｆａｕｌｔ ｄｅｔｅｃｔｉｏｎ」に開示され
る定理１によって与えられる。これは、Ｆ（ｔ_ｋ－１）＝Ｉ_３、Ｂ（ｔ_ｋ－１）＝Ｉ_３、Ｚ_Δｘ（ｔ_ｋ）、Ｚ_ｑ（ｔ_ｋ－１）＝０を設定することによってフィルタ伝播に使用される。

【0152】

６．２．２．速度ノルム伝播（Velocity-norm propagation)
この適用例では、車両間測距観測はＧＮＳＳ観測と同期されず、より高率で到着する。測距観測を用いて更新する前、必要な伝播を実行するために、ノルム有界相対速度に基づく粗いモデルを使用する。伝播ホライゾンｔ_ｋ－ｔ_ｋ－１では、ノルム有界ν_ｍａｘは波形ｒ＝ν_ｍａｘ（ｔ_ｋ－ｔ_ｋ－１）を有する球を定義する。ユーザ定義の複雑性の対称ＣゾノトープＺ^Ｓによってこの球を外近似する。次に、伝播が、Ｆ（ｔ_ｋ－１）＝Ｉ_３、Ｂ（ｔ_ｋ－１）＝Ｉ_３、Ｚ_ｕ（ｔ_ｋ－１）＝０、Ｚ_ｑ（ｔ_ｋ－１）＝Ｚ^Ｓを設定することによって実行される。

【0153】

６．３．更新（Update）
６．３．３．ＧＮＳＳ更新（GNSS update）
ＧＮＳＳコード位相単一差観測は、

【数123】

によって支配され、式中、ε_Δρは、ハードウェアノイズに起因した残余誤差及び大気残余誤差を捕捉し、ΔＴは差動受信機クロック誤差であり、ｃは光速である。

【0154】

ｄ_ｓ＝ｐ_Ａ－ｐ_Ｓ及びｐ_Ｂ＝ｐ_Ａ＋ｘを使用して、よりコンパクトに

【数124】

と書く。

【0155】

差動受信機クロック誤差は、基準衛星Ｒに関して二重差を形成することによってなくなる。

【数125】

【0156】

したがって、搬送波位相二重差と全く同じように、二重差コード位相観測は、Ａ及びＢの相対位置並びに絶対車両位置と衛星軌道位置との間のベクトルの非線形関数である。搬送波位相観測と同じ方策に従い、コード位相二重差の線形包含を形成して、

【数126】

を生成し、状態包囲の更新を可能にする。

【0157】

６．３．２．測距更新（Ranging Update）
車両間測距観測は、

【数127】

によって与えられ、式中、各体座標系におけるＧＮＳＳアンテナ中心に対するＵＷＢ測距アンテナ中心の位置ｒ_Ａ、ｒ_Ｂであり、デバイスハードウェア測定誤差

【数128】

である。式１０６は、両ＵＡＳの属性Ｒ_Ａ、Ｒ_Ｂの不確実性及び不確かなアンテナ位置ｒ_Ａ、ｒ_Ｂに起因して、パラメトリック不確実性を受ける。

【0158】

以下、ｒ_Ａ、ｒ_Ｂは小さく、すなわち、測距アンテナはＧＮＳＳアンテナに近く、測距観測に対する相加効果の区間境界が特定され、結合レンジ誤差ｗ_ｒにまとめられたと仮定する：

【数129】

【0159】

式１０７を線形化するために、まず、レンジ観測を

【数130】

として書く。

【0160】

ＧＮＳＳ搬送波位相観測と同じ趣旨で式１０８の一次区間テイラー展開が続き、

【数131】

式中、ヤコビアン及びヘシアン

【数132】

である。

【0161】

幾何学的に、この線形化レンジ観測は勾配ベクトル

【数133】

に射影された相対位置ベクトルに対応し、一方、非線形勾配レンジは球に対応する。大きな距離では近似は良好であるが、球の曲率が大きくなるほど、車両は互いに近くなるため、導入される線形化誤差は大きくなり、緊密な包囲を更に一層重要にする。

【0162】

区間算術を使用して式１０９のラグランジェ剰余項を評価する一般的な手法を適用することは、ｘにおける多重従属性に起因して大きな過剰近似を提供する。

【0163】

線形化誤差区間が広いほど、線形化観測が保証付き状態集合の収縮に寄与することができる可能性は低くなる。

【0164】

以下、非線形観測式１０７とその線形化との間の誤差を明示的に評価することにより、単調性を利用して、より緊密な境界を得ることができることを示す。

【0165】

明示的な線形化誤差区間は、

【数134】

によって与えられる。

【0166】

［ｘ］に関する式１１１の単調性を確認するために、区間値勾配

【数135】

を考える。

【0167】

幾何学的に、これは２つの単位ベクトル間の差である。したがって、（１１２）の要素の符号は［ｘ］の部分区間にわたり一定であるため、式（１１１）は区分単調である。

【数136】

であるため、これらの部分区間の数は多くとも２^３であり、多くとも２７の独特な頂点（ｕｎｉｑｕｅ ｖｅｒｔｅｘ）に対応する。これらの頂点のそれぞれで式（１１１）を評価することにより、［ε_ｒ］の厳密な包囲を計算することができる。区分単調性を利用することの利点を定量化するために、両手法を用いて区間を計算し、それぞれの半径を比較する。図１６は、事前状態Ｃゾノトープ区間包囲から引き出された多数のサンプルについて式（１１１）を評価することにより得られた区間の内近似とともに、ベンチマークミッション（セクション６．４参照）の過程でのこの測定を表示する。

【0168】

なお、一般的な方式は、区分単調を利用する方式の最高で４倍大きな区間半径をもたらし、区分単調を利用する方式は、サンプリングにより得られた区間をより密に近似する。

【0169】

６．４．シミュレーション（Simulations）
ベンチマークＵＡＳ操作でフィルタを評価した。これは２つのＵＡＳ：リーダ及びフォロワからなる。フォロワは、円軌道上でｘ軸の右から左にその相対位置を変える。図１７参照。両エージェントの相対位置を１Ｈｚで追跡する。シミュレートされたＧＮＳＳ及び測距観測には記録観測誤差が混入している。ＧＮＳＳ観測では、これらは２つの低コスト単一周波数ＧＰＳ受信機のゼロベースライン測定から取得される。測距観測誤差は、既知の距離に配置された２つの静的ＵＷＢ測距ビーコンから取得される。

【0170】

６．４．１．集合ボリューム（Set Volumes）
３つ全ての状況において、ＣＺＥＳＭＦフィルタは、ＩＥＳＭＦ（反復拡張集合メンバシップフィルタ）及びＶＳＩＶＩＡ（ＳＩＶＩＡアルゴリズムのベクトル実装）の両方よりも小さな集合ボリュームを提供する。粗い速度ノルム境界モデルと比較した、伝播にＴ３Ｄ観測を使用することの利点も評価した。図１８が示すように、Ｃゾノトープボリュームは、伝播ステップに搬送波位相情報を使用することにより、略１桁小さくすることができる。

【0171】

６．５．衛星コンステレーションの影響の受けやすさ（Sensitivity to satellite constellations）
相対位置包囲の達成可能なサイズ及び形状は、馴染みのある精度低下率（ＤＯＰ）精度尺度によって捕捉されるベースライン衛星ジオメトリに依存する。したがって、提示され
た集合メンバシップフィルタの相対保守性が衛星コンステレーションの変動に関してロバストであるか否かを知ることに興味がある。ランダムに選択された１組のＧＰＳ軌道記録にわたりベンチマーク操作のシミュレーションを実行した。衛星位置は、書く軌道記録の最初の２分のスプライン近似から計算した。図１９は、その結果生成された、ＶＳＩＶＩＡによって得られたサブペイビングボリュームを用いて正規化されたＣＺＥＳＭＦ及びＩＥＳＭＦの相対位置集合ボリュームを示す。初期遷移後、ＣＺＥＳＭＦは、全てのシミュレーションにわたりＶＳＩＶＩＡよりも低い保守性を示すとともに、ＧＮＳＳ軌道に対して非常に小さな変動を示す。ＩＥＳＭＦははるかに保守的であり、衛星軌道に対して同様に低い影響の受けやすさを有する。これらのシミュレーションは、ＣＺＥＳＭＦによって計算された包囲の観測された優れた緊密性が、フィルタのオンライン実装に持ち越されると予期することができることのいくらかの信頼を提供する。

【0172】

７．まとめ（Conclusion）
Ｃゾノトープロバスト性に基づく、非線形システムに関して上述した集合状態推定方法は、シミュレーション状況において、ＶＳＩＶＩＡに基づく非線形オフライン集合推定器及び最近の高速楕円体集合状態推定器に勝る（集合ボリュームに関して）。

【0173】

提案されたフィルタは、道路車両の安全な視覚的位置特定への現在の関心に鑑みて心躍る視点である同時位置特定及び地図作成（ＳＬＡＭ：Ｓｉｍｕｌｔａｎｅｏｕｓ Ｌｏｃａｌｉｚａｔｉｏｎ Ａｎｄ Ｍａｐｐｉｎｇ）を含め、多種多様な推定問題に適用することができる。

【0174】

付録Ａ：Ｃゾノトープを用いた二次元又は三次元球の近似
ｃを中心とした半径ｒの球を近似するために、Ｈ－ｒｅｐにおける近似を計算し、それをＣゾノトープに変換する。

【0175】

θは天頂角であり、φは方位角であり、ｌは複雑性パラメータである、区間

【数137】

に均等に離間された角度ベクトルθ＝（θ_１，θ_２，・・・，θ_ｌ）及び［０，２π］に均等に離間されたφ＝（φ_１，φ_２，・・・，φ_ｌ）を所与として、角度の各組合せ（θ_ｉ，φ_ｊ）について、対応する法線ベクトル

【数138】

を計算する。

【0176】

これらの法線ベクトルの各々は半空間を定義する。次に、Ｊ．Ｋ．Ｓｃｏｔｔらによる「Ｃｏｎｓｔｒａｉｎｅｄ ｚｏｎｏｔｏｐｅｓ：Ａ ｎｅｗ ｔｏｏｌ ｆｏｒ ｓｅｔ－ｂａｓｅｄ ｅｓｔｉｍａｔｉｏｎ ａｎｄ ｆａｕｌｔ ｄｅｔｅｃｔｉｏｎ」に与えられる方法を使用して、Ｈ－ｒｅｐにおける対応するポリトープを容易に書き、

【数139】

それをＧ－ｒｅｐに変換することができる。

【0177】

二次元では、同じ手法が、［０，２π］に均等に離間された角度の１つの集合θ＝（θ_１，θ_２，・・・，θ_ｌ）にわたり実行される。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】