二次元ハールウェーブレット解析による重力異常値の線状構造の定量評価

─ ─29 （ ） 数々の分析手法が適用されている。例えば萩原・小竹 （1992）や山本・吉田（2009）はWalsh変換を用いて， 断層に関連するリニアメントを検出した。またプリマほ か（2012）は重力異常分布に空間的なバンドパスフィル ターを施した後，流域分割を適用してカルデラリムを抽 出している。 重力異常分布から断層やカルデラリムなど地殻構造に 関連する線状構造を検出できることは，従来の研究に よって示されてきたが，検出の評価は主に二次元的に表 示された分布図のパターン認識を基にしている。そこで 本研究では，断層のように直線的な線状構造の検出を定 量的に評価する手法を開発するため，二次元ハール ウェーブレット解析を重力異常データに適用した。二次 元ハールウェーブレット解析では，空間分布データを平 均値，X軸成分，Y軸成分，及びエッジ成分に展開す る。本研究で注目する線状構造はエッジ成分として検出 1．はじめに 重力異常分布は断層やカルデラ構造などの地殻構造と 密接に関連していることが知られている（例えば坪井， 1979）。近年，日本では稠密な重力測定が行われ，さら に測定値のデータベース化が地質調査所（現産業技術総 合研究所）（駒澤ほか，2000）や名古屋大学を中心とし た研究グループ（Gravity Research Group in Southwest-ern Japan， 2001）によって進められたことから，重力異 常データベースを用いて特徴的な地殻構造を検出する研 究が行われている（山本・吉田，2009；工藤ほか，2010 など）。重力異常分布は地表に現れない伏在構造を検出 することができるという利点がある。 重力異常分布から断層やカルデラのような構造を抽出 するためには，構造を形成する物質の密度差によって生 じる線状構造を認識することが重要であり，そのために

大渕 一樹

＊

_{・鵜川 元雄}

＊＊

Two dimensional (2D) Haar wavelet transformation is applied to gravity data for the purpose of evaluating predomi-nant direction of lineaments appearing on the gravity anomaly map. In order to detect the predomipredomi-nant direction, we rotate coordinate system with one degree interval, using the Lanczos resampling technique, and then calculate the sum of edge components of the 2D Haar wavelet transformation. The predominant direction is obtained as the minimum of the sum of the edge components as a function of rotation angle. Application of this procedure to synthetic data sets with linear discontinuities gives us the directions of the discontinuities, proving the method and procedure. Application to the actual gravity anomaly data around the Median Tectonic Line shows satisfactory results indicating the predominant direction of the lineaments of the gravity anomaly distribution.

Keywords: gravity anomaly, Haar wavelet transformation, lineament

二次元ハールウェーブレット解析による

重力異常値の線状構造の定量評価

Quantitative Evaluation of Gravity Anomaly Lineaments Using 2D Haar Wavelet Transformation

Kazuki OFUCHI

＊

_{and Motoo UKAWA}

＊＊

（Accepted November 30, 2017） 日本大学文理学部自然科学研究所研究紀要

No.53 （2018） pp.29−36

13 ＊ _{Graduate School of Integrated Basic Sciences, Nihon University: 3-25-40,}

Sakurajosui, Setagaya-ku, Tokyo, 156-8550, Japan

＊＊ _{Department of Earth and Environmental Sciences, College of Humanities} and Sciences, Nihon University: 3-25-40, Sakurajosui, Setagaya-ku, Tokyo, 156-8550, Japan

＊ _{日本大学大学院総合基礎科学研究科：}

〒156-8550 東京都世田谷区桜上水3-25-40 ＊＊ _{日本大学文理学部地球科学科：}

─ ─30 （ ）14 異常データに対して，座標軸を回転し，二次元ハール ウェーブレット変換を適用して，線状構造を検出した結 果について報告する。なお，座標軸の回転に関して，サ ンプルデータの間隔を一定にしなければならないため， 回転後にリサンプリングする必要がある。今回はLanc-zos4法 を 用 い た リ サ ン プ リ ン グ を 行 っ た（Pascal， 2011: 高木・下田，2004）。 2．解析方法 2-1 ハールウェーブレット変換 ここでは，イヴェス(2004)に従って，ハールウェー ブレット変換を行う。イヴェス(2004)によれば，ハー ルウェーブレット変換はステップ関数

φ

とウェーブ レット_{𝜓により下記のように表すことができる。} ，

𝜓𝜓

， ・・・① ・・・② これらの式でuとwはステップの始点と終点を表す。 ハールウェーブレット変換は，与えられたデータ列を ステップ関数

φ

で近似した関数

�

_{��，��}

��� � � � if � � � � �

_{� ���������}

_①

�Ψ

��，��

��� � �

� if � � � � �

�� if � � � � �

� ���������

� � �� � ����

_②

Φ

��

Ψ

Φ

_��ℓ���

� �

_����

��

_ℓ���

③

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

④

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑤

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑥

� � �

K�t� � ��t������t�

��t� � ������� �

_{� ���������}

� � if |�| � �

で表し，この関数

�

_{��，��}

��� � � � if � � � � �

_{� ���������}

_①

�Ψ

��，��

��� � �

� if � � � � �

�� if � � � � �

� ���������

� � �� � ����

_②

Φ

��

Ψ

Φ

_��ℓ���

� �

_����

��

_ℓ���

③

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

④

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑤

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑥

� � �

K�t� � ��t������t�

��t� � ������� �

_{� ���������}

� � if |�| � �

を

φ

と_{𝜓の組み合わせにより置き換えることでウェー} 数を導入し，一次元ウェーブレットを二次元ウェーブ レットに拡張する。このとき下記の第一次元（X座標） の 基 本 ウ ェ ー ブ レ ッ ト と 第 二 次 元（Y座標）の基本 ウェーブレットのテンソル積を用いる。jは周波数，k， ℓは平方ステップの始点の位置を示す（イヴェス， 2004）。

Ψ

Φ

_��ℓ���

� �

_����

��

_ℓ���

③

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

④

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑤

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑥

� � �

K�t� � ��t������t�

��t� � ������� �

_{� ���������}

� � if |�| � �

・・・③ ・・・④ ・・・⑤ ・・・⑥ 二次元ハールウェーブレットにおけるデータ列の近似 をFig.1(a)に， 変 換 の 計 算 を Fig.1(b)に 示 す。Fig.1(a) において座標と値によって構成されるデータは各正方形 の左上に位置している。各正方形内の値は全て左上の値 としてデータ列を近似し，平方ステップ関数とする。平 方ステップ関数を ③ ∼ ⑥ 式に従いウェーブレット変換 するが，具体的にはFig.1(b)に示す式による計算であ る。この操作により二次元ハールウェーブレット変換で は正方形の頂点に位置する４つのデータを平均，X軸成 分，Y軸成分，エッジ成分に展開することができる。 二次元ハールウェーブレットを用いた各成分の計算例

─ ─31 （ ） 二次元ハールウェーブレット解析による重力異常値の線状構造の定量評価 15 るには，重力異常分布を回転して，エッジ成分が最小と なる方向を見出すことが考えられる。あるいは座標軸を 回転しても同じ結果が得られる。今回，我々は座標軸を 回転させる操作を行った。回転角は0°∼90°まで1°ず つとし，方向は反時計回りである。座標軸のX軸に活断 層のトレースが平行な場合と，Y軸に平行な場合でエッ ジ成分の合計は同じとなるので，回転角は0°∼90°の範 をFig.2に示す。仮想データとしてFig.2(a)に示すもの を用いる。Fig.2(a)は60×60個のデータで構成され，色 が同じ箇所には同じ値が位置している。Fig.2(b)は計算 した各成分を成分ごとにまとめたものである。今後デー タを画像化するとき，カラーバーは全て同じものを使用 し，色の表示は場合によって二つの方法を使用する。 「値が2 なら赤」「値が 1 なら黒」というように絶対的に 色を表示する方法と，「最大値なら赤」「最小値なら黒」 というようにデータ内の値を用いた相対的に表示する方 法である。この仮想データを画像化する目的はウェーブ レット変換の結果を把握することであるため，色は相対 的に表示している。X軸成分は画像のX軸に向かう変化 のことであり，X軸方向に値（色）が変化しているとこ ろを抽出している。Y軸成分はY軸，エッジ成分は斜め 方向の変化を指す。 2-2 卓越方向の算出方法 二次元ウェーブレット変換の数値計算例をFig.3に示 す。計算例として4×4マスの数値で構成される二つの 場合を仮定した。変換する値はFig.3(a-1)(a-2)に示し， 変換後の値をFig.3(b-1)(b-2)に示す。 Fig.3(a-1)はグリッドデータ内で設定されている座標 軸と直線的な変化が平行な場合，Fig.3(b-1)は平行でな い場合の原データで，それぞれの変換後をFig.3(a-2)， Fig.3(b-2)に示す。Fig.3のように直線的な変化が一つで かつ軸に平行な場合，エッジ成分の合計が0 になる。原 データの0 のラインを活断層のトレースとすると，活断 層のトレースが座標軸に平行な場合，エッジ成分の合計 が0 に最も近くなる。 重力異常分布から直線状の活断層のトレースを検出す

Fig. 2　Example of the 2D Haar wavelet transformation. (a) Original data, and (b) the transformed data. Top left part in (b) is average data, top right is horizontal component, bottom left is vertical component and bottom right is edge component.

Fig. 3　Relationship between original data and the transformed result. (a-1) and (a-2) show the data with variation parallel to vertical axis, and the result, respectively. (b-1) and (b-2) show the data with variation obliquely crossing to vertical axis with 45 degrees, and the result, respectively.

45°傾いた段差のある分布である。Fig.5ではデータの変 化の性質をわかりやすくするため図の下部に示した相対 的なカラーバーで表示した。 解析によって得られるエッジ成分をFig.6(a)(b)(c)に 示す。エッジ成分を画像化する目的はデータ間の比較で あるため，色はデータ毎の結果の違いを比べやすくする ため絶対的に表示している。Fig.5(a)(b)はY軸に平行， Fig.5(c)では加えて45°傾いた方向に段差もしくは不連 続があるため，0°と45°回転した結果を示したいが， 0°の場合仮想データが離散的なものであるためエッジ成 分が完全に0となるので，Fig.5(a)(b)では10°と45°， Fig.5(c)では1°と45°回転した図を用いる。Fig.6を見る とエッジ成分は段差もしくは不連続がXもしくはY軸と 平行なとき弱くなっている。これは段差もしくは不連続 とXもしくはY軸が平行なとき，エッジ成分が0に近づ くことを示す。さらに仮想データを回転させながらハー ルウェーブレット解析を行うとき，各角度でエッジ成分 を合計するが，エッジ成分の絶対値を足し合わせる。ま た今回Lanczos4法において64点確保できない点は除外 しており，その数が異なる。エッジ成分の合計値を正確 に比較するため，合計値をエッジ成分の数で除算した結 果をそれぞれFig.7(a)(b)(c)に示す。 Fig.5(a)では回転角が0°，90°で直線状の変化と座標 軸が平行となるため，エッジ成分の合計が最も0に近づ く。また直線状の変化がX軸と平行でも Y軸と平行でも 同じ結果となるため，回転角0°と90°でエッジ成分の合 計値が同じになると考えられる。これらからFig.7(a)の 結果は予想通りである。0°と90°でエッジ成分の合計値 は一致していないが，これはリサンプリングによる影響 と考えられる。また45°付近で値が大きく変動している が，これ仮想データの不連続な部分が急変をしているた めと考えられる。 Fig.5(b)はPP 以外の場所ではX軸に一定の割合で傾 斜しているが，結果ではFig.6 (b)のように不連続部分を 囲とする。例えば0°でX軸に平行な場合，90°回転する とY軸に平行になるが，0°と90°でエッジ成分の合計は 変わらない。 2-3 リサンプリング 我々は座標を回転させるが，その結果データを収納す る等間隔のグリッドも回転し，元のグリッド値から新し い グ リ ッ ド 値 を リ サ ン プ リ ン グ す る 必 要 が 生 じ る （Fig.4）。ハールウェーブレット変換では，周囲のデー タの差をとる演算が含まれるため，リサンプリングの手 法が重要となる。我々は，リサンプリング手法として， ニアレストネイバー，バイキュービック，Lanczos4法 を試した。その結果Lanczos4法を除き，リサンプリン グ後にデータの差をとる演算結果に人為的な模様が生じ ることがわかった。Lanczos4法ではX -Y 座標（座標 回転後の座標）に近いX-Y座標（座標回転前の座標）を 64点選び，それぞれの距離（X -Y 座標からX-Y座標 への）を用いて重み付け平均した値を得る。重み付け係 数には以下のKを用いる。tは距離を示し，nによって 計算に用いられるデータ数が変化する。今回はn = 4 と し計算する。

� � �� � ����

Φ

��

Ψ

Φ

_��ℓ���

� �

_����

��

_ℓ���

③

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

④

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑤

Ψ

_��ℓ�，���

� �

_����

��

_ℓ���

⑥

� � �

K�t� � ��t������t�

��t� � ������� �

_{� ���������}

� � if |�| � �

Lanczos4法を用いることで不連続の発生を完全では ないが解決できる。Lanczos4法の場合64点の座標が必 要となるため，64点確保できないデータの縁を補間す ることができない。データに対して縁の補間をすること Fig. 4　Relationship between the original grids (black) and the

─ ─33 （ ） 二次元ハールウェーブレット解析による重力異常値の線状構造の定量評価 17 Fig.7の全ての結果で，回転角 45°で値が抜きん出てい て，かつ周辺に振動が見られる。値が抜きん出るのは， 仮想データの不連続が人為的なものであるためと考えら れる。また振動については，リサンプリングによる影響 の可能性が考えられる。 3-2 重力データ 重力異常値は，産業技術総合研究所地質調査総合セン ターが公表しているブーゲー異常値のグリッドデータを 抽出できた。Fig.6 (a)と(b)を比較すると (b)ではうっす らと斑模様があるが，仮想データの傾斜は階段状のもの であるため，補間方法の影響によると考えられる。 Fig.7 (c)では45°付近で最も0に近づいている。追加 した不連続の影響が強いのは，追加分のほうが約 だ け長いため影響が強くなっていると考えられる。また 0°と90°でも値が小さくなっているが，Y軸に平行な不 連続が回転角0°と90°で軸に平行となっているため，影 響が出ているためと考えられる。

Fig. 5　Synthetic data sets. Color bars indicate relative values.

Fig. 6　The results of the 2D Haar wavelet transformation of the synthetic data sets shown in Fig.5. The top image in (a) and (b) are the transformed results with 10 degrees rotation, the top image in (c) is that with 1 degree rotation, and the bottom image in (a), (b) and (c) are those with 45 degrees rotation. The color bars indicate absolute values.

示す。Fig.9では73°付近でエッジ成分の合計値が最小と なっているため，73°回転した座標軸と直線構造が最も 平 行 に 近 く な っ て い る と 考 え ら れ る。Fig.10(a)では 73°における4成分と軸に平行な補助線を示し，色は分 布を見やすくするため相対的に表示する。Fig.10(b)は 解析に使った元の重力異常データで，73°回転したもの である。平均成分は実際のデータを圧縮したものであり 実際のデータが粗くなったものとなるため，実際のデー タと傾向は同じものとなる。Fig.10(b)を見ると，全体 用いる（駒澤，2013）。グリッド間隔は南北，東西とも に500mで あ る。 こ こ で は 北 緯33 °∼ 北 緯34 °， 東 経 132°∼東経133°の部分を長方形に切り出した四国西部 から九州東部にかけての中央構造線付近のデータを利用 する。切り出したデータをFig.8に示し，色はブーゲー 異常分布を見やすくするため相対的に表示する。 今回の解析はデータの全体的な卓越方向を推定するも のなので，結果と画像全体の整合が取れているかについ て確認した。Fig.9に使用データのエッジ成分の合計を

Fig. 7　The sum of the edge components of the 2D Haar wavelet transformation as a function of rotation angle. Figures (a), (b) and (c) correspond to the synthetic data sets shown in Fig.5.

Fig. 8　The gravity anomaly data around the Median Tectonic Line. The color bar indicates relative values.

Fig. 9　The sum of the edge components of the 2D Haar wavelet transformation of the gravity anomaly data shown in Fig.8 as a function of rotation angle.

─ ─35 （ ） 二次元ハールウェーブレット解析による重力異常値の線状構造の定量評価 19 る可能性がある。 ウェーブレット解析は，平均成分に対して同様の解析 を繰り返すことでより低周波な成分の解析を行うことが できる。今回は周波数を1としているため，ウェーブ レットで解析できる最も高周波な成分のみ使用してい る。使用するのに適切な周波の目的に応じた検討が必要 と考える。 謝辞 重力異常データは産業技術総合研究所地質調査総合セン ターにより公表されているデータを使用させていただいた。 また日本大学の竹村教授には論文について有益な助言をいた だいた。ここに記して感謝いたします。 としてY 軸に平行な成分が卓越していることがわか る。特に青から緑へ大きく変化している部分は，中央構 造線に対応しているが，Y 軸にほぼ平行している。 4．まとめ 重力異常に対するX-Y軸の回転とウェーブレット解析 により，全体的な卓越方向を得ることを試みた。データ 全体を対象にした解析で，仮想データでは予想通りの結 果となった。実際のデータにおいても，全体の傾向と対 応する結果となった。今回はエッジ成分を利用して卓越 方向を得たが，X-Yどちらの軸に平行になっているのか を見るにはX軸成分，Y軸成分も利用する必要がある。 X-Y軸の回転により解析結果が変化することを逆に利 用して卓越方向を導いた。今後重力異常を解析する際 も，X-Y軸の設定が結果に影響する可能性がある。結果 への影響が予想される場合，回転しながら解析データ全 体の傾向を導くことにより，適切な角度で解析するため に活用できると考える。 ウェーブレット解析の結果を，全体の卓越方向の検出 に適用した。これはエッジ成分を主に利用したものだ が，ほかの三成分（特にX-Y軸成分）を十分に利用して いない。四成分を効率よく利用して，より情報を得られ

Fig.10　The results of the 2D Haar wavelet transformation applied to the gravity anomaly data with 73 degrees rotation of the coordinate system. Top left part in (a) is average data, top right is horizontal component, bottom left is vertical component and bottom right is edge component. The original gravity anomaly data are shown in (b). The color bar indicates relative values.

（2010）：重力異常からみた東北本州弧近く構造の特徴， 月刊地球，32（6），373-382. 駒 澤 正 夫・ 広 島 俊 男・ 村 田 泰 章・ 森 尻 理 恵・ 牧 野 雅 彦 （2000）：地質調査所の重力基本図・日本重力 CD-ROM， 情報地質，11（2），86-87. 駒澤正夫（2013）：日本重力異常グリッドデータ， 日本重力 72，11-22 イヴェス ニイベルゲルト（Yves Nievergelt）（2004）：ウェー ブレット解析の基礎（松本 忠・雛元孝夫・茂呂征一郎 訳），森北出版株式会社.