第一工業大学研究報告第33号（2021）pp.51-56

ローカライゼーションプログラムの作成
田中1

龍児 1・外山

第一工業大学（〒899-4332

泉 2・石澤 直樹 3

鹿児島県霧島市国分中央 1-10-2）

2 砂防エンジニアリング株式会社（〒350-0033
3RTK 研究会（〒899-2513

埼玉県川越市富士見町 31-9）

日置市伊集院町麦生田 813-24）

Development and Evaluation of Localization Algorithm
Ryoji TANAKA1, Izumi TOYAMA2, and Naoki ISHIZAWA3
1 Department of Regional revitalization design and Engineering, Daiichi Institute of Technology
2 SABO Engineering Co., Ltd.
3 RTK Research and Study GroupIn recent years, low-cost and highly precise Global Navigation Satellite Systems (GNSSs) receiver has beenmarketed. In addition the receiver supports the Quasi-Zenith Satellite System (Michibiki), GPS, GLONASS, BeiDou.
We can easily get longitude, latitude and altitude values with the system as Realtime Kinematic (RTK) and open source
software (FOSS) package called RTKLIB. However, there is a difference of several centimeters between the GNSS
survey and the TS survey due to the difference in the coordinate system. We have developed an algorithm for affine
transformation using the least squares method and created a program.
Key Words: RTK-GNSS, GCP, Affine

1. まえがき
国土交通省 ICT 建設機械精度確認要領(案)では，精度検証にトータルステーション（以下 TS）を用いることを前提としている．計測の元データとなる３次元座標は，TS を設置した基準点からの相対位置として算出されるため，出来形値の誤差精度は，基準点の残差に依存する．このため，構築物の施工精度を確実に確保するには，設計照査の段階で，出来形管理用 TS を設置する可能性がある各基準点の 3 次元座標と，TS またはRTK-GNSS を用いて計測される 3 次元位置座標との残差，あるいは各基準点に対し GNSS 座標系上で算定された 3 次元位置座標との残差を確認するとともに，基準点の位置座標を包括する面との残差（回転，移動，大きさ，水平ゆがみ，垂直ゆがみ）を確認することが必要である 1）．この残差が比較的大きいと判断する場
合は，基準点の確認を行う必要があると判断され，この残差が比較的小さい場合は，残差の影響を最小限に留める対応として，GNSS 座標系と現場座標系にローカライゼーションする（図 1）．
一方，ローカライゼーションの計算方法については，明確に決められているわけではないが，国土地理院の公共測量座標変換マニュアル 2）においては，変換の数学モデルとしてアフィン変換を採用している．アフィン変換の計算式については，多くの文献に見られるが，そのほとんどが，画像処理などに用いられるような 2次元の拡大縮小，回転，平行移動で，その変換パラメータが既知である場合の解説に終始している．土木設計など実務的な処理は，3 次元座標の変換であり，過剰観測をともなうため最小二乗法の適用も考慮しなければならない．しかしながら，その解説はほとんど見当たらず，既存のシステムもブラックボックス化され，新たなシステムの開発を困難にしている．
本研究では，RTK-GNSS データの現地座標へのローカライゼーションプログラムの開発を行うため，その計算式を明らかにし，数値データを代入して確認した．

図 1 ローカライゼーションのイメージ

第一工業大学研究報告 第33号 （2021）

𝑎𝑑𝑔𝑗

2. ローカライゼーションの計算式
2.1 ヘルマート変換とアフィン変換
点 𝑥, 𝑦 を点 𝑋, 𝑌 に移すことを座標変換というが，拡大・縮小，移動，回転を組み合わせることで，様々な変換を行うことができる．
2 次元座標の場合で示すと，式（１）になる．式（１）をアフィン変換という．また，アフィン変換において，𝑎 𝑑, 𝑐𝑏の特別な場合をヘルマート変換という．
𝑋𝑎𝑥𝑐𝑦𝑒𝑌𝑏𝑥

𝑑𝑦𝑓

1𝑏𝑒ℎ𝑘𝑐𝑓𝑖
𝑙0001𝑥𝑥𝑥𝑥𝑦𝑦𝑦𝑦𝑧𝑧𝑧𝑧

1
111𝑋𝑌𝑍𝑑𝑦𝑔𝑧

𝑗𝑌𝑏𝑥

𝑒𝑦ℎ𝑧𝑘𝑍𝑐𝑥

𝑓𝑦𝑖𝑧

𝑙𝑥

1𝑦𝑧

1𝑎𝑑𝑔𝑗𝑏𝑒ℎ𝑘0001𝑌𝑌𝑌𝑌𝑍𝑍𝑍𝑍𝑧𝑧𝑧𝑧

1
1
11𝑎𝑑𝑔𝑗𝑏𝑒ℎ𝑘𝑐𝑓𝑖𝑙00016𝑥𝑋

2299.949
1118.903
238.159292.7695

0.967
0 1863.691
13.202 16066.745 2181.702 14877.184 5227.604 1

1111

変換パラメータ A の計算
4818.084
834.309
7.94863680.119
9.145
2.070
2373.236 6664.633 2043.360
2360.798 4846.740 2009.967

𝑥4

式(4)より，変換パラメータを求めると，式(5)になる．

1111より，

𝐴𝑥𝑋

1.265682
0.322720
1.197714
3538.480.297768
0.623819
1.551843
1968.430
1.273584
1.828210
7.838460
4673.250

𝑦𝑦𝑦𝑦

5𝑋

834.3085
7.94869.14532.0706664.633 2043.3604846.74 2009.9671111𝑥𝑥𝑥𝑥𝑥

4818.0843680.1192373.2362360.798

3

式(3)における未知数（変換パラメータ）は，𝑎~𝑙の12 個であるが，3 次元座標を用いて求めるには，少なくとも４点の座標が必要である．４点より多い場合，すなわち過剰観測の場合は，最小二乗法を適用するが，ここでは簡単のために，４点の座標より変換パラメータを求める．
４点の座標の変換は，式(4)で表される．
𝑋𝑋𝑋𝑋

1111

変換前座標𝑥と変換後座標𝑋を次のように定める．

2𝑐𝑓𝑖𝑙𝑍𝑍𝑍𝑍

2.3 計算例１
説明のために，式(5)を式(6)のように表す．
𝐴𝑎𝑥𝑌𝑌𝑌𝑌

もちろん，式(5)における逆行列が存在する必要がある．

2.2 3 次元のアフィン変換
式(1)と同様に点 𝑥, 𝑦, 𝑧 を点 𝑋, 𝑌, 𝑍 に移す場合の式は，式(2)になる．これを行列表記すると，式(3)になる．
𝑋𝑋
𝑋
𝑋𝑋

1.0842𝐸
2.1684𝐸
01

計算結果を検算すると，

𝑥𝐴

4818.084 834.3085
7.94863680.1199.14532.070
2373.236 6664.633 2043.360
2360.798
4846.74 2009.967

11
1
1

19
19

田中・外山・石澤：ローカライゼーションプログラムの作成

1.265682
0.322720
1.197714
3538.480.297768
0.623819
1.551843
1968.430
1.273584
1.828210
7.838460
4673.250

2299.949
1118.903
238.159292.7695

1.0842𝐸
2.1684𝐸
01

𝜕𝑓
𝜕𝑎⎧𝑓
⎪
⎪⎪⎪
⎪
⎪𝑓
⎪

19
19

𝜕𝑓
𝜕𝑏

⎨
⎪⎪𝑓
⎪
⎪⎪
⎪
⎪
⎩

0.967
0 1863.691
13.202 1
6066.745 2181.702 1
4877.184 5227.604 1

𝑋

𝜕𝑓
𝜕𝑐

834.3085

1.265682
0.322720
1.197714
3538.48

7.9486

0.297768
0.623819
1.551843
1968.430

2299.949

1.0842𝐸
2.1684𝐸
01

∆𝑑

𝜕𝑓
𝜕𝑗

∆𝑗

𝑋

∆𝑏

𝜕𝑓
𝜕𝑒

∆𝑒

𝜕𝑓
𝜕𝑘

∆𝑘

𝑌

∆𝑐

𝜕𝑓
𝜕𝑓

∆𝑓

𝜕𝑓
𝜕𝑙

∆𝑙

𝑍

𝜕𝑓
⎧
⎪ 𝜕𝑎
⎪ 𝜕𝑓
⎪
⎪ 𝜕𝑑
𝜕𝑓
⎨ 𝜕𝑔
⎪
⎪ 𝜕𝑓
⎪ 𝜕𝑗
⎪
⎩⋮

1

1.273584
1.828210
7.838460
4673.250𝜕𝑓𝜕𝑑𝜕𝑓𝜕𝑔∆𝑔

𝜕𝑓𝜕ℎ∆ℎ𝑣

7𝑣
𝜕𝑓𝜕𝑖∆𝑖𝑣⋮となり，変換が確認できた．
次に𝑥の 1 行目の座標から𝑋の 1 行目の変換を確認す
る．
4818.084∆𝑎1919𝑥𝑦8𝑧
1

したがって，式(7)は式(9)の観測方程式になる．

0.967 0 1𝑥 ∆𝑎⎧𝑥 ∆𝑏⎪⎪𝑥 ∆𝑐𝑥 ∆𝑎⎨𝑥 ∆𝑏⎪⎪𝑥 ∆𝑐⎩となり，変換が確認できた．
式(2)あるいは，式(3)における未知数（変換パラメータ）は，𝑎~𝑙の 12 個であり，3 次元座標を用いて求めるには，少なくとも４点の座標が必要である．４点より多い場合，すなわち過剰観測の場合は，最小二乗法を適用することになる．

𝑦 ∆𝑑𝑦 ∆𝑒𝑦 ∆𝑓𝑦 ∆𝑑𝑦 ∆𝑒𝑦 ∆𝑓𝑧 ∆𝑔 ∆𝑗 𝑋𝑧 ∆ℎ ∆𝑘 𝑌𝑧 ∆𝑖 ∆𝑙 𝑍𝑧 ∆𝑔 ∆𝑗 𝑋𝑧 ∆ℎ ∆𝑘 𝑌𝑧 ∆𝑖 ∆𝑙 𝑍
⋮𝑓𝑓𝑓𝑓𝑓𝑓𝑣𝑣𝑣𝑣𝑣𝑣

式(9)を，2.4 観測方程式
式(2)あるいは，式(3)を観測方程式という．この場合の未知数は，𝑎, 𝑏, ⋯ , 𝑙であることに注意する．
式(2)の左辺をそれぞれ𝑓 , 𝑓 , 𝑓 とおくと，式(7)のように線形化される．
ここで， 𝜕𝑓 ⁄𝜕𝑎 , 𝜕𝑓 ⁄𝜕𝑑 , ⋯は，𝑓 を𝑎, 𝑑, ⋯ で偏微分して得られる関数の𝑎 𝑎 , 𝑑 𝑑 , ⋯ における値，𝑓 は式(2)の𝑎 𝑎 , 𝑑 𝑑 , ⋯における値，∆𝑎, ∆𝑑, ⋯は，それぞれ，𝑎, 𝑑, ⋯の修正量， 𝑋, 𝑌, 𝑍 は変換後の座標そして𝑣 はそれらの残差である．
𝜕𝑓 ⁄𝜕𝑎 , 𝜕𝑓 ⁄𝜕𝑑 , ⋯の偏微分を実行すると，式(8)になる．

𝐴𝐿𝑥𝑥
⎛𝑥
⎜𝑥
⎜
⎜𝑥𝑥
⎝⋮𝑦𝑦𝑦𝑦𝑦𝑦⋮𝑧𝑧𝑧𝑧𝑧𝑧⋮11⎞1⎟1⎟, 𝑋
1⎟
1
⋮⎠

∆𝑎∆𝑑
∆𝑔
∆𝑗𝑋
𝑋𝑓𝑓

𝑌𝑌𝑓𝑓𝑓𝑓

⋮

𝑉𝑣𝑣
𝑣𝑣𝑣𝑣
⋮𝑣𝑣
𝑣

𝑍𝑍

∆𝑏
∆𝑒
∆ℎ∆𝑘,

∆𝑐
∆𝑓,
∆𝑖
∆𝑙

9

第一工業大学研究報告 第33号 （2021）と置いて行列表記すると，
𝐴𝑋𝐿𝑉

10

2.5 正規方程式
式(10)より，正規方程式を作成すると，式(11)になる．
𝐴 𝐴𝑋

1.273584
1.828217.838464673.25𝐴𝐿

1.0842𝐸
2.1684𝐸
011919f i 0 の計算
f10 ～ f12 0 （４点の変換後の座標）については，計11

算例１の結果を用いる．
ただし，重量（重み）は 1 とする．
2.6 解 𝑋式(11)を１次方程式の解法で解けば，式(12)になる．
𝑋𝐴𝐴𝐴𝐿

５点目の座標𝑓12

このとき，残差𝑉は，
𝑉𝐴𝑋𝐿

𝐴 𝐴𝐴𝐴

𝐼 𝐿, 𝑓,𝑓

2360.790

4846.700

1.265682
0.32272
1.197714
3538.480.297768
0.623819
1.551843
1968.4313

1.273584
1.828217.838464673.25となる．
ただし，𝐴′は𝐴の転置行列，𝐼は単位行列である．
292.772
2.7 計算例２（計算例１とは文字が異なっていることに注意する）
変換前座標𝐴と変換後座標を次のように置く．
𝐿4818.084
834.3093680.119
9.145⎛⎜2373.236 6664.633
2360.798 4846.740⎝2360.790 4846.700

𝐴

𝑋
𝑋

𝑌𝑌
⋮

7.94862.070
2043.360
2009.967
2009.970

11⎞
1⎟
1
1⎠

0.967
0 1863.691
13.202 1⎞
6066.745 2181.702 1⎟
4877.184 5227.604 1
4877.180 5227.600 1⎠

𝑋
𝑋
⎛
⎜𝑋
𝑋
⎝𝑋

𝑓𝑓𝑓𝑓
𝑓

0
⎛ 0
⎜ 00
⎝ 0.003

𝑍𝑍

2299.949⎛ 1118.903⎜ 238.159
292.7695⎝ 292.769

𝐴

𝑌
𝑌
𝑌
𝑌
𝑌𝑓𝑓𝑓𝑓𝑓00
000.0210001𝐴𝐴
1.265682
0.32272
1.197714
3538.480.297768
0.623819
1.5518431968.43𝑍𝑍𝑍𝑍𝑍𝑓𝑓𝑓𝑓𝑓

1
1
⎞1⎟
11⎠

00 ⎞0 ⎟00.098⎠

2360.798 2360.790
4846.740 4846.700
2009.967 2009.970
11

計算例１で求めた値（４点のみの変換）を用いる．
𝑐𝑓𝑖
𝑙

19
194818.084 3680.119 2373.236
834.3099.145 6664.633
7.949
2.070 2043.360
111

変換パラメータの初期値（ a0 , b 0 , c 0 ,  など）は，

𝑏𝑒ℎ𝑘

1.0842𝐸
2.1684𝐸
014877.159 5227.698 1

ただし，行列計算のために最後の列に１を追加してある．

𝑎𝑑𝑔
𝑗2009.970 153536155.018
42754386.181
14385533.53215593.027

42754386.181
92094876.852
33108404.059
17201.527

田中・外山・石澤：ローカライゼーションプログラムの作成

14385533.532 15593.027
33108404.059 17201.527
12255334.297
6073.316
6073.316
5.000

1.273563
1.82818
7.838352
4673.17292.770

1.84537𝐸 06
4.06256𝐸 07
2.74016𝐸 06
0.00768569

𝐴′𝐴

4.06256𝐸 07
5.17892𝐸 07
1.67601𝐸 06
0.001521037

2.74016𝐸 06
1.67601𝐸 06
6.95982𝐸 06
0.011233308

0.000888274
1.7578𝐸 05
0.000626006
3.270088832

𝐴𝐴′𝐴

0.000414622
0.000583756
0.001678863
1.523993845

𝑋

2.8 ４点以上のアフィン変換の計算の流れ
以上の計算式により作成したアフィン変換プログラ
ムの流れを図 2 に示す．実際のローカライゼーション
プログラムでは，基準点が３点以下の場合は，１点のみの場合は平行移動，３点まではヘルマート変換，４点以上はアフィン変換で計算されるようにプログラミングした．

0.000892568
2.72344𝐸 05
0.001150126
4.286889929

0.00020945
0.0002967
0.00110146
0.7703809

4.36966𝐸 06
6.18991𝐸 06
2.29792𝐸 05
0.016071814

0.0002095
0.000297
0.0011015
0.770414

2.06134𝐸 05
2.92002𝐸 05
0.000108402
0.075816873

変換パラメータ𝑎𝑑𝑔
𝑗𝑏𝑒ℎ𝑘𝑐𝑓𝑖
𝑙

0001

1.266
0.323
1.198
3538.478𝑎𝑑𝑔
𝑗𝑏𝑒ℎ𝑘𝑐𝑓𝑖
𝑙00010.298
0.624
1.552
1968.446

𝑋

1.274
1.828
7.838
4673.1740001

検算（５点目の座標の変換で）2360.794846.7 2009.97

1.265681
0.32272
1.197711
3538.48

19194877.170 5227.649 1

5.43213𝐸 07
7.69496𝐸 07
2.85665𝐸 06
0.001997961

𝐴𝐿𝐴′𝐴

0.0076857
0.00152104
0.0112333
32.580488

1.084𝐸
2.168𝐸
0110.297772
0.623813
1.551866
1968.45
図 2 ローカライゼーションプログラムの流れ

第一工業大学研究報告 第33号 （2021）

3. むすびアフィン変換の解説は，多くの文献でなされているが，ほとんどが式の説明に重きを置いており，実際の計算方法について記されたものはない．本研究では，ローカライゼーション計算の自動化を目的としているので，具体的な変換の方法を明らかにする必要があった．ローカライゼーションプログラム作成のために，ステップごとに数値を入力し，出力された結果を再入力し逆計算することにより，計算式の正当性を確認した．この中で測量の誤差処理で用いられる過剰観測の場合の最小二乗法の考え方は，写真測量におけるバンドル法や，CG 分野やロボット分野の SfM や SLAM などにも使われている．なお，プログラム言語は，Web上で共有できるように CGI（Common Gateway Interface）の Perl で作成した．

参考文献
1) 国土交通省：RTK-GNSS を用いた出来形管理要領（案）, 2018 年 3 月．
2) 国土地理院：UAV を用いた公共測量マニュアル（案），2017 年 3 月改正．