ヘルマート変換式とアフィン変換式の係数考察

|1
ヘルマート変換式とアフィン変換式の係数考察
-（実例）平城京・長岡京・平安京の「尺」・「条坊の方位」の計算目次
1.ヘルマート変換
1.1 ヘルマート変換式 I
1.2 ヘルマート変換式 II
2.アフィン変換
2.1 アフィン変換式 I
2.2 アフィン変換式 I の近似式
2.3 アフィン変換式 II
3.コロケーション(3 パラメータ変換)
[付録 A]ヘルマート変換式の代数的特性
[付録 B]最小二乗法によるパラメータの決定
[付録Ｃ]ヘルマート変換・アフィン変換の逆式
[付録 D]8 パラメータのアフィン変換
[付録 E]最小二乗法による幾何学的要素の算出
|2
1.ヘルマート変換
1.1 ヘルマート変換式 I
測量では、平面の座標変換において 1 つの任意座標系 XY から国家座標系 xy に変換したい場合や、解析写真
測量においてコンパレータ座標 XY から指標を用いて写真座標系 xy に直したい場合、Helmert 変換と呼ばれる、
座標軸の回転(θ)と座標原点の平行移動量(xo , yo )で表す変換式が多く利用される。
図 1 に示すとおり、座標軸の回転（θ）および座標原点の移動 ( xo , yo ) を考慮して、座標系 XY から他の座標
系 xy に変換する場合、それらの関係は次式で表される。もしも、パラメータの計算を最小二乗法で求めたい場
合には、式(1-1)では左辺の座標は観測値でなければならない。もし、その式で右辺の X,Y が観測値の場合には、
条件方程式にし、線形化した誤差方程式を用いてパラメータを計算すればよい。
(1-1)
�
x = X ∙ cosθ − Y ∙ sinθ + xo
y = X ∙ sinθ + Y ∙ cosθ + yo
これを行列に置き換えると、
x
cosθ
�y � = �
sinθ
(1-2)
ここで、
xo
−sinθ X
� � � + �y �
o
cosθ Y
a = cosθ
b = sinθ
�
c = xo
d = yo
(1-3)
とおけば、式(1-1)は次のように簡単に表せる。
(1-4)
�
又は
x = aX − bY + c
y = bX + aY + d
x
a
�y� = �
b
(1-5)
Y
c
−b X
�� � + � �
d
a Y
y
Y
y'
y
x'
X
Θ
xo
y'
Ysin
Θ
または
図 1 座標軸の回転
X
Xsin
Θ
X
yo
Ysin
Θ
Xsin
Θ
Xcos
Θ
P
Ycos
Θ
Θ
Θ
Y
Yco
sΘ
P
Y
X
Xcos
Θ
x
x'
x
図 2 座標原点の移動
|3
x
X
�y� = �
Y
(1-6)
a
−Y 1 0 b
�� �
X 0 1 c
d
式(1-1)～式(1-6)は、二次元ヘルマート変換式と呼ばれている。3 点以上の（x,ｙ）と（X,Y）が与えられたとき、
パラメータ（a,b,c,d）は簡易的に最小二乗法により解ける。そして、式(1-3)より３つの要素は次式で計算できる。
b
tanθ =
a
� xo = c
yo = d
(1-7)
ただし、2 点でも解けるが、最小ニ乗法は適用できない。ただし、４個のパラメータに対し３個の要素で表して
いるので、数学的に不完全な解といわねばならない。
[応用例]
1）ヘルマート変換は測地学者ヘルマート(Friedrich Robert Helmert, 1843-1917)によって考え出された座標系
の二次元及び三次元変換である。
2) 写真測量では、内部標定であるコンパレータ座標系（CRT 座標系）から写真座標系（指標中心座標）に変換
する場合に使用される。
3) ラスター地図において地図座標を CRT 座標系で観測した場合、それらを UTM や平面直角座標系に直す際に
用いられる。
4) 地上測量において国家基準点を含んで任意座標系で測量しておき、これを測地座標系に変換する場合に使用
する。
5) 空中三角測量ではヘルマート変換、6 個のパラメータによるアフィン変換の他、非平面性（「ねじれ」）を考慮
して変換式に XY 項を加味した 8 個のアフィン変換式が用いられる。
（ただし、パラメータ数を増加させると見
かけ上残差は小さくなるが、空中三角測量の成果の精度が実際上向上するとは限らないので注意する。付録 E に
実例を記す。
）
6) 平面座標系における回転と移動の変換に一般に使用される変換式である。
7）ヘルマート変換又は Affine 変換の利用は、座標変換における幾何学的要素が現実のモデルである座標値の相
互関係に正しくフィットしているかの判定などにより決定すべき問題である。座標軸の回転や移動量などは定量
的なものであるが、観測値にそのような定量的なもののみが含まれているかどうかの判断は非常にむずかしい。
1.2
ヘルマート変換式 II
上のヘルマート変換式において座標軸の回転（θ）、座標原点の移動 ( xo , yo ) および座標(x,y)の縮率λを考慮す
れば、式(1-1)は次式で書ける。
(1-8)
�
x = λ ∙ X ∙ cosθ − λ ∙ Y ∙ sinθ + xo
y = λ ∙ X ∙ sinθ + λ ∙ Y ∙ cosθ + yo
そこで、上の式において
(1-9)
(1-10)
x
cosθ
�y � = λ �
sinθ
a = λ ∙ cosθ
b = λ ∙ sinθ
�
c = xo
d = yo
xo
−sinθ X
� � � + �y �
o
cosθ Y
|4
とおけば、式(1-8)は次のように簡単に表せる。
(1-11)
x=a∙X−b∙Y+c
�
y=b∙X+a∙Y+d
(1-12)
x
X
�y� = �
Y
又は
−Y
X
a
1 0 b
�� �
0 1 c
d
この式は、上に示した式と同様、二次元のヘルマート変換式である。
これら 4 個のパラメータを最小二乗法で求めると，それらの 4 個の要素は式(1-10)から求めることができる。
(1-13)
b
⎧ tanθ = a
⎪
λ = √a2 + b 2
⎨ xo = c
⎪
⎩ yo = d
[注]通常、二次元ヘルマート（Helmert）変換式といえば式(1-1)～式(1-12)をいう。また、ヘルマート変換式の
別名は数学では等角写像変換、またはコンフォーマル（共形）変換(conformal transformation)ともいう。さら
に、変換前後の座標軸が常に直角であることから「直交座標変換」(orthogonal transformation)（Decartes 座標
変換）とも呼ばれている。その座標系(x、y)と（X,Y）における変換では、直交性が保たれるところに特徴があ
る。
ただし、係数 a,b,c,d は線形的に解けるが、各要素は線形的に（つまり厳密には）正しく解けていないことに
留意する。
2.アフィン変換
2.1
アフィン変換式 I
通常、二次元のアフィン変換式(affine transformation)は、次のように書く。
(2-1)
�
x = a1 X + b1 Y + c1
y = a2 X + b2 Y + c2
しかしながら、それら 6 個のパラメータ（変換式の係数）の幾何学的意味が明らかでないので、この意味を次
のように考えてみよう。
式(1-8)におけるヘルマート変換式において測定座標 X、Y 軸の縮率 ( λ x , λ y ) と、x 軸の傾き ¸ x と y 軸の傾き ¸ y
を与えると、それらは次式で表される。
(2-2)
�
x = λx ∙ X ∙ cosθx − λy ∙ Y ∙ sinθy + xo
y = λx ∙ X ∙ sinθx + λy ∙ Y ∙ cosθy + yo
そこで、上の式において
(2-3)
λx cosθx
x
�y� = �
λx sinθx
λy sinθy X
xo
� � � + �y �
λy cosθy Y
o
|5
そこで、
y
Y
y'
Θy
yo
Ysin
Xcos
Θy
Θx
P
Y
xo
X
Θx
Xsin
Θx
Ycos
Θy
X
x'
x
図 3 アフィン変換
(2-4)
a1 = λx cosθx
⎧ b = −λ sinθ
y
y
⎪ 1
c1 = xo
⎨ a2 = λx sinθx
⎪ b2 = λy cosθy
⎩ c2 = yo
とおけば、式(2-3)は次のように表せる。
(2-5)
�
x = a1 ∙ X + b1 ∙ Y + c1
y = a2 ∙ X + b2 ∙ Y + c2
したがって、最初に示した式 2-1 を導いた。通常、二次元のアフィン変換式と呼ばれている。これら 6 個の変
換係数（パラメータ）を最小二乗法で求めると，それらの６個の要素は次式で求めることができる。
a
(2-6)
2.2
tanθx = 2
a1
⎧
⎪ tanθ = − b1
y
b2
⎪
2
λx = �a1 + a22
⎨
2
2
⎪λy = �b1 + b2
xo = c1
⎪
yo = c2
⎩
アフィン変換式 I の近似化
測定座標系 XY において X,Y 軸は微小角(θx , θy )とすると、cosθx = cosθy = 1, λx sinθx = θx , λy sinθy = θy と置
くことができ、ＸY とｘｙとの関係は次のような近似式で表すことができる（図 3 参照）。
|6
(2-7)
�
x = λx ∙ X ∙ cosθx − λy ∙ Y ∙ sinθy + xo
y = λx ∙ X ∙ sinθx + λy ∙ Y ∙ cosθy + yo
x = λx ∙ X − θy ∙ Y + xo
�
y = θx ∙ X + λy ∙ Y + yo
x
X
�y� = �
0
(2-9)
0
Y
2.3 アフィン変換式 II の意味
y
λx
⎡λ ⎤
⎢ y⎥
0 −Y 1 0 ⎢θx ⎥
�
X 0 0 1 ⎢θy ⎥
⎢xo ⎥
⎣yo ⎦
Y'
Y
X
X-Ysinβ
cos
θ
Ysinβ
y'
Yco
sβ
X'(X)
β)
cos
θ
Yco
sβ
sin
θ
β)
sin
(XYsi
n
（X
-Ys
in
θ
yo
θ
xo
Y
β
Ycosβ
(2-8)
x'
x
測定座標系 XY において、Y’軸は Y 軸に対して 90°を超えて ² だけ傾いているものとすれば、(X,Y)系と(X’,Y’)
系との関係は
（2-10）
1 −sinβ X
�� �
�X′� = �
0 cosβ Y
Y′
で表され、これをさらにθだけ回転させると
x ′ = (X − Ysinβ)cosθ − Ycosβsinθ
又は
y ′ = (X − Ysinβ)sinθ + Ycosβcosθ
x ′ = Xcosθ − Y(sinβcosθ + cosβsinθ)
y ′ = Xsinθ + Y(cosβcosθ − sinβsinθ)
|7
x ′ = Xcosθ − Ysin(β + θ)
y ′ = Xsinθ + Ycos(β + θ)
となり、縮率を考慮して、原点の移動を加えると
x = λX Xcosθ − λY Ysin(β + θ) + xo
(2-11)
y = λX Xsinθ + λY Ycos(β + θ) + yo
で示される。ここで、
(2-12)
a1 = λx cosθ
⎧b = −λ sin (β + θ)
y
⎪ 1
c1 = xo
⎨ a2 = λx sinθ
⎪ b2 = λy cos (β + θ)
⎩ c2 = yo
とおくと、アフィン変換式
(2-5)
�
x = a1 ∙ X + b1 ∙ Y + c1
y = a2 ∙ X + b2 ∙ Y + c2
を得る。これを解けば、以下のように６個の要素が解ける。
a
(2-13)
tanθ = 2
a1
⎧
⎪tan (β + θ) = − b1
b2
⎪
2
2
λx = �a1 + a2
⎨
2
2
⎪ λy = �b1 + b2
xo = c1
⎪
yo = c2
⎩
したがって、数学的（幾何学的）には、式(2-6)と式(2-13)は定義が異なるだけで求めているものは同じ要素であ
ることが分かる。
[注]
ここでは、平面の変換式であるヘルマート変換式とアフィン変換式の係数の意味を明らかにした。デジタル画
像やラスター地図の CRT 画面上での計測時に考えられる誤差には（又はデジタイザやコンパレータでの計測も
ほぼ同様）、次のようなものが考えられる。
1）画像座標 xy 軸のスケールの縮率（通常は x 軸観測値を 1 とし、y 軸の縮率を考えればよい。しかし、今回
わかったことは、測定座標単位が画素でｍｍ単位に換算する場合にはｘｙ軸の縮率は両方与えなければならな
い。）
2）画像座標 xy 軸におけるｙ軸の傾き
3）その他（スキャナーのガタ、写真の伸縮における異方性誤差、デジタイザの非平面性など）
4）観測誤差（ラスター変換時の解像度も大きく影響する。
）上には、ヘルマート変換式とアフィン変換式の係
数の幾何学的意味を明らかにした。
3. コロケーション（ヘルマート変換式による 3 パラメータ変換式の誘導）
（概要）以前に用いたある一つの平面座標ＸＹと時間が経過した時点での同じ点の平面座標ｘｙがあり、前者の
座標ＸＹをすべてｘｙに変換したい場合がある。また、両座標系が任意座標系で、座標軸の回転と原点の移動量
のみが考えられる場合である。ここで、両座標系の原点の縮率は mo=1.0 とする。
|8
そのときには、４個のパラメータの Helmert 変換、又は６個のパラメータのアフィン変換により、殆ど正確
に変換が行えるものとみなしていたが、これが誤りであることが判明した。以下に 3 個の変数（3 個の要素）の
変換式を導くことにする。
たとえば、通常の Helmert 変換では、次の式で表される。
(3-1)
�
x = λ ∙ X ∙ cosθ − λ ∙ Y ∙ sinθ + xo
y = λ ∙ X ∙ sinθ + λ ∙ Y ∙ cosθ + yo
�
x = X ∙ cosθ − Y ∙ sinθ + xo
y = X ∙ sinθ + Y ∙ cosθ + yo
しかし、今回の設定はλ=1.0 なので、
(3-2)
この式(2)を通常の最小二乗法で解くには、cosθ = a, sinθ = b, xo = c, yo = dとおいて、4 個のパラメータから 3
個の要素を解くことになるが、これは正しくない。しかし今まで、変数は 3 個なので、3 個のパラメータ（又は
3 個の要素）を直接解くことを考えなければならない。そうしないと、上に mo=1.0,又はλ=1.0 と仮定して座標
を求めたことの設定に反することになる。
そのため、式(3-2)を次のような条件式にする。
(3-3)
�
fx = X ∙ cosθ − Y ∙ sinθ + xo − x = 0
fy = X ∙ sinθ + Y ∙ cosθ + yo − y = 0
回転角（θ）と平行移動量（xo , yo ）に関し、テーラー展開して誤差式を作るには
(3-4) �
(3-5)
∆x =
∆y =
�
∂fx
∂θ
∂fy
∂θ
∆θ +
∆θ +
∂fx
∂xo
∂fy
∂yo
∆xo
∆yo
vfx = −X ∙ sinθ ∙ ∆θ − Y ∙ cosθ ∙ ∆θ + ∆xo − fxo
vfy = X ∙ cosθ ∙ ∆θ − Y ∙ sinθ ∙ ∆θ + ∆yo − fyo
又は、これを行列に置き換えて、繰り返し計算により要素を解く。
(3-6)
または
(3-7)
(3-8)
ここで
vfx
∆x
fx
−(X ∙ sinθ + Y ∙ cosθ)
�vf � = �
� ∆θ + � o � − � o �
∆y
fyo
X ∙ cosθ − Y ∙ sinθ
y
o
∆θ
vfx
fx
−(X ∙ sinθ + Y ∙ cosθ) 1 0 ∆x
�vf � = �
� � o� − � o�
fyo
X ∙ cosθ − Y ∙ sinθ
0 1 ∆y
y
o
∆θ
vfx
fx
−e 1 0 ∆x
�vf � = �
� � o� − � o�
fyo
f 0 1
y
∆yo
e = −(X ∙ sinθ + Y ∙ cosθ) = −(X ∙ b + Y ∙ a)
f = X ∙ cosθ − Y ∙ sinθ = X ∙ a − Y ∙ b
fxo = (X ∙ cosθ − Y ∙ sinθ + xo )o − x oo = (X ∙ a − Y ∙ b + c) − x oo
fyo = (X ∙ sinθ + Y ∙ cosθ + yo )o − y oo = (X ∙ b + Y ∙ a + d) − y oo
(X ∙ cosθ − Y ∙ sinθ + xo )o ： X ∙ cosθ − Y ∙ sinθ + xo の近似値
(X ∙ sinθ + Y ∙ cosθ + yo )o ：X ∙ sinθ + Y ∙ cosθ + yo の近似値
x oo , y oo ：x,y の近似値
fxo , fyo ：fx , fy の近似値
|9
(3-9)正しい要素の繰り返し計算
θ = θo + ∆θ
xo = xooo + ∆xo
yo = yooo + ∆yo
[付録 A]
（ヘルマート変換式の代数的特性）
x − xo
cosθ −sinθ X
�� �
�y − y � ＝ �
o
sinθ cosθ Y
…(A-1)
または、式（A-1）を
x − xo
cosθ
℘ = �y − y � , ℝ = �
o
sinθ
X
−sinθ
� , ℧=� �
Y
cosθ
とおくと、ヘルマート変換式は次のように書ける。
℘=ℝ∙℧
（A-2）
ここで、℘は定数行列、ℝは回転行列および℧はベクトル行列と呼ばれる。
式（A-2）の回転行列ℝは次のような性質がある。まず、式(A-2)の逆式を書くと次のようになる。
℧ = ℝ−1 ∙ ℘
（A-3）
ここで、ℝの逆行列ℝ−1 は次式で解ける。
ℝ−1 =
1
|ℝ|
adjℝ
（A-4）
ここで、|ℝ|は行列ℝ の行列式、adj ℝ は行列ℝ の随伴行列である。
|ℝ| = �cosθ
sinθ
ℝT = �
−sinθ
� = cos 2 θ + sin2 θ = 1
cosθ
cosθ sinθ
�
−sinθ cosθ
adjℝ = �
cosθ sinθ
�
−sinθ cosθ
∴ ℝ−1 = �
したがって、行列ℝ は次のような性質があることが証明された。
ℝT = ℝ−1
cosθ sinθ
�
−sinθ cosθ
（A-5）
また、ℝ は行列乗法において可換なので、ℝT ℝ = ℝℝT = I（単位行列）という特別な関係が成り立つ。
cosθ
ℝT ℝ = �
sinθ
1
−sinθ cosθ sinθ
��
�=�
0
cosθ −sinθ cosθ
0
�
1
(A-6)
このような回転行列を直交行列といい、この性質を直交性 orthogonality という。
式 A-1 の逆式℧ = ℝ−1 ∙ ℘は、次式で表される。
[付録 B]
x − xo
X
cosθ sinθ x − xo
� �y − y � = ℝT �y − y �
� �＝�
Y
o
o
−sinθ cosθ
（最小二乗法によるパラメータの決定）
上に述べた式
（A-7）
| 10
または
x
X
�y� = �
Y
a
−Y 1 0 b
�� �
X 0 1 c
d
𝐲=𝚨∙𝑿
（B-1）
（B-2）
と書く式を誤差論では観測方程式 observation equation と呼ぶ。ｎ個の対の（x,y）と（X,Y）とが与えられると
き、この方程式の数は２ｎ個となる。つまり、パラメータ（未知数）の数が４個であるのに対し、方程式が２ｎ
個となることから過剰観測(redundant observation)と呼ばれる。通常の代数ではその２ｎ個の方程式から４個の
式を適当に選んで解いても解は得られる。しかしながら、
（x,y）は観測した値であり、これらには微小であるが
誤差を含んでいる。つまり、式 B-1、B-2 には accidental 誤差が含まれているから２ｎ個の式は独立 independent
である(行列のランクが２ｎを表す)。2n 個の方程式から 4 個ずつ選んで解いたいくつかのパラメータ解をみると、
すべての解は同一にはならず、少しずつ違った値になる。したがって、2n 個の方程式から合理的に 4 個の式を
作成して 4 個のパラメータを解く方法を最小二乗法(Least Squares Method)という。
式 B-2 の 2n 個の観測方程式𝐲 = 𝚨 ∙ 𝑿 から次のように正規方程式(normal equations)（B-3）を計算する。
または
または
𝚨T 𝐲 = 𝚨𝑻 𝜜 ∙ 𝑿
∑(X 2 + Y 2 )
⎡
0
⎢
∑
X
⎢
⎣
∑Y
ℕ∙𝑿=𝓕
（B-3）
0
+ Y2)
−∑Y
−∑X
∑(X 2
（B-5）
∑(Xx + Yy)
∑X ∑Y
a
⎤
⎡
⎤
∑(Yx − Xy)
− ∑ Y − ∑ X⎥ b
⎢
⎥
� �=
∑x
n 0
⎢
⎥
⎥ c
∑y
⎣
⎦
0 n
⎦ d
（B-4）
と書く。これが正規方程式である。
ここで、
ℕ = 𝜜𝑻 𝑾𝜜
,
𝓕 = 𝜜𝐓 𝐖𝐲
ただし、式（B-1）は観測値（x、y）は等精度（W=I）とおいた場合である。不等精度の場合、𝜜T 𝐖𝐲 = 𝜜𝑻 𝑾𝜜 ∙ 𝑿によ
って正規方程式をたてる。ここで、W は重量行列であり、W はその対角要素を観測値の分散の逆数にすればよ
い。
通常、観測値は等精度とし、W=I（単位行列）として、式（B-3）により計算することが多い。
[付録Ｃ]
[式(1-2)の逆式]
[式(1-9)の逆式]
𝐗
cosθ sinθ x − xo
��
� �=�
�
𝐘
−sinθ cosθ y − yo
[式(2-2)の逆式]
1 cosθ
𝐗
sinθ x − xo
� � = λ�
��
�
𝐘
−sinθ cosθ y − yo
(C-1）
(C-2)
| 11
λy cosθ λy sinθ x − xo
1
𝐗
��
�
� �=λ λ �
𝐘
x y −λx sinθ
λx cosθ y − yo
[付録 D]
(C-3)
（8 パラメータのアフィン変換）
(D-1)
（二次多項式）
�
x = a1 X + b1 Y + c1 XY + d1
y = a2 X + b2 Y + c2 XY + d2
(D-2)
�
x = a1 X + b1 Y + c1 + d1 X 2 + e1 Y 2 + f1 XY
y = a2 X + b2 Y + c2 + d2 X 2 + e2 Y 2 + f2 XY
リモートセンシング画像の幾何補正では、式(１-11)のほか、式(D-2)が採用される。たとえば、ラスター地図
をジオコード化したいとき、その画像を回転と原点の移動により偏位修正（リサンプリング）する場合、式(1-11)
を利用する。
（縦横の縮率は変わっても良い）また、LANDSAT 画像において簡易補正（バルク補正）がなされ
た画像を地図の座標（UTM）に変換するには、式 D-2 が使用される。（下図参照）
グリニジL=0度（31系）
Lo=135度（53系）
207ｋｍ
BL→UTM
207ｋｍ
精度1/2500（縮率4/10000）
図 4 UTM への変換
(2y)2
m = mo [1 +
8R2 m2o
]
m = mo [1 +
2R2 m2o
]
y2
(D-3)
投影原点の縮率: mo = 0.9996 , m = 1.0004
y = √2Rmo �
m
mo
− 1 = 207km
（D-4）
距離 160ｋｍで 64ｍの違いであるから、LANDSAT-TM の 2 画素の誤差になる。したがって、東西 200ｋｍ以
内で平面とみなして測量が可能になる。
〔付録 E〕最小二乗法による幾何学的要素の算出
| 12
上では、Helmert 変換式における 4 個の係数や Affine 変換式における 6 個の係数を最小二乗法で解いたが、
これから求めた幾何学的意味を表す要素については厳密ではない。たとえば Helmert 変換式の要素である、縮
率、回転角、移動量を厳密に求めたければ、それらの要素によって表された関数式を線形化（テーラー展開）し
なければならない。線形化された誤差方程式を用いて正規方程式を作成し、これにより微小変量パラメータを計
算し、逐次近似により最確パラメータ（幾何的要素）を求めていく。
（1） Helmert 変換における幾何要素の誤差方程式
式 1-8 から、
（E-1）
�
f1 = λ ∙ X ∙ cosθ − λ ∙ Y ∙ sinθ + xo − x = 0
f2 = λ ∙ X ∙ sinθ + λ ∙ Y ∙ cosθ + yo − y = 0
そして、式（E-1）の線形化（テーター展開）は次式で表されます。
(E-2)
vx =
�
vy =
∂f1
∂λ
∂f2
∂λ
Δλ +
Δλ +
∂f1
∂θ
∂f2
∂θ
Δθ +
Δθ +
∂f1
∂xo
∂f2
∂xo
Δxo +
Δxo +
∂f1
∂yo
∂f2
∂yo
Δyo − f1o
Δyo − f2o
ここで、Δλ, Δθ, Δxo , Δyo は微小変動量を表す変数であり、f1o , f2o の近似値はそれぞれの各要素の近似値を用い次
式より計算します。
f o = λo X o cosθo − λo Y o sinθo + xoo − x o
� 1o
f2 = λo X o sinθo + λo Y o cosθo + yoo − y o
式 E-2 の微係数を計算すると次式になる。これが観測方程式、又は補正（残差）方程式になる。
�
（E-3）
vx = (Xcosθ − Ysinθ)∆λ − λ(Xsinθ + Ycosθ)Δθ + Δxo − f1o
vy = (Xsinθ + Ycosθ)∆λ + λ(Xcosθ − Ysinθ)Δθ + Δyo − f2o
これより正規方程式をつくりパラメータ(Δλ, Δθ, Δxo , Δyo )を解く。なお、式 E-3 における意味は、f1o , f2o に対応
する(Δλ, Δθ, Δxo , Δyo )の関数が表されていて、これを解けば最確パラメータが逐次的に求められるのである。
パラメータの最確値は、次式で計算します。
λo = λoo + Δλ
θ = θo + Δθ
� oo ooo
xo = xo + Δxo
yoo = yooo + Δyo
(E-4)
各要素値は式 E-3 などを用いて繰り返し計算で算出する。
ここで、式 E-4 の左辺はⅰ番目の値、右辺第 1 項目は i-1 番目の値を示す。
また、式（Ｅ-3）の係数は最初未知量なので、近似的 Helmert 変換で求めたもので近似する。
（2） Affine 変換式において 5 要素とする場合の幾何要素の誤差方程式
式(2-2)より
(E-5)
�
x = λx ∙ X ∙ cosθ − λy ∙ Y ∙ sinθ + xo
y = λx ∙ X ∙ sinθ + λy ∙ Y ∙ cosθ + yo
これを線形化すると
vx = (Xcosθ)Δλx − [Ysin(β + θ)]Δλy − �λx Xsinθ + λy Ycos(β + θ)�Δθ − λy cos(β + θ) Δβ + Δxo − f1o
vy = (Xsinθ)Δλx + [Ycos(β + θ)]Δλy + �λx Xcosθ − λy Ysin(β + θ)�Δθ − λy sin(β + θ) Δβ + Δyo − f2o
（E-6）
を得る。この式を用いて上と同様に斬近法により 5 つの要素を求めていく。
ここで、
| 13
f1o = λox ∙ X o ∙ cosθo − λoy ∙ Y o ∙ sinθo + xoo − x o
�
f2o = λox ∙ X o ∙ sinθo + λoy ∙ Y o ∙ cosθo + yoo − y o
（3） Affine 変換式において 6 要素とする場合の幾何要素の誤差方程式
式(2-2)より
(E-7)
�
x = λx ∙ X ∙ cosθx − λy ∙ Y ∙ sinθy + xo
y = λx ∙ X ∙ sinθx + λy ∙ Y ∙ cosθy + yo
これを線形化すると
vx = (Xcosθx )Δλx + �−Ysinθy �Δλy + (−λx Xsinθx )∆θx + (−λy Ycosθy )Δθy + Δxo − f1o
vy = (Xsinθx )Δλx + (Ycosθy )Δλy + (λx Xcosθx )Δθx + (−λy Ysinθy )∆θy + Δyo − f2o
（E-8）
を得る。この式を用いて斬近法により 6 つの要素を求めていく。
ここで、
f1o = λox ∙ X o ∙ cosθxo − λoy ∙ Y o ∙ sinθyo + xoo − x o
�
f2o = λox ∙ X o ∙ sinθxo + λoy ∙ Y o ∙ cosθyo + yoo − y o
(4)
Affine 変換係数の最小二乗法
Affine 変換係数は式(2-1)を観測方程式として解いても結果的には同じであるが、次のようにしてパラメータが
導かれる。
(2-1)
�
x = a1 X + b1 Y + c1
y = a2 X + b2 Y + c2
より、誤差方程式はテーラー展開を行い、次のように表される。
(E-9)
ここで、
vx = X o Δa1 + Y o Δb1 + Δc1 − f1o
�
vy = X o Δa2 + Y o Δb2 + Δc2 − f2o
f o = a1o ∙ X o + b1o ∙ Y o + c1o − f1o
� 1o
f2 = ao2 ∙ X o + bo2 ∙ Y o + c2o − f2o
（付録Ｅでのまとめ）
1） Helmert や Affine 変換での誤解は、係数と幾何学的要素を同じものと考えるところにある。
2） 係数と要素は、これらを最小二乗法で求める際、基本的にまったく異なる変量と考えなければならない。
3） 係数を解く場合、Helmert では４個（３個の場合もある）、Affine では６個の係数であり、最小二乗法適用
のとき、それらの式は線形なので線形化の必要もなく、繰り返し計算も不要である。誤差式にして繰り返し
計算しても、最初の式を用い直接解を計算しても結果はほとんど変わらない。
4） 幾何学的要素解は係数解の場合とは、理論的には異なる、という前提で解析を行わなければならない。した
がって、Helmert や Affine における式で幾何学的要素をまず明確にし、その変量に関し線形化を行わなけれ
ばならない。これが幾何学的要素を最小二乗法で求める際の基本的考え方である。内田論文（内田賢二、長
岡京条坊復元のための平均計算、長岡京跡発掘調査所ニュース第３１号、１９８４年 3 月１日、pp.2‐8）で
| 14
は、すでにこの考えで、Helmert 変換における幾何学的要素の線形化がなされ、回転角（条坊の振れ）、縮率
（尺の定義につながる）
、原点移動量が正しく求められている。さらに、条坊直角座標と平面直角座標との関
係を Helmert 変換にあてはめ、それらの幾何学的要素を解き、長岡京の条坊の位置、すなわち長岡京の方位
を分析している。
5） 今回、その条坊座標の計算を行うためプログラムを改良し作成し直した。その際の留意点は、写真座標の解
析の際には数学座標においてｘからｙ軸への回転κは+なのに対し、測地座標の変換において西偏は-、つま
り東偏+とする。Helmert、Affine および最小二乗法で解いた、条坊座標の振れ、すなわち軸の回転角などを
求めた結果は、表 E-1～表 E-4 にまとめて示してある。
表 E-1
平城京
θ（度）
β（度）
Sx（丈）
Sy（丈）
Xo（ｍ）
Yo（ｍ）
Rmse
表 E-2
2.9550
2.9550(0.4mm)
-18585.166
-146019.056
1.576m
-18585.166(0.2)
-146019.056(0.2)
1.592
LS(rms)
-11’53.4”(58.0”)
-4’16.3”(72.0”)
2.9545(0.7mm)
2.9540(0.5mm)
-18585.431(0.2)
-146018.992(0.2)
1.698
長岡京
Helmert（22 点使用）
LS
係数
-8’57.5”(10.5”)
西偏 8’57.5”
θ（度）
β（度）
Sx（ｍ）
Sy（ｍ）
Xo（ｍ）
Yo（ｍ）
Rmse
表 E-3
Affine
係数
-11’53.4”
-4’16.3”
2.9545
2.9540
-18585.431
-146018.992
1.698
Helmert（36 点使用）
LS(rms)
係数
西偏-14’54.5” -14’54.5”(3.2”)
2.9760
2.9760(1.2mm)
-26842.104
-117535.035
2.934m
-26842.104(0.7)
-117535.035(0.7)
2.997
Affine
係数
西偏-2’52.0”
西偏-4’39.3”
2.9709
2.9818
-26840.972
-117532.351
2.379
LS
-1’55.1”(117.3”)
-2’07.0”(155.6”)
2.9655(1.7mm)
2.9708(1.5mm)
-26840.162(0.3)
-117533.779(0.3)
2.379
Affine
係数
-14’22.7”
西偏-40.1”
2.9843
2.9848
-23233.749
-110862.408
0.861
LS
-14’22.7”(25.6”)
-40.1”(33.5”)
2.9843(0.4m)
2.9848(0.3m)
-23233.749(0.2)
-110862.408(0.2)
0.861
平安京
Helmert(全点)
LS
係数
西偏-14’43.7” -14’43.7”(1.4”)
θ（度）
β（度）
Sx（ｍ）
Sy（ｍ）
Xo（ｍ）
Yo（ｍ）
Rmse
2.9847
2.9847(0.2m)
-23233.786
-110862.440
0.692m
-23233.786(0.1)
-110862.440(0.1)
0.699
（上記は各変換で LS により計算した結果を、また（
）内数値は正規方程式の逆行列から求めた幾何学的要素
の平均 2 乗誤差を示す。ただし、係数の rms は示していない。）
[F]相関係数について
平城計算における幾何学的要素の相関係数の検討
（Helmert）
λ
θ
xo
yo
| 15
λ
θ
xo
yo
1.0
-0.0
-0.1
0.4
-0.0
1.0
-0.4
0.1
-0.1
-0.4
1.0
0.0
0.4
-0.1
0.0
1.0
1）Helmert ではλ-yo 間に 0.5 の相関、θ-xo 間の-0.5 の強い相関が認められる。
(Affine)
xo
yo
λx
λy
θ
β
1.0
-0.0
-0.0
-0.2
0.1
-0.0
λx
-0.0
1.0
-0.3 0.3
-0.0
0.5
λy
0.0
0.3
1.0
-0.8
0.0
-0.1
θ
0.2
0.3
-0.8
1.00
0.3
0.1
β
xo
0.1
-0.0
0.0
0.3
1.0
0.0
yo
0.0
0.5
-0.1
0.1
0.0
1.0
2）一方、Affine では、λy-yo 間で 0.5 であり、これは Helmert の値に対応する。θ-β間では-0.8 で非常に強
い。β-xo 間では-0.3 の相関であり、これは Helmert のθ-xo に対応する。
（本検討資料は未発表。1995 年 6 月作成。2000 年 11 月 25 日修正。2001 年 2 月 14 日、2002 年 12 月 17 日修
正。2004 年 11 月 30 日修正。2010 年 3 月 23 日小林和夫）(2013 年 2 月 10 日修正追加)