(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-03-21
(45)【発行日】2024-03-29
(54)【発明の名称】乱気流センシングシステム、航空機及び乱気流センシング方法
(51)【国際特許分類】
B64C 13/16 20060101AFI20240322BHJP
B64D 45/00 20060101ALI20240322BHJP
G01W 1/00 20060101ALI20240322BHJP
G01S 17/95 20060101ALI20240322BHJP
【ＦＩ】
B64C13/16 A
B64D45/00 A
G01W1/00 Z
G01S17/95
【請求項の数】
6
(21)【出願番号】P 2019206279
(22)【出願日】2019-11-14
(65)【公開番号】P2021079731
(43)【公開日】2021-05-27
【審査請求日】2022-09-21
【新規性喪失の例外の表示】特許法第３０条第２項適用 第５６回飛行機シンポジウム（山形テレサ、平成３０年１１月１４日開催） ［刊行物等］ 第５６回飛行機シンポジウム講演集ＣＤ－Ｒの写し ［刊行物等］ 第５６回飛行機シンポジウム開催会告 ［刊行物等］ 第５６回飛行機シンポジウムプログラム ［刊行物等］ ＩＥ０１ 地形性乱気流環境での無人機の飛行安全性に関する研究（第二報）
(73)【特許権者】
【識別番号】503361400
【氏名又は名称】国立研究開発法人宇宙航空研究開発機構
(74)【代理人】
【識別番号】110003339
【氏名又は名称】弁理士法人南青山国際特許事務所
(74)【代理人】
【識別番号】100104215
【弁理士】
【氏名又は名称】大森 純一
(74)【代理人】
【識別番号】100196575
【弁理士】
【氏名又は名称】高橋 満
(74)【代理人】
【識別番号】100168181
【弁理士】
【氏名又は名称】中村 哲平
(74)【代理人】
【識別番号】100117330
【弁理士】
【氏名又は名称】折居 章
(74)【代理人】
【識別番号】100160989
【弁理士】
【氏名又は名称】関根 正好
(74)【代理人】
【識別番号】100168745
【弁理士】
【氏名又は名称】金子 彩子
(74)【代理人】
【識別番号】100176131
【弁理士】
【氏名又は名称】金山 慎太郎
(74)【代理人】
【識別番号】100197398
【弁理士】
【氏名又は名称】千葉 絢子
(74)【代理人】
【識別番号】100197619
【弁理士】
【氏名又は名称】白鹿 智久
(72)【発明者】
【氏名】牧 緑
(72)【発明者】
【氏名】濱田 吉郎
(72)【発明者】
【氏名】加藤 博司
(72)【発明者】
【氏名】菊地 亮太
【審査官】山本 賢明
(56)【参考文献】
【文献】特開２０１９－１８２０７５（ＪＰ，Ａ）
【文献】特開２００７－３１７６５６（ＪＰ，Ａ）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｂ６４Ｃ １３／１６
Ｂ６４Ｄ ４５／００
Ｇ０１Ｗ １／００
Ｇ０１Ｓ １７／９５
(57)【特許請求の範囲】
【請求項1】
航空機に搭載されたドップラーライダーと、前記ドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測する乱気流予測部と、
前記航空機の主翼表面の複数点に配置された複数の静圧計測部とを具備し、
前記複数の静圧計測部により計測された各前記点の計測値に基づいて、前記予測される乱気流の情報を補償する
乱気流センシングシステムであって、
前記ドップラーライダーによる計測値に基づき、前記予測される航空機前方の乱気流に応じた直線渦糸の回転軸位置と強度を推定する渦推定部と、
前記推定された各直線渦糸の回転軸位置と強度に基づき前記航空機の主翼表面の第１の圧力場を推定する第１の圧力場推定部と、
前記複数の静圧計測部により計測された各前記点の計測値から、前記主翼表面の第２の圧力場を推定する第２の圧力場推定部とを具備し、
前記第１の圧力場と前記第２の圧力場との一致度を求め、当該一致度に基づき、前記予測される乱気流の情報を補償する
乱気流センシングシステム。
【請求項2】
航空機に搭載されたドップラーライダーと、前記ドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測する乱気流予測部と、
前記航空機の主翼表面の複数点に配置された複数の静圧計測部とを具備し、
前記複数の静圧計測部により計測された各前記点の計測値に基づいて、前記予測される乱気流の情報を補償する
乱気流センシングシステムであって、
前記ドップラーライダーによる計測値から、Δｔ秒先で機体に誘起する未来の第１の鉛直方向風速度を予測する第１の鉛直面風速場予測部と、
前記ドップラーライダーによる計測値及び前記複数の静圧計測部により計測された各前記点の計測値から、前方渦がΔｔ秒先で機体に誘起する未来の第２の鉛直方向風速度を予測する第２の鉛直面風速場予測部と、
前記第１の鉛直方向風速度と前記第２の鉛直方向風速度との相関度を求め、当該相関度に基づき、前記予測される乱気流の情報を補償する
乱気流センシングシステム。
【請求項3】
請求項２に記載の乱気流センシングシステムであって、既通過領域における相関度を評価し、当該評価結果に応じて、前記予測される乱気流の情報の補償の是非を決める
乱気流センシングシステム。
【請求項4】
請求項１乃至３に記載の乱気流センシングシステムを搭載した航空機。
【請求項5】
航空機に搭載されたドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測し、前記航空機の主翼表面の複数点の静圧を計測し、
前記計測した各前記点の静圧の計測値に基づいて、前記予測される乱気流の情報を補償する
乱気流センシング方法であって、
前記ドップラーライダーによる計測値に基づき、前記予測される航空機前方の乱気流に応じた直線渦糸の回転軸位置と強度を推定し、
前記推定された各直線渦糸の回転軸位置と強度に基づき前記航空機の主翼表面の第１の圧力場を推定し、
前記計測した各前記点の静圧の計測値から、前記主翼表面の第２の圧力場を推定し、
前記第１の圧力場と前記第２の圧力場との一致度を求め、当該一致度に基づき、前記予測される乱気流の情報を補償する
乱気流センシング方法。
【請求項6】
航空機に搭載されたドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測し、前記航空機の主翼表面の複数点の静圧を計測し、
前記計測した各前記点の静圧の計測値に基づいて、前記予測される乱気流の情報を補償する
乱気流センシング方法であって、
前記ドップラーライダーによる計測値から、Δｔ秒先で機体に誘起する未来の第１の鉛直方向風速度を予測し、
前記ドップラーライダーによる計測値及び前記計測した各前記点の静圧の計測値から、前方渦がΔｔ秒先で機体に誘起する未来の第２の鉛直方向風速度を予測し、
前記第１の鉛直方向風速度と前記第２の鉛直方向風速度との相関度を求め、当該相関度に基づき、前記予測される乱気流の情報を補償する
乱気流センシング方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、航空機に搭載される乱気流センシングシステム及び乱気流センシング方法に関する。本発明は、このような乱気流センシングシステムを搭載する航空機に関する。
【背景技術】
【0002】
航空機事故の半数以上は乱気流に起因しており、耐乱気流安全技術の必要性は論をまたない。特に晴天乱気流（ＣＡＴ）は気象レーダでは検知できず、寿命が短いため予報も困難である。この問題に対して、本出願人は、飛行機に搭載可能な晴天乱気流の事前検知用小型ドップラーライダーの開発を進め、危険空域回避、警報発出に役立つことを実証した。加えて、回避飛行が間に合わない場合は、乱気流空域通過時の機体動揺を可能な限り抑え、被害を最小化する飛行制御技術も必要である。
【0003】
これに関連する技術として、特許文献１や特許文献２、非特許文献１、非特許文献２、非特許文献３などが開示されている。
【先行技術文献】
【特許文献】
【0004】
【文献】特開２０１７－６７６８０号公報
【文献】特開２０１７－２１５５８号公報
【非特許文献】
【0005】
【文献】http://www.aero.jaxa.jp/publication/magazine/pdf/fp_no13.pdf
【文献】http://www.aero.jaxa.jp/research/star/safeavio/news120831.html
【文献】https://www.jstage.jst.go.jp/article/jjsass1969/34/390/34_390_369/_pdf/-char/ja
【発明の概要】
【発明が解決しようとする課題】
【0006】
晴天乱気流は、対流圏の任意の高度で起こり得るが、特にジェット気流（高度8000～13000m）の周辺で発生する晴天乱気流は、巡行する航空機の安全運行に脅威となっている。乱気流遭遇時の機体動揺を顕著に低減するためには、慣性センサ情報のフィードバック制御だけでなく、ドップラーライダーを用いた前方風速情報の先読みにもとづく予見フィードフォワード制御が有効である。ただし、達成性能は予見情報の精度に依存し、特に予見風速変動の位相が実際と大きくずれていた場合、動揺を増幅、不安定化につながる危険性がある。
【0007】
以上のような事情に鑑み、本発明の目的は、より高い安全性及び信頼性で乱気流遭遇時の機体動揺を抑えることができる乱気流センシングシステム、航空機及び乱気流センシング方法を提供することにある。
【課題を解決するための手段】
【0008】
本発明者らは、鋭利検討し、動揺低減技術の航空機実装を可能にするためには、ドップラーライダーとは別系統で安全性や信頼性を補強するための、乱気流センシング補強システムを導入した新しいシステム構成を工夫する必要がある、との結論に至った。
【0009】
本発明に係る乱気流センシングシステムは、航空機に搭載されたドップラーライダーと、前記ドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測する乱気流予測部と、前記航空機の主翼表面の複数点に配置された複数の静圧計測部とを具備し、前記複数の静圧計測部により計測された各前記点の計測値に基づいて、前記予測される乱気流の情報を補償する。
【0010】
本発明では、ドップラーライダーとは別系統である航空機の主翼表面の静圧を複数点計測し、この計測結果からドップラーライダー側で予測した乱気流の情報を補償しているので、ドップラーライダーの観測誤差に対するロバスト性を補強できる。これにより、より高い安全性及び信頼性で乱気流遭遇時の機体動揺を抑えることができる。
【0011】
本発明の一形態に係る乱気流センシングシステムは、前記ドップラーライダーによる計測値に基づき、前記予測される航空機前方の乱気流に応じた直線渦糸の回転軸位置と強度を推定する渦推定部と、前記推定された各直線渦糸の回転軸位置と強度に基づき前記航空機の主翼表面の第１の圧力場を推定する第１の圧力場推定部と、前記複数の静圧計測部により計測された各前記点の計測値から、前記主翼表面の第２の圧力場を推定する第２の圧力場推定部とを具備し、前記第１の圧力場と前記第２の圧力場との一致度を求め、当該一致度に基づき、前記予測される乱気流の情報を補償する。
【0012】
本発明の別形態に係る乱気流センシングシステムは、前記ドップラーライダーによる計測値から、Δｔ秒先で機体に誘起する未来の第１の鉛直方向風速度を予測する第１の鉛直面風速場予測部と、前記ドップラーライダーによる計測値及び前記複数の静圧計測部により計測された各前記点の計測値から、前方渦がΔｔ秒先で機体に誘起する未来の第２の鉛直方向風速度を予測する第２の鉛直面風速場予測部と、前記第１の鉛直方向風速度と前記第２の鉛直方向風速度との相関度を求め、当該相関度に基づき、前記予測される乱気流の情報を補償する。このシステムにおいて、既通過領域における相関度を評価し、当該評価結果に応じて、前記予測される乱気流の情報の補償の是非を決めてもよい。
【0013】
本発明の一形態に係る航空機は、ドップラーライダーと、エアデータセンサとを搭載した航空機であって、前記ドップラーライダーによる計測値と、前記エアデータセンサによる計測値と、前記航空機に関する拡大機体モデルのパネルデータとを用い、擾乱ポテンシャルを算出することによって前記航空機の圧力場を推定する圧力推定部を有する。この場合、前記圧力推定部は、前記推定された圧力場から当該航空機に働く空気力変動を得るものであり、前記空気力変動をパラメタの１つとして機体の姿勢安定化フィードバック制御を行う姿勢安定化制御部をさらに有してもよい。
本発明の一形態に係る航空機は、複数の静圧計測部を搭載した航空機であって、前記静圧計測部による計測値と、前記航空機に関する基本圧力場データとを用いて、前記航空機の迎角変化と後縁形状効果とを最小二乗推定により算出し、前記航空機を構成する機体の圧力場を推定する圧力推定部を有する。本発明の一形態に係る航空機は、さらに、ドップラーライダーと、エアデータセンサと、前記圧力推定部が推定する圧力場と、前記ドップラーライダーによる計測値と、前記エアデータセンサによる計測値とを用いて、前記ドップラーライダー視野内の渦の渦強度及び渦誘起迎角を算出し、前記航空機周囲の空間における鉛直面風速場を予測する鉛直面風速場予測部とを有してもよい。
本発明に係る乱気流センシング方法は、航空機に搭載されたドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測し、前記航空機の主翼表面の複数点に配置された静圧を計測し、前記計測した各前記点の静圧の計測値に基づいて、前記予測される乱気流の情報を補償する。
本発明の一形態に係る航空機は、複数の静圧計測部を搭載した航空機であって、前記静圧計測部による計測値と、前記航空機に関する基本圧力場データとを用いて、前記航空機の迎角変化と後縁形状効果とを最小二乗推定により算出し、前記航空機を構成する機体の圧力場を推定する圧力推定部を有する。本発明の一形態に係る航空機は、さらに、ドップラーライダーと、エアデータセンサと、前記圧力推定部が推定する圧力場と、前記ドップラーライダーによる計測値と、前記エアデータセンサによる計測値とを用いて、前記ドップラーライダー視野内の渦の渦強度及び渦誘起迎角を算出し、前記航空機周囲の空間における鉛直面風速場を予測する鉛直面風速場予測部とを有してもよい。
本発明に係る乱気流センシング方法は、航空機に搭載されたドップラーライダーによる計測値を用い、当該航空機前方の乱気流を予測し、前記航空機の主翼表面の複数点に配置された静圧を計測し、前記計測した各前記点の静圧の計測値に基づいて、前記予測される乱気流の情報を補償する。
【発明の効果】
【0014】
本発明によれば、より高い安全性及び信頼性で乱気流遭遇時の機体動揺を抑えることができる。
【図面の簡単な説明】
【0015】
【図1】本発明の第１の実施形態に係るドップラーライダー２を使ったシステムの説明図である。
【図3】渦推定部２２及び第１の圧力場推定部２３の構成及び動作を説明するためのフロー図である。
【図4】渦の回転軸をｚ軸と一致するようにシフトした場合の、動点Ｑとの位置関係を示す説明図である。
【図5】第２の圧力場推定部２４の構成及び動作を説明するためのフロー図である。
【図6】主翼３を概略的に示した斜視図である。
【図7】本発明の第２の実施形態に係る乱気流センシングシステム１１０の構成を示すブロック図である。
【図8】本発明の第３の実施形態に係る乱気流センシングシステムの説明図である。
【発明を実施するための形態】
【0016】
以下、図面を参照しながら、本発明の実施形態を説明する。
【0017】
【0018】
図１に示すように、例えば高度方向の幅３００ｍで水平方向距離１０ｋｍの範囲の晴天乱気流を時速２００ｍ/ｓ程度の航空機１が通過するのに要する時間は５０秒程度である。このような乱気流を回避する飛行をしたいが、回避飛行が間に合わない場合は、乱気流空域通過時の機体動揺を可能な限り抑え、被害を最小化する飛行制御技術も必要である。
本実施形態に係る航空機１は、図２に示す乱気流センシングシステム１０を搭載し、慣性センサ情報のフィードバック制御だけでなく、ドップラーライダーを用いた前方風速情報の先読みにもとづく予見フィードフォワード制御を行い、乱気流遭遇時の機体動揺を顕著に低減し、更により高い安全性及び信頼性で乱気流遭遇時の機体動揺を抑えている。
本実施形態に係る航空機１は、図２に示す乱気流センシングシステム１０を搭載し、慣性センサ情報のフィードバック制御だけでなく、ドップラーライダーを用いた前方風速情報の先読みにもとづく予見フィードフォワード制御を行い、乱気流遭遇時の機体動揺を顕著に低減し、更により高い安全性及び信頼性で乱気流遭遇時の機体動揺を抑えている。
【0019】
図２に示すように、乱気流センシングシステム１０は、機体１１の機首付近に搭載されたドップラーライダー２と、主翼表面に複数分散配置された静圧計測部４と、乱気流予測部２１と、渦推定部２２と、第１の圧力場推定部２３と、第２の圧力場推定部２４と、一致度検定部２５とを有する。
ドップラーライダー２としては、典型的には、晴天乱気流の事前検知用小型ドップラーライダーを用いることができる。乱気流予測部２１は、ドップラーライダー２による計測値を用い、航空機１の前方の乱気流を予測する。このドップラーライダー２を用いて乱気流を予測する技術は例えば特許文献１や非特許文献１、非特許文献２に開示されている。
静圧計測部４は、典型的には、主翼３の表面に設けられた圧力孔と、圧力トランスデューサとを有する。圧力孔を構成する技術は例えば特許文献２に開示されている。しかし、本発明に係る静圧計測部４は圧力孔による構成に限定されない。
ドップラーライダー２としては、典型的には、晴天乱気流の事前検知用小型ドップラーライダーを用いることができる。乱気流予測部２１は、ドップラーライダー２による計測値を用い、航空機１の前方の乱気流を予測する。このドップラーライダー２を用いて乱気流を予測する技術は例えば特許文献１や非特許文献１、非特許文献２に開示されている。
静圧計測部４は、典型的には、主翼３の表面に設けられた圧力孔と、圧力トランスデューサとを有する。圧力孔を構成する技術は例えば特許文献２に開示されている。しかし、本発明に係る静圧計測部４は圧力孔による構成に限定されない。
【0020】
この航空機１では、姿勢安定化制御部３１が慣性センサ３２による情報から機体１１のフィーバック操舵を行うとともに、フィードフォワード部３３が乱気流予測部２１の計測値から機体１１のフィードフォワード操舵を行う。フィードフォワード部３３の後段には、複数の静圧計測部４により計測された各点の計測値に基づいて、予測される乱気流の情報を補償する乱気流情報補償部３４が設けられている。乱気流情報補償部３４については後述する。
【0021】
渦推定部２２は、ドップラーライダー２による計測値に基づき、予測される航空機１の前方の乱気流に応じた単数または複数の直線渦糸の回転軸位置と強度を推定する。
第１の圧力場推定部２３は、推定された各直線渦糸の回転軸位置と強度に基づき航空機１の主翼表面の第１の圧力場を推定する。
本実施形態では、晴天乱気流に内在する最も危険な渦構造として直線渦糸モデルを仮定する。この渦糸モデルは、気象モデルではない。晴天乱気流に内在する主要現象を近似表現し、以下に述べるロジック設計に活用するための工学モデルである。
第１の圧力場推定部２３は、推定された各直線渦糸の回転軸位置と強度に基づき航空機１の主翼表面の第１の圧力場を推定する。
本実施形態では、晴天乱気流に内在する最も危険な渦構造として直線渦糸モデルを仮定する。この渦糸モデルは、気象モデルではない。晴天乱気流に内在する主要現象を近似表現し、以下に述べるロジック設計に活用するための工学モデルである。
【0022】
このような渦推定部２２及び第１の圧力場推定部２３の構成及び動作を、図３を参照しながら説明する。
【0023】
渦推定部２２は、典型的には、ドップラーライダー２が出力する風速分布から、大気において渦形状をなす構造と、渦構造の速度をデータ抽出する（ステップ３０１）。大気において渦形状をなす構造のデータ抽出の一例としては、ＳＮ画像を用い、渦形状をなす構造の抽出には公知の特徴検出を用いることができる。この時、同時に渦の位置、渦の回転軸も併せて抽出する。渦構造の速度をデータ抽出する一例としては、ドップラー画像を用い、抽出した渦構造についての回転速度を抽出する。
【0024】
渦推定部２２は、ステップ３０１の出力に対し、直線渦糸が誘起する風速場をフィッテイングさせるように、直線渦糸の回転軸配置及び渦の強度を決定する（ステップ３０２）。このとき、ステップ３０２の出力を、ステップ３０５で使用するポテンシャル方程式における機体の座標系で扱える形に座標変換、単位変換等により補正する。
【0025】
航空機１は典型的には対気速度、迎角、横滑り角などの計測するエアデータセンサ（図示を省略）を備えており、第１の圧力場推定部２３は、このエアデータセンサの計測値（対気速度、迎角、横滑り角）を取り込む（ステップ３０３）。
【0026】
第１の圧力場推定部２３は、互いに接する複数の四辺形パネル（平板）が閉空間を構成するように予め作られた拡大機体モデルから、全パネルデータを抽出する（ステップ３０４）。ここで、拡大機体モデルとは、機体モデルと、その後方に仮想的に設ける後流モデルとから構成される。後流モデルとは、機体の後縁に対し、その後方に細長く伸びた薄い後流シート（ウエイク）のモデルであり、流れ場を解析するために仮想的に設けたものである。機体モデルと後流モデルとは機体１１の後縁において連続的に接続しており、拡大機体モデルは一体をなしている。各パネルデータは、パネル番号（例えば１～Ｚ）、パネルの空間座標、パネルの法線ベクトルを含み、これら以外にパネルごとの物理パラメタを含んでもよい。
【0027】
第１の圧力場推定部２３は、機体１１本体及び機体１１の後方空間に仮想的に設けたウエイクからなる拡大機体が一様流を受ける際の速度ポテンシャル方程式(連続の式あるいは質量保存則)に、機体とは別に配置された外部渦に関するデータを加味して、機体表面の速度ポテンシャルを算出する。この算出に際し、ステップ３０２から入力する渦の回転軸方向、渦の座標及び渦の強度と、機体の座標、ステップ３０３から入力する上記一様流に相当する対気速度、機体１１の迎角及び横滑り角と、ステップ３０４から入力するパネルデータとを、速度ポテンシャル方程式に代入する。
上記ポテンシャル方程式は、公知のグリーンの公式とパネルごとの離散化により、典型的には以下の式を用い、各パネルについて擾乱ポテンシャルを算出する。
ここで、機体１１に向かってｘ方向から一様流が流れているとし、鉛直方向をｚとする。このとき、
Ｓ：拡大機体における境界面
φ：擾乱ポテンシャル
ｉ：着目パネルを指定する添字（機体表面上、あるいはそれ以外の任意の点、ただし拡大機体内部の空間は除く）
ｊ：機体を構成する全パネルを示す添字（ただし、ｊ≠ｉ）
ｋ：後流を構成する全パネルを示す添字（ただし、ｋ≠ｉ）
ｕ：後流パネルのうち上方（upper）のパネル
ｌ：後流パネルのうち下方（lower）のパネル
∫ｐ：各パネルにおける面積分
ｒ：着目位置ｉと拡大機体表面を覆うパネルの図心との距離
φν：渦誘起ポテンシャル
ｎ：境界面Ｓ上の法線ベクトル
ｎｘ：境界面Ｓ上の法線ベクトルのｘ成分
また、φν（渦誘起ポテンシャル）を含む、上式の右辺２番目の積分記号は、以下のように展開できる。
ここで、機体パネルと渦の回転軸をひとまとめにして原点周りに回転させ、回転軸がｚ軸と平行になるようにする。その平面（便宜上ここではｘｙ平面と記す）において、
ｘＱとｙＱ：機体表面上の動点Ｑのｘｙ座標
ｘ０とｙ０：渦の回転軸がｚ軸と平行になるように回転後の渦のｘｙ座標
Γ：渦強度（渦の循環）
また、ａ、ｂは機体表面上の動点Ｑと渦の回転軸との相対位置により決まる角度であり（図４参照）、渦の始点Ａと終点Ｂの距離が十分大きい場合はａ＝０、ｂ=０であるから、（cosａ＋cosｂ）の項は２としてよい。
上記ポテンシャル方程式は、公知のグリーンの公式とパネルごとの離散化により、典型的には以下の式を用い、各パネルについて擾乱ポテンシャルを算出する。
Ｓ：拡大機体における境界面
φ：擾乱ポテンシャル
ｉ：着目パネルを指定する添字（機体表面上、あるいはそれ以外の任意の点、ただし拡大機体内部の空間は除く）
ｊ：機体を構成する全パネルを示す添字（ただし、ｊ≠ｉ）
ｋ：後流を構成する全パネルを示す添字（ただし、ｋ≠ｉ）
ｕ：後流パネルのうち上方（upper）のパネル
ｌ：後流パネルのうち下方（lower）のパネル
∫ｐ：各パネルにおける面積分
ｒ：着目位置ｉと拡大機体表面を覆うパネルの図心との距離
φν：渦誘起ポテンシャル
ｎ：境界面Ｓ上の法線ベクトル
ｎｘ：境界面Ｓ上の法線ベクトルのｘ成分
また、φν（渦誘起ポテンシャル）を含む、上式の右辺２番目の積分記号は、以下のように展開できる。
ここで、機体パネルと渦の回転軸をひとまとめにして原点周りに回転させ、回転軸がｚ軸と平行になるようにする。その平面（便宜上ここではｘｙ平面と記す）において、
ｘＱとｙＱ：機体表面上の動点Ｑのｘｙ座標
ｘ０とｙ０：渦の回転軸がｚ軸と平行になるように回転後の渦のｘｙ座標
Γ：渦強度（渦の循環）
また、ａ、ｂは機体表面上の動点Ｑと渦の回転軸との相対位置により決まる角度であり（図４参照）、渦の始点Ａと終点Ｂの距離が十分大きい場合はａ＝０、ｂ=０であるから、（cosａ＋cosｂ）の項は２としてよい。
【0028】
第１の圧力場推定部２３は、各パネルにおける擾乱ポテンシャルと一様流の速度ポテンシャルを加えて速度場を算出する（ステップ３０６）。この速度場の算出法は、非特許文献３に記載されている。
【0029】
第１の圧力場推定部２３は、速度場から、ベルヌーイの定理より圧力場を算出する（ステップ３０７）。
【0030】
以上のステップ３０１～ステップ３０７によって、機体１１にかかる圧力場を算出することができる。例えば非特許文献３の第７図に示されるように、翼にかかる圧力分布を算出することができる。
【0031】
一方、第２の圧力場推定部２４は、複数の静圧計測部４により計測された各点の計測値から、最小二乗推定により主翼３の表面全体の第２の圧力場を推定する。
【0032】
本実施形態に係る第２の圧力場推定部２４は、過渡応答が穏やかな準非定常圧力場の推定である簡易圧力場推定を行うもので、図５を参照しながら説明する。
【0033】
図６は一方の主翼３を概略的に示した斜視図であり、主翼３の表面に点在する符号４１は静圧計測部４の圧力孔を示している。Δαは迎角変動、Ｓｊは後縁形状変化に対応する翼素の面積、δｊは後縁形状変化を示す。
【0034】
圧力場をＮ次元状態変数Ｘ（Ｎ～１０２）で近似する。Ｘは表面圧力場（無限次元）を離散化した多次元状態ベクトル（有限次元）である。
【0035】
第２の圧力場推定部２４は、静圧計測部４の多点分散配置に応じた観測行列Ｈを選択する（ステップ５０１）。
観測行列Ｈは状態ベクトルＸの中で、圧力孔４１により実際に直接観測できる変数を指定する行列であり、通常０あるいは１である。圧力孔４１を主翼３の表面のどこに設置するかは自由度があるので、推定系を構成する際の設計パラメタ（最適化変数）とみなすことができる。
観測行列Ｈは状態ベクトルＸの中で、圧力孔４１により実際に直接観測できる変数を指定する行列であり、通常０あるいは１である。圧力孔４１を主翼３の表面のどこに設置するかは自由度があるので、推定系を構成する際の設計パラメタ（最適化変数）とみなすことができる。
【0036】
第２の圧力場推定部２４は、ＣＦＤ（Computational Fluid Dynamics：数値流体解析）や風洞試験により得られる基本圧力場をデータベース化したデータを取り込む。
第２の圧力場推定部２４では、基本圧力場からの観測誤差Ｙεを、迎角変化Δαと後縁形状効果δの重ね合わせ（係数行列Ｂα、Ｂδとの掛け算）として以下に示す線形確率過程モデルで表現する。
第２の圧力場推定部２４では、基本圧力場からの観測誤差Ｙεを、迎角変化Δαと後縁形状効果δの重ね合わせ（係数行列Ｂα、Ｂδとの掛け算）として以下に示す線形確率過程モデルで表現する。
【0037】
ここで、Δαは迎角変動分の真値であり、表面圧力場全体を変化させるためのパラメタである。δは翼後縁形状変化分の真値であり、表面圧力場をスパン方向に変化させるパラメタである。図６では８舵面に分割されているので８次元ベクトルであるが、もっと粗くても細かくてもよい。Ｙεは飛行中の迎角計測値αとデータベース参照により得られる圧力場のノミナル値と、実観測圧力の誤差である。Ｂαは圧力場の迎角変動に対する感度行列であり、同定作業が必要である。Ｂδは圧力場の後縁形状効果δに対する感度行列であり、同定作業が必要である。Ｗは迎角変動と後縁形状の変化では説明できない圧力場の変化であり、正規分布に従う確率変数と考える。
【0038】
第２の圧力場推定部２４は、Ｙεを説明できる迎角変化Δαと後縁形状効果δを最小二乗推定により決定する（ステップ５０２）。
第２の圧力場推定部２４は、下記の式により推定圧力場を決定する（ステップ５０３）。
ここで、α（オーバライン）は飛行中の迎角基準値であり、飛行機に標準搭載される迎角計による測定値又は何らかの推定値である。Ｘ（α）（いずれもオーバライン）は迎角基準値にもとづくノミナル圧力場であり、あらかじめテーブル化されたデータベースである。Δα（オーバハットライン）とδ（オーバーハットライン）は最小二乗推定値である。
また、第２の圧力場推定部２４は、下記の式により推定誤差共分散Ｐを計算する（ステップ５０４）。
ここで、E（・）は期待値である。Δα（オーバウェーブライン）は真値との誤差である。δ（オーバウェーブライン）は真値との誤差ベクトルである。
つまり、第２の圧力場推定部２４は、この行列P
の対角成分により圧力場推定の精度を評価する。
第２の圧力場推定部２４は、下記の式により推定圧力場を決定する（ステップ５０３）。
また、第２の圧力場推定部２４は、下記の式により推定誤差共分散Ｐを計算する（ステップ５０４）。
ここで、E（・）は期待値である。Δα（オーバウェーブライン）は真値との誤差である。δ（オーバウェーブライン）は真値との誤差ベクトルである。
つまり、第２の圧力場推定部２４は、この行列P
【0039】
なお、推定圧力場は揚力へ換算できる。
また、カルマンフィルタ型への拡張により非定常圧力場推定も可能である。この推定は、典型的には、カルマンフィルタ型逐次推定アルゴリズムにより行うことができる。
まず、圧力場の非定常性が無視できない場合のカルマンフィルタ型動的逐次アルゴリズムは以下のように定式化できる。 非定常圧力変動分を以下のように離散時間状態方程式として表現する。
ΔXは非定常成分(過渡応答分)を示す変数、A，Γ1，Γ2は過渡応答特性を表現するパラメタであり、風洞試験入出力データ等から同定する。V，Wはモデル化誤差を表現し、適当な正規分布に従うノイズを仮定する。ｋは時系列を示す添字である。
拡大ベクトル、拡大行列をまとめ、上記の数式を以下のように書き直す。
上式に対して以下のように線形カルマンフィルタを構成して、
これにより、
を求める。
この結果を用いて第２の圧力場の推定値を以下のように更新する。
ここで、
は基本空力特性を格納したデータベースである。
また、カルマンフィルタ型への拡張により非定常圧力場推定も可能である。この推定は、典型的には、カルマンフィルタ型逐次推定アルゴリズムにより行うことができる。
まず、圧力場の非定常性が無視できない場合のカルマンフィルタ型動的逐次アルゴリズムは以下のように定式化できる。 非定常圧力変動分を以下のように離散時間状態方程式として表現する。
拡大ベクトル、拡大行列をまとめ、上記の数式を以下のように書き直す。
この結果を用いて第２の圧力場の推定値を以下のように更新する。
【0040】
一致度検定部２５は、上記のように推定された第１の圧力場と第２の圧力場との一致度（γ：例えば０≦γ≦１）を求める。一致度検定部２５は、例えば部分静圧情報を、翼断面圧力係数分布、さらにスパン方向も考慮した主翼表面圧力分布という物理的特徴量に変換し、その一致性を評価すればよい。
【0041】
乱気流情報補償部３４は、一致度検定部２５により求められた一致度γに応じて、乱気流予測部２１による予測値を補償する。例えば、乱気流情報補償部３４は、第１の圧力場と第２の圧力場との一致度γが高くなると、乱気流予測部２１による予測値による機体１１のフィードフォワード操舵量を多くし、一致度γが低くなると、フィードフォワード操舵量を少なくする。このように本実施形態に係る乱気流センシングシステムでは、ドップラーライダー２とは別系統である航空機１の主翼３の表面の静圧を複数点計測し、この計測結果からドップラーライダー２側で予測した乱気流の情報を補償しているので、ドップラーライダー２の観測誤差に対するロバスト性を補強できる。これにより、より高い安全性及び信頼性で乱気流遭遇時の機体動揺を抑えることができる。
なお、上記実施形態では、乱気流情報補償部３４が一致度γに応じてフィードフォワード制御を行っていたが、ドップラーライダー２、渦推定部２２、第１の圧力場推定部２３の出力を利用し、乱気流情報補償部３４（フィードフォワード部）の出力をゼロとして、姿勢安定化制御部３１のゲインを適応的に調整してもよい（図示せず）。その場合、第１の圧力場推定部２３の出力より機体１１に働く空気力変動に関する情報が得られるので、この機体１１に働く空気力変動に関する情報を姿勢安定化制御部３１に入力し、姿勢安定化制御部３１はその情報もパラメタとして再設計してフィードバック制御を行う。再設計は、例えば飛行機搭載状態で、実時間で行う。また、実時間動作が不安視される場合は、ゲインスケジューリング化してもよい。この形態では、フィードバックゲインのみの制御系となるが、制御系内部のゲインが乱気流遭遇時用に自動調整されることとなる。以上の制御は、乱気流情報補償部３４（フィードフォワード部）の出力をゼロに切り替えるように構成しているが、乱気流センシングシステム１０を、ドップラーライダー２と、乱気流予測部２１と、渦推定部２２と、第１の圧力場推定部２３とから構成し、第１の圧力場推定部２３の出力である機体１１に働く空気力変動に関する情報を姿勢安定化制御部３１に入力し、フィードフォワード制御をせずにフィードバック制御を行うようにしてもよい。また、乱気流センシングシステム１０を、ドップラーライダー２と、渦推定部２２と、第１の圧力場推定部２３とから構成し、第１の圧力場推定部２３の出力である機体１１に働く空気力変動に関する情報を姿勢安定化制御部３１に入力し、フィードフォワード制御をせずにフィードバック制御を行うようにしてもよい（図示せず）。
なお、上記実施形態では、乱気流情報補償部３４が一致度γに応じてフィードフォワード制御を行っていたが、ドップラーライダー２、渦推定部２２、第１の圧力場推定部２３の出力を利用し、乱気流情報補償部３４（フィードフォワード部）の出力をゼロとして、姿勢安定化制御部３１のゲインを適応的に調整してもよい（図示せず）。その場合、第１の圧力場推定部２３の出力より機体１１に働く空気力変動に関する情報が得られるので、この機体１１に働く空気力変動に関する情報を姿勢安定化制御部３１に入力し、姿勢安定化制御部３１はその情報もパラメタとして再設計してフィードバック制御を行う。再設計は、例えば飛行機搭載状態で、実時間で行う。また、実時間動作が不安視される場合は、ゲインスケジューリング化してもよい。この形態では、フィードバックゲインのみの制御系となるが、制御系内部のゲインが乱気流遭遇時用に自動調整されることとなる。以上の制御は、乱気流情報補償部３４（フィードフォワード部）の出力をゼロに切り替えるように構成しているが、乱気流センシングシステム１０を、ドップラーライダー２と、乱気流予測部２１と、渦推定部２２と、第１の圧力場推定部２３とから構成し、第１の圧力場推定部２３の出力である機体１１に働く空気力変動に関する情報を姿勢安定化制御部３１に入力し、フィードフォワード制御をせずにフィードバック制御を行うようにしてもよい。また、乱気流センシングシステム１０を、ドップラーライダー２と、渦推定部２２と、第１の圧力場推定部２３とから構成し、第１の圧力場推定部２３の出力である機体１１に働く空気力変動に関する情報を姿勢安定化制御部３１に入力し、フィードフォワード制御をせずにフィードバック制御を行うようにしてもよい（図示せず）。
【0042】
＜第２の実施形態＞
図７は、本発明の第２の実施形態に係る乱気流センシングシステムの構成を示すブロック図である。なお、第１の実施形態と同一の要素には同一の符号を付す。
図７は、本発明の第２の実施形態に係る乱気流センシングシステムの構成を示すブロック図である。なお、第１の実施形態と同一の要素には同一の符号を付す。
【0043】
図７に示すように、乱気流センシングシステム１１０は、２軸タイプ又はスキャンタイプのドップラーライダー２と、乱気流予測部２１と、第１の鉛直面風速場予測部１１１と、静圧計測部４と、渦強度及び渦誘起迎角推定部１１２と、第２の鉛直面風速場予測部１１３と、相関部１１４とを有する。
２軸タイプ又はスキャンタイプのドップラーライダー２の構成は、例えば特許文献１に記載されている。
２軸タイプ又はスキャンタイプのドップラーライダー２の構成は、例えば特許文献１に記載されている。
【0044】
ドップラーライダー２は、晴天乱気流コア位置の自機に対する相対位置Ｒｔ（時刻ｔでの距離）を算出する。 第１の鉛直面風速場予測部１１１は、乱気流予測部２１による航空機１の前方の乱気流を予測値から、Δｔ秒先で機体に誘起する未来の鉛直方向風速度を予測する。
【0045】
渦糸までの距離Ｒｔはドップラーライダー２より得られるので、直線渦糸の未知パラメタは回転軸方向ψ（ベクトル）と強度Γ（スカラー）である。ただし、以下の理由により、直線渦糸の回転軸方向を、航空機１の進行方向と直交する水平面内に限定する。第１の理由として、晴天乱気流の正体は高度方向のシアによるヘルムホルツ不安定であり、竜巻のように回転軸が有意な鉛直成分をもたない。第２の理由として、航空機１のスパン長（主翼の左端から右端までの長さ）は、距離Ｒｔや渦の空間スケールに比較して圧倒的に小さいので、ψの影響は無視できる。この仮定により、未知パラメタは強度Γ（短時間では一定値）のみとなる。
【0046】
渦強度及び渦誘起迎角推定部１１２は、静圧計測部４による主翼３の表面の多点圧力計測結果に基づき、最小二乗推定により、第１の実施形態と同様のアルゴリズムにより主翼表面圧力場を推定し、その圧力場を積分することにより揚力Ｌを計算し、その時間変化ΔＬｔ＝Ｌｔ－Ｌｔ－Δｔと、揚力傾斜ＣＬα（機体形状に応じて決まる揚力曲線の傾き）より、迎角変動換算値Δαｔ＝ΔＬｔ／ＣＬαを計算する。飛行機の標準的なエアデータセンサである迎角計を用いて渦誘起分の迎角変動Δαｔを算出することも考えられるが、従来迎角計の計測点が擾乱大気中の１点のみで、定常風はともかく、渦誘起迎角分を十分な精度で抽出することは困難である。これに対して、主翼表面流れの多点圧力観測から主翼発生揚力変化を算出し、渦強度及び渦誘起迎角相当値の算出に結び付けることで、機体上下動揺に直結する渦の影響をよりロバストに峻別できる。なお、ここで第１の実施形態と同様のアルゴリズムとは、定常圧力場推定のアルゴリズム（ステップ５０１～５０４）と、上述した非定常圧力場推定のアルゴリズム（典型的には、カルマンフィルタ型逐次推定アルゴリズム）とが含まれる。
【0047】
渦強度及び渦誘起迎角推定部１１２は、算出した渦誘起効果の迎角換算値Δαｔを用いて、次式より、渦糸モデルの強度を計算する。
第２の鉛直面風速場予測部１１３は、前方渦がΔｔ秒先で機体に誘起する未来の鉛直方向風速度を以下のように予測する。
【0048】
相関部１１４は、第１の鉛直面風速場予測部１１１で予測されたΔｔ秒後の未来の風速時系列予測値と、第２の鉛直面風速場予測部１１３で予測されたΔｔ秒後の未来の風速時系列予測値との相関係数（０≦γ≦１）を求める。
【0049】
乱気流情報補償部３４は、求められた相関係数に応じてフィードフォワード操舵量を補正する。例えば、相関γが弱いと判断された場合はドップラーライダーの信号処理のみによる鉛直方向風速の誤差が大きい可能性ありと判断しゲインを絞る。相関γが強い場合はドップラーライダー計測の信頼性は高いと判断し、補正なしでフィードフォワード操舵量をそのまま加える。
【0050】
この航空機１では、慣性センサ３２による情報のフォードバックループにより機体１１の姿勢安定制御部３６がフィーバック操舵を行うとともに、空力角センサ３５で計測されたα、β（迎角、横滑り角）は渦強度及び渦誘起迎角推定部１１２に入力される。
本実施形態に係る乱気流センシングシステム１１０では、Δｔ秒後の未来のフィードフォワード制御をより高い安全性及び信頼性で行うことができる。加えて、第２の実施形態では、第１の圧力場推定部におけるようなパネル法計算を用いないので、計算負荷が小さく、オンボード演算に向いている。
本実施形態に係る乱気流センシングシステム１１０では、Δｔ秒後の未来のフィードフォワード制御をより高い安全性及び信頼性で行うことができる。加えて、第２の実施形態では、第１の圧力場推定部におけるようなパネル法計算を用いないので、計算負荷が小さく、オンボード演算に向いている。
【0051】
＜第３の実施形態＞
図８は本発明の第３の実施形態に係る乱気流センシングシステムの説明図である。
図８に示す航空機１は、例えば図７に示した乱気流センシングシステム１１０を搭載する。現在と過去のセンシング情報を用いて、ドップラーライダー２と後段の乱気流予測部２１の信頼性を評価しながら、予見情報を活用するか否かを判定する。既通過領域（図８のＡ）である過去の評価に基づき、遭遇予定（図８のＢ）である未来の予見制御の適用の是非を判定する。その判定指標として、既通過領域（図８のＡ）での静圧計測部４～第２の鉛直面風速場予測部１１３による主翼表面圧力場センシング情報を活用する。例えば、既通過領域（図８のＡ）において、相関部１１４における相関係数の累積値が所定以上の値の場合には、遭遇予定（図８のＢ）である未来の予見情報を用いたフィードフォワード制御を適用し、そうでない場合には予見制御を適用しない。
これにより、予見制御のリスクを低減する最も注意深い方法で、乱気流遭遇時の機体動揺を抑えることができる。
図８は本発明の第３の実施形態に係る乱気流センシングシステムの説明図である。
図８に示す航空機１は、例えば図７に示した乱気流センシングシステム１１０を搭載する。現在と過去のセンシング情報を用いて、ドップラーライダー２と後段の乱気流予測部２１の信頼性を評価しながら、予見情報を活用するか否かを判定する。既通過領域（図８のＡ）である過去の評価に基づき、遭遇予定（図８のＢ）である未来の予見制御の適用の是非を判定する。その判定指標として、既通過領域（図８のＡ）での静圧計測部４～第２の鉛直面風速場予測部１１３による主翼表面圧力場センシング情報を活用する。例えば、既通過領域（図８のＡ）において、相関部１１４における相関係数の累積値が所定以上の値の場合には、遭遇予定（図８のＢ）である未来の予見情報を用いたフィードフォワード制御を適用し、そうでない場合には予見制御を適用しない。
これにより、予見制御のリスクを低減する最も注意深い方法で、乱気流遭遇時の機体動揺を抑えることができる。
【0052】
＜その他＞
本発明は、上記の実施形態には限定されず、本発明の技術思想の範囲内で変形や応用が可能であり、その変形や応用による実施の範囲も本発明の技術的範囲に属する。
本発明は、上記の実施形態には限定されず、本発明の技術思想の範囲内で変形や応用が可能であり、その変形や応用による実施の範囲も本発明の技術的範囲に属する。
【0053】
例えば、上記の実施形態では、乱気流の検出にドップラーライダー２を用いたが、これに限るものではない。ドップラーライダー２と同様に、大気の流れを動的に測定し、渦の形状、渦の位置、渦の強度（風速）を検出できるものであればよく、航空機に搭載可能なドップラーレーダであってもよい。
【符号の説明】
【0054】
１ ：航空機
２ ：ドップラーライダー
３ ：主翼
４ ：静圧計測部
１０ ：乱気流センシングシステム
１１ ：機体
２１ ：乱気流予測部
２２ ：渦推定部
２３ ：第１の圧力場推定部
２４ ：第２の圧力場推定部
２５ ：一致度検定部
３１ ：姿勢安定化制御部
３２ ：慣性センサ
３３ ：フィードフォワード部
３４ ：乱気流情報補償部
３５ ：空力角センサ
４１ ：圧力孔
１１０ ：乱気流センシングシステム
１１１ ：第１の鉛直面風速場予測部
１１２ ：渦強度及び渦誘起迎角推定部
１１３ ：第２の鉛直面風速場予測部
１１４ ：相関部
２ ：ドップラーライダー
３ ：主翼
４ ：静圧計測部
１０ ：乱気流センシングシステム
１１ ：機体
２１ ：乱気流予測部
２２ ：渦推定部
２３ ：第１の圧力場推定部
２４ ：第２の圧力場推定部
２５ ：一致度検定部
３１ ：姿勢安定化制御部
３２ ：慣性センサ
３３ ：フィードフォワード部
３４ ：乱気流情報補償部
３５ ：空力角センサ
４１ ：圧力孔
１１０ ：乱気流センシングシステム
１１１ ：第１の鉛直面風速場予測部
１１２ ：渦強度及び渦誘起迎角推定部
１１３ ：第２の鉛直面風速場予測部
１１４ ：相関部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】