特許 7597596 谷地形の抽出方法及び水系の抽出方法並びにプログラム

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

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

(11)【特許番号】

(24)【登録日】2024-12-02

(45)【発行日】2024-12-10

(54)【発明の名称】谷地形の抽出方法及び水系の抽出方法並びにプログラム

(51)【国際特許分類】

   G06T 17/05 20110101AFI20241203BHJP

   G01C 7/02 20060101ALI20241203BHJP

   G09B 29/00 20060101ALI20241203BHJP

【ＦＩ】

G06T17/05

G01C7/02

G09B29/00 Z

【請求項の数】 10

(21)【出願番号】P 2021011819

(22)【出願日】2021-01-28

(65)【公開番号】P2022115289

(43)【公開日】2022-08-09

【審査請求日】2024-01-17

(73)【特許権者】

【識別番号】000173809

【氏名又は名称】一般財団法人電力中央研究所

(74)【代理人】

【識別番号】100101236

【弁理士】

【氏名又は名称】栗原 浩之

(74)【代理人】

【識別番号】100166914

【弁理士】

【氏名又は名称】山▲崎▼ 雄一郎

(72)【発明者】

【氏名】平田 康人

【審査官】渡部 幸和

(56)【参考文献】

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

【文献】特開２０１２－００３４００（ＪＰ，Ａ）

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

Ｇ０６Ｔ １７／００

Ｇ０１Ｃ ７／００

Ｇ０９Ｂ ２９／００

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

【請求項1】

地表面を一辺がＬｍｉｎである正方形のセルで区切り、各セルに標高が割り当てられた数値標高モデル（Ｌｍｉｎ）を用いた谷地形の抽出方法であって、
数値標高モデル（Ｌｍｉｎ）に含まれるｎ×ｎ個のセルについて平均標高を計算し、一辺の長さがＬ（Ｌｍｉｎ×ｎ）である正方形のセルの標高として前記平均標高を割り当てることで低解像度の数値標高モデル（Ｌ）を作成する集約ステップと、
数値標高モデル（Ｌ）に含まれる３×３のセル（以下セル群）を抽出し、抽出した全てのセル群について中央のセルの平面曲率を計算する工程を、前記集約ステップで作成した複数の数値標高モデル（Ｌ）のそれぞれについて実行する平面曲率計算ステップと、
前記平面曲率にセルの辺の長さＬを乗じた正規平面曲率を計算する工程を、前記集約ステップで作成した複数の数値標高モデル（Ｌ）のそれぞれについて実行する正規平面曲率計算ステップと、
前記セル群の中央のセルが凹地領域と判定されるための正規平面曲率に対する第１閾値であって複数の数値標高モデル（Ｌ）に共通の第１閾値を設定し、前記第１閾値以下となる前記正規平面曲率のセルを凹地領域と判定する凹地判定ステップと、
複数の数値標高モデル（Ｌ）ごとに抽出された凹地領域を重ね合わせ、連続する凹地領域からなるポリゴンを作成し、前記ポリゴンの谷線方向の長さが第２閾値以上となる前記ポリゴンを谷領域として抽出する谷領域抽出ステップと、を備える
ことを特徴とする谷地形の抽出方法。

【請求項2】

請求項１に記載する谷地形の抽出方法であって、
前記谷領域抽出ステップは、複数の前記数値標高モデルごとに抽出した凹地領域を重ね合わせて、それらの前記数値標高モデルの中で最小の傾斜角が第３閾値以上であるセルを重ね合わせた凹地領域から抽出し、抽出した連続するセルからなる前記ポリゴンを作成する
ことを特徴とする谷地形の抽出方法。

【請求項3】

請求項１又は請求項２に記載する谷地形の抽出方法であって、
前記谷領域抽出ステップは、前記ポリゴンの面積を前記ポリゴンに含まれる最大のセルの長さの二乗で割った値が前記第２閾値である２．４以上であるならば前記ポリゴンを谷領域として抽出する
ことを特徴とする谷地形の抽出方法。

【請求項4】

請求項１から請求項３の何れか一項に記載する谷地形の抽出方法であって、
前記凹地判定ステップでは、第１閾値の絶対値として式（１）を用いる
ことを特徴とする谷地形の抽出方法。

【数1】

ただし、Ｋｐｌａｎ，ｎｏｒｍａｌは第１閾値、ｍは整数、ｋは谷線方向の高低差ａと、谷に直交する方向の標高差の倍率である。

【請求項5】

請求項１から請求項３の何れか一項に記載する谷地形の抽出方法であって、
前記集約ステップでは、前記凹地判定ステップにより得られた凹地領域が尾根を超えないように、集約に用いるセルの一辺の長さＬを定める
ことを特徴とする谷地形の抽出方法。

【請求項6】

請求項１から請求項５の何れか一項に記載する谷地形の抽出方法によって得られた谷領域を用いた水系の抽出方法であって、
前記谷領域を一辺がｌの正方形の解析用セルで分割し、各解析用セルに一定値を与えて累積流量を計算し、累積流量が第４閾値を越える解析用セルからなる水系を抽出する水系抽出ステップを備え、
前記第４閾値は、Ｖ字谷地形の谷頭における累積流量である
ことを特徴とする水系の抽出方法。

【請求項7】

請求項６に記載の水系の抽出方法であって、
前記水系抽出ステップは、谷次数で区切られた谷の長さが第５閾値より長い解析用セルからなる水系を抽出する
ことを特徴とする水系の抽出方法。

【請求項8】

請求項７に記載の水系の抽出方法であって、
前記水系抽出ステップは、水系を細線化し、水系の最も上流から再度累積流量を計算し、累積流量が前記第４閾値を越える解析用セルからなる水系を抽出する
ことを特徴とする水系の抽出方法。

【請求項9】

請求項１から請求項５の何れか一項に記載の谷地形の抽出方法の各ステップをコンピュータに実行させる
ことを特徴とするプログラム。

【請求項10】

請求項６から請求項８の何れか一項に記載の水系の抽出方法の各ステップをコンピュータに実行させる
ことを特徴とするプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、数値標高モデル（ＤＥＭ；Ｄｉｇｉｔａｌ Ｅｌｅｖａｔｉｏｎ Ｍｏｄｅｌ）から谷地形及び水系を抽出する方法並びにプログラムに関する。

【背景技術】

【0002】

従来、山地の地形表現図を得たり、斜面災害のリスク評価のために谷地形や水系図を数値標高モデルから抽出することが行われている。

【0003】

水系図は、地図上の等高線の曲がり具合を人の目で判読し、谷地形と判断される最も上流側の位置から下流に向かって谷線を連続させることにより作成されている。谷地形と判断する基準は、等高線の折入角が９０度又は５３度とすることが一般的である。

【0004】

このような方法においては、谷地形の抽出結果は作成者に依存してしまうため、作成者によって異なる谷地形や水系図が得られることになる。また、広範囲の水系図の作成には相応の時間を要してしまう。

【0005】

このような作成方法とは別に、数値標高モデルを用いて上流側の集水面積を計算し、ある閾値以上の集水面積を谷地形と、水系図を得る方法もある。この閾値は作成者により任意に決定される。

【0006】

このような方法においては、深い谷でなくとも集水量が増えれば水系とみなすことになるので、厳密には谷地形の基準に従った水系図ではない。また、人工構造物がなす水平面も水系図に含まれてしまう。

【0007】

さらに、上述した２つの方法の何れについても、地図の縮尺によって作成される水系図が異なるという問題がある。

【0008】

特許文献１に記載の技術では、集水地形とされる斜面などを数値標高モデルから抽出するが、上述したような問題を解決できるものではなかった。

【先行技術文献】

【特許文献】

【0009】

【文献】特開２０１３－２２１８８６号公報

【非特許文献】

【0010】

【文献】河田 佳樹, 仁木 登. 2001. 3次元曲率特徴の抽出アルゴリズム. ＭEDICAL IMAGING TECHNOLOGY 19: 142-153

【発明の概要】

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

【0011】

本発明は、上記事情に鑑み、地図の縮尺や作成者の主観によらず客観的に谷地形及び水系を作成することができる谷地形の抽出方法及び水系の抽出方法並びにプログラムを提供することを目的とする。

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

【0012】

上記目的を達成するための本発明の態様は、地表面を一辺がＬｍｉｎである正方形のセルで区切り、各セルに標高が割り当てられた数値標高モデル（Ｌｍｉｎ）を用いた谷地形の抽出方法であって、数値標高モデル（Ｌｍｉｎ）に含まれるｎ×ｎ個のセルについて平均標高を計算し、一辺の長さがＬ（Ｌｍｉｎ×ｎ）である正方形のセルの標高として前記平均標高を割り当てることで低解像度の数値標高モデル（Ｌ）を作成する集約ステップと、数値標高モデル（Ｌ）に含まれる３×３のセル（以下セル群）を抽出し、抽出した全てのセル群について中央のセルの平面曲率を計算する工程を、前記集約ステップで作成した複数の数値標高モデル（Ｌ）のそれぞれについて実行する平面曲率計算ステップと、前記平面曲率にセルの辺の長さＬを乗じた正規平面曲率を計算する工程を、前記集約ステップで作成した複数の数値標高モデル（Ｌ）のそれぞれについて実行する正規平面曲率計算ステップと、前記セル群の中央のセルが凹地領域と判定されるための正規平面曲率に対する第１閾値であって複数の数値標高モデル（Ｌ）に共通の第１閾値を設定し、前記第１閾値以下となる前記正規平面曲率のセルを凹地領域と判定する凹地判定ステップと、複数の数値標高モデル（Ｌ）ごとに抽出された凹地領域を重ね合わせ、連続する凹地領域からなるポリゴンを作成し、前記ポリゴンの谷線方向の長さが第２閾値以上となるポリゴンを谷領域として抽出する谷領域抽出ステップと、備えることを特徴とする谷地形の抽出方法にある。

【0013】

上記目的を達成するための本発明の他の態様は、上記谷地形の抽出方法によって得られた谷領域を用いた水系の抽出方法であって、前記谷領域を一辺がｌの正方形の解析用セルで分割し、各解析用セルに一定値を与えて累積流量を計算し、累積流量が第４閾値を越える解析用セルからなる水系を抽出する水系抽出ステップを備え、前記第４閾値は、Ｖ字谷地形の谷頭における累積流量であることを特徴とする水系の抽出方法にある。

【0014】

上記目的を達成するための本発明の他の態様は、上記谷地形の抽出方法の各ステップをコンピュータに実行させることを特徴とするプログラムにある。

【0015】

上記目的を達成するための本発明の他の態様は、上記水系の抽出方法の各ステップをコンピュータに実行させることを特徴とするプログラムにある。

【発明の効果】

【0016】

本発明によれば、地図の縮尺や作成者の主観によらず客観的に谷地形及び水系を作成することができる谷地形の抽出方法及び水系の抽出方法並びにプログラムが提供される。

【図面の簡単な説明】

【0017】

【図1】谷地形の抽出方法及び水系の抽出方法のフロー図である。

【図2】ＤＥＭ（１）及びＤＥＭ（２）の一部を示す図である。

【図3】ＤＥＭ（１）及びＤＥＭ（２）の一部を示す図である。

【図4】ＤＥＭの３×３のセルの標高を示す図である。

【図5】Ｖ字谷地形における平面曲率及び正規平面曲率と谷地形の関係を示す図である。

【図6】ｎが奇数である場合におけるＤＥＭの集約を説明するための図である。

【図7】ｎが偶数である場合におけるＤＥＭの集約を説明するための図である。

【図8】ｎが偶数である場合におけるＤＥＭの集約を説明するための図である。

【図9】ｎが偶数である場合におけるＤＥＭの集約を説明するための図である。

【図10】凹地領域の実施例を示す図である。

【図11】谷領域と傾斜を示す図である。

【図12】長さに基づいて谷領域を抽出する方法を説明するための図である。

【図13】Ｖ字谷地形における累積流量を示す図である。

【図14】水系から短い谷を除去する方法を説明するための図である。

【図15】水系を細線化する方法を説明するための図である。

【図16】ＤＥＭから谷地形及び水系の抽出方法により得られた実施例を示す図である。

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

【0018】

本発明の実施形態に係る谷地形の抽出方法及び水系の抽出方法について説明する。図１は谷地形の抽出方法及び水系の抽出方法のフロー図である。まず、谷地形の抽出方法では、基準となる最も高解像度の数値標高モデルから、それよりも低解像度の数値標高モデルを作成する。

【0019】

数値標高モデルとは、地表面を一辺がＬｍｉｎである正方形のセルで区切り、各セルに標高が割り当てられたデータである。本実施形態ではＬｍｉｎが１ｍである場合について説明する。以後、セルの長さＬ（メートル）の数値標高モデルをＤＥＭ（Ｌ）と表記する。

【0020】

図２（ａ）に、ＤＥＭ（１）の一部を例示する。各セルには標高が割り当てられており、各セルの標高をｈｉｊ（ｉ＝１，２，・・・；ｊ＝１，２，・・・）とする。このようなＤＥＭ（１）に含まれるｎ×ｎこのセルについて平均標高を計算し、一辺の長さがＬ（Ｌｍｉｎ×ｎ）である正方形のセルの標高としてその平均標高を割り当てる。このようにして得られる長さＬのセルからなる低解像度の数値標高モデル（Ｌ）を作成する（集約ステップ；図１のステップＳ１）。

【0021】

例えば、図２（ｂ）に、ｎを２とし、長さＬが２であるＤＥＭ（２）を示す。まず、ＤＥＭ（１）に含まれる２×２個のセルについて平均標高を計算する。例えば、左上の２×２個のセルの標高ｈ１１、ｈ１２、ｈ２１、ｈ２２の平均標高をＨ１１とする。

【0022】

そして、２×２個のセルと同じ辺の長さである２ｍを有する一つのセルの標高として上述した平均標高Ｈ１１を割り当てる。同様にして、ＤＥＭ（１）に含まれる２×２個のセルについて平均標高を計算し、長さ２ｍのセルの標高とすることでＤＥＭ（２）が作成される。

【0023】

図２の例では、ＤＥＭ（２）のセルの境界（辺）は、ＤＥＭ（１）の境界と一致する場合について説明したが、このような低解像度化に限定されない。

【0024】

図３（ａ）はＤＥＭ（２）のセル（点線で示されている）の境界がＤＥＭ（１）のセルの境界に対して半分だけズレて配置された状態を表している。図３（ｂ）は図３（ａ）のＤＥＭ（２）を抜き出して表したものである。

【0025】

ＤＥＭ（２）は、ＤＥＭ（１）の各辺から０．５（ＤＥＭ（１）の辺の長さの半分）だけズレている。このように低解像度、高解像度のセルの境界が一致しない場合は、ＤＥＭ（２）の一つのセルに重なるＤＥＭ（１）の複数のセルの標高から平均標高を計算する。

【0026】

例えば、ＤＥＭ（２）の左上のセルの標高Ｈ１１は、そのセルが重なっているＤＥＭ（１）の９個のセルの標高ｈ１１、ｈ１２、ｈ１３、ｈ２１、ｈ２２、ｈ２３、ｈ３１、ｈ３２、ｈ３３から求めれば良い。詳細な計算は、図９で説明するので後述する。

【0027】

なお、特に図示しないが、境界のズレがＤＥＭ（１）のセル長の半分以外である場合は、ＤＥＭ（１）とＤＥＭ（２）の境界のズレた位置に応じて適宜加重平均を計算すればよい。

【0028】

このようにして、境界の位置に応じて図２（ｂ）や図３に示したように異なる複数のＤＥＭ（２）がＤＥＭ（１）から得られることになる。本実施形態では、ＤＥＭ（１）から、境界が一致するＤＥＭ（２）、及び境界がＤＥＭ（１）のセルの長さの半分だけズレたＤＥＭ（２）を作成することとする。

【0029】

低解像度のＤＥＭは、１つに限らず、複数個作成することが好ましい。後述する実施例では、ＤＥＭ（１）を元にして、１ｍごとにＤＥＭ（２）～ＤＥＭ（５００）を作成した。また、それぞれのＤＥＭ（Ｌ）は、上述したように境界が一致するものと、境界がＤＥＭ（１）のセルの長さの半分だけズレた２パターンを作成した。

【0030】

次に、複数のＤＥＭ（Ｌ）について平面曲率の計算を行う（平面曲率計算ステップ；図１のステップＳ２）。平面曲率をＫｐｌａｎ、横断面曲率をＫｃｒｏｓｓ、谷線上の最大傾斜角をθ、図４に示す一辺の長さがｄの３×３のセル（以下「セル群」とも称する）の各標高をＺ１からＺ９とすると、次のようにして求めることができる。Ｋｐｌａｎ、Ｋｃｒｏｓｓの計算は、非特許文献１に記載の公知の方法であるので詳細な説明は省略する。

【数1】

【0031】

平面曲率及び横断面曲率は、何れも斜面方位に対して垂直な方位に延びる断面線の曲率の一つである。平面曲率は、最大傾斜角法項に対して直角方向の曲率である。横断面曲率は、曲率を調べる断面を標高方向と平行にとった最大傾斜角方向の曲率である。平面曲率及び横断面曲率は、負の値であれば上方に凹型の地形であり、正の値であれば上方に凸型の地形であることを表す。

【0032】

このように、平面曲率計算ステップでは、あるＤＥＭに含まれるセル群を抽出し、抽出した全てのセル群について上述の式（１）を用いてセル群の中央のセルの平面曲率を計算する。そしてこの計算を集約ステップで作成した複数のＤＥＭの全てについて実行する。

【0033】

次に、正規平面曲率を計算する（正規平面曲率計算ステップ；図１のステップＳ３）。正規平面曲率は、平面曲率にセルの辺の長さを乗じたものである。具体的には、ＤＥＭ（Ｌ）に含まれるセルの平面曲率にセルの辺の長さＬを乗じて正規平面曲率を計算し、この計算を集約ステップで作成したＤＥＭの全てについて実行する。

【0034】

ここで、図５を用いて、平面曲率及び正規平面曲率と谷地形との関係について説明する。図５（ａ）～（ｃ）は、ある解像度のＤＥＭに含まれる一つのセル群とその標高を表している。図５（ｄ）～（ｆ）は図５（ａ）～（ｃ）の標高差及び高低差を表している。図５（ｇ）は図５（ａ）～（ｃ）のＸ方向の断面を表している。図５（ｈ）は図５（ａ）～（ｃ）の谷線におけるＹ方向の断面を表している。図５（ｉ）は図５（ａ）を水平方向及び鉛直方向を２倍にしたときのＸ方向の断面を表している。

【0035】

最も単純なＶ字谷地形は、同じ傾斜角、同じ高さを持つ２つの斜面が谷線で接合する形態である。このようなＶ字谷地形は、図５（ａ）から図５（ｃ）に示すセル群において標高の特有のパターンで表現される。図５（ａ）から図５（ｃ）では、点線は等高線を示しており、各セル内の数字及び記号は標高を示している。ｉ行ｊ列のセルをＣｉｊとも表記する。また、セル群のうち中央に位置してＹ方向に沿った３つのセル（Ｃ１２、Ｃ２２、Ｃ３２）が谷線を含むセルとなっている。

【0036】

セル内の全ての等高線の形状が一様になるためには、第一に、谷線となる中央のセルとそれに隣接するセルの標高差は一定であることが必要である。第二に、等高線の間隔が一定になるためには、セル間の流下方向の高低差は一定となる必要がある。

【0037】

第一の条件を図示すると、図５（ａ）のパターンは図５（ｄ）に示すように、中央のセルであるＣ１２は、両隣のＣ１１及びＣ１３との標高差は２ａである。中央のセルＣ２２、Ｃ２３についても両隣のセルとの標高差は２ａである。つまり、標高差は一定である。同様に、図５（ｂ）のパターンは図５（ｅ）に示すように標高差はａで一定であり、図５（ｃ）のパターンは図５（ｆ）に示すように標高差は１／２ａで一定である。

【0038】

第二の条件を図示すると、図５（ａ）のパターンは図５（ｄ）に示すように、上流側のセルであるＣ２１と、下流側のセルであるＣ１１との高低差はａである。その他の上流側と下流側の高低差についても高低差はａである。つまり、高低差は一定である。同様に、図５（ｂ）のパターンは図５（ｅ）に示すように高低差はａで一定であり、図５（ｃ）のパターンは図５（ｆ）に示すように高低差はａで一定である。

【0039】

図５（ａ）から図５（ｃ）に示したＶ字谷形状は、図５（ｇ）及び図５（ｈ）に示すように谷線方向（Ｙ方向）の高低差ａと、谷に直交する方向（Ｘ方向）の標高差の倍率ｋによって決定される。ここで、ｋは斜面の標高が谷線に向かって減少する場合を正とする。つまり、Ｖ字谷形状の場合はｋは正である。

【0040】

図５（ａ）～図５（ｃ）に示すＶ字谷形状のｋはそれぞれ２、１、１／２である。また、谷線上の最大傾斜角θは、ｔａｎ^－１（ａ／Ｌ）であり、側壁の傾斜角Φはｔａｎ^－１（ｋａ／Ｌ）である。図５（ａ）～図５（ｃ）の何れのＶ字谷地形は、中央のセルにおける最大傾斜方向の曲率は０であり、横断面曲率Ｋｃｒｏｓｓは、中心のセルと左右のセルとの標高差、及びセルの幅Ｌから計算される。したがって、図５（ａ）～図５（ｃ）に示すようなパターンにおける平面曲率は式（２）のように表せる。

【数2】

【0041】

例えば、図５（ａ）のように等高線の開口幅と奥行きが等しいとき（等高線の折入れ角が約５３度であるとき）、平面曲率Ｋｐｌａｎはｋ＝２として－４／Ｌとなる。図５（ｂ）のように等高線の折入れ角が９０度であるとき、平面曲率Ｋｐｌａｎはｋ＝１として－２／Ｌとなる。これらの平面曲率はＤＥＭを平面上で回転させても同じ値を持つ。

【0042】

上述の式（２）の通り、平面曲率Ｋｐｌａｎはセルの長さＬに依存している。例えば、図５（ｉ）に示すように谷の形状を水平と鉛直方向へ２倍に拡大すると、平面曲率は２分の１となる。一方、平面曲率にセルの幅を掛けた値は一定となる。この値を正規平面曲率と称し、Ｋｐｌａｎ，ｎｏｒｍａｌとも表記する。

【0043】

正規平面曲率の絶対値は、谷の幅やセルサイズに依存せず、セルの中心がＶ字谷地形の中央にあるときに大きくなり、谷線から離れるほどに小さくなる。この性質を利用すると、セル群を用いて正規平面曲率を異なるセルの長さで計算することにより、異なる縮尺のＶ字谷地形を一つの閾値を用いて抽出することができる。

【0044】

例えば、図５（ａ）、図５（ｂ）のパターンのＶ字谷地形において、異なるセルの長さＬ（異なる縮尺）ごとにＫｐｌａｎ，ｎｏｒｍａｌを計算すると次のようになる。
Ｌ Kplan,normal
１ －４ －２
２ －４ －２
３ －４ －２
・・・

【0045】

あるセルについてのＫｐｌａｎ，ｎｏｒｍａｌが－４以下であれば、どの様なセルの長さ、縮尺であっても、当該セルは図５（ａ）のようなパターンのＶ字谷地形、又はそれよりも折入角が小さいＶ字谷地形（厳密にはセル１つであるので凹地である）として抽出することができる。同様に、あるセルについてのＫｐｌａｎ，ｎｏｒｍａｌが－２以下であれば、どの様なセルの長さ、縮尺であっても、当該セルは図５（ｂ）のようなパターンのＶ字谷地形、又はそれよりも折入角が小さいＶ字谷地形として抽出することができる。

【0046】

このような性質の正規平面曲率に基づいて凹地判定を行う（凹地判定ステップ；図１のステップＳ４）。具体的には上述の「２」や「４」のようにセル群の中央のセルが谷地形と判定されるための正規表面曲率に対する第１閾値を適宜設定する。第１閾値は、上述したように複数のＤＥＭ（Ｌ）に共通なものである。ＤＥＭに含まれる各セルの正規平面曲率が第１閾値以下となるものを谷地形（凹地）として判定する。そして、同じ第１閾値を用いて、このような谷地形の抽出を全てのＤＥＭ（Ｌ）について行う。

【0047】

ただし、基準となるＤＥＭ（本実施形態ではＤＥＭ（１））から低解像度のＤＥＭを得る場合、平均化した際に用いたセルの辺の長さ（平均化した際に用いたｎ×ｎのセルにおけるｎで決まる長さ）によって正規平面曲率は変化する。なぜならば、Ｖ字谷地形の中央のセルと、周囲の壁面との比高とセルの長さとの比が小さくなるからである。

【0048】

Ｌｍｉｎ＝１としたＤＥＭ（１）について、ｎ×ｎ個のセルについて平均標高を計算して新たな低解像度のＤＥＭを作成したとき、Ｖ字谷地形を示すセル群の正規平面曲率の分布は、ｎが偶数の場合と奇数の場合とで異なる。

【0049】

図６を用いて、ｎが奇数である場合について説明する。図６（ａ）は、セルＣ５５を中心とし、谷線が第５列（ｊ＝５）に沿って延びるＶ字谷地形のＤＥＭ（１）を表している。図６（ｂ）は、ＤＥＭ（１）をｎ×ｎ個のセルで集約したＬ＝ｎ×ＬｍｉｎのＤＥＭ（Ｌ）を表している。例えば、ＤＥＭ（Ｌ）のＣ１１の標高は、（ｋ－１）ａｎとなり、ｎを３とすれば３（ｋ－１）ａとなる。なおＤＥＭ（Ｌ）のＣ１１は、ＤＥＭ（１）の太線で囲まれた３×３個のＣ１１、Ｃ１２、Ｃ１３、Ｃ２１、Ｃ２２、Ｃ２３、Ｃ３１、Ｃ３２、Ｃ３３の標高を平均化したものである。ＤＥＭ（Ｌ）の他のセルについても同様である。

【0050】

図７から図９を用いて、ｎが偶数である場合について説明する。図７の実線で示す９×９のセルは図６（ａ）と同じであるが、標高は一部のみを表示していある。点線で示す３×３のセルは、ＤＥＭ（１）を２×２のセルで集約したＬ＝２（２×Ｌｍｉｎ）のＤＥＭ（２）を表しており、そのＤＥＭ（２）は谷線が第２列（ｊ＝２）に沿って延びるＶ字谷地形を表している。ＤＥＭ（２）は、ＤＥＭ（１）に対して、セルの長さＬ＝１の半分ズレて配置されている。

【0051】

図８（ａ）は図７に示したＤＥＭ（２）を抜き出したものであり、ＤＥＭ（１）のセルから求めた標高が表記されている。図８（ｂ）は、図８（ａ）を一般化したもの、すなわちｎが偶数である場合についての低解像度化したＤＥＭ（Ｌ）を表している。

【0052】

なお、図８（ａ）の標高の計算は、次のように行っている。図９（ａ）の実線で示す９×９のセルは図６（ａ）と同じであるが、標高は一部のみを表示している。図９（ａ）の点線は、図（７）に示したＤＥＭ（２）に重なる位置に配置した６×６のセルからなるＤＥＭ（１’）を表している。図９（ａ）の点線で示すＤＥＭ（１’）は、実線のＤＥＭ（１）からセルの長さの半分だけズレて配置されている。

【0053】

図９（ｂ）に点線で示すＤＥＭ（１’）の各セルＣＩＪ（Ｉ行Ｊ列で示す）の標高は、セルＣＩＪに重なる、図９（ａ）に示す実線のＤＥＭ（１）の４個のセルの標高を平均化したものである。例えば、図９（ｂ）に示すＤＥＭ（１’）のセルＣ１１（Ｉ＝１、Ｊ＝１）は、図９（ａ）に示すＤＥＭ（１）の４個のセルＣ２２（ｉ＝２、ｊ＝２）、Ｃ２３（ｉ＝２、ｊ＝３）、Ｃ３２（ｉ＝３、ｊ＝２）、Ｃ３３（ｉ＝３、ｊ＝３）に重なっている。ＤＥＭ（１’）のセルＣ１１の標高は、ＤＥＭ（１）の４個のセルＣ２２、Ｃ２３、Ｃ３２、Ｃ３３の標高を平均化した（ｋ－１）５ａ／２となる。

【0054】

そして、このような点線で示すＤＥＭ（１’）に含まれる２×２のセルについて平均標高を求めることで、図８（ａ）に示したＤＥＭ（２）の標高が得られる。例えば、図９（ｂ）のＤＥＭ（１’）のセルＣ１１、Ｃ１２、Ｃ２１、Ｃ２２の標高の平均を求めると、２（ｋ－１）ａとなり、図８（ａ）のＤＥＭ（２）のセルＣ１１の標高となる。

【0055】

なお、図９に示したようなＤＥＭ（１）から、セルの長さの半分ズレて配置されたＤＥＭ（２）を計算する方法は、図３の標高Ｈｉｊを求める場合についても同様に適用できる。

【0056】

上述した図６（ｂ）や図８（ｂ）に示したように、低解像度化する際に用いるｎ×ｎのセルのｎが偶数か奇数かによって標高は異なる。このような偶数奇数の相違を考慮にいれて正規平面曲率を補正すると、Ｖ字谷地形の中央のセルにおける正規平面曲率は次のように表される。

【数3】

【0057】

ｍは整数である。ｎが奇数の場合、正規平面曲率の絶対値は集約する個数に依存し、最小値１．５ｋにすぐに収束する。例えば、Ｌ＝１であり、ｋ＝１であるＶ字谷地形は、Ｌ＝３に集約すると約１．５５のＫｐｌａｎの値を示す。一方、ｎが偶数の場合には、最小セルサイズＬｍｉｎの２つのセルの間に谷線が通り、谷線を挟む２つの標高が一致し、正規平面曲率は常に１．５ｋとなる。

【0058】

例えば、折入角が９０度であるとき、つまりｋ＝１としたとき、式（３）から得られる正規平面曲率の絶対値はｎが偶数の場合１．５であり、奇数の場合も最小値１．５×１に収束する。したがって、折入角が９０度であるＶ字谷地形の正規平面曲率である－１．５を第１閾値とすれば、その第１閾値以下の正規平面曲率となるセルは折入角が９０度以下の凹地領域として抽出される。

【0059】

セル群の中央のセルが凹地領域であるとみなすための正規平面曲率に対する第１閾値を、低解像度化する際に用いたｎ×ｎのセルのｎに応じて変化させれば、高解像度ＤＥＭを平均化した低解像度のＤＥＭであっても、平均化の影響を加味して凹地領域を抽出することができる。

【0060】

このようにして定めた第１閾値は、低解像度化する際に用いたｎ×ｎのセルのｎに応じて異なるものであるが、同じｎ×ｎのセルで低解像度化した複数のＤＥＭ（Ｌ）について共通して用いることができる。

【0061】

以上に説明した集約ステップ、平面曲率計算ステップ、正規平面曲率計算ステップ、凹地判定ステップをまとめると次のようになる。まず、基準となる高解像度のＤＥＭ（１）の解像度を多段階で低下させたＤＥＭ（Ｌ）を作成する。次に各ＤＥＭ（Ｌ）に含まれる各セルについて正規平面曲率を計算する。

【0062】

次に、理想的なＶ字谷地形と判断される、等高線の折入角が９０度と５３度であるときの正規平面曲率を上記の式（３）から求めて第１閾値とする。例えば、折入角が９０度であればｋ＝１とし（図５（ｂ）参照）、式（３）より得られた値を第１閾値とする。同様に折入角が５３度であればｋ＝２とし（図５（ａ）参照）、式（３）より得られた値を第１閾値とする。

【0063】

次に、ＤＥＭ（Ｌ）ごとに計算した正規平面曲率が上記の第１閾値以下となるセルを凹地領域として抽出する。図１０に抽出された凹地領域を例示する。

【0064】

図１０（ａ）～（ｆ）は、第１閾値として－１．５を用い、その第１閾値以下となる正規平面曲率のセル、すなわち凹地領域として判定されたセルを黒い四角で表したものである。図１０（ａ）の背景は１０ｍ間隔の等高線を表し、図１０（ｂ）～（ｆ）の背景は２ｍメッシュの傾斜量図である。

【0065】

図１０（ａ）は、セルの辺の長さＬが１０ｍ、つまりＤＥＭ（１０）において凹地領域として判定されたセルが黒い四角形で表されている。図１０（ｂ）は、１ｍ間隔のＤＥＭ（２）～ＤＥＭ（１０）を一つのグループとし、各ＤＥＭにおいて凹地領域として判定されたセルが黒い四角形で重ねて表されている。図１０（ｃ）～（ｆ）は、それぞれＤＥＭ（１１）～ＤＥＭ（２０）、ＤＥＭ（４１）～（５０）、ＤＥＭ（９１）～ＤＥＭ（１００）、ＤＥＭ（１９１）～ＤＥＭ（２００）の各グループについて、図１０（ｂ）と同様に凹地領域として判定されたセルが黒い四角形で重ねて表されている。

【0066】

図１０（ａ）では図１０（ｂ）～（ｃ）と比べて凹地領域が点在しているが、図１０（ｂ）～（ｃ）では凹地領域が連なって抽出されていることが分かる。図１０（ｄ）～（ｆ）では、セルが大きくなり尾根を超えてしまっている。したがって、セルの長さの上限としては尾根を超えないように適宜設定する（図１のステップＳ５）。同図では、例えば、セルの長さの上限としてＬ＝２０ｍとした。なお、凹地領域が尾根を越えているか否かは、凹地領域と等高線との位置関係に基づいて人が判断してもよいし、画像処理技術を用いてコンピューターにより判断してもよい。

【0067】

次に、上述のようにして得られた凹地領域から谷領域を判定する。谷地形は凹地領域が落水線の方向に連続する地形であることを踏まえると、傾斜角と長さに基づいて凹地領域から谷領域を抽出することができる。

【0068】

まず、１つのグループに属する複数のＤＥＭごとに抽出した凹地領域を重ね合わせる（図１０参照）。次に、１つのグループに属する複数のＤＥＭの中で１ｍのセルごとに傾斜角を計算する。つまり、一つの１ｍのセルには、ＤＥＭ（Ｌ）ごとに計算した複数の傾斜角が得られる。次に、重ね合わせて得られた凹地領域から、最小の傾斜角が第３閾値以上であるセルを抽出し（図１のステップＳ４）、抽出された連続するセルからなるポリゴンを作成する。

【0069】

傾斜角の計算は、図５（ａ）のセル群を例にとると、セルＣ２２が凹地領域として判定されたセルであるが、その周囲のセルＣ１１～セルＣ３３の全てのセルを元にセルＣ２２の傾斜角を計算する。このような傾斜角の計算は公知の方法であるので詳細な説明は省略する。

【0070】

図１１は、凹地領域とその傾斜角を示す図である。図１１（ａ）はＤＥＭ（２）～ＤＥＭ（１０）のそれぞれで抽出された凹地領域を重ね合わせたものである。凹地領域に含まれる１ｍのセルごとに、ＤＥＭ（２）を用いて計算した傾斜角、ＤＥＭ（３）を用いて計算した傾斜角、・・・、ＤＥＭ（１０）を用いて計算した傾斜角の中から最小の傾斜角を第３閾値に比較される傾斜角として採用する。同図の例では、セルは、最小の傾斜角が０度～２度、３度～４度、・・・など、どの範囲に属するかに応じて色分けされている。同様に、図１１（ｂ）は、ＤＥＭ（１１）～ＤＥＭ（２０）についてセルが最小の傾斜角に応じて色分けされている。

【0071】

セルの傾斜角と比較される第３閾値は、谷地形とみなす傾斜角である。この第３閾値は、図１１のように地図上に谷地形とその傾斜角を表示し、視覚的に判断して設定する。例えば、ＤＥＭと重ね合わせて表された等高線などを見て谷地形と判断できるセルがあったとき、そのセルにおける傾斜角が採用されるように第３閾値を定める。他にも、ＤＥＭと重ね合わせて表された等高線などを見て尾根や林道と判断できるセルがあったときに、そのセルにおける傾斜角が採用されないように第３閾値を定める。図１１の例では、７度を傾斜角に対する第３閾値として定めた。

【0072】

さらに、谷地形の傾斜角は、谷地形以外の尾根や林道の傾斜角と比べて急であると考えられる。ある傾斜角を第３閾値としたときには、尾根や林道などが除外される結果、谷地形として抽出される広さ（セルの数）は大きく変化すると考えられる。したがって、第３閾値を変化させたときに、抽出される谷地形の広さの変化に基づいて第３閾値を決定することができる。

【0073】

なお、地域によって谷地形と判断される傾斜角は異なるので、傾斜角に対する第３閾値は地域毎に設定することが好ましい。

【0074】

次に、図１２を用いて、長さに基づいて谷領域を抽出する方法について説明する（谷領域抽出ステップ；図１のステップＳ６）。まず、低解像度化した複数のＤＥＭについて、セルの長さに応じてグループ分けする。例えば、セルの長さが２～１０ｍ、１１～２０ｍ、２１～３０ｍ、・・・であるＤＥＭをグループ化する。そして、グループ内の何れかのＤＥＭにおいて、正規平面曲率が第１閾値以下かつ傾斜角が第３閾値以上となるセルを抽出したものをポリゴンとする。

【0075】

図１２には、セルの長さが２～１０ｍのグループ（ＤＥＭ（２）～ＤＥＭ（１０））について抽出された６個のセルからなるポリゴンＰが例示されている。セルの大きさが同じであれば大きな面積のポリゴンは細帯形であるが、小さいポリゴンは楕円形に近似される。また、ポリゴンＰに含まれるセルのうち最大のセルの長さは４ｍである。ポリゴンＰの面積は、楕円形の面積Ｓｐである（４／２）×（１７／２）×π＝１７πに近似できる。

【0076】

一方、一つのセルの長さをポリゴンＰの最大のセルの長さとした３×３のセル群においては、最短の谷の長さは３個のセルであると考えられる。３個のセルに内接する楕円形の面積Ｓｔｈは、３πＬ^２／４である。

【0077】

面積Ｓｔｈをセルの長さＬの二乗で割った値は約２．４となる。ポリゴンＰの近似の面積ＳｐをＬ^２で割った値が２．４以上であるならば、ポリゴンＰは最短の谷よりも長い谷であると判別できる。同図の例では、ＳｐをＬ^２で割った値は約３．３であり、第２閾値の２．４以上であるので、ポリゴンＰは十分な長さを有する谷領域であると判定することができる。

【0078】

凹地領域から谷領域を抽出するステップを纏めると次のようになる。まず、グループ化した複数のＤＥＭ（Ｌ）の凹地領域についてセル毎に最小の傾斜角が第３閾値以上となるものを抽出する。そして、複数のＤＥＭ（Ｌ）ごとに抽出された傾斜角が第３閾値以上の凹地領域（セル）を重ね合わせ、連続した凹地領域をポリゴンＰとして抽出する。ポリゴンＰの谷線方向の長さが第２閾値以上であれば当該ポリゴンＰを谷領域として抽出する。ポリゴンＰの長さが第２閾値以上であるかを判定する具体的な方法としては次のように行う。ポリゴンＰを近似した楕円形の面積Ｓｐを求めるとともに、その面積ＳｐをポリゴンＰに含まれる最大のセルの長さの二乗で割った値を計算する。その値が第２閾値である２．４以上であればポリゴンＰは谷領域であり、２．４未満であれば谷領域ではなく局所的な凹地として除外する。

【0079】

次に、水系を作成する（水系抽出ステップ；図１ステップＳ７）。具体的には、上述のようにして作成した谷領域を解析用に設定した大きさのセル（以下解析用セル）で分割し、各解析用セルに一定値を与え、上流から下流へ累積流量を計算する。累積流量が第４閾値を越える解析用セルを水系とする。この水系は、第４閾値を始めて越えた解析用セルを基点とする。累積流量の計算については公知の方法を用いることができるので詳細な説明は省略する。

【0080】

図１３は、Ｖ字谷地形の開始点付近に重なる谷領域から生じる累積流量を示す図である。解析用セルは一辺がｌの正方形のセルとする。Ｌは谷の幅を表している。これらの比であるＬ／ｌをｗとする。点線で示すＶａｌｌｅｙ ｌｉｎｅは谷線を表し、灰色のセルで示すＶａｌｌｅｙ ｚｏｎｅは谷領域を表している。同図には比が３、５、７のそれぞれについて累積流量が示されている。谷領域に含まれる解析用セルに流量として１を与え、累積流量を計算すると図に記載したとおりとなる。なお、谷線における累積流量を下側に抜き出してある。

【0081】

図１３に示すようなＶ字谷地形の開始点付近に重なる谷領域に流量として１を与えた場合において、Ｖ字谷地形の開始点における累積流量ｖは式（４）で表される。

【数4】

【0082】

図１３に示すように、Ｖ字谷地形の開始点における累積流量は、ｗ＝３のときは４、ｗ＝５のときは１４、ｗ＝７のときは３０となっている。したがって、式（４）で表される累積流量を第４閾値とすることで、累積流量が第４閾値以上であるセルは（同図では下側に抜き出した連続するセル）、最初に第４閾値を越えたセル（同図の例では累積流量が「４」「１４」「３０」のセル）を基点（谷頭）とする水系となる。

【0083】

例えば、ＤＥＭ（２）～ＤＥＭ（１０）のグループから得られた谷領域について、ｌ＝５、Ｌ＝１０、ｗ＝２として累積流量の計算を行う場合は、式（４）より第４閾値として１（式の計算結果の小数点以下を切り捨ててある）を用いる。同様に、ＤＥＭ（１１）～ＤＥＭ（２０）のグループから得られた谷領域についてｌ＝５、Ｌ＝２０、ｗ＝２、第４閾値として８（同様に小数点以下を切り捨ててある）を用い、ＤＥＭ（２１）～ＤＥＭ（３０）のグループから得られた谷領域についてｌ＝５、Ｌ＝３０、ｗ＝６、第４閾値として２１（同様に小数点以下を切り捨ててある）を用いる。このようにしてＤＥＭのグループごとに、谷領域の最上流の位置から開始する水系を得ることができる。

【0084】

図１４に示すように、上述のようにして得られた水系は、屈曲部や合流点の近傍において第４閾値以上の累積流量を示す可能性がある。同図のセルは谷領域を表し、そのうち斜線のセルは水系を表している。「×」が付されたセルは、屈曲部や合流点の近傍において第４閾値以上の累積流量を示すセルを示している。このように一つの谷領域内であるにも関わらず、主要な谷から分岐した短い谷（「×」のセル）が多数生じる。

【0085】

そこで、谷幅のセルの長さｌを第５閾値に設定し、谷次数で区切られた谷の長さがその第５閾値よりも長いものを取り出す（水系抽出ステップ；図１のステップＳ８）。谷次数の区分はＳｔｒａｈｌｅｒ方式をとり、同次の谷が合流するときに、合流後の谷の次数は谷の次数＋１となる。図１４の例では、Ｌ１の谷は、第５閾値ｌより長いので水系として残る。Ｌ１、Ｌ３、Ｌ４及びＬｂは一次谷、Ｌａは二次谷となり、これらは個別に長さの判定を行う。Ｌ１、Ｌａ及びＬｂは第５閾値ｌよりも長いので水系として残る。Ｌ３及びＬ４の「×」が付された谷は第５閾値以下であるので水系から除外される。ＬａとＬｂは最終的に一つの一次谷を成す水系であるＬ２として残る。第５閾値は、例えば、ＤＥＭ（２）～ＤＥＭ（１０）では５ｍ、ＤＥＭ（１１）～ＤＥＭ（２０）では１５ｍ、ＤＥＭ（２１）～ＤＥＭ（３０）では２５ｍとする。

【0086】

次に、水系を細線化する（水系抽出ステップ；図１のステップＳ８）。図１５を用いて細線化について説明する。同図のセルは谷領域を表し、そのうち斜線のセルは水系を表している。図１５（ａ）に示すように、水系は、一方向（同図の上下方向）に延び、セル１個分の幅を有する２本の流路を含んでいる。このうちの一本の流路を残し、その余の流路を除去する。その結果、図１５（ｂ）に示すように、一本の流路からなる細線化された水系を得ることができる。

【0087】

図１４及び図１５で説明した処理によって、本流の近傍に発生する谷の開始点が除去されることになる。換言すれば、谷の最も上流側に特定された谷の開始点のみを取り出すことができる。このようにして取り出された谷の開始点を元に、再度、累積流量の計算を行うことで、谷幅分布図に調和的な水系を得ることができる（水系抽出ステップ；図１のステップＳ９）。

【0088】

図１６を用いてＤＥＭから谷地形及び水系の抽出方法により得られた結果を示す。
図１６（ａ）には、黒で塗られたセルからなるセルが凹地領域として抽出されている。図１６（ｂ）には、凹地領域について傾斜角及び長さの第２閾値に基づいて抽出して得られた谷領域（ポリゴン）が示されている。赤で塗られたセルは傾斜角が０度～７度であり、青で塗られたセルは傾斜角が７度～４０度である。傾斜角が７度～４０度の谷領域は、実際の谷地形に適合していることが分かる。

【0089】

図１６（ｃ）には、傾斜角が７度～４０度の谷領域について累積流量を計算し、第４閾値以上である累積流量を示す緑色のセルからなる水系が示されている。図１６（ｄ）には、水系について短い谷を除去したものが赤色で示されている。図１６（ｅ）には、細線化された水系が黄色で示されている。図１６（ｆ）には、青色のセルで示される谷頭から累積流量を計算して得られた水色の水系が示されている。

【0090】

以上に説明した谷地形の抽出方法及びそのプログラムは、低解像度化した複数のＤＥＭ（Ｌ）について正規平面曲率を計算する。正規平面曲率は、ＤＥＭのセルの長さＬに依存しない値である。したがって、正規平面曲率に対して設定した第１閾値は、異なる解像度の複数のＤＥＭ（Ｌ）に共通して用いることができ、客観的にＶ字谷地形を凹地領域として抽出することができる。
また、谷地形の抽出方法は、凹地領域の長さと第２閾値に基づいて谷領域を抽出する。これにより、孤立した凹地領域を除外し、Ｖ字谷地形として適切な長さを有する谷領域を客観的に抽出することができる。
さらに、谷地形の抽出方法は、複数のＤＥＭ（Ｌ）ごとに得られた凹地領域を重ね合わせたポリゴンを谷領域として抽出する。例えば、横断面の傾斜が緩やかなＶ字谷地形は、短いＬのＤＥＭ（Ｌ）では凹地領域として抽出できないが、長いＬのＤＥＭ（Ｌ）では凹地領域として抽出できる場合がある。本発明の谷地形の抽出方法では、複数の解像度のＤＥＭ（Ｌ）から得られた凹地領域を重ね合わせてポリゴン（谷領域）とするので、様々な形状のＶ字谷地形が漏れにくく、網羅的に谷領域を抽出することができる。したがって、谷地形の抽出方法により抽出された谷領域は、現実の谷地形に適合したものとなる。

【0091】

また、谷地形の抽出方法は、傾斜角が第３閾値以上である凹地領域について谷領域（ポリゴン）を作成する。これにより、現実の谷地形に、より一層適合した谷領域を得ることができる。

【0092】

また、谷地形の抽出方法は、ポリゴンの面積と第２閾値に基づいて谷領域を抽出する。これにより、現実の谷地形に、より一層適合した谷領域を得ることができる。

【0093】

また、谷地形の抽出方法は、第１閾値の絶対値として本明細書に記載の式（３）を用いる。これにより、より一層客観的に凹地領域を抽出することができる。

【0094】

また、谷地形の抽出方法は、凹地領域が尾根を越えないように、集約に用いるセルの一辺の長さを定める。これにより、現実の谷地形により一層適合した凹地領域を得ることができる。

【0095】

水系の抽出方法及びそのプログラムは、谷地形の抽出方法によって得られた谷領域から水系を抽出する。水系の抽出は累積流量と第４閾値に基づいて行われる。これにより、Ｖ字谷地形の谷頭から始まる水系を得ることができる。また、谷地形の抽出方法によって得られた谷領域の範囲内で累積流量の計算を行えばよいので、ＤＥＭ全体から累積流量を計算するよりも効率的に水系を得ることができる。
また、水系の抽出は、解析用セルの長さｌと谷領域の幅Ｌの比ｗとし、図１３に示したように比ｗごとに累積流量を計算する。ある比ｗでは谷頭が検出できなくても、他の比ｗでは谷頭を検出することができる場合もある。したがって、様々な形状のＶ字谷地形であっても谷頭が漏れにくく、網羅的に谷頭を抽出し、その谷頭から開始する水系を抽出することができる。
また、谷地形に適合した谷領域を対象として水系を作成するので、従来技術のように谷地形ではないが集水面積が大きい箇所を水系として抽出することを回避でき、かつ人工構造物がなす水平面を含まない水系を抽出することができる。
また、異なる解像度のＤＥＭ（Ｌ）から作成された谷領域に基づいて水系を作成するので、地図の縮尺によって異なる水系が作成されることはない。

【0096】

また、水系の抽出方法は、谷の長さ及び第５閾値に基づいて短い谷を除外する。これにより、その短い谷を起点とする水系を除外することができる。さらに、水系を細線化し、水系の最も上流から再度累積流量を計算し、第４閾値を越える解析用セルからなる水系を抽出する。これにより、地域で主要な谷地形に基づく水系を得ることができる。

【0097】

なお、谷地形の抽出方法及び水系の抽出方法における各ステップをプログラムとし、このプログラムをコンピュータに実行させてもよい。これにより各ステップがコンピュータ上で実現される。

【0098】

また、 上記プログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等が挙げられる。また、上記プログラムは、インターネット等の通信手段を介してダウンロード可能としてもよい。

【産業上の利用可能性】

【0099】

山地の地形表現図を利用したり、斜面災害のリスクを評価する産業分野などで利用したりすることができる。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】