1 序論
これまでの境界層での熱収支の観測では一点での観測が行われてきた。しかし近年、一点での観測は熱輸送などが正確に測定できていないのではないかとの指摘がある。
2002年5月の水文・水資源学会誌の神田らの論文「LESによる熱収支インバランス問題に対する検討」によれば、LES(Large Eddy Simulation)を使ったシミュレーションによって大気境界層内に形成される「大規模組織構造対流」が熱輸送において重要であり、しかもその対流は準定常的に存在することが述べられている。
図1-1はその解析により求められた、ある条件における縦4km、横4kmの高度100mにおける時間平均鉛直流速の空間分布である。この図では白色部分が下降流、等高線が密に描かれている部分が上昇流を表している。この図を見ると上昇流が下降流に比べて小さく、地表付近の場の不均一性を形成している恐れがある。つまりこの大規模組織構造対流が現実に存在すると、従来の点観測では境界層における熱収支の全貌を捕らえることが困難である。
例えば図1-1A点で観測していたとすると移動速度の遅い大規模組織構造対流が移動するまで下降流のみ観測していることとなる。そして上昇気流が下降流域に比べ非常に小さいため下降流域で観測している可能性が高くなる。これはあくまでシミュレーション結果であり、その存在が確認されたわけではないが、その存在の可能性は否定することはできない。したがって大規模組織構造対流の観測を正確に行うのであれば点ではなく線での観測が必要であり、その観測を実現するのに地上からのリモートセンシングが有効である。
図1-1 大規模組織構造対流
現在存在する風向風速の観測機器としてレーダ(Radar)とソーダ(sodar)がある。レーダは電波を使った観測機器であり、降雨の観測をする降雨レーダ、機器上空の風向風速を観測するウィンドプロファイラーなどがある。ソーダとは音波を発信して観測機器上空の風向風速を観測するものであるが現在ソーダには鉛直上向きに音波を発信するものしかない。その理由は、水平方向に音波を発信すると地表面からの反射波など、多くのノイズが入ってしまい測定が難しいからである。
レーダとソーダの違いは発する波が電波であるか音波であるかの違いである。電波の速度は光の速度と同じく約30万km/sであり、音波の速度は約300m/sである。同じパルス幅では音波のほうが移動距離が小さいためソーダのほうが距離分解能が良い。神田らのLESによるシミュレーションで顕された図1は100mオーダーであり、その大規模組織構造対流の観測を行うにはソーダによる地上からのリモートセンシングが適している。
2002年度岐阜大学同研究室の加藤氏により低コストでの水平方向観測用ドップラーソーダの試作が行われ、これによりパルス音波発信と録音まで可能になった。そこで、このシステムに解析プログラムの導入を行った後、ハードウェアの改良を行い、この装置の完成を目指す。
2. ドップラーソーダ
2.1 概要
ソーダ(Sodar)とはSound detection and ranging の略で、ドップラーソーダとは観測機からある周波数の音波をパルスで発信し、観測点の空気からの反射波を捉えてそのドップラーシフトを観測することによって、発信源から離れた場所の風速を観測するリモートセンシング機器である。
パラボラを利用したドップラーソーダは、アンテナの内側のカーブが放物線の形をしている。放物線というのは準線と焦点から等しい距離にある点の集合であり双曲線となるが焦点の位置をパラボラの正面となる位置からずらして設計、使用されることが多い。正面からの波はパラボラに反射し位相がそろった状態で焦点に集まるため、パワーゲインを期待できる。また、正面以外からの波は焦点において位相がずれるため、ノイズは互いに打ち消しあいパワーは小さくなる。
2.2 ドップラー効果の関係式
発信音波の周波数を(Hz)、音速を(m/s)、観測点の風速(音波の発信方向を正とする)を(m/s)、発信音波の周波数と受信音波の周波数の差をとするとその関係の1次近似式は以下で表される。(1)
(1)
またパラボラの特性としての半値幅は、発信パルス波長を、パラボラの直径をとして
(2)
と表される(1)。また、気温から音速を算出するときは音速を(m/s)、気温を(摂氏温度)として以下の関係式を用いる。(6)
(3)
観測の手順は次のとおりである。
@スピーカーからある周波数の音波を発信すると同時に受信側PCで録音を開始する。
Aその反射波を受信する。
B観測時の気温から音速を算出し、音速と観測点までの距離から何秒後のデータを見るのか判断する。
C録音した反射波の受信波形からスペクトル解析により周波数を求め、ドップラーシフトにより観測点の風速を求める。
3 音波発信、受信システム
3.1 システム概要
3.1.1 発信側システム
@発信側PCからパルス波を発信する。
APCから発信された電気信号をパワーアンプで増幅する。
B増幅された電気信号がスピーカより出力され音波となる。
Cパラボラにより指向性が高いパルス波がパラボラ正面に発信される。
3.1.2 受信側システム
@パラボラにより正面に垂直に入ってきた音波のみが位相がそろってマイクに到達する。
Aマイクで取り込んだ音波を電気信号にし、マイクアンプに送られる。
B送れられてきた電気信号をマイクアンプにより増幅する。
C増幅された電気信号は受信側PCに取り込まれ、録音される。
3.2 発信パルス音波の設計
発信パルスは加藤氏(2002岐阜大学同研究室卒業研究)により作成された「周波数、パルス長、再生周期」を指定してリニアPCMデータを作成するプログラムにより用意する。使用する発信パルス音波は周波数を小さくしすぎると波長が大きくなり、距離分解能が悪くなってしまう。逆に周波数を大きくしすぎるとAD変換した際の一波分のデータ数が少なくなってしまうと同時に、人間の可聴域にあわせて作られているマイク、スピーカで正常に発信できなくなる恐れがある。また、反射波を受信するにあたり、観測点からの反射音と発信パルスを区別するには、発信パルスの終点より後に反射波の始点が受信できなければならない。よって距離分解能は音速を、発信パルス幅をとして次の式で表される。
(4)
測定可能な風速の最小単位は周波数分解能を、風速分解能をとおくと(1)式より次式が導かれる。
(5)
発信パルス周波数5000Hz、発信パルス幅0.1秒の場合、周波数分解能10Hz、音速=337㎧、発信パルス周波数=5000を式に代入すると=0.337㎧となり、風速分解能は約0.3㎧である。
3.3 使用機器
今回パラボラおよびシステムのハードウェアは加藤氏(2002)により作成されたものを使用した。ソフトウェアについては自作のものを用意した。
3.4 WAVEデータの再生と録音、解析プログラム
パルス音波の発信と受信にあたりPCで波形データの再生・録音が必要となる。今回はwindowsで一般に使用されているWAVE形式で再生と録音を行うプログラムを用意した。
3.4.1 WAVEデータの再生
WAVEデータの再生には、指定したWAVEファイルをplayボタンにより再生するというプログラムを用意した。(3)
WAVE再生は波形データを再生デバイスに書き込むことであり、時間の経過を考えなければ、WAVE再生もデータ転送の一種と考えることができる。しかしファイルへの書き込みなど通常のデータ転送とは異なり、WAVE再生には時間の経過というサウンドプログラム特有の問題がある。実際このアルゴリズムでは、再生が終わるまでストップボタンが押せないという問題があるため、デバイスが波形データを再生し終えるまで待つのではなく、デバイスからの終了メッセージに応答するというスタイルで実装する。
3.4.2 WAVEデータの録音
WAVEデータの録音には、指定したファイルにWAVEデータをRecボタンにより録音開始し、Stopボタンで録音終了するプログラムを用意した。(3)
録音は、オーディオ入力に入ったオーディオ信号を最終的にアプリケーション側のメモリに取り込む作業であり、録音デバイスはオーディオ信号を決められたサンプリング周波数でAD変換する。AD変換されたWAVEデータはバッファへ蓄積され、アプリケーションはデバイスドライバーと連携しながらそのデータをアプリケーション側に転送する。
3.4.3 解析プログラム
録音データから各周波数のパワーの表記を行う為にFFT(Fast Fourier Transform)プログラムを用意した(4)。今回使用するスピーカーから発信される音波の音圧レベルやパラボラの音の反射率や透過率など明らかではなくFFT実行により得られた数値そのものに意味はないので、そのまま振幅の二乗スペクトルをパワースペクトルと定義する。FFTは2n個のデータにより実行するが、1ファイル分のデータ数が多くなると時間の区切りが大きくなり、細かい距離の測定が不可能になる。逆にデータ数を少なくすればFFTを実行する為のデータが少なくなるので細かい周波数の測定が不可能になる。よって、実験目的に応じたデータ数に変えて実行する。
4 実験
4.1 壁面反射実験
4.1.1 概要
このシステムによりパルス音波の発信と受信が可能であるか試験する。
また、実際の水平方向観測では反射音の録音された時間から距離を算定することが必要となる。そこで、壁に向けての試験による録音データから発信から受信までの時間を計測することで、このシステムの距離を測る能力を試験する。
4.1.2 試験方法
岐阜大学農学部A棟外壁の南側を利用する。そして壁から巻尺で50mを測定し、その位置にパラボラ2台を据え付け、壁に垂直に向ける。次にパルス音波発信側のPCで発信音波を5000Hz、パルス幅0.01秒に設定する。受信側のPCで録音を開始し、音波を発信する。
実験により得られた録音データの解析を行う。
・録音データの解析について
発信パルス音波が0.01秒(約3.4m)であり、FFTを実行する範囲よりも小さいと発信パルス音波を特定できない可能性があるのでデータ数を256(28)個づつに分ける。また、この際の1ファイル分の時間は0.0058秒である。
この録音データに対してFFTを実行すると172.27Hzごとにパワーを示すがこの試験は発信音波をそのまま録音するのが目的であるので4995.83Hzのパワーが極端に大きくなる時刻を発見できれば良い。
録音データの発信が始まった時間と受信した時間の差と音速から実際の距離50m(往復距離100m)となっていることを確認する。
図4-1に概略図を示す。なお、マイクをパラボラに取り付けたものをパラボラ-m、スピーカーをパラボラに取り付けたものをパラボラ-sと表す。
図4-1 壁面反射試験概略図
4.1.3 測定データ
録音したデータの発信および受信したと思われる部分を抜き出したものを図4-2に示す。縦軸はWAVEデータをAD変換した際の数値で、横軸が時間である。
@ A B C D
図4-2 壁面反射実験録音データ
図4-2中のBに見られる波形は録音データ中ずっとみられ、電源等ノイズの周波数と思われる。そこで、図4-2中AとCに見られる大きな波形の乱れに注目する。発信パルス音波の録音された部分と思われるので図4-2中5箇所のパワースペクトルをFFTにより表す。注目しているのは発信パルス音波であるのでFFTの結果の中からその周波数である5000Hz前後だけを以下の図(図4-3〜4-7)に示す。
図4-3 @でのFFTの結果
図4-4 AでのFFTの結果
図4-5 BでのFFTの結果
図4-6 CでのFFTの結果
図4-7 DでのFFTの結果
図4-3から図4-7の5つのFFTの結果より、AとCに見られる波形の大きな振れが5000Hzであり、発信パルス音波であることが確認できる。
4.1.4 考察
FFTの結果より、受信波形に見られる振れが発信パルス音波であると確認できた。ここで、反射波を録音した図4-2のC付近だけを抜き出すと図4-8となる。
E F
図4-8 録音データのC付近
図4-8の波形において反射音が受信され始めた時刻を特定するのは難しいので、波形に変化が見られる図4-8中EとFで反射音の時刻を採った場合を考える。
発信パルス音波が直接受信されてから反射音が受信されるまでの時間は受信音波を図4-8のEで採った場合は約0.293秒間で、図4-8のFで採った場合約0.298秒間となる。試験当日の気温の12.8℃と第2章の式(3)より空気中の音速を339.18㎧とすると往復距離が99.38〜101.1mとなる。これは本来の往復距離100mを値の範囲内に含んでおり、距離の測定が可能であると言える。
しかし、発信音が直接録音されたとする図4-2中A部分では発信パルス幅を0.01秒に設定したのにかかわらず、約0.05秒間パルス波が受信されている。
ここで、図4-2のA付近だけを抜き出すと図4-9になる。
図4-9 録音データのA付近
図4-9より、受信波形の振幅が徐々にBやDに見られる波形に近づいていくのが確認できるので、要因の一つとしてパラボラが振動していると考えられる。
また、もう一つの要因として地表面からの反射音の受信が考えられる。ここで、録音開始から各時間の5000Hzのパワーを抜き出すと図4-10のようになる。
G
図4-10 録音開始からの5000Hzの受信パワー
発信音波を直接録音した部分と壁面からの反射音を録音した部分だけの5000Hzのパワーが大きくなることを期待したが、図4-10より、直接音を録音している(図4-10中G)時刻から約0.04秒後に再び5000Hzのパワーが大きくなっていることが確認できる。これは地表面からの反射音を受信しているとすると、前方約6.8m先の地表面からの反射音を録音しており、受信側パラボラの指向性が悪く、10°方向の反射音も録音してしまっているということになる。
これらの問題を解決するにはパラボラの振動防止と指向性の向上が必要となる。
4.2 大気反射実験
4.2.1 概要
水平方向観測用パラボラにおいてもっとも大きな障害となると思われるのが水平方向にパラボラを向けたときに入るノイズである。特に受信機器近くに存在する構造物からの反射や、受信機から最も近くに存在する地面からの反射が問題となる。そこでまず水平方向よりもノイズの少ない鉛直方向斜め45°にパラボラを傾けて試験する。
4.2.2実験方法
50m以内に音を反射させると思われる障害物のないところに2台のパラボラを上向きに45°傾けて並べる。次にパルス音波発信側のPCで発信音波を5000Hz、パルス幅は0.05秒に設定する。受信側のPCで録音を開始し、発信側PCでパルス音波を発信する。実験により得られた録音データの解析を行う。
・ 録音データの解析について
パルス幅が0.05秒(17m)であり、FFTを実行する範囲よりも小さいと発信パルス音波を特定できない可能性があるのでデータ数を2,048(211)個づつに分ける。また、この際の1ファイル分の時間は0.0464秒である。
この録音データにFFTを実行すると21.53Hzごとにパワーを示す。これと第2章の式(1)より風速分解能は0.7m/sである。
ノイズを除去するために30回行い、直接発信音波を受信した時刻を0秒とし、同じ時刻のデータを平均してパワーを比較する。受信パワーの大きくなっている周波数と録音時刻から、位置と風速の算定を行う。
4.2.3 測定データ
録音データの発信パルス部分を図4-11に示す。縦軸がWAVEデータをAD変換した際の数値で横軸が時間である。
@
図4-11 大気反射実験録音データ
録音データの中で他に見られない波形をした図4-11中@が直接音波を受信した部分であると考えられ、FFTの結果5000Hzであった。よってその時刻を直接音波を受信した時刻とする。
4.2.4 考察
録音データの直接音波を受信した図4-11中@の時刻から0.0464〜0.0928秒後のデータを取り出し、FFTを実行した結果を図4-12に示す。縦軸は対数目盛であり、一目盛100倍である。
図4-12 0.0464〜0.0928秒後のFFTの結果
ドップラーシフトにより受信パワーの最大値を示す周波数が変化していることを期待したが、図4-12より発信音波の周波数に1番近い4995Hzにピーク値を示していることがわかる。
原因として、パラボラの指向性が悪くビームが絞れていないために大気からの反射音だけでなく、地面や建物など固定物体からの反射音も録音していることが考えられる。
ここで、5000Hzのsin波にデータ数2048個でFFTを実行した場合の結果を図4-13に示す。縦軸は対数目盛である。
図4-13 5000Hzのsin波のFFTの結果
図4-13より、5000Hzのsin波にデータ数2048個でFFTを実行すると4995Hz以外の周波数にもパワーを示すことが分かる。ここで、5000Hzのsin波の4995Hzのパワーを、周波数のパワーをとし、を得られたデータの周波数に各々掛けるとドップラーシフトを起こした大気からの反射音が見えやすくなる。
実験により得られた録音データにFFTを実行して得られた各周波数の振幅の2乗の値と距離の関係は図4-14となる。縦軸は対数目盛である。ここで、発信パルス音波の影響を受けていないと思われる3014Hzから4026Hzのパワー平均値をノイズ平均とした。
図4-14 各周波数のパワーと距離
固定物体からの反射音を受信してしまっていると考えられるので、大気からの反射音を見やすくするためを得られたデータの周波数に各々掛けると図4-15の結果が得られる。
図4-15 振幅の2乗にを掛けた値と距離
この結果より、約8m先(上空約5.5m)では4866Hzの受信パワーが大きく、式(1)に代入すると西南西の風、風速4.5m/sとなる。これは視線方向の風速であり、鉛直方向の風速は水平方向の風速と比較して十分小さいと考えられるので鉛直方向の風速を0m/sとすると水平方向に6.4m/sである。また約16m先(上空約11m)では4844Hzの受信パワーが大きく、西南西の風、風速5.3m(水平方向7.5m)となる。
岐阜市加納にある岐阜地方気象台のデータによると実験を行った時刻は西南西4m/sの風であり、今回得られた実験結果と矛盾していない。
しかし、録音データは4995Hzのパワーが大きくを得られたデータに掛けなければ大気からの反射音を確認することができなかった。これは、パラボラの指向性が悪く、動かない地表面からの反射音を録音してしまっていると思われ、大気からの反射音をはっきりと観測するにはパラボラの指向性を向上させる必要があると考えられる。
5 パラボラ2号機の製作と実験
5.1 パラボラの設計と製作
水平方向観測を行うにあたって、非常に微弱な空気からの反射音が地表面からの反射音やノイズに埋もれてしまわないようにするには加藤氏により製作されたパラボラ(以下パラボラ1号機)の実験結果より、パラボラ本体の改良が必要となった。そこで、パラボラ2号機を製作する。
パラボラの放物線の式は焦点の座標をとして
(6)
で表される。第4章の大気反射実験の結果から、30mのオーダーで観測を行うにはパラボラ1号機の約129倍以上のパワーゲインが必要となる。しかし実験を行うにあたって、持ち運びができる大きさのものが良い。それに加え、使用しているマイクは90°方向からの音波も受信していると考えられ、パラボラ1号機はマイクが放物線の外側に設置されているので90°方向からの音波が直接マイクに録音されている可能性がある。よって、今回パラボラ2号機は焦点を放物線の内側に設定することで高い指向性を目指し、直径120cm、焦点距離20cmで設計した。マイク受信部の直径が3.5cm程度であり、またパラボラ1号機の直径が80cmであるから、マイク単体での受信の約1000倍、パラボラ1号機の約2倍のパワーで受信できることになる。さらに、使用していたスピーカーを一般のスピーカーと比較して出力音圧レベルが約128倍のスピーカーに取り替えることでパラボラ1号機の合計200倍以上のパワーで受信可能となる。図5-1は製作するパラボラの概略図である。
図5-1 製作するパラボラの概略図
受信音波を5000Hzとすると1波長は6cmであり、パラボラの表面の凹凸は波長より十分に小さい6mmを目標とした。またこの場合、式(2)にあてはめると半値幅は3.5°となるが、この値は発信音波周波数の値を変えることで調整が可能である。
また、安価で作り上げることを目標とし、パラボラの材料はホームセンターで手に入るもので自作した。
パラボラ2号機の製作に使用した材料は以下のとおりである。
・アルミ棒材
・アルミ板
・木板
・針金
・亜鉛製工作ネット
製作手順
@ 薄い木板を切り放物線の形にする。
A パラボラ本体の骨組みとなるアルミ棒材をその放物線に沿って曲げ、固定するための穴を空けておく。
写真 1
B パラボラの中央部分となるアルミ板を正十二角形に切り、Aで作ったアルミ棒材を取り付ける。
写真 2
C 中央板に取り付けたアルミ棒材同士を針金でつなぎ、固定する。
写真 3
D 亜鉛製工作ネットをCで作った骨組みに固定し、その内側表面に薄いアルミ板をエポキシ樹脂系接着剤で接着する。
写真 4
E 接着剤が乾いたら用意しておいた土台に取り付ける。
写真 5 写真 6
完成したパラボラ2号機に@で作った木板をあてて放物線との誤差を測定したところ、最大1cmの誤差が生じていた。原因として、手順Aで作ったアルミ棒材に後の工程で力が加わり変形してしまったことと、手順Dの際に亜鉛製工作ネットとアルミ板がうまく接着できておらず浮いてしまった部分があるということが考えられる。
目標誤差に対して二倍近くの誤差が生じてしまったが、発信音波を5000Hzに設定した場合の1波長の1/6の大きさであったのでこのまま実験を行う。
5.2 特性試験
5.2.1 概要
今回作成したパラボラ2号機の特性を知るために音波の送受信テストを行う。
5.2.2 試験方法
二台のパラボラを50mの間隔をおき向かい合わせに並べる。マイクをパラボラに取り付けたものをパラボラ‐m、スピーカーをパラボラに取り付けたものをパラボラ‐sとし、お互いを正面に向けた状態をとしてパラボラ‐mの正面からの角度を0°から180°まで5°づつずらして録音を行う。180°方向までの観測が終了したらパラボラ‐sとパラボラ‐mを交換して同様の観測を行う。発信音波は5000Hz、パルス幅は0.5秒に設定し音波を発信、受信する。録音データの発信した音波の受信部分についてFFTを実行し4995Hzのパワーを比較することで解析を行う。
5.2.3 測定データ
録音したデータを見ると発信したと思われるパルス音波を確認することができた。図5-2はにおける録音データである。
@ A B
図5-2 特性試験の録音データ
図5-2中@やAに見られる波形は録音開始から終了まで見られ、電源等ノイズの周波数と思われるのでBに見られる部分が発信したパルス音波を受信したと思われる。その部分を取り出し、FFTを実行することでパワーの比較が可能で、容易に解析可能である。
5.2.4 考察
以下にFFT実行により得られたパラボラ2号機の特性を示す。
図5-3 パラボラ‐mの特性図
図5-4 パラボラ‐sの特性図
パラボラ‐mでは、=0°のパワーを1として60°方向が0.01、90°方向が0.02、180°方向が0.3となった。 パラボラ‐sでは、=0°のパワーを1として60°方向が0.12、90°方向が0.01、180°方向が0.4となった。
しかし図5-3と図5-4から、60°方向まで受信パワーが下降せず設計値ほど良い指向性が出ていないことがわかる。
原因として、録音のAD変換の際に受信音量が大きいとデータが‐32768から32768までの数値しか取れないために矩形波になってしまい、FFTが正しく実行ができなかったことが考えられる。実際に60°方向までそれぞれの受信音波に矩形波になってしまった部分が多く見られた。さらに、パラボラ本体に歪みが生じてしまっていることや、パラボラ本体を音が透過してしまっていることが考えられる。
5.3 壁面反射実験
5.3.1 概要
パラボラ2号機の距離を測る能力を試験する。パラボラ1号機でも距離の測定は可能であったが、4.1.4で指摘したパラボラ本体の振動や指向性の課題が残ったので、その点も評価する。
5.3.2 実験方法
第4章の壁面反射実験と同じ方法で行う。
5.3.3 測定データ
録音したデータの発信および受信したと思われる部分を抜き出したものを下図に示す。縦軸がWAVEデータをAD変換した際の数値で、横軸が時間である。
@ A B C E D
図5-5 壁面反射実験の録音データ
図5-5中のBに見られる波形は録音データ中ずっとみられ、電源等ノイズの周波数と思われる。そこで、図4-2中AとCに見られる大きな波形の乱れに注目する。発信パルス音波の録音された部分と思われるので図5-5中@〜Dの各周波数のパワースペクトルをFFTにより表す。注目しているのは発信パルス音波であるのでFFTの結果の中からその周波数である5000Hz前後だけを以下の図5-6〜5-10に示す。
図5-6 @でのFFTの結果
図5-7 AでのFFTの結果
図5-8 BでのFFTの結果
図5-9 CでのFFTの結果
図5-10 DでのFFTの結果
図5-6から図5-10の5つのFFTの結果より、AとCに見られる波形の大きな振れが5000Hzであり、発信パルス音波であることが確認できる。
5.3.4 考察
FFTの結果より、受信波形に見られる振れが発信パルス音波であると確認できた。また、波形からこの発信パルス音波が直接受信されてから反射音が受信されるまでの時間を特定するには難しいので4.1.4章と同様にして波形に変化が見られる地点での時間から距離の測定が可能か判断すると、実験当日の気温が10.2℃より空気中の音速が337.62㎧であり、往復距離が98.59〜100.6mとなる。これは本来の往復距離100mを値の範囲に含んでおり、パラボラ2号機での距離の測定も可能であると言える。
しかし、図5-5中Eの部分に注目すると波の重ねあわせにより波形が徐々にBやDに見られる波形に近づいていくのが見られるので、パラボラが振動してしまっている可能性がある。
次に、録音開始から各時間の5000Hzのパワーを抜き出したものを図5-11に示す。
F H G
図5-11 録音開始からの5000Hzの受信パワー
図5-11より、パルス音波の反射音を録音した部分(図中F)のパワーが直接音を録音した部分(図中G)のパワーの約5倍の大きさで録音できている。直接音は録音側のパラボラの横90度方向からの音を録音しており、反射音は録音側のパラボラの正面からの音を録音しているので、この結果よりパラボラ2号機はパラボラ1号機よりも高い指向性があると確認できる。
しかし、第1章の式(2)をそのまま使うと直接音と反射音のパワーの比は5×106倍であり、それには及ばない結果になってしまった。さらに、直接音を録音している時刻から約0.047秒後(図5-11中H)に再び5000Hzのパワーの上昇がみられる。これは地表面からの反射音を録音しているとすると前方約8m先(パラボラの10°方向)の地表面からの反射音を録音してしまっていることとなる。
ここで、今回使用したマイクの指向特性が明らかではないので一般的なマイクの指向特性図を図5-12に示す。
図5-12より、一般のマイクは正面のパワーを1とすると90°方向が約0.5、180°方向が約0.016である。また150°方向から180°方向に向かって受信パワーが上昇していることが分かる。
今回使用したマイクがこのような指向特性をもっているとしたら、図5-11中のHではマイクの170°方向(パラボラの10°方向)からの音波がパラボラにより位相がずらされることなく、直接受信されてしまっている可能性があるということが分かった。
今回の壁面反射試験により得られたデータでは、パラボラ90°方向からの音波のパワー(図5-11中のF)がパラボラ正面からの音波のパワー(図5-11中のG)よりも十分小さく、焦点を放物線の中に設定したことで90°方向からの音波の受信パワーをパラボラによって低下させることができたことが分かる。
今回は残念ながらできなかったが、次に指向性の高いマイクに取り替えての実験が必要となる。今後取り替えるマイクとして図5-13のような指向特性をもったマイクがあげられる。
図5-12 一般的なマイクの指向特性 図 5-13 高指向性マイクの特性
壁面反射実験の結果より、地表面からの反射音を受信してしまっていることが確認されたが、パラボラ1号機よりも2号機が高い指向性を持っていることが確認できた。固定物体からの反射波はドップラーシフトを起こすことなく発信音波を反射するので発信音波の周波数以外のパワーに注目することでドップラーシフトをおこした大気からの反射音を確認することが可能であると考えられる。よって大気反射実験を行った。
5.4 大気反射実験
5.4.1 概要
水平方向観測よりもノイズが入りにくい鉛直方向斜め45°方向にパラボラを傾けて実験することで、パラボラ2号機で大気からの反射音を録音可能であるか試験する。
5.4.2 実験方法
第4章の大気反射実験と同じ方法で行う。
5.4.3 測定データ
録音データの発信パルス部分を図5-14に示す。縦軸がWAVEデータをAD変換した際の数値で横軸が時間である。
@
図5-14 大気反射実験の録音データ
録音データの中で他に見られない波形をした図5-14中@が直接音波を受信した部分であると考えられ、FFTの結果5000Hzであった。よってその時刻を直接音波を受信した時刻とする。
5.4.4 考察
録音データの直接音波を受信した時刻から0.0464〜0.0928秒後のデータを取り出し、FFTを実行した結果を図5-15に示す。縦軸は対数目盛である。
図5-15 0.0464〜0.0928秒後のFFTの結果
次に、各周波数の振幅の2乗の値と距離の関係を図5-16に示す。
図5-16 各周波数のパワーと距離
図5-3からもわかるように、指向性が悪いので地表面からの反射音を録音してしまっていると考えられる。
ここで、第4章の大気反射実験と同様にして大気からの反射音を見やすくするために各周波数にそれぞれを掛けると図5-17が得られる。
図5-17 振幅の2乗にを掛けた値と距離
図5-17より、約8m先(上空約5.5m)では4888Hzの受信パワーが大きく、式(1)より西南西の風、風速3.8m/s(水平方向5.3m/s)であり、約16m先(上空約11m)では4909Hzの受信パワーが大きく西南西の風、風速3.1m(水平方向4.3m/s)となる。
岐阜市加納にある岐阜地方気象台のデータによると実験を行った時刻は西南西2.5m/sの風であり、矛盾していない。
しかし、今回の実験はパラボラ1号機の大気反射実験時と同じ場所、同じ風向で風速も同程度であったので完全に風が観測できたと断言できない。例えば、パラボラの固有振動数が不明であるので発信パルス音波の周波数から少しずれた固有振動数をもつパラボラの共振現象とも考えられる。そこで、発信パルス音波の周波数や音波の発信方向を変化させた実験が必要となるが今回は残念ながらこれ以降の実験ができなかった。
2回の大気反射実験のデータ解析において5000Hzのsin波にFFTを実行した際の結果を利用して大気からの反射音を見え易くするという方法を採ったが、水平方向へパラボラを向けて観測を行った場合パラボラの指向性を向上させても地表面からの反射音を録音してしまう事が予想され、大気からの反射音を観測するには今回行ったような方法が必要であると考えられる。
6 結論
今回、大規模組織構造対流の観測を目指して水平方向観測用ドップラーソーダの試作を行った加藤氏のシステムに改良を加え、この装置の完成を目指して研究した。
解析(FFT)プログラムの導入を行ったことにより、録音データから各周波数のパワースペクトルの表示が可能となった。これにより、受信データに見られる波形の大きな振れが発信パルス音波であるという確認ができるようになった。
パラボラ1号機の壁面反射実験と大気反射実験の結果から、反射音からの距離の測定が可能であり、大気からの反射音の録音ができそうであるとわかったが、指向性と振動の問題が残ったため水平方向観測に移る前にパラボラ2号機を製作した。
製品のドップラーシステムの価格は千数百万円するが、パラボラ2号機の製作費はノートパソコンの値段を除外すると合計15万円以下で製作できた。また、パラボラ1号機は土台が非常に不安定であり実験が困難であったが、パラボラ2号機は風によってパラボラが転倒したり傾いてしまう心配が無くなった。
パラボラ2号機の壁面反射実験の結果から、反射音からの距離の測定が可能であり、さらにパラボラ1号機よりも指向性が向上されていることが確認できた。大気反射実験は大気からの反射音の録音が可能であると思われる結果であった。しかし、大気反射実験の結果はパラボラ2号機にも5000Hzの大きなパワーの音波が受信されており、水平方向にパラボラを向けて観測を行うには更なる指向性の向上が必要であった。
パラボラ2号機の指向性はパラボラ‐mでは正面のパワーを1として90°方向が0.02、180°方向が0.3となった。パラボラ‐sでは正面のパワーを1として90°方向が0.01、180°方向が0.4となった。しかし、60°方向まで受信パワーが下がらなかったので今後追加実験として1°刻みで受信音波が矩形波にならない音量での特性試験を行った上で設計値との比較・検討が必要である。
水平方向の観測を行うには、高指向性マイクに取り替えた上で発信パルス音波の周波数や音波の発信方向を変化させた実験を行い、高指向性のパラボラ3号機の製作を行う必要がある。
謝辞
本研究を遂行するにあたり、ご多忙であるにも関わらず、終始適切なご指導・ご援助を頂きました岐阜大学工学部玉川一郎助教授には心から感謝の意を表します。また、当研究室の学部生、院生の方にも数々のご協力・ご支援を頂きました。
ありがとうございました。
参考文献
1)加藤聡明(2003):卒業論文 水平方向観測用ドップラーソーダの試作
2)神田ら(2002):水文・水資源学会誌 LESによる熱収支インバランス問題に関する検討
3)田辺義和(2003):windowsサウンドプログラミング
4)小池慎一(1994):Cによる科学技術計算
5)William H.Press et al著 丹慶勝市ら訳(1993):NUMERICAL RECIPES in C [日本語版]
6)原康夫(2000)物理学通論T
7)小林隆久(2004)ウィンドプロファイラー‐電波で探る大気の流れ‐
8)株式会社カイジョー wwwページ (http://www.kaijo.co.jp)
9)株式会社オーディオテクニカ wwwページ
(http://www.audio-technica.co.jp/index.html)
10)JR4MDA wwwページ (http://www.chukai.ne.jp/~jr4mda/index.html)
11)JA9NEG wwwページ (http://www2.nsknet.or.jp/~okino/index.htm)
付録1 リニアPCMデータ出力プログラムソース
#include<stdio.h>
#include<math.h>
#define PI 3.141529
int main(void)
{
FILE *fp;
float playtime;
float puls;
int hz;
unsigned short data2;
int data4;
long i;
playtime=0;
puls=0;
data2=0;
data4=0;
i=0;
fp=fopen("output.wav","wb");
fprintf(stdout,"何秒おきに発信しますか?(秒)\n");
fscanf(stdin,"%f",&playtime);
fprintf(stdout,"パルス幅を入力してください(秒)\n");
fscanf(stdin,"%f",&puls);
data4=playtime*44100.0*2+38;
fprintf(stdout,"周波数を入力してください(Hz)\n");
fscanf(stdin,"%d",&hz);
fwrite("RIFF",4,1,fp);/*RIFFヘッダ*/
fwrite(&data4,4,1,fp);/*これ以降のファイルサイズ(ファイルサイズ-8)*/
fwrite("WAVE",4,1,fp);/*WAVEヘッダ*/
fwrite("fmt ",4,1,fp);/*fmtチャンク*/
data4=16;/*fmtチャンクのバイト数*/
fwrite(&data4,4,1,fp);
data2=1;/*フォーマットID*/
fwrite(&data2,2,1,fp);
data2=1;/*チャンネル数*/
fwrite(&data2,2,1,fp);
data4=44100;/*サンプリングレート*/
fwrite(&data4,4,1,fp);
data4=88200;/*データ速度*/
fwrite(&data4,4,1,fp);
data2=2;/*ブロックサイズ*/
fwrite(&data2,2,1,fp);
data2=16;/*サンプルあたりのビット数*/
fwrite(&data2,2,1,fp);
fwrite("data",4,1,fp);/*dataチャンク*/
/*バイト数*/
data4=44100.0*playtime*2;
fwrite(&data4,4,1,fp);
/*波形データ*/
for(i=0;i<=44100.0*playtime;i++)
{
if(i<=44100.0*puls)
{
data2=32700.0*sin(2*PI*hz*i/44100);
fwrite(&data2,2,1,fp);
}
else
{
data2=0;
fwrite(&data2,2,1,fp);
}
}
fclose(fp);
}
付録2 AD変換プログラムソース
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define SKIP_CH 22
#define SKIP_TO_DATANUM 42
#define SKIP_DATA 46
#define OUTPUTNUM 2048 /* 2048 データずつファイルに分ける */
int main( int argc, char *argv[] ) {
FILE *fp, *fp2;
int i, j, count=0, count2=0;
char outfilename[1024];
int datanum;
short data2[2],ch;
fprintf( stdout, "%d %s %s\n", argc, argv[0], argv[1] );
if( argc < 3) {
fprintf( stdout, "command datafile.wav outputfile\n");
exit(EXIT_FAILURE );
}
fp = fopen( argv[1], "rb" );
if( fp == NULL ) {
fprintf( stdout, "Can not open %s\n");
}
fprintf(stdout, "open %s\n", argv[1]);
fseek( fp, SKIP_CH, SEEK_SET );
fread( &ch, 2, 1, fp);
fprintf( stdout, "ch = %d\n", ch );
/*fseek( fp, SKIP_TO_DATANUM, SEEK_SET );
fread( &datanum, 4, 1, fp);
fprintf(stdout, "data number = %d\n", datanum );*/
fseek( fp, SKIP_DATA, SEEK_SET );
for(i=0; /*i<datanum*/; i++){
for(j=0; j<ch; j++){
if(0==fread( &data2[j], 2, 1, fp)) goto close;
/* fprintf( stdout, "%d, ", data2[j]); */
}
if( count % OUTPUTNUM == 0 ) {
if( count != 0 ) {
fclose( fp2);
}
sprintf( outfilename, "%s%0d.csv", argv[2], count2);
/*
fprintf( stdout, "open %s for output\n", outfilename);
*/
fp2 = fopen( outfilename, "w" );
count2++;
}
if( ch == 1 ) {
fprintf( fp2, "%d\n", data2[0] );
fprintf( stdout, "%d %d %d\n",count,count2, data2[0] );
} else {
fprintf( fp2, "%d, %d\n", data2[0], data2[1] );
fprintf( stdout, "%d, %d\n", data2[0], data2[1] );
}
count++;
}
close: fclose(fp);
}
付録3 FFTプログラムソース
#include <math.h>
#include <stdio.h>
#define M_PI 3.14159265358979323846
void fft(double* x, double* y, int l, double f);
#define DATANUM 2048 /* データ数2048(N=2048)でfftを行う */
int main( int argc, char *argv[] ) {
double x[DATANUM], y[DATANUM], hz[DATANUM];
int i, count=0, count2=0;
FILE *fp, *fp2;
char outfilename[1024];
fprintf( stdout, "%d, %s, %s\n", argc, argv[0], argv[1] );
fp = fopen( argv[1],"rb");
if(fp == NULL) {
fprintf(stdout, "Can not open %s\n");
}
for(i=0; i<DATANUM; i++){
fscanf( fp,"%lf", &x[i]);
y[i] = 0;
}
sprintf( outfilename, "%s.csv", argv[2]);
fp2 = fopen( outfilename, "w");
fft( x, y, 11, -1.0 );
for( i=0; i<DATANUM; i++){
hz[i] = 44100*i/DATANUM;
/*fprintf ( stdout, "%d, %lf, %lf, %lf\n", i, x[i], y[i], x[i]*x[i]+y[i]*y[i] );*/
fprintf ( fp2, "%lf, %lf\n", hz[i], x[i]*x[i]+y[i]*y[i] );
}
}
void fft(double* x, double* y, int l, double f) /* FFT 1 (c) */
{
int i, i0, i1, j, l1, n, ns, k;
double s, c, s1, c1, sc, x1, y1, t;
/*n = ipow2(l);*/ /* n = 2^l*/
n=1;
for(i=0; i<l; i++){
n = 2*n;
}
sc = M_PI;
j = 0;
for( i = 0; i < n-1; i++ ) { /* bit reversi counter */
if( i <= j) {
t = x[i]; x[i] = x[j]; x[j] = t;
t = y[i]; y[i] = y[j]; y[j] = t;
}
k = n / 2;
while( k <= j ) {
j = j-k; k /= 2;
}
j = j + k;
}
ns = 1;
while( ns <= n/2 ) { /* main loop */
c1 = cos(sc); s1 = sin(f * sc);
c = 1.0; s = 0.0;
for( l1 = 0; l1 < ns; l1++ ) {
for( i0 = l1; i0 < n; i0 += (2*ns) ) {
i1 = i0 + ns;
x1 = x[i1] * c - y[i1] * s; y1 = y[i1] * c + x[i1] * s;
x[i1] = x[i0] - x1; y[i1] = y[i0] - y1;
x[i0] = x[i0] + x1; y[i0] = y[i0] + y1;
}
t = c1 * c - s1 * s; s = s1 * c + c1 * s; c = t;
}
ns = ns * 2; sc = sc / 2.0;
}
if( f < 0.0 ) /* test fft or oft */
for( i = 0; i < n; i++ ) {
x[i] /= (double)n; y[i0] /= (double)n;
}
}