Vol. 22, 2004

海洋情報部技報

海底地殻変動観測における局位置解析ソフトウェアの開発
藤田雅之、佐藤まりこ：航法測地室
矢吹哲一朗：海洋研究室
Development of seafloor positioning software using inverse method
Masayuki FUJITA, Mariko SATO: Geodesy and Geophysics OfficeTetsuichiro YABUKI: Ocean Research Laboratory

1

2

はじめに

局位置解析ソフトウェア SGOBS

2.1

海洋情報部では，東京大学生産技術研究所との概要

技術協力の下，GPS 音響結合方式による海底地殻SGOBS は，KGPS 解析結果と音響解析結果を変動観測の技術開発を行っている（浅田及び矢吹，結合して，海底局位置を決定するソフトウェアで2001；矢吹，2002）．観測システムの概念図を第ある．全体の解析の流れを第 2 図に示す．

1 図に示す．

SGOBS への具体的な入力データは，GPS アン

このシステムにより取得されたデータの解析過

テナ位置，音波走時，動揺計測パラメータ（ヘデ程は，大きく３つのプロセスに分けられる．一つ

ィング，ロール，ピッチ），及び海中の音速度構造

目 は ， 船 の 位 置 を 求 め る キ ネ マ テ ィ ッ ク GPS

プロファイルである．

（KGPS）解析，二つ目は，船と海底局間の音波

まず，動揺計測パラメータを用いて，KGPS 解の往復走時を求める音響解析，そして三つ目は，析により求められた GPS アンテナの位置を，観これら二つの解析結果を結合し，海底局の位置を測支柱のもう一方の端に取り付けた海中のトラン求める局位置解析である．

スデューサ位置へと変換する．このトランスデュ

現在我々は，この最後の段階にあたる局位置解ーサ位置と音波走時から，音速度構造を介して，

析 の た め の ソ フ ト ウ ェ ア 「 SGOBS (Seafloor

幾何学的原理に基づき海底局位置を決定する．そGeodetic OBServation)」の開発を行っている．の際，局位置と共に，音速度の時間変化係数も推

本稿では，現行バージョン（Ver. 2.5）について，定パラメータとしている．これらのパラメータのその基本構造を紹介すると共に，この中で行って解法としては，ベイジアン最小自乗法に基づくイ

動揺計測装置
測量船

Fig.1

Schematic Figure of the observation
system.

第2図
Fig.2

−50−

解析の流れ．
Flow chart of data analysis.

音速度構造

観測システム．

トランスデューサ位置

第1図

動揺補正

海底基準局
（ミラートランスポンダー）

アンテナ位置CTD

ＫＧＰＳ解析ＫＧＰＳデータ

音響トランスデューサー
XBT

音響解析

局位置解析

海底局位置

精密測位用GPS受信アンテナ

海底局までの往復走時

音響測距データ

いる音速度の補正方法について考察を行う．

Vol. 22, 2004

海洋情報部技報
ンバージョン法(松浦，1994)を適用している．
初期値

2.2

(複数局)

(複数時間ウィンドウ)局位置

音速度時間係数

ベイジアンインバージョン法について

この手法は，正規分布をもつ誤差 e を含む観測解析終了

収束？

データ y0 とモデルパラメータ x の間に，A を偏微

第3図局位置解析ソフトウェアのアルゴリズム．

分行列とした線形の観測方程式

Fig.3

Algorism for seafloor positioning software.

y 0 = Ax + e（１）が成り立つと仮定し，モデルパラメータの初期値からの誤差分布がアプリオリに与えられる場合に，Nモデルパラメータ解 x̂ が次式によって計算されるというものである(松浦，1994)．

Wxˆ = x 0 + DA t ( E + ADA t ) −1 ( y 0 − Ax 0 ) （２）
ここで，x0 はモデルパラメータ初期値，D はモデE
S

ルパラメータ解の初期値からのずれを特徴づける共分散行列，E は e の誤差分布を特徴づける共分散行列である．
一般には，観測データとモデルパラメータの間第4図の関係は非線形であることがほとんどであるが，Fig.4

その場合でも真値近傍における微小な変化に対し

海底基準局配置図と測線例．
Distribution of the seafloor stations withan example of one-day observation lines.ては線形関係が成り立つと仮定し，与えた初期値と真値との差を解として求める．これは式（1）

れ同一の音速度構造で局位置解を求めた後，その及び（2）において， x0=0 として、y0 と x をそれ残差データを用いて，一定の時間ウィンドウ毎にぞれ初期値に対するΔ値と考えることに相当する．

音速度の時間変化係数を求める．このサイクルをそして，実際の非線形の影響については、繰り返局位置が収束するまで繰り返し，最終的な局位置し計算（iteration）により真値に収束させる．

解を求める．音速度係数の時間ウィンドウ幅は，ベイジアンの手法，すなわちモデルパラメータ測線単位から最大 1 日単位まで選択できる．の共分散行列 D を導入する意味は，初期値の信頼このように局位置と音速度構造を分けて推定す性に応じて拘束をかけられることである．その最る最も大きな理由は，それぞれの推定パラメータも簡単な適用例として，パラメータの固定，推定毎に，用いる走時データの取捨あるいは重み付けの選択を D の対角成分により自由に設定できるに違いを持たせたいことである．例えば，後に述ことがある．

べる伝播距離による重み付けは，海底局毎に違うため，全局を同時に推定することとは両立しない．

2.3

アルゴリズムSGOBS の簡易なアルゴリズムを第 3 図に示す．

2.4

我々が展開している海底基準点では，1 点につ次に，具体的な観測データとモデルパラメータにき原則 4 局の海底局を東西南北に配置し，測線観ついて述べる．上記アルゴリズムに示したように，測を行っている．第 4 図に，海底局配置と 1 日の観測方程式は局位置解を求める部分と音速度構造測線例を示す．まずこれら複数局についてそれぞを求める部分との 2 つに分かれている．

−51−

モデルパラメータと観測方程式

Vol. 22, 2004

海洋情報部技報(a) 海底局位置の推定

その伝播経路の速度構造の情報を含んでいる．そまず局位置解について，式（１）の観測データ

こで SGOBS では，音速度構造を推定パラメータy0 は，例えば海底局 A について，以下のように記とすることにより，誤差の補正を試みている．こ述される．

れは，GPS 測位で通常行われている大気遅延量の ∆t1( A) t∆
2 ( A) y0 = M ∆t 
 n ( A) 

推定などと同様の考え方である．
音速度構造を求める部分の観測データ y0 は，海底局を 4 局（A〜D）とすると，例えばタイムウ（３）ィンドウ[tk(A),tl(D)]について，以下のように記述される．

 ∆tk ( A)  ∆ tk ( B ) 
 ∆tk ( C ) y 0 =  ∆tk ( D )  M 
 ∆tl ( D ) 

右辺は，KGPS 解析と動揺補正により求められたトランスデューサの位置と海底局 A の初期位置
座標から音速度構造を用いて計算された往復走時と，観測された往復走時の差（O-C）の時系列を表している．
また式（１）のモデルパラメータ x は，海底局位置座標（ローカル座標系の 3 成分）の初期値からの補正量として次のように表現される．

（５）

音速度構造推定におけるモデルパラメータは， ∆x A 
x =  ∆y A  ∆z  A

現行バージョンでは，音速度の時間変化を二次式で表したときの係数としている．すなわち，平均（４）音速度の時間関数 V(t)をV (t ) = V (t0 ) + a0 + a1 (t − t0 ) + a 2 (t − t0 ) 2（６）(b) 音速度構造の推定（誤差の補正手法）と表すと，モデルパラメータ x は，その係数の補海底局位置を精密に求めるためには，要求精度内正値として次式のように表現される．

で正確な海中の音速度構造が必要である．我々の ∆a0 
x =  ∆a1  ∆a  2

通常の観測では，原則として一日の測距観測前と後に CTD 観測を行い，測距観測中に XBT による水温観測を行っている．しかしながら，これらの（７）

観測機器の精度は十分ではなく，また，観測値（水式（６）で t0 は基準時刻，V(t0)は音速度初期値の温，電気伝導度，水圧）から音速度を求める経験

t0 における値である．

式 （ 例 え ば ， Del Grosso, 1974; Wilson, 1962;

なお，ここで時間変化係数を定義するにあたり，Chen and Milleo, 1977）の選択によっても有意な深さ方向のプロファイルの形状は時間的に一定と違いが生ずる．さらに，音速度は時間空間で変化仮定している．現実には浅い部分ほど時間変化がするため，これらを全てカバーする観測は実際上大きいなど，この仮定は厳密には成り立っていな不可能である．したがって観測値を確定値として，いと考えられるが，数値実験によると，局位置推

これのみからセンチメートルレベルの測位を行う

定結果はほぼ平均音速に依存していると言ってよ

ことは大変困難であるといわざるを得ない．

いため（佐藤及び藤田，2004），現行精度におい

しかしながら音波走時は距離と速度の関数なのては問題ないと判断される．

で，走時データには距離の情報を含むと同時に，−52−

また，入力データとしての音速度構造は，CTD

Vol. 22, 2004

海洋情報部技報
等の観測値から 1m 層で与えているが，プログラ定パラメータ間の分離，すなわち局位置解の一意ム内で，これを 200m，400m，800m，1600m，性の確保である．既に述べたように，SGOBS V2.5それ以深という層毎に勾配を平均化して計算してでは，音速度変化の時間二次式の係数をパラメーいる．この場合，局位置推定結果は，現在の数 cmタとして補正する機能を有しているが，我々の取という精度においては，1m 層をそのまま用いる得するデータの中に，局位置と音速度構造を要求場合と比較してほとんど差異がないことを確認し精度内で分離できるだけの情報が含まれているかており，これにより計算時間が大幅に短縮されるどうか，またそれが含まれていたとして，その情メリットがある．

報を適切に取り出すための最善の手法はどういうただ将来さらに高精度の議論になった場合，層ものかが問題となる．

厚を含めたこれらの仮定を見直す必要が出てくるかもしれない．

その一つの指標に，初期値に対する依存性がある．ある範囲で異なる初期値を与えた場合に，そこから求められる局位置解が一意的であれば，そ2.5

データの重み付けによる誤差の軽減の手法は一意性が確保されているといえる．逆に，

音速度誤差については，上記手法で正確に補正異なる初期値から求められた局位置解に有意な差しようと試みてはいるが，時間関数近似の限界やがあれば，音速度と局位置のパラメータの分離が空間不均質の影響及びパラメータ分離の不十分さうまくできていないことを意味する．

等の理由により，補正後も有意な誤差を含む場合

このような観点から，ここでは音速度の初期値が少なくない．この影響は，伝播距離が長くなるを変えて局位置解を比較すると共に，現行の機能ほど大きくなるため，伝播距離の長いデータの重で行える範囲で補正ストラテジーを変え，より良みを適切に下げることが，結果の改善につながる

い方法について検討する．

可能性がある．
この観点から，SGOBS では，局位置推定の際，3.1それぞれの海底局毎に，観測データの誤差行列 Eの対角成分に重み付けを行えるようにしている．

比較方法

まず，音速度構造の初期値として次の二通りを与えた．
(C) 朝夕の CTD 観測値を時間で線形補間した

具体的には，海底局への音波の入射角を θ として，cosnθ の重みを用いている．

もの(X) (C)について，水温のみ XBT 観測値で置きこのように遠くのデータの重みを下げることにより、データの幾何学的配置が実効上改善される換えて計算したもの
また，音速度構造の補正は，以下の 4 通りの方効果もある．佐藤及び藤田（2004）は，音速度構造に誤差を含む場合でも，海底局の周りのデータ法で行った．の幾何学的配置が均等であれば，水平成分につい(A)補正しなかった場合ては，誤差が相殺され，精度劣化をある程度防ぐ(B)測線毎に補正した場合ことができることを示している．現在の観測方法(C)1 日毎に補正した場合

では，第 4 図に示したように，それぞれの局に対(D)まず 1 日毎に補正した後，これを初期値としするデータ配置は一方に偏ってしまうため，このてさらに測線毎に補正した場合重み付けが有効である場合が多い．

以上の初期値と補正手法を組み合わせた 8 通りの場合について，局位置解を求め，その結果を比3較した．

音速度初期値依存性と補正ストラテジー

局位置推定を行う際に重要なことの一つは，推−53−

Vol. 22, 2004

海洋情報部技報

結果

Difference (m)3.2

ここでは比較例として，当庁が 2002 年 5 月 13日〜16 日の 4 日間，宮城県沖海底基準点で取得し0.3(a)Off‑Miyagi 0205 (EW)

0.2
0.1
0‑0.1

た観測データ（藤田及び佐藤，2004）について行‑0.2

った結果を例示する．

‑0.3AX AC

BX BC

CX CC

DX DC

第 5 図に比較結果を成分毎に示す．またそれぞDifference (m)れの場合の解析後の走時残差を第１表に示す．なお，図のプロットは，海底局 4 局の座標値の平均について，CTD の初期値（C）から（D）の手法で求めた解（DC と表記，他も同様）を基準とし0.3(b)
Off‑Miyagi 0205 (NS)0.2
0.1
0‑0.1‑0.2て，それぞれから求められた平均座標値との差を‑0.3AX AC

BX BC

CX CC

DX DC示している．また図のエラーバーは，局毎の差のDifference (m)自乗平均値を示している．
図を見ると，補正しなかった場合（A）には，上下成分に 50cm 程度，水平成分でも 30〜40cm程度の差が見られる．これだけを見ても，観測値0.6(c)
Off‑Miyagi 0205 (UD)0.40.20‑0.2‑0.4のみからセンチメートルレベルの測位を行うこと‑0.6
AX ACは困難であることがわかる．
次に，測線単位に補正した場合（B）には，差第5図

CX CCDX DC

音速度初期値及び補正ストラテジーの違いによる局位置推定結果の比較：(a) 東西成

はかなり小さくはなるものの，局毎の大きなばら分、(b) 南北成分、(c) 上下成分．

つきを含めまだ有意な差が残っている．これに対して，1 日単位で補正した場合（C）を見ると，

BX BC

Fig.5

その差は 1cm レベルとなっており，どちらの初期値を用いても，結果はほぼ変わらない，すなわちComparison of positioning results due todifferencesof

initialvaluesstrategies

forestimationofandsoundvelocity: (a) EW component, (b) NS

初期値依存性がほとんどなくなっていることがわ

component, (c) UD component for the Off

かる．

Miyagi reference point.

さらに，（C）と（B）の組み合わせ，すなわち一旦 1 日単位の補正を行ってから，これを改めて初期値として測線単位で補正した場合（D）にも，第1表

初期値依存性は見られない．これは（C）の補正を行った時点で，初期値がほぼ同じになるため当音速度初期値及び補正ストラテジーの組み合わせによる解析走時残差(msec)．

Table 1

Travel time residuals (msec) in analyses

然である．（C）と（D）の結果を比べると，東西

for the combination of initial values and

成分に約 10cm の差があるが，どちらが真値に近

strategies

いかはこれだけでは判断できない．ただ第 1 表よcorrection.

り，最後に測線毎の補正を行うことにより，走時

（A）

（B）

（C）

（D）

（C）

0.397

0.1190.166

0.120（X）

0.277

0.118

0.1640.119

残差が小さくなっていることがわかる．
（C）と（D）の手法をさらに詳細に比較するため，4 日間のデータを 1 日毎のサブセットに分け，得られた解を比較した．第 6 図，第 7 図は，13−54−forthesoundvelocity

Vol. 22, 2004(a)Difference (m)ALL(4D)

0.4
0.3
0.2
0.1
0‑0.1
‑0.2
‑0.3
‑0.4(b)ALL(4D)

第 6 図

Fig.6

Strategy (C) ‑ EW

1314

Difference (m)

0.4
0.3
0.2
0.1
0‑0.1
‑0.2
‑0.3
‑0.4

151415Strategy (D) ‑ EW(a)

ALL(4D)

Strategy (C) ‑ NS

130.40.3
0.2
0.1
0‑0.1
‑0.2
‑0.3
‑0.4

16

Difference (m)Difference (m)海洋情報部技報

0.4
0.3
0.2
0.1
0‑0.1
‑0.2
‑0.3
‑0.4

16(b)

第 7 図

1415161516

Strategy (D) ‑ NS

ALL(4D)

音速度補正ストラテジー（C）による全日

131314

音速度補正ストラテジー（D）による全日

解と 1 日解の比較：(a) 東西成分、(b) 南北

解と 1 日解の比較：(a) 東西成分、(b) 南北成分．

成分．

Comparison of positioning results between

Fig.7

Comparison of positioning results between

solutions from all-day data and one-day

solutions from all-day data and one-day

data in the case of the strategy (C) fordata in the case of the strategy (D) for

estimation of sound velocity: (a) EW

estimation of sound velocity: (a) EW

component, (b) NS component for the Offcomponent, (b) NS component for the Off

Miyagi reference point.

Miyagi reference point.
Average sound velocity (May 16, 2003)

水平成分について示したものである．図から明ら

1480

SV[m/s]日〜16 日までの 1 日解の，4 日解に対する差を，1481

かに，（D）の手法により 1 日解の再現性が改善されていることがわかる．

1477
9:00:00

期値依存性はなくなるが，さらに測線単位の補正判断される．
3.3

1479

1478

以上の検討結果から，（C）の手法によりほぼ初を行う（D）の手法が，より良い結果を与えると

SV0[m/s]
SV1[m/s]
SV2[m/s]

10:00:00 11:00:00 12:00:00

13:00:00 14:00:00

Time

第 8 図

音速度補正ストラテジー（D）による解析過程における音速度時間変化（2003 年 5月 16 日）：SV0 は初期値、SV1 は 1 日補正議論及び問題点

第 8 図は，5 月 16 日の音速度補正に関して，（ D）の手法をとった場合の初期値，一日推定値，及び値、SV2 は最終の測線毎補正値．
Fig. 8

最終の測線単位の推定値をプロットしている．これを見ると，まず一日単位でバイアス的な音速度の補正がなされ，さらに測線毎にさらに細かい変動を補正しているように見える．しかしながら，最後の測線単位の補正曲線では，部分的には現実−55−

Temporal variation of sound velocity inthe analysis with the strategy (D) for
May 16, 2003; SV0 is the initial values,
SV1 and SV2 are ones corrected with aday through and ones finally correctedwith every line, respectively.

Vol. 22, 2004

海洋情報部技報のものとは考えにくい急激な変化も見られる．他

い．

例では，こういった現象がさらに顕著に見られる
参 考 文 献ものもある．
現実の音速度変化は複雑であり，これをある特浅田昭，矢吹哲一朗：熊野トラフにおける長期地定の関数にあてはめて補正することには限界があ

殻変動観測技術の高度化，地学雑誌,110（4）,る．
（D）の手法の考え方は，一日の変化について二次式で近似できない部分を測線単位の推定で補おうとするものであり，ある程度の妥当性と有効529−543,（2001）．
Chen, C.T., and F.J. Milleo: Speed of sound in性をもっていると考えられるが，部分的にはデー

seawater at height pressures, The Journal of

タに含まれるなんらかの別の要因を，音速度構造

the Acoustical Society of America, 62, No.5,に無理矢理押し付けている可能性があることには1129-1135, (1977).

留意しなければならない．

Del Grosso, V.A.: New Equation for the Speed of

他方，我々の目的は音速度構造を正確に求める

Sound in Natural Water (with Comparison to

ことではなく，局位置を正確に求めることにある．

other

例えば KGPS から求められた高さ成分に時間ドEquations),

The

Journalof

theAcoustical Society of America, 56, No.4,

リフト誤差が含まれているとすれば，数値解析的

1084-1091, (1974).には音速度補正値によって吸収されてしまい，結果を改善する可能性もある．その場合，たとえ推藤田雅之，佐藤まりこ：海底地殻変動観測データ定された音速度が正しくなくても，手法としてはによる局位置再現性の評価，海洋情報部研究報有効と考えるべきであろう．

告 ， 40，72-78，(2004)．

音速度補正手法についての本稿の検討結果は，

松浦充宏：インバージョン解析法，「現代測地学」

かなりの事例で確認したが，十分な一般性をもつ

日本測地学会，477−482, (1994)．

かどうかは，さらに検証しなければならない．また音速度の補正関数については，物理的合理性を佐藤まりこ，藤田雅之：海底地殻変動観測におけ含め，より適切なものがないか検討する必要がある海中音速誤差の局位置への影響について，海る．いずれにしても完全な補正を行うことは不可

洋情報部技報，22，42-49，（2004）．

能であるが，方法の工夫により局位置推定結果のWilson, W.D.: Extrapolation of the Equation for

精度，信頼性を高める努力を行っていきたい．

the Speed of Sound in Sea Water, The4

Journal of the Acoustical Society of America,
おわりに34, No.6, 866, (1962).

本稿では現在開発中の海底地殻変動観測におけ矢吹哲一朗：海底地殻変動観測を目指した音響技

る局位置解析ソフトウェアについて紹介すると共に，音速度の補正手法に焦点を当て，その妥当性等について議論した．現行バージョンは，現在取得しているデータからある程度目標精度を満たす事例が確認されており，今後の基盤となりうる段階には達していると考えている．しかしながら，まだまだ課題は多く，さらに精度のよい安定した結果を目指して改良を重ねていかなければならな−56−

術開発，水路部研究報告，38，47-58，（2002）．