5.ちょっと考えること5.1 「傾斜量の値」 と 「傾斜量図の色区分」 国土地理院のホームページでは50mDEMを使って全国の傾斜量図を作成し、公開している。その際には傾斜量を次の式により0〜249のデジタルナンバーに換算し、デジタルナンバーにより作図している。すなわち傾斜量を角度に変換し、50°の角度が最大値の250となるように計算している(正確には49.8°以上を全部249としている)。 D=250×{arctan(S)×180/π/50} 作図するときに角度にしちゃってるので、結局、計算上は斜面勾配図とはXY方向の重み付けだけの違いとなるのだが、実はさらに傾斜量の小さい部分が良く分かるように色調補正も施されているとのこと。どういう方法・係数で色調補正しているのかわからないが、色調補正するということはすなわち計算された値をさらに加工しているということ。 ※) 後になって知ったが、差分式を使った演算による手法をオペレータと呼ぶらしい。そしてこれらは画像データの基本的な演算手法で、手法の開発者の名前を取って傾斜量として知った方はPrewittのオペレータ、斜面勾配として知った方はSobelのオペレータと呼ばれる誰でも知っているものだったようだ。 タンジェント曲線と4乗根の違いは下図のごとし。このグラフではいずれも傾斜量 S=1.2 くらいが250になるように係数をつけているが、係数にそれ以上の意味はない。 ![]() 実際どうなのか、比較してみた。 まず、左が傾斜量図。右が斜面勾配図(4近傍重み付け)。いずれもアークタンジェントで角度の換算しているが、見た目、ほとんど違いがわからない。 ![]()
国土地理院基盤地図情報(10mDEM)および同閲覧コンバートソフトを使用
次に、利尻のように単調な地形だと面白くなさそうなので、わりと低平な丘陵地の例と並べてみた。同一の色調で表現するため、上記3つの換算式で計算したものをひとつのファイルに入れた。左2列が低平丘陵地、右2列が利尻島西斜面で、同じ縮尺。それぞれの左右はポジとネガの関係。上からタンジェント曲線、平方根、4乗根による換算。 利尻は急傾斜地との印象だったのだが、急なのは山ひだの発達する上部だけで、中腹から下方は意外と勾配が緩いようだ。緩くても単調に続くとそびえたつように見えるからなのかもしれない。逆に低平丘陵地はなだらかで急斜面のイメージがなかったが、沢で削られた谷壁は結構な急斜面になっていたようだ。 両者とも4乗根はやりすぎで、平方根で換算した後さらに色調を調整するのがよさそうに思えるが、このあとは解析する人が何を見るか、ということだろう。 ![]()
国土地理院基盤地図情報(10mDEM)および同閲覧コンバートソフトを使用
* 井上誠、伊計秀明(2001):「傾斜量図の利用法について」、情報地質、vol.12、No.2、pp.72-73 * 神谷泉、黒木貴一、田中耕平(2000):「傾斜量図を用いた地形・地質の判読」、情報地質、Vol.11、No.1、pp.11-24 * 岩橋純子(1994):「数値地形モデルを用いた地形分類手法の開発」、京都大学防災研究所年報37号B-1、pp.141-156 ===== 【追記】 2013.10.10 ![]() この間に知り合った人の中に、緩傾斜面は角度で区分して検討したい、という人がいた。その場合は単純にアークタンジェントで角度換算して離散表示すればよいのだが、傾斜変化を視覚で感覚的にとらえるなら、やはり偏った強調が必要となる。ArcGISソフトではヒストグラムや標準偏差で平坦化するなどができるが、それをすると異なる地域を同じ指標で見ることができない。だからとりあえずは決まった計算式で換算しておく必要がある。 どういう計算式を使うかであるが、上記を参考に考えると、たぶん、勾配45°以上になると急崖の中の微地形はあまり問題にならないように思う。緩傾斜はやっぱり強調したいと思う。4乗根だとやりすぎの感がある。平方根、もしくは2.5乗根くらいかな? 5.2 断面図作成例 高度分散量異常図と地形断面を緯度経度から自動的に作成し、それに別な図のデータと地質断面を重ねたエクセルグラフ。地質断面測線とDEMデータからとった断面測線とがややずれているのと、縦横スケールが違うので断面図はぴったり合わないが、地質図から読み取った断層位置も入れている。 こんなこともエクセルだけでできちゃう、という例。 ![]()
高度分散量異常図は国土地理院基盤地図情報(10mDEM)および同閲覧コンバートソフトを使用
5.3 斜面曲率 (1) 西田ほか(1997)の式 斜面曲率についてネット検索してみた。 高精度DEMを使った研究例が出てきたが、偏微分方程式を使って算出せよ、とある(国総研資料、No.204、国総研資料、No.511など)。簡易計算プログラムのマニュアルがついているがプログラムそのものは見つけられなかった。
国総研資料 No.511に載っている平均曲率の読み取り標高の一覧(表2.10〜2.12)のデータを使って比較してみた。なお、この表には明らかな入力ミス(例えば、標高4505m、6639mなど)がいくつかあるため、私の計算結果と平均曲率の基本統計(表2.15)とは一致しない。
上段は国総研資料 No.511の表2.15。 下段はそれを取り出して計算した結果。 1桁違ってたり、『30.0 311.5 31.5』『131.0 131.5 313.0』のように並んでいる入力ミスは勝手に数値を推定した。こういうミスがあるということは、気づきにくいミスもほかにあると思うが、計算結果の基本統計の数値はほぼ似ているので計算は合っていると思う。 計算式は、西田ほか(1997)ではの縦横方向と斜め方向のそれぞれから算出される平均曲率を平均しているが、国総研資料 No.511の計算結果は縦横方向のみを使っているらしい。 縦横方向から算出される平均曲率(以下、縦横平均曲率)を横軸にとり、斜め方向から算出される平均曲率(以下、斜め平均曲率)と下記に記したArcGISの曲率ツール出力値の1/100(以下、ArcGIS曲率)を下左図のようにプロットした。縦横方向と斜め方向の平均曲率は、あまり相関性がよくないようだ。これは、曲面に方向性があるからかな? 例えば痩せ尾根の鞍部は、稜線方向には凹状となるがその直交方向は凸状となり、平均曲率はゼロに近くなる。稜線と45°方向を考えると、上に凸状の断面同士になって平均曲率はマイナスの値が出てくる。 ArcGISの出力値は縦横方向の4近傍ラプラシアンなので、縦横の4近傍を主に計算対象とする縦横平均曲率とよく相関するんだろうと思う。 同じく、下右図に縦横平均曲率と凸凹度をプロットした。ここで使われているデータは10mメッシュなので凸凹度はArcGIS曲率の2.5倍となるだけであるため、相関性は同じになる。 ![]() ![]() ![]() 斜面勾配を指標に縦横平均曲率との関係を見ると、勾配が急になるほど曲率のばらつきが減って小さくなる傾向が見える。これは、急斜面で曲率が大きいということは、その勾配よりさらに急な斜面と緩い斜面との組み合わせになるということであり、それはメッシュデータでは拾いにくいからと思う。 凸凹度と縦横平均曲率の比率をとり、斜面勾配を横軸にプロットしてみた。凸凹度は、斜面勾配が急になるにつれて縦横平均曲率との比が大きくなる。この傾向は、西田ほか(1997)が示しているラプラシアンと平均曲率の斜度による変化と同じことを示しているものと思う。 ![]() ここでもう一度、斜面のデコボコ程度の意味を考えてみた。 西田ほか(1997)は、地震と斜面崩壊との関連性で、出っ張っているところは地震動による揺れを増幅するのではないか、という観点で見ている。おそらく出っ張っている方向と斜面方向、すなわち重力方向とは無関係なので、デコボコの程度が斜面勾配で変わってくると尺度が変わってしまうのでラプラシアンは使えない、ということなんだと思う。 しかし、例えば表流水の流れを考えるとどうだろう? 同じ曲率であっても、急勾配な方が流速が速く、侵食力も大きくなる。そうすると斜度によって変化する方が都合良いような気がする。それは勾配で補正するからやはり変化しない方がいいのかな? (2) ArcGISの式 ArcGISの曲率計算は以下のようになっている。 [曲率(Curvature)] ツールは、入力サーフェスの二次導関数値をセルごとに計算します。 セルごとに、次の四次多項式 Z = A*x^2*y^2 + B*x^2*y + C*x*y^2 + D*x^2 + E*y^2 + F*x*y + G*x + H*y + I ![]() 図に示した、番号の付いた各セルの 9 個の標高値と係数との関係は次のとおりです。 A = [(Z1 + Z3 + Z7 + Z9) / 4 - (Z2 + Z4 + Z6 + Z8) / 2 + Z5] / L^4 B = [(Z1 + Z3 - Z7 - Z9) /4 - (Z2 - Z8) /2] / L^3 C = [(-Z1 + Z3 - Z7 + Z9) /4 + (Z4 - Z6)] /2] / L^3 D = [(Z4 + Z6) /2 - Z5] / L^2 E = [(Z2 + Z8) /2 - Z5] / L^2 F = (-Z1 + Z3 + Z7 - Z9) / 4L^2 G = (-Z4 + Z6) / 2 / L H = (Z2 - Z8) / 2 / L I = Z5 [曲率(Curvature)] ツールの出力は、サーフェスの二次導関数(たとえば、傾斜角の傾斜角)で、次のようになります。 曲率 = - 2 * (D + E) * 100 たくさん式が載っているが、ようするに曲面の一般式として四次多項式を示しているだけで、そのうち曲率に関係するのは D と E だけと書いてある。ようするに4近傍ラプラシアンである。 10mDEMと50mDEMについて、計算してみた。曲率半径の考え方もよくわからなかったので、とりあえず隣接する3点を通る円を考えた。円の半径に対し、中心から50m、10mの標高差なら簡単に計算できる。その標高差から D = E として求めた。 ![]() ![]() その結果、50mDEMでは入力半径50mで計算される D + E の逆数は25となるが、入力半径200m以上であれば算出される D + E の逆数とほぼ一致する (10mDEMだと入力半径30m以上)。なるほど、そうするとマイナスをつけるのは尾根をプラス、谷をマイナスにするためとわかる。100倍している意味が分かればそのまま曲率として使えるということになると思い、ずいぶん悩んでしまった。 あれこれ見てみたが、『曲率とは、曲率半径の逆数である』 『曲率半径とは、近似される微小円の半径である』 としか出てこない。隣接する点を通る円が近似される微小円ではないので、 D + E の逆数が曲率半径ではないと思うのだが、離散近似としては曲率半径となるのだろうか? それにしても D + E を100倍すると曲率になる理由はなんだろう? とずいぶん悩んでしまったが、実は簡単なことだった。 ArcGISのヘルプは英文を翻訳ソフトで日本語表記にしているだけなので日本語がわかりにくいと思っているのだが、 『曲率出力ラスタの単位は、オプションの出力断面曲率ラスタと出力平面曲率ラスタの単位も、Z 単位の 1/100 です。3 つの出力ラスタすべてに合理的に期待される値は、丘陵地(中程度の起伏)で -0.5 〜 0.5、急峻で険しい山岳地帯(極端な起伏)で -4 〜 4 の範囲になります。』 との記述がある。「 Z 単位の 1/100 」 の意味が分からなかった。しかし、単位が書いてないから気づかなかったが、なんのことはない。 m^-1 単位なのを cm^-1 単位にしているだけのことのようだ。 と思ったのも早とちり。 逆だ。1000倍すれば km^-1 になる。曲率0.01m^-1は曲率半径100m = 0.1km。この逆数は10km^-1。とすれば、100倍するのは単に数値を見やすくするだけのためとしか考えられない。曲率を扱う文献をいくつか見てみたが、100倍した値を使っているものはない。ArcGIS独自の単位なんだろうか? 次項のデコボコ図において、地形変換点として指標になるだろうと考えてた凸凹度は0.5くらいである。これはこの曲率ツールの出力では2.0に相当する。2.0は中程度の起伏より少し起伏がある斜面に相当するので、地形イメージとは合致する。 5.4 デコボコ図とラプラシアン 曲率について考え始めた当初、上記ArcGISのヘルプを見ても理解できなかったので、以下のように考えた。 曲率は曲率半径の逆数だということなので、両隣との3点を通る円の半径を求めようかと思った。しかし、円だと実際の斜面を近似しているとは到底思えない。そこで次に考えたのは、この D 、E という係数はようするに両隣セルとの標高差、つまりデコボコの程度だということ。曲率だって、結局は凹凸の程度を示す指標である。ならば数学的な意味はわからないが、両隣セルとの標高差をそのまま図示したらどうなるだろうか、と思って作成したのが、デコボコ図である。基本的な考え方は下記に示すようにラプラシアンと同じだ。 デコボコ図の地形学的な意味は、凸状を示すプラスの値は 尾根 および 遷急ゾーン、凹状を示すマイナス値は 谷 および 遷緩ゾーンとなる。遷急ゾーン、遷緩ゾーンという単語は聞いたことないので私の造語だと思うが、遷急線・遷緩線というほどくっきりしたラインにならないのと、斜面傾斜方向だけでなく水平方向の変化も含むのでゾーンとした。まだ確認していないが、閾値を設けることでゾーンを線に近づけられるような気がしている。さらに傾斜方向をファクターに入れれば、尾根と遷急線も区別できるかもしれない。 さて、実際に計算し作図してみて驚いた。尾根谷度図と傾斜量図を重ねた偽赤色立体地図に重ねてみたところ、緩斜面の微地形がはっきりし、段丘のような地形がくっきり浮かび上がった。ラプラシアンをこのように重ねたことがないが、ラプラシアンでも同じ効果が期待できることを知った。 段丘崖付近の地形を傾斜量図でみると、緩−急−緩 の繰り返しとなるが傾斜変換線部分は同じ色で表現される。デコボコ図だと斜面下から 平−遷緩−平−遷急−平 となり、傾斜変換線が区別される。白地図に傾斜変換線を入れるだけでも地形がイメージできることを考えれば、デコボコ図は地形判読に使えそうだ。 【その後、岩橋(1992)や岩橋(1995)で、地形分類には傾斜量とラプラシアンが重要だと書いていることを知った。】 -------- ラプラシアンとはなんなのか、これをするまで知らなかった。 ラプラシアンとは、目標の微小点における勾配の程度である。数値計算をするために離散近似している。画像処理で使われているラプラシアンフィルタでは、セルは正方形だから縦横比が 1:1 だという前提があるし、セル値そのものにもメートルのような単位がないのでセル間隔で除するなどということをしない。しかしDEMデータで扱うラプラシアンの計算ではセル値に単位があるしグリッドの大きさとで勾配の程度に大きく関連する。正方グリッドなら構わないだろうが、国土地理院の10mDEMは正方グリッドではないので、東西と南北のセル間隔 (グリッドの大きさ) で補正が必要となる。したがって4近傍のラプラシアンは、このようになる (記号はArcGISに準ずる)。 lap4 = [(Z4 + Z6) - 2 * Z5] / Lx^2 + [(Z2 + Z8) - 2 * Z5] / Ly^2 ここでラプラシアンとは別に、斜面勾配に着目して考えてみた。 凸凹度_x = - [(Z4 - Z5) / Lx + (Z6 - Z5) / Lx] / 2 X方向だけを見ると、これはZ5を基準とした、Z4との勾配と、Z6との勾配の平均値の2倍である。言い換えると、グリッド2つ分の平均標高差になっている。さらにX方向とY方向とを加えているので勾配の平均値の4倍となる。これでは実際の地形がイメージしにくいので、平均値として定義しちゃおう。さらに尾根谷度とイメージを合わせ、目標点のZ5が周囲より高い場合がプラスとなるように正負を逆にする。 勾配平均値と書くと斜面勾配と紛らわしいので、ここではこの勾配平均値のことを凸凹(デコボコ)度と称することにする。 凸凹度という観点から作図したものを見ると、どうやら絶対値が0.05より大きい個所が勾配変化点として認識するのによさそうだ。0.05という数値は、10m隣で50cm高さが変化している程度を示す。例えば、両隣が50cm低い、片側は水平だがもう一方が1m高い、という状態だ。 0.05を指標にして着色すると、尾根地形、谷地形はやや幅を持って表現され、遷急線、遷緩線も識別できる。特に等高線幅がまばらになる緩傾斜面は明瞭に分かる。 緩斜面がわかりやすいというのは、起伏量とも関係するのだろうか? 低地は起伏量が小さく、山地では大きくなることは知られている。凸凹度にも同様の傾向が推測される。ならば、山地の凸凹度を起伏量で低減させることで山地も低地と同様な濃淡表現にならないのかな。 いずれ、考えてみたい。 これまで、ラプラシアンは画像処理のためのフィルターとしての使い方しか見ていなかった。そもそもラプラス方程式なんてものは、その名称だけは聞いたことあったが何を示すものかなど知らなかった。 デコボコ図はラプラシアンの亜流になるようだが、実は隣接セルとの変化量の差を表現するのが目的なため、二乗でも単純平均でも、どっちでもよいことらしい。傾斜量の算出に4近傍で重みをつけるかつけないかというのと同じことらしい。しかしこうやっていろいろ考えてみて、実は15年前にすでに論文化されていることを後で知ったが、勉強になったと思っている。 5.5 微分 って何さ? 地形解析図のことを調べるにつれ、空間微分とか、偏微分とか、2階微分とか、まぁ、今までも耳にしてはいたけどなんであるかは調べたくもなかった言葉が何度も出てきた。しかし、仕方なく主にネット検索を使って覚えようとしてみた。結局、今でもよくはわからないのだが、何となくイメージできたこと。 1階微分 ⇒ 地形面の変化度合 = 傾き 2階微分 ⇒ 傾きの変化度合 ≒ 曲率 ≒ ラプラシアン 3階微分 ⇒ 曲率の変化度合 微分 なんて簡潔だけど捉えどころのない単語を使うからなんだかわからなくなる、というか、イメージすることを避けたくなる(私だけ?)のだが、”変化の度合い” と置き換えると、なんとイメージしやすくなることか。これまで微分がイメージできてなかったんだから、ラプラシアンがイメージできるわけがない。 これは、当たり前にみんな知ってることなんだろうか? ・・・なんだろうなぁ、、、 とある文献に重力異常の水平微分図というのがあった。じゃぁ、と高度分散量図の1階微分図 すなわち 高度分散量の変化度合図 を作ってみた。すると、値が変化するエッジ部分が強調されて浮かび上がってきた。知ってる人には当たり前なんだろうけど、私にしては驚きの新発見だった。 いろいろ検索していくと、大学の数学や情報数理の先生のページにたどり着く。そこには、”微分” のように聞いたことはあるけど意味がよくつかめない単語が羅列してある。専門分野の専門用語と言われればその通りだが、それが数学が嫌いになる理由なんじゃないかと思えてきた。 ちなみに私は、受験数学は得意で得点科目だったが、大学に入って最初の数学の授業で 「1+1=2を証明せよ」 という課題で挫折し、数学がいやになった。そしてその時に初めて気づいた。受験数学は答えが用意されたクイズであり、学問ではなかったのだということに。 5.6 高度分散量異常図とは (1)算出法 ![]() 同様に、高度分散量異常を相対値で表現しようとする場合は、平均高度分散量の標準偏差の近似式を求め、近似式を代入することで求められる。 というのが、三箇、藤原(2000)の算出法なんだろうと理解した。しかし、、、 ⇒ 計算法は(3)に続ける。 (2)高度分散量異常の意味 もう一度立ち返って、高度分散量異常の意味を考えてみた。 高度分散量、すなわち山地の起伏の程度を標高別に並べてみると、だいたい相関しそうだ、というのが出発点になっている。これは、標高50m地点では標高差50m以上の谷壁はないが、標高500m地点では100m程度の崖があることからも想像できる。このため、低標高と高標高の起伏の程度を単純に比べてもよくわからない。それで、同標高の平均的な起伏程度との差をとることで、その地点の起伏程度の評価をしようというのが高度分散量異常であると、私は理解した。 ということは、標高ごとの高度分散量のプロットから近似式を求めるのが重要ということになる。 また、どういうスケールでの異常を知りたいかという点で、適切な計算エリアの設定が必要である。 (3)標高−高度分散量曲線の近似式 近似式そのものは、エクセルのグラフ機能でクリックするだけで求められるので特に考えることはないのだが、問題はどの近似式を使うかということかと思う。 三箇、藤原(2000)では、日本全域の高度分散量と標高との関係を1km四方の計算エリアで求めると、高度分散量もその標準偏差も、平均標高とは累乗式によって近似され高い相関性を持つことを示している。そしてその近似曲線を日本全域の標準地形曲線とし、山地別の標高−高度分散量のプロットを示し、山地別に求められる近似曲線との傾きの違いを山地の隆起と関連付けて考察している。 ![]() 三箇、藤原(2000)では、日本全域の高度分散量と標高との関係は累乗式で近似されるとしているが、地形学的には累乗式で近似されることの意味はおそらくないものと思う。つまり、高度分散量異常を求めるために使う近似式は、相関係数の大きなものを使ったほうがよく、場合によっては低標高は対数式、高標高は多項式と使い分けてもいいものと思う (例えば、右下図の赤破線)。 そういう意味では、私が当初勘違いした平均標高ごとに差し引きするのもあながち間違いでもなかったのかな、と思ったり。 ※正確には天塩山地のデータではなく、国土地理院の『数値地図50mメッシュ(標高)I』の北緯44°30'より北のデータを使用して求めた。 5.7 微分図の使い道 微分 って何さ? のその2 微分は変化しているところを示すものだと知った。試しにいろんな微分図を作って、そして微分していない図に重ねてみた。 色相の配分にもよるのだろうが、変化しているところが強調されて、視覚的に立体感のある画像となることを知った。 画像処理の本では、エッジの検出に1次微分を使うとどれにも書いてある。ずっと疑問に思っていたのだが、先日、ようやく気付いた。地形データは断崖絶壁を除くと基本的になめらかな連続値。一方、通常の画像データは背景の中に物体が置かれていて、そこは不連続値となっている。だから1次微分でもエッジが検出できるんだ、と。(2014.8.7追記) 5.8 赤色立体地図画像に関連して 2013.9.6記 微分 って何さ? のその3 赤色立体地図画像とは、『地形の凹凸に応じて地上開度と地下開度を用いて明るさを決定し、傾斜に応じて彩度を変化させる』(千葉ほか、2007)ことで地形を表現したものだそうだ。地形の凹凸を表現するんだったら曲率でいいじゃん! と思ったら、『傾斜とパラメータの相関をとると、無相関とはならない』(千葉氏のブログ)から使えないと書いてある (それで尾根谷度図を傾斜の値と独立した指標として考案したらしい)。しかし、例えばこの図のように、曲率と傾斜を合成した図でもよく見えるし、むしろ微地形はこっちのほうがわかるんじゃないのかな? (300×400pixel以下の画像にするために解像度落としているので、少々わかりにくいかも) ![]() で、そんな風にして作成したのが右図。この地区は、パッと見て地形が大きく二分されるのはこれ以外の図でもわかる。この図をよく見ると、沢が細かく入っている方もいない方も、尾根がとがり気味なところと平頂なところに分かれる。これはこの図を見るまで気づかなかったが、この地区の4つの地層の分布に一致する。 いずれにしても <微分って何さ? その2> に書いたように、ポイントは「一次微分と二次微分の合成」 あるいは 「微分していない図と一次微分の合成」 ということらしい。 画像処理における鮮鋭化にラプラシアンを使うのと同じ理屈なのかな。 ![]() なんだかんだと言っても、やはり赤色立体地図は優れているので、ぜひこれに似た画像を自分で作りたいと思っていたのだが、同じことを考える人はたくさんいるようだ。ネット検索していたら、大きなヒントが出てきた(参考8 画像の合成)ので、早速やってみた。 |