JP2008304267A - 懸濁物質濃度の分布解析装置及びプログラム - Google Patents

Abstract

【課題】簡易に水中の横断面の懸濁物質濃度の分布を解析する。
【解決手段】浮遊土砂濃度の分布解析対象となる河川の横断面に計算格子を設定し（ステップ１００）、河川で実際に計測された浮遊土砂濃度の計測地点に相当する計算格子に該測定値を設定し、他の計算格子には初期値を設定した状態で、浮遊土砂の沈降量を表す項と、浮遊土砂の水平方向の拡散フラックスを表わす項、浮遊土砂の鉛直方向の拡散フラックスを表わす項、及び付加項の和との関係を表わした基礎式を用いて、付加項の値を演算する第１の演算（ステップ１０４〜１０６）と、第１の演算で演算された付加項の値を前記基礎式に設定した状態で計測地点に相当する計算格子以外の他の計算格子の浮遊土砂濃度を演算する第２の演算（ステップ１０８〜１１２）とを、付加項と浮遊土砂濃度の値が収束するまで繰り返し行なって（ステップ１１４）、浮遊土砂濃度の分布を求める。
【選択図】図７

Description

本発明は、二次元平面の懸濁物質の濃度分布を解析する懸濁物質濃度の分布解析装置及びプログラムに関する。

水系一貫土砂管理を行う上では、源流域から河口、沿岸にわたる流域圈全体の土砂動態を把握することは不可欠であり、掃流砂・浮遊砂量を自動連続的にモニタリングすることは必須である。このうち浮遊土砂輸送量（＝浮遊土砂濃度（ＳＳ）×流量）の計測には、河川の横断面において大きく変化している流速・ＳＳ分布を把握する必要があるが、これらの量をほぼ瞬時に直接計測できる観測システムは皆無である。

なお、流速を計測する装置として、水平設置型超音波ドップラー流速分布計Ｈ−ＡＤＣＰ（Horizonta1 Acoustic Doppler Current Profiler）が知られている。また、Ｈ−ＡＤＣＰは、流速のみならず、濁度と関連性の高い超音波の反射強度の計測ができるため、Ｈ−ＡＤＣＰは上記の要請に対して有望な機器である。

しかしながら、Ｈ−ＡＤＣＰはある高さの「線」データしか計測できないため、横断面全体のＳＳ分布を取得することができない、という問題がある。また、光学式センサを用いて濃度データを取得する場合も同様である。

なお、超音波ドップラー流速計を用いて水中の二次元平面の濁度分布を観測する装置としては、ドップラー流速計を移動させながら水位、流速分布、河床位を同時観測する装置が提案されている（例えば、特許文献１参照。）。また、観測船にＡＤＣＰを設置し、観測船を海上で移動しながら濁度を観測する装置も提案されている（例えば、特許文献２参照。）。
特開２０００−１１１３６５号公報 特開２００７−７１８８１号公報

しかしながら、上記従来の技術では、ＡＤＣＰを移動させながら濃度分布を観測するため、時間やコストがかかる。少ない観測値で河川の横断面の濃度分布を求める技術が望まれている。

本発明は、上述した問題を解決するためになされたものであり、簡易に水中の横断面の懸濁物質濃度の分布を解析することができる懸濁物質濃度の分布解析装置及びプログラムを提供することを目的とする。

上記目的を達成するために、本発明の懸濁物質濃度の分布解析装置は、懸濁物質濃度の分布を解析する水中の横断面内に複数の格子点を配列することにより、隣接する格子点で囲まれた区域を複数備えた仮想横断面を想定する想定手段と、測定地点に対応する前記区域に懸濁物質濃度の測定値を設定すると共に、前記測定値を設定した区域以外の区域の各々に懸濁物質濃度の初期値を設定する設定手段と、懸濁物質の沈降量を表す項と、懸濁物質の水平方向の拡散フラックスを表わす項、懸濁物質の鉛直方向の拡散フラックスを表わす項、及び付加項の和との関係を前記区域の各々について表わした基礎式、前記測定地点に対応する区域に設定された懸濁物質濃度の測定値、及び前記測定地点に対応する区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記測定地点に対応する区域の付加項の値を演算する第１の演算、並びに前記測定地点に対応する区域以外の区域の１つを注目区域として、前記基礎式、前記演算された付加項の値、前記注目区域に設定された懸濁物質濃度の値、及び前記注目区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記注目区域の懸濁物質濃度を演算し、演算により得られた演算値を前記注目区域に既に設定された値に代えて設定することを、前記注目区域を代えながら前記仮想横断面全体にわたって繰返し行う第２の演算を、前記懸濁物質濃度の演算値及び前記付加項の演算値の各々が収束するまで繰返し行う繰返演算手段と、を含んで構成されている。

このように、本発明は、分析対象となる横断面を複数の区域に分割した仮想横断面において、第１の演算及び第２の演算を、基礎式の付加項と各区域の懸濁物質の濃度が収束するまで繰り返し行なって懸濁物質濃度の分布を求めるようにしたため、懸濁物質濃度の実際の測定値が仮想横断面において少なくとも１箇所でもあれば、簡易に懸濁物質濃度の分布を求めることができる。

なお、前記基礎式を以下の式で表すことができる。

ここで、ｙは水平方向、σは鉛直方向、Ａ_Ｈは水平渦動粘性係数、Ａ_Ｖは鉛直渦動粘性係数、Ｄは水深、Ｃは懸濁物質濃度、ｗ_０は懸濁物質の沈降速度、Ｆa,cは、付加項である。

また、所定位置に設置された超音波ドップラー流速分布計により水中に発信され該水中の懸濁物質により反射され計測された超音波の反射強度のうち、前記超音波ドップラー流速分布計近傍の懸濁物質群で反射された超音波の反射強度を示す第１反射強度と、前記超音波ドップラー流速分布計から離れた懸濁物質群で反射された超音波の反射強度を示す第２反射強度と、に基づいて前記懸濁物質濃度の測定値を算出する濃度測定手段を更に含んで構成することもできる。

反射強度の濁度依存性は、超音波ドップラー流速分布計からの距離に応じて異なる。例えば、超音波ドップラー流速分布計近傍の反射強度は低濁度時に濁度依存性が高く、高濁度時には濁度依存性が低いが、超音波ドップラー流速分布計から遠い位置の反射強度は高濁度時に濁度依存性が高い、という特徴がある。従って、この特徴を利用して、超音波ドップラー流速分布計近傍の反射強度及び遠く離れた位置の反射強度に基づいて懸濁物質濃度の測定値を算出することによって、精度高く測定値を算出することができる。

そこで本発明の懸濁物質濃度の分布解析装置における前記濃度測定手段は、前記第１反射強度が所定値以下の場合には、前記第１反射強度に基づいて前記懸濁物質濃度の測定値を算出し、前記第１反射強度が所定値を超える場合には、前記第１反射強度と前記第２反射強度との比に基づいて前記懸濁物質濃度の測定値を算出するようにしてもよい。

このように、反射強度の特徴を活かして濁度に応じて懸濁物質濃度の測定値を算出することによって、精度高く懸濁物質の濃度を求めることができる。

また、本発明のプログラムは、コンピュータを、懸濁物質濃度の分布を解析する水中の横断面内に複数の格子点を配列することにより、隣接する格子点で囲まれた区域を複数備えた仮想横断面を想定する想定手段と、測定地点に対応する前記区域に懸濁物質濃度の測定値を設定すると共に、前記測定値を設定した区域以外の区域の各々に懸濁物質濃度の初期値を設定する設定手段と、懸濁物質の沈降量を表す項と、懸濁物質の水平方向の拡散フラックスを表わす項、懸濁物質の鉛直方向の拡散フラックスを表わす項、及び付加項の和との関係を前記区域の各々について表わした基礎式、前記測定地点に対応する区域に設定された懸濁物質濃度の測定値、及び前記測定地点に対応する区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記測定地点に対応する区域の付加項の値を演算する第１の演算、並びに前記測定地点に対応する区域以外の区域の１つを注目区域として、前記基礎式、前記演算された付加項の値、前記注目区域に設定された懸濁物質濃度の値、及び前記注目区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記注目区域の懸濁物質濃度を演算し、演算により得られた演算値を前記注目区域に既に設定された値に代えて設定することを、前記注目区域を代えながら前記仮想横断面全体にわたって繰返し行う第２の演算を、前記懸濁物質濃度の演算値及び前記付加項の演算値の各々が収束するまで繰返し行う繰返演算手段と、して機能させるためのプログラムである。

このように、本発明は、上記懸濁物質濃度の分布解析装置と同様に、分析対象となる横断面を複数の区域に分割した仮想横断面において、第１の演算及び第２の演算を、基礎式の付加項と各区域の懸濁物質の濃度が収束するまで繰り返し行なって懸濁物質濃度の分布を求めるようにしたため、懸濁物質濃度の実際の測定値が仮想横断面において少なくとも１箇所でもあれば、簡易に懸濁物質濃度の分布を求めることができる。

なお、前記基礎式を以下の式で表すことができる。

ここで、ｙは水平方向、σは鉛直方向、Ａ_Ｈは水平渦動粘性係数、Ａ_Ｖは鉛直渦動粘性係数、Ｄは水深、Ｃは懸濁物質濃度、ｗ_０は懸濁物質の沈降速度、Ｆa,cは、付加項である。

更に、本発明のプログラムは、更に、前記コンピュータを、所定位置に設置された超音波ドップラー流速分布計により水中に発信され該水中の懸濁物質により反射され計測された超音波の反射強度のうち、前記超音波ドップラー流速分布計近傍の懸濁物質群で反射された超音波の反射強度を示す第１反射強度と、前記超音波ドップラー流速分布計から離れた懸濁物質群で反射された超音波の反射強度を示す第２反射強度と、に基づいて前記懸濁物質濃度の測定値を算出する濃度測定手段として機能させることもできる。

なお、前記濃度測定手段は、前記第１反射強度が所定値以下の場合には、前記第１反射強度に基づいて前記懸濁物質濃度の測定値を算出し、前記第１反射強度が所定値を超える場合には、前記第１反射強度と前記第２反射強度との比に基づいて前記懸濁物質濃度の測定値を算出するようにしてもよい。

以上説明したように、本発明の懸濁物質濃度の分布解析装置及びプログラムによれば、簡易に水中の横断面の懸濁物質濃度の分布を解析することができる、という優れた効果が得られる。

以下、図面を参照して、本発明の実施の形態を詳細に説明する。

本実施の形態では、河川において水中の懸濁物質（浮遊土砂濃度）を測定し、該土砂濃度に基づいて河川の横断面（水の主流方向に対してほぼ直交する断面であって、一方の岸から他方の岸までの横断面）における浮遊土砂濃度の分布を解析して出力する浮遊土砂濃度モニタシステムについて説明する。

本実施の形態の浮遊土砂濃度モニタシステムは、図１に示すコンピュータ１０と、図２に示す水平設置型超音波ドップラー流速分布計Ｈ−ＡＤＣＰ（Horizonta1 Acoustic Doppler Current Profiler）４０とを含んで構成されている。コンピュータ１０は、浮遊土砂濃度の分布を解析する。Ｈ−ＡＤＣＰ４０は、通常は河川の流速の計測に用いられるが、本実施の形態では、浮遊土砂濃度と相関のある超音波の反射強度を計測するために用いることとする。

図１（Ａ）は、コンピュータ１０の外観図であり、図１（Ｂ）は、コンピュータ１０のハードウェア構成を模式的に示すブロック図である。

コンピュータ１０は、ＣＰＵ１２、ＲＡＭ１４、ＲＯＭ１６、通信インタフェース１８、及びＩ／Ｏ（入出力）インタフェース２０を備え、それらはバス３８を介して相互に接続されている。

Ｉ／Ｏインタフェース２０には、ＨＤＤ（ハードディスクドライブ）２２、ＣＤ−ＲＯＭドライブ２４、ＦＤ（フレキシブルディスク）ドライブ２６が接続されている。ＨＤＤ２２は、内蔵されたハードディスクに対してデータを読み書きする。ＣＤ−ＲＯＭドライブ２４は、可搬型記録媒体であるＣＤ−ＲＯＭ３４に対してデータを読み書きする。ＦＤドライブ２６は、可搬型記録媒体であるＦＤ３６に対してデータを読み書きする。

上記ＨＤＤ２２やＲＯＭ１６等の記録媒体には、後述する浮遊土砂濃度の測定値算出処理や土砂濃度の分布解析処理等、様々な処理を行なうためのプログラムや処理に必要なデータが記録され、ＲＡＭ１４等のワークメモリを利用してＣＰＵ１２により実行される。

なお、ＣＰＵ１２が実行するプログラムが記憶される記録媒体は、上記ＣＤ−ＲＯＭ３４やＦＤ３６、ＨＤＤ２２等に限定されず、例えば、図示は省略するが、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの可搬型記録媒体やコンピュータ１０の外部に備えられたＨＤＤ等の記憶装置等であってもよく、さらにまた通信回線を介して接続されたデータベース、或いは他のコンピュータシステム並びにそのデータベースや、更に通信回線上の伝送媒体であってもよい。従って、後述する本発明の実施の形態に係るプログラム等を通信媒体により提供することはもちろん、ＣＤ−ＲＯＭ等の記録媒体に格納して提供することも可能である。

また、Ｉ／Ｏインタフェース２０には、画像やデータを表示する表示装置２８や、各種情報を入力したり指定したりするための入力装置３０が接続されている。さらにまた、プリンタ等の外部機器と接続するためのコネクタ３２等も設けられている。

通信インタフェース１８は、ネットワークに接続されている。この通信インタフェース１８を介して、ネットワークに接続された他の装置と相互に情報をやりとりすることができる。

このコンピュータ１０で河川の横断面の土砂濃度分布を解析するにあたり、まず、Ｈ−ＡＤＣＰ４０を用いて水中の浮遊土砂濃度と相関のある超音波の反射強度を計測する。そして、超音波の反射強度の測定結果から水中の濁度（水の濁りの程度を示す）を求める。

なお、濁度を測定するには、超音波を用いた測定の他、光学式濁度計を用いて測定することもできるが、本実施の形態ではＨ−ＡＤＣＰ４０を用いて、濁度を測定するものとする。

Ｈ−ＡＤＣＰ４０により、水中の懸濁物質で反射した超音波の反射強度から濁度を推定（算出）するには、厳密には水中の超音波吸収等を考慮する必要があるが、超音波吸収性を正確に設定することが困難であるため、以下に説明するように、本実施の形態では、反射強度値や反射強度比を利用して濃度を簡易的に求める。このように簡易的に濃度を求めたとしても、実測濁度と良好に一致し、後述する浮遊土砂濃度の分布を解析するために問題ない精度の結果が得られる。詳細は後述する。

ここで、Ｈ−ＡＤＣＰ４０の流速及び反射強度の測定方法について説明する。

図２に示されるように、低水路側岸部の河床に設置した支持体４２（例えば水位標）に対して、予め定められた高さ（ただし水中）にＨ−ＡＤＣＰ４０を固定設置する。このＨ−ＡＤＣＰ４０により、超音波を発信し、水中の懸濁物質群で反射して戻ってきた超音波を受信してその強さ（反射強度）を計測する。

なお、本実施の形態では、Ｈ−ＡＤＣＰ４０を例えば以下のように設定した。層厚（距離分解能）０．５ｍ、不感帯幅（Ｈ−ＡＤＣＰ４０から計測第１層目までの距離）２．０ｍ、層数１２８、サンプリング間隔１０分とする。このように設定することにより、Ｈ−ＡＤＣＰ４０では、１〜１２８の各層毎に、超音波を発信してから反射波が受信されるまでの時間に応じて各層毎に超音波の反射強度を計測することができる。１層から１２８層にいくに従って、Ｈ−ＡＤＣＰ４０の設置箇所から遠くなる。

ここで、図３を参照して、超音波の反射強度に関する基本的な特徴を説明する。図３（Ａ）は、ある河川における観測期間中の水位と実測濁度との関係を示すグラフであり、図３（Ｂ）は、同じ観測期間中の、第３層の反射強度Ｉ_３及び第８０層の反射強度Ｉ_８０を示すグラフである。また、図３（Ａ）において、水位の単位は、y.P.mを用い、実測濁度の単位はFTU(Formazin Turbidity Unit)を用い、図３（Ｂ）において、反射強度の単位は、独自に基準化したcountという単位を用いている。また、グラフに示した反射強度Ｉ_３、Ｉ_８０は、水中の超音波減衰等の補正を施していない生データである。

前述した設定により、第３層目の反射強度Ｉ_３の観測値は、Ｈ−ＡＤＣＰ４０からの横断距離３．２５ｍにおける反射強度の観測値を示し、第８０層目の反射強度Ｉ_８０の観測値は、Ｈ−ＡＤＣＰ４０からの横断距離４１．７５ｍにおける反射強度の観測値を示す。

図３から明らかなように、Ｈ−ＡＤＣＰ４０近傍の第３層目の反射強度Ｉ_３は、低濁度時には実測濁度に追随し濁度依存性が高いが、高濁度時には反射強度Ｉ_３が上限値（＝２２３count）が存在しており、濁度依存性が低いことがわかる。

一方、Ｈ−ＡＤＣＰ４０から遠い位置の第８０層目の反射強度Ｉ_８０に関しては、超音波の距離減衰の影響で、反射強度レベルが第３層目の反射強度Ｉ_３よりも大幅に小さくなっている。また、実測濁度が高くなると超音波減衰が顕著になり、Ｉ_３とは逆に反射強度Ｉ_８０が減少する傾向がある。しかしながら、Ｈ−ＡＤＣＰ４０に近い第３層目の反射強度Ｉ_３のように、上限値（下限値）は生じておらず、高濁度側ではＨ−ＡＤＣＰ４０から遠い第８０層目の反射強度Ｉ_８０の方が濁度依存性が高くなる。

次に、反射強度データと実測濁度の相関性を見るために、実測濁度Turbと第３層目の反射強度Ｉ_３の相関関係を調べ、図４に示すグラフにした。

図４から明らかなように、低濁度時には実測濁度と反射強度Ｉ_３には正の相関関係が認められるものの、高濁度時には反射強度は上限値（＝２２３count）付近に留まる。このことから、単純に、反射強度Ｉ_３から濁度を推定する場合、低濁度時には概ね良好な精度が期待されるが、高濁度時では推定そのものが困難となることがわかる。

なお、反射強度Ｉ_３≦２１０countの場合、低濁度時における実測濁度Turbと反射強度Ｉ_３の相関式として次式が得られた。

Turb = 0.0034^*exp(0.0397^*I₃) ・・・（１）
ここで、Turbの単位としてはＦＴＵを用い、反射強度Ｉ_３の単位としてはcountを用いる。

次に、第８０層の反射強度Ｉ_８０と第３層の反射強度Ｉ_３との比（Ｉ_８０／Ｉ_３）と、実測濁度Turbとの相関関係を調べ、図５に示すグラフにした。

図５に示すように、低濁度の場合には同一濁度に対して幅広く反射強度比Ｉ_８０／Ｉ_３が分布するものの、高濁度の場合には実測濁度と反射強度比Ｉ_８０／Ｉ_３には良好な相関関係が見られる。これは、高濁度時には、水中の超音波減衰が濁度レベルに大きく依存するためである。

なお、ここで反射強度Ｉ_８０を単独で用いずに、反射強度の水平方向勾配を示す反射強度比Ｉ_８０／Ｉ_３を用いたのは、浮遊土砂濃度ＳＳが横断方向に一様である保証がなく、濁度推定にはＨ−ＡＤＣＰ４０近傍の反射強度値を用いる必要があるためである。

高濁度条件における実測濁度Turbと反射強度比Ｉ_８０／Ｉ_３との相関式として、次式が得られた。

Ｉ_８０／Ｉ_３≦0.56の場合、
Turb = -5516(I₈₀/I₃)²+4973(I₈₀/I₃)-1001 ・・・（２）
0.56＜Ｉ_８０／Ｉ_３≦0.64の場合、
Turb = 7334(I₈₀/I₃)²-9333(I₈₀/I₃)+2979 ・・・（３）
これらの結果から、低濁度時（反射強度Ｉ_３≦２１０countの場合）には、式（１）を用い、高濁度時（反射強度Ｉ_３＞２１０count）の場合には、反射強度比I₈₀/I₃の値に応じて式（２）または（３）を用いて、反射強度から濁度の測定値としてTurbを求めればよいことがわかる。

コンピュータ１０では、式（１）〜（３）を用いて濁度の測定値を求めた後、後述する濃度分布の分析に用いるため、更にＳＳ濃度（単位容積あたりの浮遊物質量を示す、以下では浮遊土砂濃度ＳＳと呼称）に変換する。なお、実験により実測濁度Turbと浮遊土砂濃度ＳＳとの相関関係を調べた結果、以下に示す相関式（４）が得られた。

SS = 2.77Turb ・・・（４）
ＳＳの単位は、mg/Ｌである。コンピュータ１０は、上式（４）を用いて濁度（Turb)の測定値から浮遊土砂濃度(SS)の測定値に変換する。

以上説明した、計測した反射強度から濁度の測定値を求め、濁度の測定値から浮遊土砂濃度の測定値を求める処理を実行するためのプログラムや、上記式（１）〜（４）、及び各種パラメータの設定値（ここでは反射強度値(210)、反射強度比値(0.56、0.64)）は、予めＨＤＤ２２に記憶され設定されている。ＣＰＵ２０は、記憶されているプログラムを実行し、入力された反射強度データから最適な計算式を選択して用いて最終的に浮遊土砂濃度ＳＳを求める。

なお、本実施の形態で、Ｈ−ＡＤＣＰ４０から近い位置の観測値として第３層の反射強度を採用したのは、第１層や第２層のように、Ｈ−ＡＤＣＰ４０にあまりにも近すぎると、Ｈ−ＡＤＣＰ４０の動作の影響を受けやすいためである。また、Ｈ−ＡＤＣＰ４０から遠い位置の観測値として第８０層の反射強度を採用したのは、実験に採用した河川の幅がＨ−ＡＤＣＰ４０の観測可能距離よりも短いことから、対岸の岸の反射波を観測しないように調整したためである。従って、本発明において、Ｈ−ＡＤＣＰ４０から近い位置の観測距離は、良好な観測結果が得られる距離であれば第３層に限定されず、また、Ｈ−ＡＤＣＰ４０から遠い位置の観測距離は、良好な観測結果が得られる距離であれば第８０層に限定されない。

次に、コンピュータ１０では、上記求めた浮遊土砂濃度ＳＳの測定値に基づいて、観測対象である河川の横断面における浮遊土砂濃度の分布を解析して出力する。

ここでは、図６に示すように、水が主に流れる方向（主流方向）をｘ方向とし、ｘ方向に直交し河川を横断する方向をｙ方向とし、鉛直方向をσ方向とする座標系を定義し、該座標系の横断面（ｙーσ平面）の浮遊土砂濃度の分布を解析する。

この解析にあたり、コンピュータ１０は、土砂力学を満足した形でその分布を算出する。土砂力学を考える際に基礎となる三次元土砂輸送方程式は、浮遊土砂濃度ＳＳをＣとすると、以下の式（５）のように与えられる。

ここで、u、v、wは、x方向、y方向、σ方向の流速であり、w₀は浮遊土砂の沈降速度、Ａ_ＨとＡ_Ｖは水平、鉛直渦動粘性係数、Ｄは水深である。

しかしながら、式（５）中のいくつかの項は次の理由により省略する。

左辺第１項：大きさ自体が、他の項（左辺第４項など）よりも２、３オーダー小さいので、省略する。

左辺第２項：∂Ｃ／∂ｘを求めるには、ｘ方向の浮遊土砂濃度Ｃの変化を知る必要があるが、一断面でしか計測していないの∂で、∂Ｃ／∂ｘを求めることができない。従って、この項も省略する。

左辺第３項：本実施の形態のＨ−ＡＤＣＰ４０では、横断面全体の横断方向流速ｖを計測していない。よって、この第３項の評価は困難であり省略する。

左辺第４項の一部ｗ∂Ｃ／∂σ：左辺第３項と同じく、実施の形態のＨ−ＡＤＣＰ４０では、鉛直方向流速ｗを計測していないため、評価は困難であり省略する。

右辺第１項：左辺第２項と同じ同じ理由により算定困難であり、省略する。

これらの項を式（５）中から省略すると、簡略化された次式（６）が得られる。

このように、上式（６）では、１つの横断面内の計測では推定が困難な移流項の一部や主流方向拡散項を省略しているが、式（６）のままでは式（５）で表されるような土砂動態を再現できない。そこで、これらの省略された項の代わりに、土砂力学条件を満たすために、浮遊土砂用付加項Ｆa,cを式（６）に導入する。

式（７）の左辺は、単位時間あたりの沈降土砂量を表す項、右辺は、浮遊土砂の水平方向における拡散フラックスを表す項（右辺第１項）、浮遊土砂の鉛直方向の拡散フラックスを表わす項（右辺第２項）、及び付加項Fa,c（右辺第３項）の和を表している。

コンピュータ１０は、上記式（７）を用いて水中の土砂濃度の分布解析処理を行なう。

図７は、コンピュータ１０が実行する土砂濃度の分布解析処理の流れを示すフローチャートである。

ステップ１００では、水中の鉛直面内に計算格子の設定を行なう。具体的には、図８に示すように、解析対象となる河川の横断面（ｙーσ平面）のうちＨ−ＡＤＣＰ４０が設置されている横断面を計算領域とし、該計算領域に対して鉛直方向に等間隔でかつ水平方向に等間隔に複数の格子点を配列し、隣接する格子点で囲まれた矩形区域（以下計算格子と呼称）を複数備えた仮想横断面を設定する（図９も参照）。言い換えれば、横断面を格子状に分割し、分割して得られた矩形区域を計算格子として設定し、該計算格子を設定した横断面を仮想横断面とする。なお、以下では、ここで設定された計算格子のうち、前述したようにＨ−ＡＤＣＰ４０で実際に反射強度が計測されて浮遊土砂濃度ＳＳが算出された計測地点に対応する計算格子を、計測地点の計算格子と呼称する。

図９は、計算格子を模式的に示す図であり、格子点を通る格子線により横断面を格子状に分割して計算格子が設定されている。

また、以下では、演算対象の計算格子（以下、注目格子）を順にずらしながら全計算格子についての浮遊土砂濃度の演算を行なうが、注目格子の浮遊土砂濃度ＳＳを示すＣの値だけでなく周囲の計算格子のＣの値も用いて演算を行なう。なお、ここで、用いる周囲の計算格子は、注目格子のｙ−σ平面の座標値を（ｉ、ｊ）とすると、座標値（ｉ−１、ｊ）の計算格子（注目格子の直左の計算格子）、座標値（ｉ＋１、ｊ）の計算格子（注目格子の直右の計算格子）、座標値（ｉ、ｊ−１）の計算格子（注目格子の直下の計算格子）、及び座標値（ｉ、ｊ＋１）の計算格子（注目格子の直上の計算格子）である。

このように、計算領域及び計算格子を設定した後は、ステップ１０２で、流れの計算パラメータ（w₀,D,A_H,A_V）の設定を行なう。なお、沈降速度w₀は、実測された浮遊土砂の粒径とストークス近似等を用いて得られる値を設定する。また、Ａ_ＨとＡ_Ｖは、最も一般的な０方程式モデルなどの乱流モデルを用いて、摩擦速度や水深の関数として設定する。これらにより、計算格子毎のw₀,A_H,A_Vが定まる。またＤには、観測対象の河川の水深を設定する。

ステップ１０４では、計算格子のうち、計測地点の計算格子にはＨ−ＡＤＣＰ４０で計測されたＳＳの測定値をＣの値として設定し、他の計算格子には予め定められた初期値をＣの値として設定する。初期値としては、予め定められた値を用いることができるが、計測地点のＳＳの測定値をそのまま計算領域内の各計算格子の初期値として一様に与えるようにしてもよい。

次に、ステップ１０６では、計測地点の計算格子の付加項Ｆa,cの値を演算する。

付加項を演算するには上記式（７）を用いればよいが、式（７）で数値計算するに当たり、ここでは差分法を適用する。式（７）の差分方程式（８）は以下のようになる。

なお、各パラメータに与えられた添え字i,jは、ｙ−σ平面の座標値を表す。また、Ｃ_i,jは、注目格子のＣを示しており、Ｃ_i,j+1は、注目格子の直上の計算格子のＣを示し、Ｃ_i+1,jは注目格子の直右の計算格子のＣを示し、Ｃ_i-1,jは注目格子の直左の計算格子のＣを示し、Ｃ_i,j-1は注目格子の直下の計算格子のＣを示している（図９も参照。）。

ステップ１０８では、上記演算により得られた付加項Ｆ_a,cの値を用いて、他の計算格子の浮遊土砂濃度SSを示すＣを上記式（８）を用いて演算する。なお、本実施の形態では、付加項F_a,cは計算領域内で一様としている。これは、本実施の形態では、観測データが一点しか得られていないためである。

ステップ１１０では、注目格子のＣの値を、設定した初期値に代えて今回の演算で得られたＣの値を設定する。すなわち、注目格子の浮遊土砂濃度ＳＳの演算が１回目である場合には、初期値に代えて上記演算で得られたＣの値を設定する。また、当該注目格子の浮遊土砂濃度ＳＳの演算が２回目以降である場合には、前回の演算で得られたＣの値に代えて今回の演算により得られたＣの値を設定する。なお、計測地点の計算格子のＣの値は、変更せずに測定値のままとする。

ステップ１１２では、仮想横断面の全計算格子について浮遊土砂濃度ＳＳを示すＣの演算が終了したか否かを判断する。

ステップ１１２で、仮想横断面の全計算格子についてＣの演算が終了していないと判断した場合には、ステップ１０８に戻り、次の計算格子を注目格子として、Ｃの演算を行ない、ステップ１１０で、演算により得られたＣの値を設定する。

また、ステップ１１２で、仮想横断面の全計算格子についてＣの演算が終了したと判断した場合には、ステップ１１４で付加項F_a,c及びＣの値が変化したか否かを判断することにより収束したか否かを判断する。

ステップ１１４で、付加項F_a,c及びＣの値が変化しており収束していないと判断した場合には、ステップ１０６に戻り、計測地点の計算格子の付加項Fa,cを再度演算する。なお、２回目以降の演算は、全計算格子について前回演算されたＣの値が設定されているため、その状態で付加項F_a,cの演算を行なうと共に各計算格子について新たにＣを算出する。

ステップ１０８以降は、上記と同様に再度Ｃを算出して設定する。

一方、ステップ１１４で、付加項F_a,c及びＣの値が変化しておらず収束したと判断した場合には、ステップ１１６で、上記得られた浮遊土砂濃度ＳＳの分布を表示装置２８に表示する。

このように、付加項Ｆ_a,cと浮遊土砂濃度ＳＳを示すＣの解が収束するまで、ステップ１０６〜１１４までの処理が繰り返し行なわれ、最終的に浮遊土砂濃度ＳＳの分布が得られる。

すなわち、少なくとも１箇所で測定された測定結果さえあれば、上述のように力学条件を満たすように横断面全体に内外挿して濃度分布が解析されるため、容易に横断面全体の浮遊土砂濃度の分布を導出することができる。

なお、測定値が複数個得られている場合には、測定地点に対応する計算格子の各々にＳＳの測定値の各々を設定し、測定値が設定された計算格子以外の各計算格子のＣを上記のように演算すればよい。

また、Ｈ−ＡＤＤＰ４０の反射強度特性に基づいて、低濁度時にはＨ−ＡＤＣＰ４０近傍の反射強度値を、高濁度時にはＨ−ＡＤＣＰ４０近傍の反射強度値とＨ−ＡＤＣＰ４０から遠い反射強度値との比である反射強度比を各々用いることにより、反射強度データから精度良く濁度、ＳＳを計算することができる。上記ではＨ−ＡＤＣＰ４０の測定値を用いたが、濁度計を用いて測定した値を用いてもよい。

なお、江戸川における浮遊土砂輸送量計測に本システムを適用した結果、本システムによる推定結果は、別途行われた観測結果と良好に一致しており、本システムの有効性が検証された。

また、本発明は、上述した実施の形態に限定されるものではなく、特許請求の範囲に記載された発明の範囲内で様々な設計上の変更を行うことができる。

例えば、上記では、１つの注目格子のＣを演算した後、該演算結果を設定して次の注目格子のＣの演算を行ない、これを繰り返して仮想横断面の全計算格子（計測地点の計算格子を除く）についてＣを演算する例について説明したが（ステップ１０８〜１１２）、仮想横断面の計算格子についてＣを演算する手法はこれに限定されるものではなく、Ｃの演算が全計算格子について行なわれた後に該演算結果の各々を全計算格子に設定する（すなわち、Ｃの演算及び設定をまとめて行なう）手法を用いてもよい。このような手法を用いても、付加項Ｆ_a,cとＣの解が収束するまで、繰り返し演算が行なわれるため、上記と同様に浮遊土砂濃度分布が得られる。

（Ａ）は、本実施の形態に係る懸濁物質濃度の分布解析装置として機能するコンピュータの外観図であり、（Ｂ）は、該コンピュータのハードウェア構成を模式的に示すブロック図である。 河川にＨ−ＡＤＣＰを固定設置した様子を示す図である。 （Ａ）は、ある河川における観測期間中の水位と実測濁度との関係を示すグラフであり、（Ｂ）は、同じ観測期間中の、第３層の反射強度Ｉ_３及び第８０層の反射強度Ｉ_８０を示すグラフである。 実測濁度Turbと第３層目の反射強度Ｉ_３の相関関係の一例を示すグラフである。 第８０層の反射強度Ｉ_８０と第３層の反射強度Ｉ_３との比（Ｉ_８０／Ｉ_３）と、実測濁度Turbとの相関関係を示すグラフである。 浮遊土砂濃度の分布解析対象となる河川の座標系を説明する説明図である。 コンピュータが実行する土砂濃度の分布解析処理の流れを示すフローチャートである。 計算格子が設定された仮想横断面を示す図である。 計算格子を模式的に示す図である。

符号の説明

１０ コンピュータ
１２ ＣＰＵ
１４ ＲＡＭ
１６ ＲＯＭ
２２ ＨＤＤ
２８ 表示装置
３０ 入力装置
４０ Ｈ−ＡＤＣＰ

Claims (8)

1. 懸濁物質濃度の分布を解析する水中の横断面内に複数の格子点を配列することにより、隣接する格子点で囲まれた区域を複数備えた仮想横断面を想定する想定手段と、
   測定地点に対応する前記区域に懸濁物質濃度の測定値を設定すると共に、前記測定値を設定した区域以外の区域の各々に懸濁物質濃度の初期値を設定する設定手段と、
   懸濁物質の沈降量を表す項と、懸濁物質の水平方向の拡散フラックスを表わす項、懸濁物質の鉛直方向の拡散フラックスを表わす項、及び付加項の和との関係を前記区域の各々について表わした基礎式、前記測定地点に対応する区域に設定された懸濁物質濃度の測定値、及び前記測定地点に対応する区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記測定地点に対応する区域の付加項の値を演算する第１の演算、並びに前記測定地点に対応する区域以外の区域の１つを注目区域として、前記基礎式、前記演算された付加項の値、前記注目区域に設定された懸濁物質濃度の値、及び前記注目区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記注目区域の懸濁物質濃度を演算し、演算により得られた演算値を前記注目区域に既に設定された値に代えて設定することを、前記注目区域を代えながら前記仮想横断面全体にわたって繰返し行う第２の演算を、前記懸濁物質濃度の演算値及び前記付加項の演算値の各々が収束するまで繰返し行う繰返演算手段と、
   を含む懸濁物質濃度の分布解析装置。
2. 前記基礎式を以下の式で表した請求項１記載の懸濁物質濃度の分布解析装置。
   

   

   なお、ｙは水平方向、σは鉛直方向、Ａ_Ｈは水平渦動粘性係数、Ａ_Ｖは鉛直渦動粘性係数、Ｄは水深、Ｃは懸濁物質濃度、ｗ_０は懸濁物質の沈降速度、Ｆa,cは、付加項である。
3. 所定位置に設置された超音波ドップラー流速分布計により水中に発信され該水中の懸濁物質により反射され計測された超音波の反射強度のうち、前記超音波ドップラー流速分布計近傍の懸濁物質群で反射された超音波の反射強度を示す第１反射強度と、前記超音波ドップラー流速分布計から離れた懸濁物質群で反射された超音波の反射強度を示す第２反射強度と、に基づいて前記懸濁物質濃度の測定値を算出する濃度測定手段を更に含む
   請求項１または請求項２に記載の懸濁物質濃度の分布解析装置。
4. 前記濃度測定手段は、前記第１反射強度が所定値以下の場合には、前記第１反射強度に基づいて前記懸濁物質濃度の測定値を算出し、前記第１反射強度が所定値を超える場合には、前記第１反射強度と前記第２反射強度との比に基づいて前記懸濁物質濃度の測定値を算出する
   請求項３記載の懸濁物質濃度の分布解析装置。
5. コンピュータを、
   懸濁物質濃度の分布を解析する水中の横断面内に複数の格子点を配列することにより、隣接する格子点で囲まれた区域を複数備えた仮想横断面を想定する想定手段と、
   測定地点に対応する前記区域に懸濁物質濃度の測定値を設定すると共に、前記測定値を設定した区域以外の区域の各々に懸濁物質濃度の初期値を設定する設定手段と、
   懸濁物質の沈降量を表す項と、懸濁物質の水平方向の拡散フラックスを表わす項、懸濁物質の鉛直方向の拡散フラックスを表わす項、及び付加項の和との関係を前記区域の各々について表わした基礎式、前記測定地点に対応する区域に設定された懸濁物質濃度の測定値、及び前記測定地点に対応する区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記測定地点に対応する区域の付加項の値を演算する第１の演算、並びに前記測定地点に対応する区域以外の区域の１つを注目区域として、前記基礎式、前記演算された付加項の値、前記注目区域に設定された懸濁物質濃度の値、及び前記注目区域の周囲に存在する複数の区域に設定された懸濁物質濃度の値を用いて、前記注目区域の懸濁物質濃度を演算し、演算により得られた演算値を前記注目区域に既に設定された値に代えて設定することを、前記注目区域を代えながら前記仮想横断面全体にわたって繰返し行う第２の演算を、前記懸濁物質濃度の演算値及び前記付加項の演算値の各々が収束するまで繰返し行う繰返演算手段と、
   して機能させるためのプログラム。
6. 前記基礎式を以下の式で表した請求項５記載のプログラム。
   

   

   なお、ｙは水平方向、σは鉛直方向、Ａ_Ｈは水平渦動粘性係数、Ａ_Ｖは鉛直渦動粘性係数、Ｄは水深、Ｃは懸濁物質濃度、ｗ_０は懸濁物質の沈降速度、Ｆa,cは、付加項である。
7. 更に、前記コンピュータを、所定位置に設置された超音波ドップラー流速分布計により水中に発信され該水中の懸濁物質により反射され計測された超音波の反射強度のうち、前記超音波ドップラー流速分布計近傍の懸濁物質群で反射された超音波の反射強度を示す第１反射強度と、前記超音波ドップラー流速分布計から離れた懸濁物質群で反射された超音波の反射強度を示す第２反射強度と、に基づいて前記懸濁物質濃度の測定値を算出する濃度測定手段として機能させるための請求項５または請求項６に記載のプログラム。
8. 前記濃度測定手段は、前記第１反射強度が所定値以下の場合には、前記第１反射強度に基づいて前記懸濁物質濃度の測定値を算出し、前記第１反射強度が所定値を超える場合には、前記第１反射強度と前記第２反射強度との比に基づいて前記懸濁物質濃度の測定値を算出する
   請求項７記載のプログラム。

Images