河川流計測における画像解析技術の実用化

河川流計測における画像解析技術の実用化 1. 市民工学専攻：本田 将人 指導教員：藤田 一郎 はじめに 流量観測において，接触型流速計は経済性と安全性 図-1 斜め画像（左），無歪（幾何補正）画像 （右） に不安が残る．そこで，安全で安価な方法として画像 による非接触型の流速計即方法 LSPIV（Large-Scale Particle Image Velocimetry）や STIV（Space-Time Image Velocimetory）が開発されてきた 1)．これらの手法は河 岸や橋梁に設置したビデオカメラで撮影された映像 をパーソナルコンピュータで解析する手法である．こ れらの方法は近年問題になっているゲリラ豪雨など の災害にも役立つ．また，流量解析が困難な山間部や 図-2 簡易型キャリブレーションの実験 中小河川でも流量観測が可能になる．本研究では，実 河川の画像流速解析の実用化を目的とし，これまでの 画像解析技術の応用と利便性を重視したソフトウェ アの開発と精度について述べる． 2. カメラキャリブレーション 画像記録によって，流速情報を抽出するためには， 撮影座標（ビデオ画像上の座標）から実空間座標へと 図-3 簡易型キャリブレーションによる結果（実 験） 斜め画像（左），幾何補正画像（右） 座標変換する必要がある．そのため，カメラのパラメ ータをビデオカメラ画像上に写っている標定点から 逆推定する，カメラキャリブレーションが必要となる． このキャリブレーションによって PC 内部では，ビデ オカメラからの斜め画像は無歪画像に変換される（図 -1）．ところで，通常のカメラキャリブレーションでは 画像上に物理座標を既知とする数点の標定点が必要 である．しかし，山間部に存在する河川では，任意の 図-4 簡易型キャリブレーションによる結果 （CG） 標定点を設置することが困難である．そのため，あら 斜め画像（左），幾何補正画像（右） における，ビデオカメラからの斜め画像と簡易型キャ かじめカメラのレンズ位置，俯角を既知とし，一点の リブレーションによって作成された幾何補正画像を 物理座標を与えることによって簡易的にキャリブレ 示している．斜め画像では歪んで映っているチェスボ ーションを成功させる技術が必要となる．実験によっ ードが，幾何補正画像では正方形に補正されている． てカメラの位置，俯角，床に設置したチェスボードの この結果より，キャリブレーションが簡易的に行われ 一点の座標によって，無歪画像を本開発ソフトで作成 た事がわかる．次に仮想空間を想定し，仮想空間にお することが可能であることを確認した．カメラ，三脚， けるカメラ座標と標定点から簡易型キャリブレーシ チェスボード，角度計を用意し，図-2 のように設置す ョンが可能かを実験した．仮想空間の描画にはオープ る．この時，カメラのレンズの真下を原点とする任意 ンソースの OpenGL を用いた．図 4 に OpenGL によっ の物理座標軸を考える. 図-3 には実際の実験 て描画された斜め画像と幾何補正画像を示す． 画像解析技術，ソフトウェア開発，カメラキャリブレーション，浮子追跡 PTV，STIV．

元フーリエ変換と Gaussian フィルターを用いる． 図-7 STIV（左；斜め画像，右；幾何補正画像） 図-5 Float-PTV 図-8 STI（ヒストグラム均等化） ① 図-6 Float-PTV の多重合成 ⑩ 標定点の緑の板と格子の床が正方形に幾何補正され 図-9 STIV 解析結果 ていることがわかる． （緑；自動解析結果，赤；マニュアル解析） 3. Float-PTV さらに，ヒストグラムの均等化処理を施すことによっ 土木研究所（ICHARM）との共同研究において， て，精度を向上させる工夫を行っている（図-8）．鬼怒 Float-PTV という河川画像上に写る浮子を追跡するこ 川を対象に新しい STIV による解析結果を図-9 に示す． とによって容易に流速を計測できるソフトウェアを 左岸（画面奥）に近づくほど流速が低下している．ま 開発した（図-5）．また，河川管理において専門家でな た，岸から遠い画面手間に近い検査線上①～⑤では流 くとも容易に操作できる UI を兼ね備えている．その 速値は 2.2m/s～2.6m/s となっている．自動解析による 機能の一つが多重合成画像の作成である（図-6）．画像 結果はマニュアル解析によるものと一致しており，精 上に映し出された浮子の集まりをマウスクリックす 度に問題はないと考えられる． るだけで，表面流速値が得られる．しかしながら，浮 子は河川表面下に潜ってしまい視認できないことが ５．結論 ある．このことは，水面波紋を用いて流速を測定する 技術の必要性を再認識することとなった． 今回の研究では，画像流速解析技術の実用化を目的 としていくつかのソフトウェアを作成してきた．カメ ラキャリブレーションでは，カメラの俯角と対象まで 4. STIV の距離を測定することで，簡易的に幾何補正画像を作 STIV は，検査線上の輝度を時間方向に並べた時空間 成することが可能となった．また，容易に表面流速を 画像における勾配によって流速を測定する手法であ 測定できるソフトウェア Float-PTV，新しい試みを取 る（図-7）．本研究では，先のキャリブレーションによ り入れることで精度を向上させた STIV を開発した． って得られたカメラパラメータを用いる，新たな STIV を作成した．STIV の自動解析では輝度勾配テン 6．参考文献 ソル法と統計処理によって，勾配線の平均を算出する． 1) さらに，時空間画像のエッジ検出とノイズ除去に 2 次 藤田一郎：非接触型流速計測法を用いた実河川流の計 測と問題点，ながれ，第 26 巻，pp.5-12, 2007.

河川流計測における画像解析技術の実用化 平 成 25 年 2 月 1 日 神戸大学大学院工学研究科市民工学専攻 氏名 学籍番号 本田 将人 102T132T

Development and Practical Applications of Image-Based Techniques for River Velocity Measurements Basin Masato HONDA 102T132T Department of Architecture and Civil Engineering Graduate School of Science and Technology Kobe University February 01, 2013. ABSTRACT Today, the river flow measurements are usually performed by point measure methods (ex. propeller meter) at low-water stages, and a float method in high-water stages. However, these methods are time-consuming in measurements, economically disadvantageous, and hazardous during the measurement at the huge flood. On the other hand, non-contact techniques have been developed for a stream flow observation with lower risk and lower cost. LSPIV (Large Scale Particle Image Velocimetry) and STIV (Space-Time Image Velocimetry) are representative non-contact techniques for this purpose. These methods simply use a personal computer to analyze images taken by a camcorder. Therefore, these methods can be used for measuring flow discharge when the flow was monitored continuously. In this study, I developed several Image based velocimetry softwares for practical purposes. First, I developed a camera calibration software, which enable river managers to calculate camera parameter for a oblige angle image with ease. Secondly, I developed a software ‘Float-PTV’, that tracks floats in the image screen. Although river managers can easily obtain surface velocity data with this software, floating object can not always be seen when heavy floods occurs. For solving this problem, we paid attention to surface patterns. STIV uses orientation angles on a STI image containing time series of intensity on the search line drawn on the surface of the river. For this reason, I developed a new STIV software that can be easily used by a video manager . image based velocimetry，sofware development，camera calibration，tracking PTV，STIV．

要旨 現在行われている実河川の流量観測では，低水時にはプロペラ流速計や電磁流速計による 計測，高水時には浮子による計測が行われている．しかし，接触型流速計を用いた計測の場 合，流量観測を行うときに多大な時間がかかる点や，経済的に不利である点，安全性に不安 が残る点などが問題点として挙げられる．そこで，危険度が低く，より安価な流量観測の方 法 と し て 非 接 触 型 の 流 速 計 が 開 発 さ れ て き た ． そ の 非 接 触 型 の 流 速 計 と し て LSPIV （Large-Scale Particle Image Velocimetry）や STIV（Space-Time Image Velocimetory）が存在す る．これらの手法は河岸や橋梁に設置したビデオカメラで撮影された映像をパーソナルコン ピュータで解析する手法である．これらの方法はビデオカメラで撮影した映像を用いて解析 を行う手法であるため，ビデオカメラの設置をしておけば近年問題になっているゲリラ豪雨 などの災害にも役立てることができる．また安価なため，流量解析が困難な山間部や中小河 川でも流量観測が可能になる． 今回の研究では，画像流速解析の実用化と改良を視野に入れたソフトウェア開発を行った． まず，河川画像流速解析で最重要となるカメラキャリブレーションの実用化を目的としたソ フト開発を行った．さらに，土木研究所(ICHARM)との共同研究中に，画像浮子追跡ソフト Float-PTV を開発した．このソフトは浮子だけでなく河川漂流物を追跡することで流速を簡単 に測定できる．しかしながら，洪水時における漂流物は必ずしも，映像に移るわけではなく， また，流線上を流れるわけでもない．そこで，河川表面を流れる，波紋パターンに注目して 河川流速を測定する必要がある．STIV では，河川表面に引いた検査線上の輝度を時間方向に 並べる事で作り出される，時空間画像（STI）に描かれる勾配から流速を測定する．このこと から，STIV を実際の河川管理で用いられる事を主な目的として新たに構築しなおした． 画像解析技術，ソフトウェア開発，カメラキャリブレーション，浮子追跡 PTV，STIV．

目次 1. 序論 .................................................................... 1 1.1 我が国の気候と治水対策 ............................................... 1 1.2 我が国における治水事業 ............................................... 2 1.3 我が国における流量観測の実態と問題点 ................................. 2 1.4 近年における高精度な非接触の流量測定方法 ............................. 3 1.5 近年の洪水災害の増加と河川情報システムの現状 ......................... 4 1.6 本研究による治水への貢献 ............................................. 7 2. 画像流速解析ソフトの開発 ................................................ 8 2.1 概要 ................................................................. 8 2.2 構成 ................................................................. 8 2.2.1 画像取得機器 ...................................................... 8 2.2.2 開発環境 ......................................................... 10 3. カメラキャリブレーション ............................................... 12 3.1 画像の座標変換 ...................................................... 12 3.2 アルゴリズム ........................................................ 15 3.3 精度についての考察 .................................................. 30 3.4 最小条件からのキャリブレーション .................................... 35 3.5 まとめ .............................................................. 42 4. Float-PTV .............................................................. 43 4.1 計算の仕組み ........................................................ 43 4.2 画像の多重合成 ...................................................... 45 4.3 さまざまな状況での多重合成と計算結果 ................................ 46 4.4 まとめ .............................................................. 50 5. STIV(Space-Time Image Velocimetry) ..................................... 51 5.1 理論 ................................................................ 51 5.2 STIV 解析 ............................................................ 66

5.3 まとめ .............................................................. 71 6 各々の開発ソフトの操作方法 .............................................. 73 6.1 キャリブレーションソフトの操作 ...................................... 73 6.2 Float-PTV の使い方 ................................................... 78 6.3 STIV の操作方法 ...................................................... 84 3. 結論 ................................................................... 90 謝辞 参考文献 ............................................................. 92

第1章 序論 1.1 我が国の気候と治水対策の必要性 四季の変化に富む気候を持つ日本は，その気象学的，地学的条件から自然災害の発生 しやすい状況に置かれた国である．さらに，財産および人口の都市への集中といった社 会学的な条件も加わって，数多くの甚大な自然災害を経験してきた国でもある．そのた め，現在までに進められてきた治水事業にもかかわらず，このような自然条件は全国各 地での水害を誘発し，日本社会に甚大な損失を与えてきた．以上の事から，未だに我が 国の治水対策については問題点や取り組むべき課題が多く残されている． 図 1.1 全国の水害・土砂災害の発生状況（1994 年～2003 年）1) 図 1.2 洪水氾濫域に集中する資産と人口 1) 1

1994 年から 2003 年の 10 年間に日本の市区町村の 97％以上が水害・土砂災害に襲われ ている（図 1.1）1)．洪水氾濫域（洪水時の河川水位より低い地域）は日本全土の 10%に しかすぎないが，その地域に人口の 51％，総資産の 75％が集中するため，いったん洪 水氾濫が発生すれば，その被害は深刻なものになり，実際に過去 10 年間に 10 回以上被 災した市区町村は 50%以上，1 回以上の被災市区町村は 90%以上に達している（図 1.2） 1) ． 1.2 我が国における治水事業 我が国における河川基本方針の設定や河川整備計画は利水,環境に加えて洪水防御を 主要な目的とする.また、河道計画の策定は、河道内の流出流量の安全な流下を目的とし ている.そのため、降雨と流量の関係性に関するデータが必要となる。これらのデータは、 過去の降雨実績や降雨流出解析をもとに算出される.そして、解析モデルやパラメータの 妥当性の検証と精度向上のための唯一の方法は、現実に流下している流量の観測と対応 する降雨量データの関連を解析することである.さらに，行政による計画のみでは災害か らの防御は完遂できない．そのため，災害時には現地での消防団や自衛隊による救助活 動が重要となる．この救助活動においても，災害の発生している箇所の情報が必要とな る．かつて，古代中国の戦略家である孫子はその著作のなかで「敵を知り，己を知れば 百戦危うからず」と述べている 2)．われわれ土木技術者にとって，敵とは災害，引いて は自然そのものであり，己とは我々の生活している日本の国土そのもの，戦とは自然災 害に対する挑戦，つまり，防災・減災といえる．したがって,河川の流量を把握すること は,防災や水利用,河川環境保全という様々な側面において非常に重要である. 1.3 我が国における治水事業 2013 年現在では、一級河川における定期的な流量計測が、低水から高水に至るまで,業 務化されている.古くから，低水時にはプロペラ流速計による低水観測（２点法）,高水 時には浮子による観測が行われてきた（図 1.3)3）．しかし，国土交通省の管轄外の河道 区間や二級河川以下の河道区間では，流量の観測は定期的なものではなく，情報量も極 端に少ない．こういった事情は，突発的な集中豪雨によって発生した中小河川での洪水 災害，例えば，2008 年 7 月 28 日の都賀川水難事故や 2010 年 8 月 10 日の佐用川豪雨災 害など，における流量データの把握を困難なものとしている．また，上記の浮子観測法 は，①機敏性の欠如（高確率のピーク流量の欠測），②連続観測の不可，③浮子の追随 2

性の問題といった問題から，精度を期待できない上，観測技能者の不足と言った問題を 持つ． 図 1.3. 浮子（左），浮子観測法（右）3) 1.4 近年における高精度な非接触の流量測定方法 ところで，近年では様々な高精度で高効率の流量計測方法が開発されている 4)．非接触 型の計測方法には大きく分けて２つある．一つは，超音波やマイクロ波（電波）を洪水 中の河川に発射し，ドップラー変調を受けた反射波の周波数を解析することで流速を得 る，ドップラー型の流速解析手法である（図 1.4）5)～8)． 図 1.4 ドップラータイプの測定原理８) そして，もうひとつは，画像処理によって洪水流の波紋や流下物の時間変化を解析する， 画像流速解析手法（image based velocimetry）とよばれる流速手法である．画像を用いる メリットとしては，①解析技術者の安全性，②任意の過去時点における流速解析，③解 析に用いる機材が安価である，といったものがある．本論文では，この画像流速解析手 3

法を中心とする． 河川分野における画像流速解析手法には２つある．まず，一つは藤田らによって開発 された LSPIV(Large-Scale Particle Image Velocimetry)である 9)～23)．この手法は画像のテン プレートパターンマッチングを利用し，流速を測定する PIV を実スケールの河川に適用 したものである．パターンマッチングとは画像上における輝度分布の面的な移動ベクト ルを追跡する方法であり，LSPIV においては河川表面の濃淡分布によって表面流速分布 が求められる．輝度の濃淡を考慮するため，トレーサー（例えば，浮子や木片）を特に 用意せずとも，表面流速を測定することが可能である．基本的に LSPIV では，ビデオカ メラによって撮影される斜め画像を，カメラキャリブレーションと幾何補正によって無 歪画像として変換し，解析に用いる． しかし，LSPIV においては，無歪画像を作成する過程に時間がかかる点，幾何補正変 換の限界（例えばカメラから遠地点は少ない画素数で描画され，カメラの俯角が小さけ れば特に顕著になる）といった点が，解析上の障害となる場合がある．さらに加え，画 像に映り込むノイズ，遮蔽物などを考慮しなければならないといった問題点を解決する ため，藤田らは時空間画像を利用した STIV(Space-Time Image Velocimetry)を開発した 24)-36) ．STIV では計測対象を主流方向の一次元表面流速に限定することによって，流量観 測法としての実用性がより高められている．さらに，STIV は LSPIV と比べ，計算の高 速化による短時間での計測が可能となっており，悪天候時の解析にも強くなっている． また，原によって，準リアルタイム流速計測の STIV が開発された 37)． 1.5 近年の洪水災害の増加と河川情報システムの現状 近年では台風・集中豪雨災害が急激に増加している(図 1.5)．国土交通省によれば， 10 年間で水害被害 2.6 倍と言われる．また，気象庁の統計によれば，アメダス 1000 地 点あたりにおける 1 時間降水量 50mm 以上の年間観測回数は統計期間 1976～2011 年で 増加傾向が明瞭に現れており，1 時間降水量 80mm 以上の年間観測回数についても同期 間で増加傾向が現れている(図 1.6)38)．一方、日降水量 200mm 以上の年間観測回数につ いては同期間で変化傾向は見られないが、日降水量 400mm 以上の年間観測回数について は増加傾向が明瞭に現れている(図 1.7)．また，気象庁によって「平成 21 年 7 月中国・ 九州北部豪雨」と命名された，山口，福岡県を中心に襲った豪雨災害では大規模な土砂 崩れや増水が発生し．死者は 32 人にも至った．このように，都市地域での災害は，甚 大な被害を我々の社会に与える．また仮に行政による災害に対する備えが十全に用意さ れ，救助活動が迅速に行われたとしても，住民の理解が無ければ，人的被害の減少は困 4

難なものとなる．そこで，この様な都市部での災害，特に水災害，による被害を最小限 に抑えるためには，より一層高度で河川管理と住民の正確な状況把握が必要である． 図 1.5 水害被害額および水害密度等の推移（過去 5 年平均）1) 2012 年現在，日本の災害現場では，通信ネットワークの整備が注目され，災害に強い 通信網の構築が進められている．一例として国土交通省では、多重無線（マイクロ）回 線と光ファイバ回線による通信基盤と機動性の高い移動通信や衛星通信、ヘリコプター 映像伝送システム等を運用して災害時の迅速な情報収集、災害復旧活動を支援している。 2012 年 9 月の台風第 12 号による紀伊半島大水害では、大規模な土砂災害や河道閉塞に よる被災自治体の孤立や、復旧や捜索活動における二次災害の回避の必要性から河道閉 塞箇所の監視が必要となった．そのため，近畿ほか 4 地方整備局から TEC-FORCE とし て衛星通信車等を派遣して通信の確保、監視カメラによる映像監視を行い、専用通信網 やインターネットで映像のリアルタイム配信を行っている 39)．河川整備においても，ネ ットワークシステムの発達により，洪水時の雨量・水位等の情報や画像をリアルタイム で入手できるようになった．ただし，流量に関するリアルタイムでの公表はなされてい ない．前述した非接触型計側法や H-ADCP といった流量観測システムを組み合わせるこ とによって，リアルタイムで流量を取得することが可能となる 40)-42)．しかしながら，こ れらの設備は現状では大河川のみを対象としている．さらに，中小河川では流量や水位 などの観測施設の不備が流量の把握を困難なものとしている． 5

図 1.6 アメダス地点で 1 時間降水量が 50mm（左）， 80mm（右）以上となった年間の回数（1,000 地点あたりの回数に換算）38) 図 1.7.アメダス地点で日降水量が 200mm（左）， 400mm（右）以上となった年間の回数（1,000 地点あたりの回数に換算）38) したがって，付近の住民への情報提供のため，中小河川に関する的確な河川情報の取 得を可能とするシステムネットワークの整備が急務とされている．そこで，河川情報の 取得のツールとして注目を集めているものに CCTV の整備がある．CCTV の映像から取 得したデータを発信することで河川管理，市民に対する情報提供など様々な用途に活用 される事が期待される．さらに，24 時間リアルタイムでの河道の監視は，先に説明した 画像データを利用する非接触型流速計側をどの時間帯でも行え，洪水時の流量データの 欠測防止を実現できる．また，河川流のみならず水質や動植物，人の動きなど様々な計 測も可能とする．しかし，以上のような CCTV を高度利用した例がほとんどなく，通信 ネットワークによる情報の共有を可能とする環境の整備は，改善の余地を必要としてい る． 6

図 1.8. 衛星通信車（左），配信映像（右）39) 1.6 本研究による治水への貢献 本研究では，画像を用いた河川流速計側の実用化を目的としている．具体的には，こ れまでの画像解析技術の応用と利便性を重視したソフトウェアの開発と精度について 述べる．まず，画像解析技術に必要なカメラキャリブレーションについて述べる．ここ では，カメラキャリブレーションの理論から，実際のソフトウェアのアルゴリズム，精 度と応用方法について述べる． 次に，神戸大学と土木研究所(ICHARM)との共同研究で開発した，画像浮子解析ソフ ト’Float-PTV’について述べる．本ソフトでは，河川漂流物を追跡することで河川表 面流速を簡単に計測できる． しかしながら，実際の河川洪水時には漂流物が常に画面上に表示されているというこ とはなく，さらに流れにそっているとは限らない．そのため水面上を流れる波紋によっ て，表面流速を測定する手法が必要となる．先に述べたように，STIV では検査線上の輝 度を時間方向に並べて作成した時空間画像を用いて流速を測定する．このような理由か ら，最後に新しく開発された新しい STIV での自動解析結果の精度向上の工夫について 述べる． 7

第2章 周辺機器と開発環境 2.1. 概要 画像による河川流速計側法は市販の PC とビデオカメラを必要とする．そのため，従 来から採られてきた浮子法などに比べ，本手法は簡便性・経済性に優れている．さらに， 洪水時において橋の上や河岸から撮影した画像を利用するために安全性にも優れてい る． この様にメリットの多い画像流速解析であるが，実際の河川管理の現場で利用では， 土木工学や画像工学の専門者以外からの“使いやすさ”が重視される．そのため，本研 究では使いやすく，理解しやすいことを重視した解析ソフトの改良，開発を主軸として いる．ここでは，河川管理現場における利用を目的として作成した解析ソフトについて 述べる． 本論における画像流速解析ソフトでは，カメラから取得された画像と物理空間を空間 幾何学的に結び付けるカメラキャリブレーション，そして二つの解析方法をサポートし ている．この章では，画像取得機器，本ソフトの開発環境について述べた後，画像流速 解析で最も基本となるカメラキャリブレーション，Float-PTV（画像浮子追跡)， STIV(Space-Time Image Velocimetry)について順に説明する． 2.2. 構成 2.2.1. 画像取得機器 画像を取得する機器（製品には）様々なものがあるが，対象とする使い方に合わせて 性質が異なっているが，画像解析における計測精度は，利用する画像の質に大きく影響 されるため，計測対象に対して適切な画像取得機器を利用する必要がある．以下に一般 的に利用可能であり，画像解析に用いることができる画像取得機器の概要を述べる． 1) CCTV(Closed-Circuid TeleVision)カメラ いわゆる監視カメラであり，CCD カメラからアナログビデオ信号が出力されるものが 一般的である．様々な製品が販売されており，その機能も様々である．カメラの向きを 隔離操作するための自動雲台やズーム・アイリス（絞り）等の修正可能なものも多く， 主に RS-232C 規格の信号として PC からこれらの操作を指示することができる．取得可 8

能な画像の解像度や感度も様々なものが多く，走査線が 480 本（水平解像度に相当）で ある事から，HDV カメラやデジタルカメラと比較と解像度は良いとはいえない．フレー ムレートはインターレース（飛び越し走査）を用いた 60FPS であるが，それぞれの走査 線の更新は 30FPS`である． 監視を目的として製品化されたものは隔離コントロール，耐候性．按手性がある程度 確保されているため，河川流の連続観測という目的に対して最も適している． 出力される画像は前述のようにアナログ信号が一般的であるため，画像解析を行うた めにはデジタル信号に変換を行う必要がある．また CCTV 画像の伝達には有線・無線法 があり，それぞれアナログ信号を用いるものと．デジタル信号を用いるものに分けるこ とができる．例えば，現在国土交通省により一級河川の監視を目的として設置されてい る CCTV カメラの画像はデジタル信号に変換された上で，同省の設置した光ファイバー 網を通して転送されている． デジタル変換の際には 640×480 ピクセルに変換される例が多いが，その他にテレビ 映像の変換などは CIF（Common Intermediate Format）形式（352×288 ピクセル（CIF）， 704×288 ピクセル(2CIF)）が用いられる場合もある． 2) DV（Digital Video）カメラ 主に民生用に利用されるデジタル化されたビデオカメラである．かつては，テープに デジタル情報を記録させていたが，昨今では大容量化された SD カードや内臓ハードデ ィスクにデジタル情報を記録させている．近年では．DV カメラの信号は，DCT（Discrete Cosine Transform：離散コサイン変換）によるフレーム単位の非可逆圧縮を行い 1/5`程度 に容量を圧縮している．解像度は 720×480 ピクセルであり，フレームレートは 30FPS であるが，アナログビデオ規格と同様のインターレースを用いて 60FPS の画像更新を行 うことが多い．DV カメラの画像は，IEEE1394 規格の信号としてそのまま PC 上に取り 込むことができる． 3) HDV（High-definition Video）カメラ DV カメラと同様の記録媒体に 1,280×720 あるいは 1,440×1,080 の解像度の画像を記録 する民生用カメラである．同じ記憶容量により高い解像度の画像を記録するために，時 系列圧縮を取り入れて圧縮率を上げている．昨今では，このカメラによって，洪水時の 動画データが取得されており，フォーマットは MTS，M2TS などである．HDV カメラ の持つ解像度の高さによって，より鮮明な洪水時の水面パターンを確認することができ る．一方，圧縮率の高さは画像の取り扱いを複雑化し，画質の低下を防ぐ適切な処理が 9

必要とするといった問題も生じている． 4) WEB カメラ 画像習得から圧縮，ネットワーク転送まで含めたおもに監視用途のカメラシステムで， パンチルトなどの操作のコントロール等も可能な機種もある．様々な性能の製品が比較 的安い価格等で提供されており，機能面では CCTV カメラと同様であるが，情報ネット ワークの利用コストの普及，低価格化，大容量化により普及が進んでいる． 2.2.2. 開発環境 43)-51) 開発環境として OS は Microsoft 社製 WindowsXP あるいは WindowsVista を使用．開発 プログラム言語に C++，IDE（統合開発環境）に Visual Studio となる Visual C++2008 を 用いた．Visual Studio は Microsoft が提供するソフトウェア開発製品群及びそれらを管理 する統合開発環境であり，Windows 上で最適化されたソフトウェアを作成することが容 易となる．また，画像処理にフリーライブラリの OpenCV2.0 を用いている．OpenCV は Intel 社が開発・公開したオープンソースのコンピュータビジョン向けライブラリである． 10

第3章 カメラキャリブレーション 3.1 画像の座標変換 52)-54) 画像記録により取得された画像を用いて，流速情報を抽出するためには，撮影座標（ビ デオ画像上の座標）から実空間座標へと座標変換する必要がある． 図-2.3.1-1 に示すように，ビデオカメラから取得された斜め画像では，画面上での物理間 隔は一定ではない．そのため，演算処理中においては斜め画像を無歪画像に変換する必 要がある．本節では本研究で用いた座標転換式について説明する． 図 3-1 斜め画像（左），無歪画像（右） 実際の出力として斜め画像に幾何補正を施すことによって得られる画像は無歪画像と 等価である．図-2.3.1-2 に示すように実際の斜め画像では撮影範囲が限られている．その ため，幾何補正を施すと画像情報の存在しない領域が発生する．今回開発したソフトで は，画像情報の存在しない領域は灰色で表現されている． 図 3-2 真上画像としての幾何補正画像（左）と実際の幾何補正画像（右） 12

x, y をビデオ上の投影座標，X, Y, Z を実空間の物理座標とするとその座標変換式は一般 的に， (3.1) のような関係となる．式中の関数 f および g に用いる式として，カメラの物理的なモデ ル化により算出される共線条件式を適用する．図 3.1 に示すように， カメラ座標系を xyz， 写真座標系を xy，地上座標系を XYZ とする．ここでは，カメラをそれぞれ z 軸の正方向 に対して左回りに，次に y 軸の正方向に対して左回りに，さらに z 軸の正方向に対し て左回りにの角だけ回転させた傾きで写真撮影を行ったとする． 図 3-3 単写真の幾何角 回転行列を用いると，地上の対象物 P( X, Y, Z )の，傾いたカメラ座標系におけるカメラ 座標( xp, yp, zp )は次の式で求められる． (3.2) ここで，X0, Y0, Z0：投影中心の地上座標であり，実際はカメラのレンズ位置である．投 影中心，写真像および地上の対象物 P( X, Y, Z )が一直線上にあるという共線条件式はつ ぎのようにして求められる． 13

(3.3) ここに，c：カメラの焦点距離，aij ( i, j = 1, 2, 3 ) はカメラと回転行列 R の要素であり， 次式で計算される． (3.4) ここに， である． また，式( 3.3 )の逆変換式は次式で表わされる． (3.5) ところで，以上の変換式はレンズの焦点距離および主点の位置が正確で，レンズディス トーションがなく，フィルム面も平面に保たれている測定用カメラに適用される共線条 件の基本式である．このように完全な測定用カメラでない場合には，式( 3.3 )または( 3.5 ) の写真座標 x および y から次の補正量x，y をそれぞれ差し引かなければならない． 14

(3.6) ここで，( x0, y0 )は主点位置（理論上は写真の真ん中）のずれを示し，( k1r2+k2r4 )を含む 項は放射方向歪曲収差によるレンズディストーションを示し，最後の p1～p6 を含む項は フィルムの平面度の歪を示す． 本解析ソフトで用いる幾何変換式は式( 2.3.3 )と( 2.3.5 )であり，出力される変数は X0, Y0, Z0, κ,x0, y0 k1, k2, p1, p2, p3, p4, p5, p6 の 16 個である．この変数の内 X0, Y0, Z0, κ, は カメラの外部標定要素，x0, y0, k1, … , p6 は内部標定要素と分けられる．これらの変数を 求めるために，空間座標値が既知である標定点を利用して，最小二乗法を適用する．ま た，キャリブレーションを成功させるためには，初期値を与える必要があり，その初期 値如何によっては非物理的な局所解に陥る場合や収束計算が発散する場合がある．よっ て，本論文における解析ソフトではこれらのパラメータの初期値を容易に設定・変更し， 計算結果の妥当性を視覚的に確認するための GUI（Graphical User Interface）を備えてい るようにした． 3.2 アルゴリズム 1). 全体の流れ 図 3.4 はカメラキャリブレーションの全体の流れを表わしている．カメラの焦点距離と 河川の水位は UI（ユーザーインターフェイス）から入力される．次に，カメラの座標位 置が既知であれば UI から入力しておき，未知であれば逆計算で再計算が行われる．一 般にカメラの座標位置が既知である方が良好なキャリブレーションを得やすい．次に， 標定位置の情報を入力する．入力は”points.csv”と名付けたファイルから読み込まれる場 合と，UI から入力する場合の二通りがある．初期値を入力した後に計算命令を与え，キ ャリブレーションを開始する．まずは，標定要素から最小二乗法を利用してカメラの外 部標定要素の初期値を計算する．次に，初期値を用いて，カメラの外部標定要素と主点 位置のずれをニュートン法による収束計算によって算出する．最後に，カメラのレンズ ディストーション係数 k1,k2 と，フィルムの平面度の歪を表わす変数 p1,p2…,p6 を最小二 乗法によって計算する．UI からのオプションによって，レンズディトーション係数と， フィルムの平面度の歪の変数を計算しないことも選択できるようにした．これは，標定 15

要素が少ない場合や最小限のカメラのパラメータを計算したい場合に用いる．最後に出 力された変数から作られた式( 3.3 )と( 3.5 )を用いて標定要素の再計算が行われる．物理 座標と写真座標の再計算と実際における誤差データと作成された幾何補正画像からキ ャリブレーションの成否を判断する． 図 3.4 キャリブレーションソフトの全体プロセス 2) カメラの外部標定要素の初期値の算出 式( 3.5 )は，測定対象物の標高 Z が与えられれば，その点の地上での位置( X, Y )が，対応 16

する単写真の写真測量( x, y )から求められることを示している．測定対象物が平面の場 合には，式( 3.5 )のかわりに，次式( 3.7 )に示すような二次の射影投影変換式をもちいて よいことが，射影幾何学の原理より誘導されている． (3.7) 一方，物理座標( X, Y, Z )と投影座標( x, y )の関係は次式の共線条件式( 3.8 )でも表される． (3.8) 式( 3.8 )の逆変換式は次式( 3.9 )となる． (3.9) 式( 3.7 )と( 3.9 )は次式群( 3.10 )と( 3.11 )によって等価な式とされる． 17

(3.10) (3.11) 後述するように，既知の標定点が 6 つ以上の時，式( 3.9 )を，6 未満の時は式( 3.7 )を最 小二乗に用いている．ただし，河川の水面は一般的に平面の方程式として次式( 3.12 )を 用いている． (3.12) 係数 D1, D2, D3 は水面上で 3 点以上の座標を与えれば求められる．ただし，水面が緩勾 配ならば，Z=D3 とおける．そのため，式( 3.10 )における初期値の計算では， としている． 式( 3.7 )と式( 3.8 )においてそれぞれ，最小二乗法を当てはめる時には次に示すような線 形の観測方程式( 3.13 )，( 3.14 )を立てる． (3.13) (3.14) 18

この時 xi,yi,Xi,Yi は i 番目の基準点の測定値の組である．式( 3.8 )及び( 3.14 )から最小二乗 法を用いて得られた係数 A1,A2…C3 は式( 3.10 )と( 3.11 )によって b1,b2…,b8 に変換される． 標定点が 6 点未満の時は b1,b2 ...,, b8 を式( 3.13 )から最小二乗法によって直接求める．そ して，b1,b2…,b8 からカメラの外部標定要素の初期値を以下の式群( 3.15 )によって算出す る．ここで，カメラの位置 X0,Y0,Z0 が既知の固定値であるとすれば，この変数に関する 初期値計算のプロセスは除外できる．さらに，カメラの回転角 κ,が既知であり固定 値である場合も UI からのオプションによって計算を除外できる．各係数間の関係式は 以下の通りである． (3.15) ここで，A1,A2,A3,A4 は であり，Zm は河川の水位である． 最後に，全体の流れを図 3.5 に示す． 19

図 3.5 外部標定要素の初期値を算出するプロセス 3) ニュートン法による外部標定要素とレンズ主点位置のずれの算出 先に求めた外部標定要素の初期値，もしくは既知の値によってニュートン法（逐次近 似代入法）による未知変量の算出を行う．まず，式( 2.3.3 )から写真座標に主点位置を差 し引いた次式( 3.16 ) (3.16) 20

を考える．このような非線形の連立方程式から未知変量を求める場合，未知変量の近似 値を与え，近似値のまわりにテーラー展開して線形化し，最小二乗法により補正量を求 めて近似値を補正し，再び同様の操作を繰り返し，収束解を求めるという手順を必要と する． 次に，式( 3.16 )を外部標定要素と主点位置のすれに関する関数として，次式( 3.17 )の ように書きなおす． (3.17) ここでの x,y は観測された写真座標を表わす．外部標定要素と主点位置のずれの近似値 をそれぞれ X00,Y00,Z00,0,xy として，その補正量をそれぞれ  X0,Y0,Z0, ,x yとすると次のようにかける (3.18) 近似値のまわりに式( 3.17 )をテーラー展開し，2 次元以上を無視すると，次式( 3.19 )の ように線形化できる． 21

(3.19) 式( 3.19 )の微係数は，つぎのように与えられる． 22

(3.20) ただし，微係数を求める式における x,y には式( 3.16 )に 6 つの外部標定要素 X0,Y0,Z,0κ,  の近似値を与えた時に得られる値を用いる．カメラの焦点距離 c は，ソフトの UI から 初期値を入力され，主点位置のずれ x0,y0 は未知であれば初期近似値を 0 で計算される． これらの値に加え，カメラの位置 X0,Y0,Z0，カメラの回転角 κ, が既知の固定値である 場合，収束計算におけるそれぞれの補正値は 0 であるとして計算される． 式( 3.19 )，式( 3.20 )からの 2×n 個の 9 元一次の観測方程式に最小二乗法を適用して，補 正項X0,Y0,Z0,,xyの最確値を求める．次に式によって，近似値 を修正し，同様の操作を繰り返す．このようにして繰り返し近似解を補正して収束解を 求める． 収束したか否かの判定は，次の式( 3.21 )で与えられる中等誤差の前後繰り返し 1 回おけ る偏差iが，10以下になったとき収束したとみなす．普通 4，5 回の繰り返しで収 束する． 23

(3.21) ニュートン法によるカメラの外部標定要素と主点の位置のずれの算出プロセスを図 3.6 に示す． 図 3.7 ニュートン法によるカメラ標定要素の算出 24

4) レンズディストーション係数の算出 式( 3.3 )，( 3.6 )に収束計算で求めた x0,y0, X0,Y0,Z0, κ, を代入し，レンズディストーシ ョン係数を最小二乗法によって算出する．ここでは先に求めた，カメラの標定要素を基 準として，より正確に標定点の写真座標を再計算できるようなレンズの歪み係数を求め ることになる．最初に接線方向のレンズの歪み係数 k1,k2 を最小二乗法によって求める． 次にレンズの歪み係数を考慮して計算された標定点の写真座標と画面から与えられて いるそれの平方誤差が 10以下に収束するかを確認する． 図 3.8 レンズディストーション係数の算出プロセス 25

誤差が 10を超えていれば，計算された k1,k2 を式(2.3.6)に代入して，カメラの標定要素 を再び収束計算によって求める．計算が収束すれば，残るレンズの法線方向 p1,p2…,p6 を最小二乗法で求め，キャリブレーションを終了する．図 3.8 にこの計算のフローチャ ートを示す． ここまでのキャリブレーションが成功していれば，図 3.9 に示す斜め画像上で，標定点 の真ん中にある青い丸に赤い丸が収まり，幾何補正画像が生成される．青い丸は既知の データとして与えられた標定点の写真座標であり，赤い丸はキャリブレーションによっ て再計算された標定点の写真座標である． 図 3.9 キャリブレーション成功時における斜め画像と幾何補正画像 また，図 3.9 の斜め画像において標定点から生じている緑の直線は，標定点の真下にあ る水面での位置を表している．これが水面と直角であるほど，キャリブレーションの精 度は高いといえる．しかしながら，レンズの歪が大きい場合，緑の直線は水面と斜めに 交わることもある． 次にカメラキャリブレーションのレンズ歪補正における効果について検証を行う．実 験では，正方形タイルが敷き詰められた床の上に 3 種類の箱と正方形を描いた紙を用意 した．3 つの箱の長さをそれぞれ測定し，箱の角における物理座標と写真座標を標定座 標とする．この時の物理座標はカメラを原点とした任意の直交座標系である．図 3.10 は， 実験での様子をビデオカメラで習得したものである． 26

図 3.10 カメラキャリブレーション実験図 図 3.11 はカメラキャリブレーションによって作成された幾何補正画像である．左がレン ズ歪を補正したカメラキャリブレーションによる幾何補正画像である．右がレンズ歪を 補正しなかったカメラキャリブレーションによる幾何補正画像である．二つの画像を比 較すると，白い方眼紙の角度が違うのが分かる． 図 3.11 幾何補正画像（左；レンズ歪補正有り，右；レンズ歪補正なし） 床のタイルに注目する．分かりやすくするために，図 3.12 に示すように，両方の画像に グレースケール処理とエッジ強調を施した画像を用意する．レンズ歪無しの幾何補正画 27

像では，床のタイルが奥（画面上部）にあるほど歪んでいるのがわかる． 図 3.12 グレースケールとエッジ強調（左；レンズ歪補正有り，右；レンズ歪補正なし） エッジ検出を用いて床のタイルを検出し，直線を引いて見やすくした図 2.3.2-8 を示す． レンズ歪無しの幾何補正画像上では，幾何補正画像の右上部にあるタイルの歪みが大き くなっている． 図 3.13 床のタイル比較図（左；レンズ歪補正有り，右；レンズ歪補正なし） この事から，レンズ歪補正は幾何補正の精度に少なからず影響を与えている．しかし 28

ながらレンズ歪の影響を補正しなかった場合でも，手前のタイルの形は正方形の形をし ている．レンズ歪なしの画像であっても，幾何補正画像としての役割を果たせる． 5) 幾何補正画像の作成 最後に幾何補正画像について述べる．幾何補正に用いる変換は透視投影変換（ホモグ ラフィー）と言われ，式( 3.9 )による変換式で表わされる．また，作成の手順として最初 に対象領域の切り出しを行う．ここでは，斜め画像の上下左右の両端から上下，左右そ れぞれの総ピクセルの 5 分の 1 に囲まれる領域を対象とする（図 3-14）． 領域が斜め画像の端であると，レンズ歪の影響を受けるためである．各点 A,B,C,D を 写真座標から，物理座標へとキャリブレーションによって求められたカメラの標定要素 を用いて式(3.6）によって算出する．次に，算出された物理座標を包括する最小限の長 方形領域を算出，これの各頂点における物理座標を写真座標に式（3.16）を用いて変換 し直す．この過程によって，カメラキャリブレーションの是非を幾何補正画像によって 確認する事ができるようにしている．変換に用いられる補間には高精度の補間法である ３次補間法（バイキュービック法）を用いている． この補間の用いられる式は， (3.21) であり，ここに，(xK,yL)：(x0, y0)の周囲の格子点座標，C(x)は次式で与えられる． (3.22) 図 3-14 幾何補正領域の算出方法 29

3.2 精度についての考察 カメラの位置が既知の時，計算に用いられる標定点の個数によってキャリブレーショ ン精度は左右される．鬼怒川の左岸から撮影されたビデオ画像を用いて，カメラキャリ ブレーションの精度について検証する．キャリブレーション計算中に最小二乗法を用い ることから，標定点は手前と奥にバランス良く設置されている事が望ましい．今回の精 度検証に用いられている標定点の位置配置はキャリブレーションの計算にもっとも適 した位置配置である． 図 3-15 に,既知の標定点 7 点を用いて，キャリブレーションを行った結果を示す．左 の斜め画像によって，キャリブレーションによって再計算された標定点の座標を視覚的 に確認できる．再計算が上手くいけば，標定点の中心から緑の垂線が引かれる．右の図 はキャリブレーションによって作成された幾何補正画像を示す．図によれば，標定点の 座標が精度よく再計算されていることを確認できる．下の図 3-16 にそれぞれの座標の実 測値と計算結果の差を示す．全ての標定点を用いてキャリブレーションを行った場合， 再計算された値は元のオリジナルの値を精度よく再現できていることが図 3-16 より読 み取れる． 図 3-15 標定点 7 点におけるキャリブレーション結果 （左；斜め画像，右；幾何補正画像） 30

図 3-16 標定点 7 点における物理座標と写真座標での 実際の値と再計算結果の誤差（上；物理座標，下；写真座標） 次に，標定点を 3 点に減らしてキャリブレーションを行った結果を示す(図 3-18)．キ ャリブレーションに用いた標定点は①，③，⑤の位置にある点である（図 3-17）．物理 座標では，奥の点ほどオリジナルと再計算の誤差が大きい．これは，カメラのレンズか ら奥に行くにしたがい，1 ピクセルあたりの物理量が大きくなるためと考えられる．ま た，計算に用いていない標定点の周囲において誤差が大きい．写真座標に関しては計算 の用いられた標定点の座標は誤差が小さい．逆に，計算に用いられていない標定点は写 真の手前，奥に関わらず誤差が大きい傾向にある． 図 3-17 標定点 3 点（①，③，⑤）におけるキャリブレーション結果 （左；斜め画像，右；幾何補正画像） 31

図 3-18 標定点 3 点（①，③，⑤）における物理座標と写真座標での 実際の値と再計算結果の誤差（上；物理座標，下；写真座標） 3 番目に 3 点の標定点の配置が奥に 2 点手前に 1 点とした例を示す（図 3-19）．キャリ ブレーションに用いた標定点は②，⑤，⑦の位置にある点である．物理座標では，奥の 点ほどオリジナルと再計算の誤差が大きい傾向にあり，先ほどの 3 点の例と比べても誤 差が大きい（図 3-20）．写真座標に関して奥にある標定点座標ほど誤差が大きくなる傾 向にある．幾何補正に関しても手前の岸が歪んでいるのが確認できる． 図 3-19 標定点 3 点（②，⑤，⑦）におけるキャリブレーション結果 （左；斜め画像，右；幾何補正画像） 32

図 3-20 標定点 3 点（②，⑤，⑦）における物理座標と写真座標での 実際の値と再計算結果の誤差（上；物理座標，下；写真座標） 4 番目に，標定点に偏らせた 3 点を用いた例を示す．キャリブレーションに用いた標 定点は②，④，⑥の位置にある点である（図 3-21）．図 3-22 に物理座標，写真座標の実 測値と再計算値の誤差を示す．物理座標，写真座標共に画面端の座標において，オリジ ナルの座標と再計算による座標の誤差が大きい．幾何補正画像も標定要素の無い領域で は歪んでいる．写真上の右側の領域における標定点を計算に用いていないため，バラン スの偏ったキャリブレーション結果になっていると考えられる．さらに，これまでの 3` 点におけるカメラキャリブレーションの結果を比較すると，カメラ手前の情報がキャリ ブレーションの情報に重要になっている傾向にある． 図 3-21 標定点 3 点（②，④，⑥）におけるキャリブレーション結果 （左；斜め画像，右；幾何補正画像） 33

図 3-22 標定点 3 点（②，④，⑥）における物理座標と写真座標での 実際の値と再計算結果の誤差（上；物理座標，下；写真座標） 最後に 2 点のみをキャリブレーションの計算に用いた結果を図 2.3.3-5 に示す．標定点 は②，⑥である．再計算による座標は誤差が大きく，画面上にも表示されていない．幾 何補正画像は非常に歪みの大きい画像が作成される．このことから，カメラキャリブレ ーションによるカメラパラメータの再現には最低でも 3 点が必要である事が考えられる． また，3 点の標定要素を直線に与えた場合，キャリブレーションは成功しなかった．そ のため，通常のカメラキャリブレーション計算で用いられる標定点では，奥行き方向の 情報も必要となる．ただ，後述のようにカメラの座標と回転角が既知の時には 1 点のみ でカメラの標定要素が計算でき，幾何補正画像を作成する事ができる． 図 3-23 標定点 2 点（②，⑥）におけるキャリブレーション結果 （左；斜め画像，右；幾何補正画像） 34

3.3 最小条件からのキャリブレーション これまでのカメラキャリブレーションでは，先に述べたように標定点を実河川に設け， 事前に物理座標を取得しておく必要があった．しかしながら，山中にある実河川などで は，対岸に標定点を設けることが難しい．そのため，カメラの俯角と河川までの距離の みという少ない情報から簡易的なキャリブレーションを成功させる技術の開発が必要 となる．ここでは，簡易型キャリブレーションに関する実空間と仮想空間における実験 と成果について述べる．カメラ，三脚，チェスボード，角度計を用意し，図 3-20 のよう に設置する．この時，カメラのレンズの真下を原点とする任意の物理座標軸を考える. 図 3-21 には実際の実験で習得した，ビデオカメラからの斜め画像を示している．ここで用 いられているチェスボードは一辺 0.2m の白と黒の正方形を 4×4 四方につなぎ合わせて 作成している．角度計は 0.1Degree までの値を測定できる Apple 社製 iPhone，iPad での 計測系アプリケーションソフトを使用した． 図 3-24 簡易キャリブレーション実験図 35

図 3-25 実験映像 表.3-1 にカメラのパラメータと実験における a 点の座標のデータを示す．ここで，標定 点となる座標は最低でも一組は必要である．最低限のカメラキャリブレーションでは， カメラの位置と俯角等の情報のほかにレンズの中心位置の座標が必要となる．レンズの 中心位置の座標のみが未知数の時，一組の標定座標によって方程式が閉じる． 表 3-1 実験におけるカメラのパラメータと標定点座標（真正面） 図 3-26 にキャリブレーションによって作成された幾何補正画像を示す．斜め画像では歪 んで映っているチェスボードが，幾何補正画像では正方形に補正されている．この結果 より，キャリブレーションが簡易的に行われた事がわかる． 36

図 3-26 実験結果（真正面） 次に，任意の座標系でも簡易型のキャリブレーションができる事を確認する．標定点の 座標に回転行列による変換 を施し，回転分の角度をカメラのパラメータに与える事でキャリブレーションが成功す るかを試した．ここで，回転させる角度は 45°を与えた．カメラのパラメータでは z 軸 周りの回転角にを与えた方向を見ていることになる（図）． 表 3-2 実験におけるカメラのパラメータと標定点座標（回転） 37

図 3-27 回転座標の与え方 図 3-28 に回転された座標系におけるキャリブレーションによって作成された，幾何補正 画像を示す．ここでもチェスボードは，正方形に幾何補正されている．幾何補正画像が 図 3-26 における幾何補正画像と比較して斜めに傾いているのは z 軸周りの回転角を与え た影響と考える事ができる． 図 3-28 実験結果（回転座標） 38

次に仮想空間を想定し，仮想空間におけるカメラ座標と標定点から簡易型キャリブレ ーションが可能かを実験した．仮想空間の描画にはオープンソースの OpenGL を用いた． 図 3-29 に OpenGL によって描画された斜め画像と内部設定での視点移動によって得られ る，真上からの画像を示す．緑の板は実験におけるチェスボードと同じ役割を持ち，標 定点でもある．ここでは，座標を( 0, -3.0, 0 )に緑の板の一辺が触れるように配置した． 斜め画像では OpenGL の仕様により緑の板の辺が座標( 0, -3.0, 0 )より離れている． 図 3-29 仮想空間（OpenGL） カメラのパラメータは以下の様に与えた．図 3-30 に幾何補正画像結果を示す．標定点の 緑の板と格子の床が正方形に幾何補正されていることがわかる． 表 3-3 仮想空間におけるカメラのパラメータと標定点座標（真正面） 39

図 3-30 幾何補正画像（OpenGL） 最後に，標定点となる緑の板を 45°回転させた方向に設置し，カメラの位置と姿勢， 標定点の一組の座標を与え，キャリブレーションを行った．表 3-4 にカメラの姿勢と標 定店に関する情報を示す．図 3-31 に斜め画像と内部設定によって真上からみた画像を示 す．原点から 45°の方向に標定点を設置している．やはり幾何補正画像は右回転に回転 しており，床のグリッドが正方形に補正されているのがわかる． 表 3-4 仮想空間におけるカメラのパラメータと標定点座標（45°回転） 40

図 3-30 仮想空間（OpenGL，回転） 図 3-31 幾何補正画像（OpenGL，回転） 最後に，図 3-23 で示したような 2 点のみの標定点座標が既知の場合に，あらかじめ成功 したカメラキャリブレーションで算出しておいたカメラの位置 X0, Y0, Z0，カメラの角 度ω,φ,κを与えてキャリブレーションを行った結果を図 3-32 に示す．2 点のみの標定 点のデータにも関わらず，良好な幾何補正画像が得られているのが分かる． 41

図 3-32 幾何補正画像 （鬼怒川，標定点 2 点に既知のカメラパラメータを用いた場合） 3.4まとめ ここまで，カメラキャリブレーションに関する内容を述べてきた．ここで，簡単にま とめておく． ① カメラキャリブレーションとは 複数の物理座標と写真座標の組によって，カメラパラメータ（標定要素）を算出する 方法．カメラパラメータを明らかにする事によって，空間幾何学に基づいた画像の校正 を行え，解析に用いる事ができる． ② カメラキャリブレーションのレンズ歪補正の効果は レンズ歪の補正を行うことによって，より正確な幾何補正を作成する事ができる．し かしながら，実用的にはレンズ歪を用いずとも実用的な幾何補正画像の作成は可能であ る． ③ カメラキャリブレーションに必要な標定点は（カメラの座標は既知である場合） 最低，三点であり，標定点の配置にも左右される． ④ 簡易型カメラキャリブレーションについて 0.1Degree の精度での角度とカメラの位置が既知であるなら，一点の標定点（手前の岸ま での長さ）によってカメラキャリブレーションを成功させる事が可能． 42

第4章 Float-PTV 4.1 計算の仕組み Float-PTV は，河川上を流れる浮子を画像上で追跡し，前節で述べた写真座標 x,y と物 理座標 X,Y の関係式 (4.1) (4.2) と二点間を流れる時間を用いることで,，簡易に表面流速を測定することができるソフト ウェアである．図 4.1 に斜め画像と幾何補正画像上における，河川を流れる浮子をあら わしている． 図 4.2`に示すように，浮子は 2 点間を流れている．左側の図に写っている浮子がt を かけて右側の図に写るような位置に移動したとする．この時，浮子の 2 点間における物 理座標(Xstart, Ysrtart)と(Xend, Yend)が斜め画像上の写真座標から式( 3.5 )を用いて求められ る．物理座標が求められば，X 方向，Y 方向`の流速と絶対値を式(4.1)から求めることが できる． 43

図 4.1 Float-PTV における斜め画像と幾何補正画像の関係 図 4.2 Floart-PTV の計算方法 44

(4.1) 4.2 画像の多重合成 以上に示すように，Float-PTV は二点間の浮子の位置から簡単に流速を求められる方 法である．しかしながら，河川管理の現場において，ソフトの使用者が，どの位置の河 川表面流速をどの程度取得すればよいのかを，必ずしも把握しているとは限らない．そ こで，本ソフトでは河川行政の専門家以外でもソフトを使用して流速を測定できるため の工夫として，浮子の流れている動画像を多重合成し，浮子の軌跡を簡単に追うことが 出来る機能を設けた．UI からの操作では合成画像によって描画された浮子の上をマウス クリックして，計算ボタンを押すだけで流速値が算出できるようにした． 図 4.3 浮子の画像合成 浮子が写っているビデオ動画を時間方向に等フレーム間隔で合成させる．合成画像は， ある基準となるフレーム上における座標に対して，他のフレームにおける同座標上の輝 45

度がもっとも高い（あるいは低い）ピクセル情報を抜き出し，それを当てはめることに よって作成される．上の図 4.3 における画像合成ではフレーム合成間隔 10F（動画のフ レームレートは 30FPS である）としている．図 4.4 に 3 通りのフレーム間隔で合成させ た斜め画像と幾何補正画像を示す．10F 間隔では浮子が多すぎて，クリックするのが大 変である．一方で 60F 間隔では浮子の数が少なすぎるので，このケースの場合では 30F 間隔で多重合成するのがもっとも使用に適していると考えられる． 図 4.4 さまざまな条件での合成画像 このマウスクリックの動作では，浮子の根元をクリックするように気をつけなければな らない．これは，浮子の頭は河川上を流れながらに大きく上下に動くからである．斜め 画像上を大きく上下に指定すると，正しく表面流速を追うことが出来ない． 4.3 さまざまな状況での多重合成と計算結果 図 4.5 に鬼怒川を左岸から撮影したビデオ画像と Float-PTV を用いて流速計測を行っ た例を示す．左に斜め画像と流速を測定した領域を示す．右に幾何補正画像の合成画像 を示す．また，合成時間間隔は 100F（= 約 3sec）としている．この例では浮子が合成画 像上に鮮明に映っているため，流速を測定できる領域が多い．表 4.1 に流速データを示 す．このソフトでは，漂流物が合成画像上に鮮明に映っているか否か流速計側において 重要となる．そのため，以降に示す例では，漂流物がどのように映るかを示していく． 46

図 4.5 鬼怒川（右岸）での合成画像（左；斜め画像，右；幾何補正画像） 表.4.1 鬼怒川の流速データ例 夜間に鬼怒川を右岸から撮影したビデオ画像に Float-PTV を使用した例を示す．夜間 では浮子を視認できないため，蛍光塗料によって発光する浮子を利用している．図 4.6 に合成画像を示す．合成間隔は 10 フレーム（= 約 0.3sec）である．高い輝度に注目して 合成画像を作成することによって，発光浮子を追跡できている． 47

図 4.6 夜間における鬼怒川（右岸）での合成画像（左；斜め画像，右；幾何補正画像） 次に鬼怒川を左岸から撮影した場合について述べる．図 4.7 に合成画像を示す．この 時の合成間隔は 10 フレームである．合成画像を見ると，画面手前の浮子が大きく上下 に移動した事が分かる．これは，水面の凹凸が大きい事を表わしており，この付近での 流速も急変する．また，合成間隔が小さいために浮子の軌跡が途切れなく続いているが， フレーム間隔を大きく取ることによって浮子が水面下に潜り，見えなくなることもある． 図 4.7 鬼怒川（左岸）での合成画像（左；斜め画像，右；幾何補正画像） 夜間での鬼怒川左岸からの撮影画像について述べる．こちらも，発光浮子を利用して， 軌跡を合成している． 流れの様子がよく表わされており，水面の凹凸の様子についても見る事ができる． 48

図 4.8 夜間における鬼怒川（左岸）での合成画像（左；斜め画像，右；幾何補正画像） 最後に鬼怒川を上流から下流に撮影した画像について述べる．図 4.9 に合成画像と幾 何補正画像を示す．斜め画像では，浮子の軌跡が見て取れる．しかしながら，画面の奥 に行くほど，浮子が小さくなり追跡が困難となっている．図 4-10 に衛星画像と幾何補正 画像の比較を示す．衛星画像と幾何補正画像上での特徴がよく一致する事から，幾何補 図 4.9 鬼怒川（上流→下流）での合成画像（左；斜め画像，右；幾何補正画像） 図 4.10 鬼怒川（上流→下流）での衛星画像（左）と幾何補正画像（右）の比較 49

正画像では，上空からの地理情報を良く表わすことが可能であると考えられる．夜間に おける画像を図 4.10 に示す．発光浮子の軌跡を視認する事ができる． 図 4.11 夜間における鬼怒川（上流→下流）での合成画像 （左；斜め画像，右；幾何補正画像） 以上から，浮子の追跡による画像解析は表面流速を測定できる．しかしながら，河川の 状況によっては浮子が水面下に隠れてしまう．また，浮子が流れ中の障害物に引っかか るなど，流れと平行に移動しない場合もある．そのため，河川表面のパターンを用いて 流速計側を可能とする技術が必要となる． 2.4.3. まとめ ① Float-PTV は浮子を追跡することで，流速を計測するソフトウェアである． ② 浮子追跡には，多重合成画像を用いる． ③ 多重合成によって浮子追跡は容易になる．しかし，河川の状況によって，浮子が見え なくなる場合もある．そのため河川表面の波紋パターンを追跡する画像解析技術が 必要となる． 50

第5章 STIV(Space-Time Image Velocimetry)55) 5.1 理論 1) 時空間勾配による流速計側法 STIV は動画像の移動追跡に時空間勾配，輝度の空間勾配を利用した手法であり，輝度 の勾配値を画像処理によって求めている．そのため STIV は輝度の勾配値を求めるため に時空間画像（STI : Space-Time Image）を作成する必要がある． 図 5-1 にビデオカメラからの斜め画像と斜め画像上における検査線から作成した時空 間画像を示す．この例では時空間画像生成に 1 秒間 30 フレームのビデオ画像を用いて いる．STIV で用いる時空間画像はまず図-5-1 に赤線で示すような検査線を主流方向に設 置し，その検査線上の輝度分布を時空間方向に積み重ねて生成する．この時，ビデオカ メラ画像上の座標(x,y)における輝度 L(x,y)を算出するために以下の式( 5.1 )を利用する． (5.1) ここで R(x,y)は座標(x,y)における赤色の輝度値，G(x,y)は座標(x,y)における緑色の輝度値， B(x,y)は座標(x,y)における青色の輝度値である． 図 5-1 ビデオカメラからの画像（左），時空間画像（右） 実際には，斜め画像から河川の主流方向に平行に線を引くことが難しい場面もある．そ 51

の時は，図 5-2 に示すように，幾何補正（真上）画像上で主流方向を確認，検査線を設 置し，前節で示した式(3.3)と(3.12) を利用して斜め画像に投影させることで適正な検査線の設置を行う事ができる． ちなみに，上の式で用いられている写真座標 x と y には式( 3.6 )に示す補正量x，y を 考慮した値を用いる． 図 5-2 真上画像（左）上に設置された検査線が 斜め画像（右）に反映されている また時空間画像を作成する際に用いられる，輝度情報は周囲の点からの補間によって 算出している．算出方法を以下に記す．物理座標における検査線長(m)を補正後の検査線 長(pixel)の数に当分割する．各点の物理座標を式(2.3.3)により写真座標に変換する．撮影 画像における輝度を時空間画像の長さ軸方向に並べることで，検査線を補正する．輝度 の補間処理には周囲 4 個の格子点における輝度値を用いて補間する線形補間を用いる ( 図 5-3 )． 52

図 5-3 輝度値の線形補間 ここでは式( 2.5.2 )を用いて格子点間の距離ではなく，面積の比を用いて輝度値に重み付 けを行う．求める輝度値を Cn とすると (5.2) で示すことができる． ところで，時空間画像では多くのノイズも確認できるが，全体的に左上から右下へ向かう縞模様が 確認できる．この縞模様はほぼ一方向を向いているため，縞パターンの平均的な傾きが求められれ ば検査線上を通過している平均流速が得られる．検査線上の平均流速は次の式( 2.5.3 )から求めるこ とができる． (5.3) ここに Sx ：検査線軸の単位長さスケール(m/pixel)，St ：時間軸の単位時間スケール (sec/pixel)である．Sx の値は場所によって変化され，個々の検査線の実長は前述のカメラ キャリブレーションの項で示した関係式 53

を用いて求められる． 2) 輝度勾配テンソル法 STIV においては，平均流速は時空間画像における縞パターンの傾きから求めているた め，STIV の精度は傾斜角の計測精度に依存する．よって傾斜角を精度よく解析すること が非常に重要となる．そこで藤田，椿 4)によって精度よく計測できることが報告されて いる輝度勾配テンソル法を用いる．以下にその手法を示す． 輝度勾配テンソル法は，時空間画像内に含まれる縞パターンの傾きを求めるために， 縞パターンの勾配ベクトル∇g，即ち (5.4) と局所的な単位方向ベクトル n の内積の二乗値 (5.5) が最大となるような方向ベクトルを局所的な検査領域 A に対して求める．すなわち (5.6) が目的関数となる．上式( 5.6)のカッコ内の積分は画像の輝度勾配に関するテンソルで (5.7) 54

とおける．ここで p, q は時空間の成分を表す指数でここでは，xp=x および xq=t とする． Jpq は対称テンソルなので座標系を適当に回転すれば対角化でき，そのとき式( 5.7 )を成 分表示すれば，回転後の変数には´をつけて式( 5.9 )のように書き表すことができる． (5.8) (5.9) この形にすれば， としたとき，単位ベクトル がこの条件を満たす．このことから，テンソル Jpq の固有値を与える回転角が求めるべ き角度ということがわかる．二次元問題の場合の固有値を λ1 および λ2 とすると固有値と 回転角の関係は次式( 5. 10)で与えられる． (5.10) この式の非対角成分の比較より 55

(5.11) が得られる．式( 5.11 )より回転角，すなわち縞パターンの勾配は次式( 5.12 )で求めるこ とができる． (5.12) また縞パターンのコヒーレンシーは次式 で与えられる． このコヒーレンシーは縞パターンの強さを示すパラメータであり，理想的な縞パターン で 1，全くパターンがない場合に 0 の値をとる指標である． 以上より，式( 5.7 )を用いて時空間各方向の輝度勾配の積を検査領域内 A で積分すれ ば，式( 5.12 )を用いて縞パターンの局所的な勾配が算出できる．輝度勾配の差分計算に は 5 次精度の中央差分を用いた．この差分式を空間軸方向について示せば， となり，また時間軸方向について示せば， となる．解析においては，検査領域 A を一辺が 10～50 画素程度の正方領域として時空 間画像上で縞模様の傾きおよびコヒーレンシーCc の値の分布を求める．ここで，ある 時空間画像に対して前述の計算を行い，縞模様の勾配ベクトルを求める．得られた勾配 ベクトルとコヒーレンシーの分布を図 5-4 に示す． 56

図 5-4 時空間上における勾配ベクトル（真ん中）とコヒーレンシー（右） 縞模様の勾配ベクトルに関するコヒーレンシーはつまり，縞パターンがはっきりと強く 現れている時にコヒーレンシーの値が 1 に近くなり，縞パターンがはっきりしておらず， 縞パターンが現れていない時にコヒーレンシーは 0 に近づく．図 5-4 のコヒーレンシー 分布図では，はっきりと見える斜め線の近くで濃い青色になっており，コヒーレンシー の値が 1 に近づいている．逆に勾配線が検出できない場所ではコヒーレンシーの値が 0 に近い値を描画しているのがわかる． STIV では時空間画像から得られる情報をもとに解析を進めていくため，精度の高い情 報を得ることが重要である．特に画像の情報が不十分な場合などはより正確な解析結果 を得ることが困難になる．そこで，コヒーレンシーの値を用いて時空間画像から得られ る情報をうまく選択して解析することにより，より精度の高い結果の導出が可能となる 3) 二次元画像信号に対するフーリエ変換 次に 2 次元画像信号に対するフーリエ変換について述べる． 1 次元離散フーリエ変換(DFT)を 2 次元に拡張することで，フーリエ係数 Fk,j を 2 次元画 像信号 fm,n，画素数 M×N を用いて (5.13) 57

と表されることが知られている．ただし，j は虚数単位である． 2 次元 DFT は 1 次元 DFT を繰り返すことによって計算する．すなわち，式(2.4.13)は次 式( 5.14 )と( 5.15 )のように分解できる． (5.14) (5.15) 各行( n =0,1,････ N -1)について式( 5.14 )の 1 次元 DFT を行い，次に各列( k =0,1,･･･ M -1)について式( 5.15 )の 1 次元 DFT を行えばよい．この過程を図 5-11 に示す．ちなみ に実際の計算は高速フーリエ変換を使用している． 図 5-11 2 次元 DFT のプロセス N 個のデータに対する 1 次元 DFT では，そのスペクトルは k =0~ N -1 に現れ， k =0 あるいは N-1 に近いほど低周波成分に，N/2 に近いほど高周波成分に対応していた．そ れが，図 5-7 の周期的な 2 次元画像の場合は，( k =0, l =0 )，( k =2, l =4 )，( k =126, l =124 ) の位置に強いスペクトルが現れる． 画像の処理は通常左上から走査するので縦軸の正方向を下向きとする．これに合わせ て，2 次元 DFT の結果も縦軸の空間周波数は下向きを正方向とする．得られた結果は図 5-12 の中央の図に示されるように，左上( k = 0, l =0 )の位置が直流分(平均濃度値)に，左 上，左下，右下，右上の 4 隅の近くが低周波成分に，中心部付近が高周波成分となる． しかし通常は，2 次元画像に対しては中心が直流分に，4 隅が高周波成分になるように スペクトルの配置を並び換えることが普通である．これは図 5-12 に示すように DFT の 結果に対し，左上と右下，右上と左下をそれぞれ交換することによって可能である． 58

図 5-12 2 次元 DFT のスペクトル位置関係 4) Band Pass Filter 最後に，時空間画像の直線成分とノイズを除くフィルタ操作である Band Pass Fileter に ついて述べる．Band Pass Filter では，フーリエ逆変換を用いている．フーリエ逆変換は， 簡単に言うと先に述べた FFT の逆の作業であり，スペクトル画像つまり周波数空間のデ ータを二次元画像にする作業である．そしてその変換の際，周波数空間のデータを操作 することで，フィルタとしての役割を果たすことができる．低周波数成分をゼロとし， 高周波成分のみを返すハイパスフィルタ，逆に高周波数成分をゼロとし，低周波成分の みを返すローパスフィルタが存在する．一般的にローパスフィルタは，平滑化フィルタ， ハイパスフィルタは輪郭強調のフィルタとして知られている．そこで本研究では，この 二つのフィルタを組合すことで任意の周波数帯を切り出し，これをフィルタとしている． 各フィルタに関する説明を図 5-13 に示す．ただ実際，解析を行う上ではなるべく元デー タをフルに利用したいといった考えがあるため，ローパスフィルタによるデータカット 及びスムージングは基本的には行わず，Gaussian Filter をフーリエ逆変換した後の画像に 適用することで画像のスムージングを行っている．一方，ハイパスフィルタによる低周 波成分のカットに関しては，雨の影響や，光の照射環境の違いによる時空間画像の色の 濃淡の除去などの効果があり解析精度に大きく影響するため，最低限の範囲(縦方向及び 横方向共に波数で 1，2 個程度)でカットしている．ただ，低周波成分にノイズがなく， 高周波成分にノイズが含まれているようなケースなどが存在すれば，状況に応じてカッ トする周波数の範囲を変化させ，ハイパスフィルタは適用せず，ローパスフィルタを適 用するといったように，時空間画像に適した周波数帯で解析を行っている． 59

図 5-13 フーリエ逆変換を用いたフィルタ 5) 空間周波数とスペクトル STIV 解析において，平均流速は時空間画像における縞パターンの傾きから求めている ため，得られる平均流速の精度は時空間画像の傾斜角の計測精度に大きく依存する．そ のためプログラムでいかに精度よく時空間画像の縞パターンの勾配を得るかが重要で ある．これまでに，時空間画像によってはマニュアルとプログラムで算出角度に差が生 じるケースが確認されている．その中で，時空間画像に対し 2 次元フーリエ変換を用い たフィルタ(以下 Band Pass Filter) を適用することで時空間画像のノイズが除去でき，こ れらの誤差を軽減できるといった報告が原ら 36)によってなされている．画像全体をフー リエ変換する 2 次元フーリエ変換及び空間周波数の概念について以下に説明する． 時間信号と同じように空間的に変化する信号に対しては，空間周波数が定義され，ま 60

たスペクトルによって表現することができる．長さ L を周期とする１次元の正弦波信号 は (5.16) で表される．ただし，画像の最大濃度値を 1 とし，振幅 A は とする． を空間周波数(spatial frequency)といい，長さの単位が cm であれば，1cm 当りの正弦波の 個数(波数)を表す．画像全体の幅を L0 としたとき，これを基準にした空間周波数L0 を基本周波数として定義し，任意の空間周波数をkのように表現する．次数 k は L0 に含まれる正弦波の数を表している．全体の画素数を M，画素間隔を d とすると，L0=Md となる．x = dm( m = 0, 1, … )とすれば，離散化された正弦信号は (5.17) のように表現できる． M =30， L =10 d ， k =3 の場合を図 5-5 に示す． 図 5-14 1 次元画像の濃度分布 式( 5.17 )にオイラーの公式を用いることで次式(5.15)が得られる． 61

(5.18) また図 5-14 の信号に対する振幅スペクトルは，横軸 k で示すと，図 5-15 のように表現 される． 0.5 A 2 A 2 -3 3 k 図 5-15 図 5-14 に対する振幅スペクトル 式( 5.15 )を 2 次元デジタル画像に拡張すると となる． M , N はそれぞれ横方向，縦方向の画素数である． k は横方向空間周波数， l は 縦方向の空間周波数に相当する．図 5-16 に M = N =128， k =2, l =4 の画像例と振幅スペ クトルを示す．横方向，縦方向の明暗の現れる回数は，それぞれ 2 回および 4 回になっ ている．なお縦方向は下向きが正方向である． 図 5-16 2 次元画像の振幅スペクトル（左：原画像，右：振幅スペクトル） また図 5-14 と同じ傾きで周期の違う画像についても，図 5-17 にスペクトル画像を示 す．図 5-17 (a)から(d)は，規則的な正弦波の濃度分布であるのに対して，図 5-17 (e)は(a) から(d)の濃度分布をすべて足し合わせた不規則な濃度分布である．規則的な(a)から(d) のスペクトル画像は，それぞれ( 0,0 )および(± k ，± l )の位置に強いスペクトルが見られ 62

る．足し合わせによってできた(e)に関しては，( 0,0 )および足し合わせた波の持つ k，l がすべてプロットされている．今回のケースであれば l =2k といった数式で表される直 線状にプロットされていることがわかる．この直線がスペクトル画像上における傾きで ある．この角度は，原画像の傾きの 90 度回転した角度となっている．これは時空間画 像に FFT を適応した場合でも同様である．フーリエ変換により，実際の時空間画像は縦 横ともに無限の正弦波の波に分解されるため，今回同様，卓越した周期のライン上によ り多く点がプロットされることとなる．よって，その点の集まりから角度を算出できれ ば，時空間画像の平均的な角度を算出したこととなる． 図 5-17 2 次元画像の周波数スペクトル まとめると，画像は特有のスペクトルを持ち，その位置や方向によって元の画像の周 期性や構成要素の大きさなどが判別できる．つまり，各 STI のスペクトルの位置及び方 向を確認することで，STI のもつ周期性や勾配を検出することができる．そこで，スペ クトルの位置情報を分布化したパワースペクトル分布を用いることで，エネルギーの強 い勾配(STI における主要な勾配)，弱い勾配(細かな波が作り出す勾配)などエネルギーを 考慮した波の周期性や勾配を確認することができる．パワースペクトル Pk,j は，2 次元画 像 fm,n を FFT した際に得られるフーリエ係数 Fk,j を用いて表現すると で表され，空間周波数(k,l)の強さを表している．スペクトルの対象性から，Fk,j の上半分 または下半分で調べることができる． パワースペクトル分布を動径方向と角度方向分布の二通りで確認している．まず動径 63

方向とは，図 5-18(a)に示すように，中心(直流成分)から の距離に存在する環状領域内のパワースペクトルの和を，角度方向分布は，図 5-18(b) に示すように，水平軸から角度 θ の線形領域内のパワースペクトルの和を表している． 動径方向における分布は． 角度方向における分布は という式によって表わされる． 図 5-18 パワースペクトル分布 図 5-19 に実際の時空間画像とそれに対応するパワースペクトル分布例を示す．時空間画 像では，はっきりとした勾配線が見える．これは，河川上を流れる浮子が検査線上を通 ったためである．時空間画像上における勾配線の角度は-15°である．この画像の角度方 向におけるパワースペクトル分布を見ると，-15°におけるエネルギーが卓越しているの が分かる．卓越した角度から算出される流速値は 2.80m/s であり，通常の STIV 解析で得 64

られる値は 2.66m/s であった．この違いは，通常の STIV 解析では卓越した勾配線以外の 角度も検出するため，流速値が若干低下したものによるものと考えられる． 図 5-19 パワースペクトル分布例 5.1.1 ヒストグラム均等化（equalization） また，今回の研究開発では，ヒストグラムの均等化を時空間画像に施すことで，輝度 勾配テンソル法における勾配検出の精度を高めることを発見し，本ソフトウェアに導入 している． ヒストグラム均等化とは，ヒストグラムの累積度数(輝度値０から画素数を累積したも の)のグラフの傾きが一定になるように変換する処理である。この処理によって，コント ラストが悪かったり、明るさが偏っている画像の全体的なバランスを改善することが可 能となる．図 5-14 に時空間画像にヒストグラム均等化を施した例を示す． 時空間画像上では勾配線がぼんやりとしているが，ヒストグラム均等化処理を施した後 は勾配線が明確に表示されている． 65

図 5-20 フーリエ逆変換を用いたフィルタ 5.2 STIV 解析 本研究における STIV 解析とは，5.1 で示した手法を取り入れたプログラムを用いての 流速算出解析である．以下解析手順について説明する．また，STIV ソフトでは動画像を 解析に用いることができる．現在対応している動画フォーマットは avi, mp4, m2ts である． まず解析の前処理として，流量観測を行いたい場所においてカメラキャリブレーショ ンを行うための標定点の測量を行う．そして，画像の標定点を確認できるビデオ画像を 用いてカメラキャリブレーションを行い，幾何補正のためのカメラパラメータ（標定要 素）を求める． 以上の準備が整うと，プログラムでの解析が可能となる．まずプログラムに動画像を 読み込む．そして流速を求めたい位置に検査線を設置する．検査線上の輝度値を動画の フレーム分読み込み時空間画像を作成する．この時，検査線は流下方向に引く必要があ る．次に作成した時空間画像から輝度勾配を求める手順へと移る．本研究ではより精度 良く勾配を検出できるように下処理として時空間画像にガウシアンフィルタ及び二次 元フーリエ変換を用いたフィルタを適用している．そうすることで，撮影画像のコマ落 ちの影響，細かな光の反射などによるノイズ除去が行える 3５), 3６)．そして次に，STI に ヒストグラム均等化処理を施した後，時空間画像全体に複数の検査領域を設置する．こ のとき検査領域の大きさ，および設置間隔はプログラム上自由に設定できるものとなっ ている．そして，その検査領域ごとに算出された角度をヒストグラムで表し，ヒストグ ラムのピークから平均角度を算出する．ヒストグラムでは明瞭な縞パターンの角度に重 みをつけるために，コヒーレンシーの値によって角度に重みをつけて積み重ねたものを 用いている．以上の手順で時空間画像のもつ平均角度を求め，式( 5.3 )より流速を算出す 66

る． ここでは鬼怒川を対象とした，STIV 解析例を示す．図 5-15 は鬼怒川の右岸から取得 したビデオ画像を用いて STIV 解析を行った例である．主流方向は画面の左から右であ り，浮子が 3 つ流れている．検査線は幾何補正画像上で主流方向に平行に引いた線を， 斜め画像上に投影することで設置している．表 5-1 に STIV 解析結果を示す．また，図 5-16 に解析によって作成された時空間画像(STI)と図 5-17 に 2 次元フーリエ変換によっ て作成された時空間画像（FFTSTI）を示す．時空間画像から，勾配線が明らかである事 が分かる．しかしながら，⑥，⑨，⑩の時空間画像では縦方向のノイズも見える．フー リエ変換をかけることによって，時空間画像から縦方向のノイズを除去できている．ま た，時空間画像①，⑥，⑦に見えている赤い斜線は浮子の軌跡を表わしている．また， 各時空間画像の下部にある，黒い領域は avi ファイルの不具合により生じたものである が，フーリエ変換によって除去されている．流速値に関しては，左岸（画面奥）に近づ くほど流速が低下しているが，岸から遠い画面手間に近い検査線上①～⑤では流速値は 2.2m/s～2.6m/s となっている．ここまでの値は，STIV の自動解析結果に基づいている． 解析結果の精度を確かめるために，マニュアル解析による結果と比較した．結果は図 5-21 に示す．自動解析による結果はマニュアル解析によるものとほぼ一致しているため，精 度に問題はないと考えられる． 図 5-21 検査線設置画像（鬼怒川右岸，左；斜め画像，右；幾何補正画像） 67

表 5-1 STIV 解析結果（左；流速値，右；動画像情報と解析条件） 図 5-22 図 5-17 時空間画像（鬼怒川右岸） FFT 変換後の時空間画像（鬼怒川右岸） 68

① ⑩ 図 5-18 自動解析結果とマニュアル解析の比較（鬼怒川右岸） 次に，鬼怒川の左岸から撮影されたビデオ画像を用いて STIV 解析を行った例を示す． 図 5-19 に検査線を設置した斜め画像と幾何補正画像を示す．流れは斜め画像の右から左 に流れており，浮子も二つほど流れている．またこの時間帯の河川表面は水面の凹凸が 激しく，局所的に早い流速となっている場所がある．表 5-2 に STIV 解析結果を示す． また，図 5-20 に解析によって作成された時空間画像(STI)と２次元フーリエ変換によっ て作成された時空間画像（FFTSTI）を示す．時空間画像上で勾配線は明確に写っている． しかしながら，⑦，⑩の時空間画像では縦方向の赤色の垂線が見える．検査線⑦で見え る赤い線は，河川上で流速値を測定している ADCP が時空間画像で映っているためであ る．検査線⑩では，検査線に標定点が触れてしまっているために直線が映っている．こ れらの直線はフーリエ変換を用いることによって除去できている．また，時空間画像全 体において，河川表面が太陽の反射によって明るくなっている．この時空間画像上の輝 度の急激な変化もフーリエ変換と輝度ヒストグラム均等化処理によって解消される．流 速値に関しては，右岸（画面奥）に近づくほど流速が低下している．また，岸から遠い 画面手間に近い検査線上流速値は 3.5m/s～4.5m/s となっている．これは，水面の凹凸に よって流速が局所的に早くなっているためである．水面の凹凸は時空間画像⑫や⑭にお いて，画像上の勾配線が急激に変化している部分にも表れている．ただ，検査線⑳にお ける流速値は，時空間画像の輝度分布が平坦であることから異常値であると考えられる． 同じように，検査線⑯～⑲においても時空間画像上の勾配分布が平坦である．実際，こ の領域における 流速値をマニュアル解析で求めると 2.8～3.5m/s である（図）．他の領域 に関しては自動解析とマニュアル解析はほぼ同じ値を記録している． 69

図-5-19 検査線設置画像（鬼怒川左岸，左；斜め画像，右；幾何補正画像） 表 5-2 STIV 解析結果（左；流速値，右；動画像情報と解析条件） 図 5-20 時空間画像（鬼怒川左岸） 70

図 5-21 FFT 変換後の時空間画像（鬼怒川左岸） ⑳ ① 図 5-22 自動解析結果とマニュアル解析の比較（鬼怒川左岸） 5.3 まとめ ① STIV`について STIV は，検査線上の輝度を時間方向に並べた時空間画像における勾配によって流速 を測定する．自動解析では輝度勾配テンソル法と統計処理によって，勾配線の平均を算 出する．STIV では，エッジ検出とノイズ除去に 2 次元フーリエ変換と Gaussian フィル ターを用いる．さらに，ヒストグラムの均等化処理を施すことによって，精度を向上さ 71

せる工夫を行った． ② STIV`自動解析の制度について 自動解析とマニュアル解析結果を比較した結果，自動解析の精度は高いことがわかっ た． 72

第6章 各々の開発ソフトの操作方法 6.1 キャリブレーションソフトの操作 以下に開発してきたキャリブレーションソフトについての使用手順を簡単に説明する． ① キャリブレーションソフト‘CameraCalib.exe’をクリックして起動する 図 6-1 CameraCalib.exe 起動 ② CameraCalib.exe が起動し図 2.61-2 に示すような UI が表示される． 図 6-2 CameraCalib.exe の表示画面 ③ 左上のメニューバーにある File(F)をクリック（Alt+F）してツールバーを展開（図 673

３），Open（静止画）か Open（動画）をクリックして，ダイアログを開き，標定点 を写した斜め画像を開く．この時，標定要素画像は sample.bmp という名前である． もしくは Ctrl+P で図 2.3.5-4 にしめすダイアログが開く． 図 6-3 標定点を写した画像（左）を読みだす為に， フォーム左上にある File（F）からバー（右）を開く ダイアログは開いたら，標定点の映っている画像を選択する．デフォルトでは sample.bmp であるが，幾何補正画像が映っている画像が別名のファイル，あるいは別形 式の静止画ファイルであっても，読み込むことができ，自動的に sample.bmp が作成され る．また，動画を選択することで，動画のフレームを取りだして標定点画像 sample.bmp を作成する事ができる．この時，ソフトの作業ディレクトリ（他のファイルを読み取っ たり，出力するディレクトリ）は読み込んだ画像ファイルのあるディレクトリとなる． 図 6-4 ダイアログを開いて標定点の映っている画像 sample.bmp を開く 74

④ 標定点画像を開くと，図 2.6.1-5 のように UI は表示される． 図 6-5 標定点画像を開いた CameraCalib.exe 標定点画像をクリックしてドラッグすると画像を平行移動させる事ができる（図 2.6.1-6）． 図 6-6 画像の平行移動 75

図 6-7 画像の拡大縮小 また，画面をクリックしてマウスホイールをする事によって，画像の拡大縮小ができる． 画面上で右久利幾をする事によって，図 2.3.5-6 にしめすようにポップアップで標定点の 物理座標を入力する UI が現れる． 図 6-8 物理座標入力 UI UI からの入力のほかに，points.csv というファイルを用意する方法がある．このファイ ルが存在していると自動的に標定点の物理座標と写真座標が読み込まれる．points.csv の 内部のデータは，次の表のような構造をしている． 表 6-1 points.csv の内部データ 76

ここで，Num は標定点の番号，Wx，Wy，Wz は物理座標，x, y は写真座標，Flag は標定 点を，-1 の時は使用し，0 の時は使用しない． 図 6-9 キャリブレーション直前のディレクトリ構成 ⑤ キャリブレーションが成功すると，図 6-10 の様に使用した標定点から水面に向かっ て緑色の垂線が描かれ，幾何補正画像が出力される．また，精度を確かめるために 標定点の中心に青丸の中心に，赤丸が入っているかを確かめる．計算に使用してい ない標定点に関しては，青丸の中心に黄色の丸が入っているかを確かめる（図 6-11）． 図 6-10 キャリブレーション成功時の UI 画面 77

図 6-11 キャリブレーションによる標定点座標の再計算結果の表示 ⑥ キャリブレーションが成功したら，左上の File(F)をクリックし，Save をクリック， もしくは Ctrl+S によってキャリブレーションデータを保存する．キャリブレーショ ンデータを保存すると savefile.dat と Area.dat が出力される．savefile.dat と Area.dat は それぞれ,キャリブレーションによって得られたカメラパラメータと幾何補正画像 の 物 理 範 囲 で あ る ． こ の 二 つ フ ァ イ ル は 後 に 説 明 す る STIV に 用 い る ． Out_oblique.bmp と out_rectified.bmp はそれぞれキャリブレーション再計算結果を表 示した，斜め画像と幾何補正画像である． 図 6-12 キャリブレーションによるファイルの出力 6.2 Float-PTV の使い方 ① 計測ソフトの使い方． 画像によって河川の流速を図るために，Float-PTV.exe を起動させる．ソフト起動させる 78

前に同じフォルダに，以下のファイルを入れておく（図 6-13）． 図 6-13 Float-PTV を起動させるフォルダ ここで savefile.dat, Area.dat は先ほどの Camera.exe（別に開発されたカメラキャリブレー ションソフト） で作成したファイル，対象動画.avi は解析対象の画像（今回は動画）を 表わしている． ② Float-PTV.exe を起動した時に savefile.dat といった先ほどの製作データが読み込まれ れば，“Load savefile.dat”と表示される． 図 6-14 Measure.exe 起動時の様子 ③ Measure.exe を起動させると，下の様なウィンドウが３つ開く． 79

図 6-15 Measure.exe 起動 ④ 次に，Composite ボタンを押す． 図 6-16 Composite ボタンを押す 80

“Frame Interval”の数字を調整して、何フレーム毎に画像を合成していくかを決定する． 例えば、ここで使用する avi ファイルはフレーム率が 29FPS なので、Frame Interval を 10 にすると、約 0.3 秒ごとに動画を重ねて合成していくことになる．ここでは Frame Interval = 10 と設定する。Frame Interval を設定したら、Composite ボタンをクリックする。成功 すると“ALL SET”と表示されるので OK ボタンをクリックする。 図 6-17 ALL SET ⑤ しばらくすると，下図のように合成された写真が作成される． Flow Flow 図 6-18 多重合成画像が作成された ⑥ 流速ベクトルの解析を開始する。”Start for Calc”をクリックする。 81

図 6-19 Start for Calc をクリック “Now, Point on FLOATS!”と表示されるので、OK ボタンをクリックする。 図 6-20 Now, Point on FLOATS ! 浮子投下ライン一本ずつに対して、画面上の浮子（画面上の赤い点）の根元をクリッ クしてトレースする。クリックしたポイントに緑色の小さい丸が追加される。これを画 面の最後（この場合は右端）まで飛ばさずにクリックする。全部の点をクリックすると 下の画面になる． Flow 図 6-21 多重合成画像上の浮子の押し方 82

⑦ Calc ボタンをクリックする． 図 6-22 Calc ボタンをクリックする 次の画像が表示される（図 6-23）． 図 6-23 ベクトル図 ⑧ 使用しているフォルダ内に“Point_VectorDATA.csv”, “VelocityImage_Geo.jpg”,”VelocityImage_Obl.jpg”の各ファイルが作成されているこ とを確認． 83

図 6-24 出力ファイルの確認 以上で終了である． 2.5.3. STIV の操作方法 ①. STIV.exe をクリックして起動． 図 6-25 STIV.exe を起動する ②. ソフトを起動し，メニューにある画像ファイル読み込み→動画像の読み込みを選択 84

し，解析対象の動画を選択する．この時ソフトの作業ディレクトリは動画が存在す るフォルダとなる．この時，動画のあるフォルダにはキャリブレーションソフトで 作成した，savefile.dat と Area.dat が含まれていなければならない．動画を開くと，図 6-25 のような UI が表示される．左に斜め画像，右に幾何補正画像が表示される．両 方の画像上でクリック後．マウスホイールで拡大縮小，左クリックを押しながらド ラッグで平行移動ができる． 図 6-26 STIV.exe の UI ③. UI の右上にあるコントロールパネルで検査線の本数と水平方向と鉛直方向における 検査線の間隔を調整できる． 図 6.-27 検査線のコントロールパネル ④. UI の右上にあるコントロールパネルで検査線の本数と水平方向と鉛直方向における 検査線の間隔を調整できる．幾何補正画像上で右クリック+マウスドラッグで検査線 を引くことが出来る．幾何補正画像上で線を引くと，斜め画像上に検査線を再現す ることができる．また，斜め画像上でも右クリック+マウスドラッグで検査線を引く ことができ，幾何補正画像上に反映される． 85

図 6-28 検査線は斜め画像幾何補正画像，両方から引ける ⑤. 検査線が引かれると，下のセルに検査線の情報が表示される．ここで野情報とは検 査線の番号，物理長さ，時間間隔，始点と終点における物理座標と写真座標である． 検査線の位置が適正であることを確認して，UI 上部にある“解析開始“ボタンをク リックする．解析が開始すると，下のセルに順番にそれぞれの検査線における流速 が出力される． 図 6-29 86 解析ボタンを押す

⑥. 解析が終了すると，二つのポップアップが順番に現れる．STI の個別解析では図 6-30 に示すように，各 STI 画像での輝度勾配テンソル法結果を確認しながら解析するこ とが出来る．STI`のフーリエ解析では各 STI を個別に 2 次元フーリエ変換を行うプロ セスを確認できる．それぞれのメッセージにおいてはい(Y)を選択すると，図 6-31， 図 6-32 の UI がそれぞれ出現する．この UI で，自動解析によって得られた流速値を 更新することが出来る．また，ここで個別の解析を選択しなくとも，動画のフォル ダに STIV の解析データと時空間画像，フーリエ変換された画像を収納したフォルダ があれば，メニューの解析(A)から，UI を呼び出すことが出来る． 図 6-30 解析終了後に出現するメッセージボックス 図 6-31 STI 個別解析 UI 87

図 6-32 STI フーリエ解析 UI ⑦. 解 析 が 完 了 す れ ば ， 自 動 的 に “ Line_Oblique.bmp ” と ”Line_GeoRectified.bmp” ， STIV_DATA.csv が 出 力 さ れ る ． Line_Oblique.bmp は 斜 め 画 像 に ， Line_GeoRectified.bmp は幾何補正画像にそれぞれ検査線を引いた図を保存したデー タである．STIV_DATA.csv は STIV 解析データを保存したデータである．また，解 析で出力された時空間画像と 2 次元フーリエ変換後の画像はそれぞれ STI フォルダ と FFTSTI フォルダに収納されている． 88

図 6-33 STIV による出力データ 89

7. 結論 本研究では，画像を用いた河川流速計側の実用化を目的とし，これまでの画像解析技 術の応用と利便性を重視したソフトウェアの開発と精度について述べてきた．まず，画 像解析技術に必要なカメラキャリブレーションについて述べた． 次に，神戸大学と土 木研究所(ICHARM)との共同研究で開発した，画像浮子解析ソフト’Float-PTV’につい て述べる．本ソフトでは，河川漂流物を追跡することで河川表面流速を簡単に計測でき る．しかしながら，実際の河川洪水時には漂流物が常に画面上に表示されているという ことはなく，さらに流れにそっているとは限らない．そのため水面上を流れる波紋によ って，表面流速を測定する手法が必要となる．先に述べたように，STIV では検査線上の 輝度を時間方向に並べて作成した時空間画像を用いて流速を測定する．このような理由 から，最後に新しく開発された新しい STIV での自動改易結果の精度向上の工夫につい て述べた． 7.1 カメラキャリブレーション ① カメラキャリブレーションとは 複数の物理座標と写真座標の組によって，カメラパラメータ（標定要素）を算出する 方法．カメラパラメータを明らかにする事によって，空間幾何学に基づいた画像の校正 を行え，解析に用いる事ができる． ② カメラキャリブレーションのレンズ歪補正の効果は レンズ歪の補正を行うことによって，より正確な幾何補正を作成する事ができる．し かしながら，実用的にはレンズ歪を用いずとも実用的な幾何補正画像の作成は可能であ る． ③ カメラキャリブレーションに必要な標定点は（カメラの座標は既知である場合） 最低，三点であり，標定点の配置にも左右される． ④ 簡易型カメラキャリブレーションについて 0.1Degree の精度での角度とカメラの位置が既知であるなら，一点の標定点（手前の岸 までの長さ）によってカメラキャリブレーションを成功させる事が可能． 7.2 Float-PTV ① Float-PTV は浮子を追跡することで，流速を計測するソフトウェアである． ② 浮子追跡には，多重合成画像を用いる． ③ 多重合成によって浮子追跡は容易になる．しかし，河川の状況によって，浮子が見 90

えなくなる場合もある．そのため河川表面の波紋パターンを追跡する画像解析技術が必 要となる． 7.3 STIV ① STIV`について STIV は，検査線上の輝度を時間方向に並べた時空間画像における勾配によって流速 を測定する．自動解析では輝度勾配テンソル法と統計処理によって，勾配線の平均を算 出する．STIV では，エッジ検出とノイズ除去に 2 次元フーリエ変換と Gaussian フィ ルターを用いる．さらに，ヒストグラムの均等化処理を施すことによって，精度を向上 させる工夫を行った． ② STIV`自動解析の制度について 自動解析とマニュアル解析結果を比較した結果，自動解析の精度は高いことがわかっ た． 91

謝辞 本研究を遂行するにあたり，終始一貫して適切な御指導および御鞭撻を賜った神戸大 学大学院工学研究科藤田一郎教授に深く感謝の意を表します．本研究，引いては，大学 院生活を通して自分の目指す理想の技術者の姿を藤田先生から学ばせて頂きました．常 に技術者としての結果を出す為に，学び続ける事，アンテナを張り続ける事，枠を狭く 取らないこと，こういったことを実際に我々学生に対して，自身の姿勢によって指導し てくださりました．私自身，学問的な知識や単なるプログラミングなどの技術のみでは なく，土木の技術者としてあるべき姿を学ぶ事が出来ました．また，画像解析技術に関 する研究・開発を通して，様々な体験をすることをさせていただいた事も生涯忘れるこ とはできません．藤田先生の指導のもとでは，技術や知識を学べば学ぶほど，それを生 かす為のチャンスを与えていただけ，応援して頂けるので，自分もやりがいがあり，ど んどん学んでいく事が出来ました． また大学院生活に際し，さらに，有益なご助言をいただきました，前田浩之技術職員， 広島大学の椿涼太准教授に厚く御礼申し上げます．前田浩之技術職員には，TA からゼ ミ旅行，飲み会まで本当にお世話になりました．ドイツ料理の御飯を一緒に食べに連れ て行って下さったこともありました．その一方で実験の時や観測の時にはてきぱきと指 示を下さる前田職員の姿は人一倍頼もしく思えました．広島大学の椿涼太先生に関しま しては，おかげさまで，簡易型カメラキャリブレーションソフトを完成する事ができた といっても過言ではなく，この場で感謝の意を述べさせていただきます． 最後に，CameraCalib や Float-PTV 開発の機会を与えてくれた土木研究所の萬矢氏と本 永氏，そして，平素より多大なご助力をいただきました研究室 の諸氏に深く感謝いたします． 参考文献 1) 国土交通省河川局．水害対策を考える http://www.mlit.go.jp/river/pamphlet_jirei/bousai/saigai/kiroku/suigai/suigai.html，2012． 2) 孫子，金谷治（訳）．新訂孫子．岩波文庫，2011． 3) 四国水門観測研究会．高水観測の用語集． http://www.skr.mlit.go.jp/kasen/mizu/yougo/sect051.htm，2012. 4) 土木研究所他．共同研究報告書第 291 号「非接触型流速計側法の開発」．2003． 92

• http://www.mlit.go.jp/river/pamphlet_jirei/bousai/saigai/kiroku/suigai/suigai.html
• http://www.skr.mlit.go.jp/kasen/mizu/yougo/sect051.htm

5) 佐藤慶太，二瓶泰雄，木水啓，飯田祐介．洪水流観測への高解像度超音波ドップラ ー流速分布計の適用 江戸川を例にして．水工学論文集，第 48 巻，pp763-768，2004． 6) 山口高志，新里邦生．電波流速計による洪水流量観測．土木学会論文集， No.497/Ⅱ-28.pp.41-50，1994. 7) Wang, C.J., Wen, B.Y., Ma, Z.G.., Yan, W.D. and Huang, X.J.: Measurement of River Surface Currents With UHF FMCW Radar Systems, J. of electromagnetic waves and applications, 21(3), 2007 , pp. 375-386, 2007. 8) 土木研究所．発表資料 非接触型流速計． http://www.pwri.go.jp/jpn/seika/pdf/newtech/ryusoku.pdf，2012． 9) 藤田一郎．実河川を対象とした画像計測技術，土木学会水工学委員会・海岸工学委 員会，水工学シリーズ0.3-A-2，2003． 10) 小林範之，金目達弥，藤田一郎．PIV による洪水時の河川流量観測装置の開発，河 川技術論文集，第8巻，pp.455-458，2002. 11) 藤田一郎，河村三郎．ビデオ画像解析による河川表面流計測の試み，水工学論文集， 第38 巻，pp.733-738，1994. 12) 藤田一郎，綾史郎，河村三郎，中島三郎：ビデオ画像を利用した洪水流量の新計測 法，土木学会第49 回年次学術講演会概要集２-A，pp.382-383，1994. 13) Fox, J.F. and Patrick, A.: Large-scale eddies measured with large scale particle image velocimetry, Flow measurement and Instrumentation, 19, pp.283-291, 2008. 14) Kantoush, S.A., De Cesare, G., Boillat, J.L. and Schleiss, A.J.: Flow field investigation in a rectangular shallow reservoir using UVP, LSPIV and numerical modelling, Flow measurement and Instrumentation, 19, pp.139-144, 2008. 15) Meselhe, E.A., Peeva, T. and Muste, M.: Large scale image velocimetry for low velocity and shallow water flows, Journal of Hydraulic Engineering,130(9), pp.937-940, 2004. 16) 藤田一郎．河川表面流のLSPIV による計測と高精度化．平成12年度土木学会関西支 部年講，2001. 17) 藤田一郎，中島丈晴．実河川流計測におけるLSPIVの汎用化と水制間流れへの適用． 水工学論文集，第44巻， pp.443-448，2000. 18) 藤田一郎．非接触型流速計測法を用いた実河川流の計測と問題点．ながれ，第26巻， 5-12, 2007. 19) Fujita, I., Muste, M. and Kruger, A.：Large-scale particle image velocimetry for flow analysis in hydraulic engineering applications ， Journal of Hydraulic Research，Vol.36, No.3, pp.397-414，1998. 93

• http://www.pwri.go.jp/jpn/seika/pdf/newtech/ryusoku.pdf

20) Muste, M., Fujita, I., and Hauet, A.: Large-scale particle image velocimetry for measurements in riverine environments, Water Resources Research, 44, W00D19, doi:10.1029/2008WR006950, 2008. 21) Hauet, A., Creutin, J.-D. and Belleudy, P.: Sensitivity Study of Large-Scale Particle Image Velocimetry Measurement of River Discharge using Numerical Simulations, Journal of Hydrology, 349(1-2), 178-190, doi:10.1016/j.jhydrol. 2007. 10.062., 2008. 22) Fujita, I., Watanabe, H. and Tsubaki, R.: Development of a non-intrusive and efficient flow monitoring technique: The space time image velocimetry (STIV), International Journal of River Basin Management, Vol.5, No.2, pp.105-114, 2007. 23) 綾史郎，藤田一郎，柳生光彦．画像解析を用いた河川洪水時の流れの観測．水工学 論文集，第39巻，pp.447-452，1995. 24) 藤田一郎．トレーサーを利用した実河川水制周辺流れのビデオ解析．水工学論文集， 第42巻，pp.505-510，1998. 25) 藤田一郎．トレーサーを利用した実河川水制周辺流れのビデオ解析，水工学論文集， 第42 巻，pp.505-510，1998. 26) 藤田一郎，椿涼太．時空間画像を利用した河川表面波紋の移流速度計測，河川技術 論文集，第9巻，2003. 27) 藤田一郎，椿涼太．小俯角のビデオ画像に対応した河川表面流計測手法の開発，河 川技術論文集，第7巻，2001. 28) 藤田一郎，椿涼太．時空間濃度勾配法による主流方向表面流速分布の現地計測．水 工学論文集，第46巻，pp.821-826，2002. 29) 冨尾恒一．準リアルタイム河川流計測手法（STIV）の実用化に関する研究．神戸大 学大学院自然科学研究科博士前期課程修士論文，2006. 30) 渡辺英樹．検査線上の時空間画像を用いた準リアルタイム河川流計測システム．神 戸大学大学院自然科学研究科博士前期課程修士論文，2005. 31) 堤志帆．STIV による夜間および悪天候時の河川洪水流計測．神戸大学工学部建設学 科卒業研究，2008.2 32) 藤田一郎，安藤敬済，堤 志帆，岡部健士．STIV による劣悪な撮影条件での河川洪 水流計測．水工学論文集, 53巻, pp.1003-1008, 2009. 33) 原浩気，藤田一郎．時空間画像を用いた河川表面流解析における二次元高速フーリ エ変換の適用．水工学論文集, 54 巻, pp.1105-1110, 2010. 34) 原浩気．STIV法を用いた準リアルタイム河川表面流計測手法の確立と流速計測精度 に関する研究．神戸大学大学院自然科学研究科博士前期課程修士論文，2011.2 94

35) 気象庁．気候変動監視レポート2011． http://www.data.kishou.go.jp/climate/cpdinfo/monitor/index.html，2012. 36) 国土交通省．平成23年度国土交通白書． http://www.mlit.go.jp/hakusyo/mlit/h23/index.html，2012． 37) 酒井雄弘，二瓶泰雄．ADCPを用いた中小河川の流量計測法に関する検討．水工学論 文集，2007. 38) 武藤 裕則，北村 耕一，馬場 康之，中川 一．ADCPを用いた水制域における流速分 布計測，水工学論文集，第49巻，pp637-642，2005. 39) 酒井雄弘，二瓶泰雄．ADCPデータに基づく大河川洪水流の更正係数に関する検討， 水工学論文集，第49巻，pp637-642，2005. 40) 矢吹 太朗．文法からはじめるプログラミング言語 Microsoft Visual C++入門 (マイク ロソフト公式解説書) ． 日経 BP ソフトプレス，2008. 41) 山田 祥寛．プログラムを作ろう！ Microsoft Visual C++ 2008 Express Edition 入門 (マイクロソフト公式解説書) ．日経 BP ソフトプレス，2008. 42) 増田 智明．ひと目でわかる Microsoft Visual C++ 2008 アプリケーション開発入門 (マイクロソフト公式解説書) ．日経 BP ソフトプレス，2008. 43) Dr. Gary Rost Bradski, Learning OpenCV: Computer Vision with the OpenCV Library, Oreilly & Associates Inc, 2008 44) Robert Laganiere, OpenCV 2 Computer Vision Application Programming Cookbook: Over 50 Recipes to Master This Library of Programming Functions for Real-time Computer Vision. Packt Publishing, 2011. 45) 谷尻 豊寿．Essential OpenCV Programming ―With Visual C++ 2008．カットシステム， 2009. 46) 北山 洋幸，OpenCV で始める簡単動画プログラミング，カットシステム，2010． 47) 昌達 慶仁，詳解 画像処理プログラミング C 言語で実装する画像処理アルゴリズム のすべて，ソフトバンククリエイティブ，2008. 48) 酒井 幸市．画像処理とパターン認識入門 - 基礎から VC#/VC++ .NET によるプロジ ェクト作成まで．森北出版，2006． 49) 日本写真測量学会．写真による三次元測定―応用写真測量．共立出版，2008. 50) K.Todd Holland, Robert A. Holman, Thomas C. Lippmann, John Stanley, Associate Member IEEE, and Nathaniel Plant. Practical Use of Video Imagery in Nearshore Oceanographic Field Stydies.. IEEE Journal of oceanic engineering, Vol22, No.1, 1997. 51) Alexandre Hauet, Jean-Dominique Creutin, Philippe Belleudy. Sensitivity study of 95

• http://www.data.kishou.go.jp/climate/cpdinfo/monitor/index.html
• http://www.mlit.go.jp/hakusyo/mlit/h23/index.html

large-scale particle image velocimetry measurement of river discharge using numerical simulation. Journal of Hydrology349, 178-190, 2008. 52) I. FUJITA, Y. KOSAKA, M.HONDA, A. YOROZUYA. TRACKING OF RIVER SURFACE FEATURES BY SPACE TIME IMAGEING.. 15th International Symposium on Flow Visualization, 2012. 96