第2章 位置の計算

第 2 章 位置の計算

測量において地球上の水平位置は経度と緯度、平面座標、又は極座標の距離と方向角(球面距離と方位角)で表 される。以下では、経緯度と平面座標の関係等について計算方法について例題（プログラム）を用いて示すこと にしよう。以下の計算での数値は、Bessel 楕円体から GRS80 楕円体に変更してあるので注意のこと。 2.1 経緯度および方位角 図 2.1.1 に示すように与点Ａの緯度と経度1 _{(φ, λ)は既知とする。いま、点Ａと新点Ｂとの球面距離（S）と点} Ａにおける点Ｂの方位角 2 _{（α）を測定したとすれば、点Ｂの緯度と経度(φ′, λ′)と点ＢにおけるＡの方位角（反} 方位角）は、次に示す式で求められる。 φ′_{= φ +} W3 1−e2u − 3W4_e2

2(1−e2₎2u2sin φ cos φ −

W4_{tan φ}

2(1−e2₎v2

− W5e2

2(1−e2₎3[1 − 2 sin2φ − e2sin2φ(5 − 6 sin2φ)]u3

− W5

6(1−e2₎2[1 + 3 tan2φ − e2tan2φ(13 − 10 sin2φ)]uv2 …(2.1.1)

�λ′− λ�cosφ = Wv + W2_{uv +}W3 3 �1 + 3 tan2φ + e2 1−e2cos2φ� u2v − 1 3(W3tan2φ)v3 …(2.1.2) 赤道 B α’ S A 9 0 -φ Nｐ λ’ λ φ φ´=φ+Δφ-δφ 90_－ φ ’ α dλ 図 2.1.1 緯度、経度と方位角 α′_{= α + (Wtanφ)v +}W2 2 (1 + 2 tan2φ + e2

1−e2cos2φ)uv

+1 6W3tanφ �5 + 6 tan2φ + e2 1−e2cos2φ − 4e2 (1−e2₎2cos4φ� u2v

−1₆W3tanφ �1 + 2 tan2φ +_1−ee2₂cos2φ� v3+ 180o …(2.1.3)

1_{latitude(緯度)、longitude(経度)：楕円体上の地点の位置}

ここで、 _{u =}S acosα , v = S asinα W = �1 − e2sin2φ e = �a2−b2 a2 [例題２．１] 新点Ｂの緯度経度(φ_B, λB)および反方位角(α′)をプログラム（ＳＯＫ２１．BAS）を用いて求めてみよう。 ただし、与点Ａの緯度経度(φ_A, λA)、AB 間の球面距離(S)および点Ａにおける点Ｂの方位角(αは次の値とする。 φA =34o41’25”.0000 , λA_＝135o_{30’19”.0000} S=14,999.930m, α=134o_29’29”.7 （解）緯度 B’＝ 34º35’ 43”.6660 経度 L’＝ 135º37’ 18”.9142 方位角 α’＝ 314º33’ 28”.404 2.2 2 点の経緯度より球面距離および方位角の計算 2 点の緯度経度が与えられている場合、前節の式 2.1.1～式 2.1.3 を用い繰り返し反復計算によって(u, v)を求 めると、これらの２点間の球面距離（Ｓ）および方位角 ( )α は次式から解ける。 _{S = a√u}2_{+ v}2 α = tan−1_�v u� …(2.2.1) （第１近似） u ≈ f1 v ≈ g1 （第２近似） u2≈ f1+ f2u12+ f3v12 v2≈ g1− g2u2v1 （第３近似） u ≈ f1+ f2u22+ f3v22+ f4u23+ f5u2v22 v ≈ g1− g2uv2− g3u2v2+ g4v23 ...(2.2.2) ここに、 f1= (φ′ − φ)1−e 2 W3 f2= 3We 2 2(1−e2₎sinφcosφ f3=Wtanφ₂ f4= W 2_e2

f5= W

2

6(1−e2₎[1 + 3tan2φ − e2tan2φ(13 − 10sin2φ)]

g1= (λ′ − λ)cosφ_W g2= Wtanφ g3=W 2 3 (1 + 3tan2φ + e2 1−e2cos2φ) g4=W 2 3 tan2φ …(2.2.3) [例題２．２] 2 点Ａ，Ｂの緯度経度が与えられているものとする。これらの値からプログラム（ＳＯＫ２２．BAS）によっ て球面距離(S)と点Ａにおける点Ｂの方位角(α)を求めてみよう。 ただし、点Ａの緯度経度(φ_A, λA)および点Ｂの緯度経度(φB, λB)は、次のとおりとする。 φA=34º41’25.”0000 λA=135º30’19.”0000 φB=34º35’43.”6660 λB=135º37’18.”9142 （解）方位 α’=134º29’ 29”.6880 球面距離S= 14,999.931m 2.3 ２点の平面座標による球面距離および方位角の計算

点A,Bの平面直角座標 1 _{はそれぞれA(xA, yA), B(xB, yA)とすれば、AB間の平面距離（s）および点Aにおける点} Bの方向角2 _{（t）は次式から計算できる。} (平面距離) _{s = �(xB− xA)}2_{+ (yB− yA)}2 _…(2.3.1) (方向角) t = tan−1_�yB−yA xB−xA� …(2.3.2) （球面距離） ∴ s ÷ (s S� ) …(2.3.3) （縮尺係数） 縮尺係数は、その点の縮率であり、日本の場合楕円体に横円筒をかぶせたときの ds/dS は次式で求められる。 m = mo�1 + y 2 2M1N1mo2+ y4 24M12N12mo4� ...(2.3.4) ここで、 平面直角座標系の場合、原点の縮尺係数_{m0= 0.9999、UTM の場合m0= 0.9996とする。}

1_{plane rectangular coordinates：地球楕円体に横円筒をかぶせて楕円体の中心を投影原点とし、地点を円筒面} に投影し、この円筒面を切り開いて平面にした座標。

（距離の補正）s/S は平面距離 s を求めたり、球面距離 S を求めたりする係数である。 S s = 1 mo�1 − 1 6rm2mo2(y1 2_{+ y} 1y2+ y22) + η 2_t 6r3_mo3(x2− x1)(y22− y12)� =_m1 o�1 − 1 8rm2mo2(y1+ y2) 2₋ 1 24rm2mo2(y2− y1) 2₊ η2t 6r3_mo3(x2− x1)(y22− y12)� …(2.3.5) (方位角) α = t + γ + T1− t1 α = t + γ + T2− t2 …(2.3.6) γ：P における子午線収差、T1− t1は方向補正 （方向補正） T1− t1=x2_6r− x1 122 (2y1+ y2) − η2_t 9r3(x2− x1)2(y2− y1) +η 2_t 6r3(y2− y1)(3y12+ 2y1y2+ y22) T2− t2= −x2_6r− x1 21 2 (y1+ 2y2) +η 2_t 9r3(x2− x1)2(y2− y1) −η 2_t 6r3(y2− y1)(y12+ 2y1y2+ 3y22) …(2.3.7) ここで、 _{R = √MN} M =a(1−e_W32) N =_Wa W = �1 − e2_sin2_φ e = �a2_a−b22 [例題２．３] 平面直角座標系における 2 点Ａ，Ｂの値が与えられているものとする。これらの値を用いてプログラム（ＳＯ Ｋ２３．BAS）により、平面距離(s)、球面距離(S)、方向角(t) および点 A における点 B の方位角( )α を求めて みよう。 ただし、点 A(-144,654.741m,107,365.335m),緯度 φA=34º41’25.”0000 点 B(-155,042.218m,118,187.713m) および点Ａにおける子午線収差はγ＝+0º43’54”.209 である。 （解） 平面距離 s= 15,000.785 方向角 t= 133°49’ 31.2016” 球面距離 S= 14,999.931ｍ 方位角 α= 133°49’ 28.272”

2.4 経緯度から平面座標の計算 地球の基準楕円体1 _{に横円筒をかぶせて各地点をその円筒上に投影しておき、その円筒を切り開けば、楕円体} から平面への写像が直接的に得ることができる。これは横メルカトール図法2_{と呼ばれている。測量で多く用い} られる横メルカトール図法としての平面直角座標系では、投影原点の縮率をmo=0.9999 として距離の精度が 1/10,000 以内に収まるように配慮されている。 地点の緯度・経度(φ, λ)から平面座標(x,y)および子午線収差(γ)を求める式は、次のとおりである。 _{x＝mo[∆φ +}∆λ2 2 Nsinφcosφ + Δλ4 24 Nsinφ cos3φ(5 − t2+ 9η2+ 4η4) +₇₂₀∆λ6Nsinφcos5_{φ(61 − 58t}2_{+ t}4_{+ 270η}2_{− 330t}2_η2_{)] …(2.4.1)} y = mo�ΔλNcosφ +Δλ 3 6 N cos3φ(1 − t2+ η2� + ∆λ5 120Ncos5φ(5 − 18t2+ t4+ 14η2− 58t2η2] …(2.4.2) γ = Δλsinφ +Δλ₃3sinφcos2_{φ(1 + 3η}2_{+ 2η}4_{) +}Δλ5 15 sinφcos4φ(2 − t2) …(2.4.3) 子 午 線 赤道 平行圏 原 子 午 線 Δλ=λ -λo Np λo λ φ φ γ ｙ ｘ δ φ Δ φ φ o P ｙ ｘ 図 2.4.1 緯度経度と平面座標 ここで、 _{Δφ = a(1 − e}2_{)[k1�φ − φo� −}1 2k2�sin2φ − 2φo� +14k3�sin4φ − 4φo� −

1₆k₄�sin6φ − 6φ_o� +1₈k₅�sin8φ − 8φ_o� −₁₀1k₆�sin10φ − 10φ_o� + ⋯ ] …(2.4.4)

k1=1.005052501813087 k2=0.005063108622224 k3=0.000010627590263 k4=0.000000020820379 1 _{reference ellipsoid：準拠楕円体}

k5=0.000000000039324 k6=0.000000000000071 Δλ＝λ-λo、λo＝座標系原点の経度 _{t=tanφ、η=e’cosφ} e’_{= �}a2−b2 b2 [例題２．４] 点Ａの緯度経度は(φ=34º41’25.”0000, λ=135º30’19.”0000)とする。この点の平面直角座標（x, y）の値および 子午線収差(γ)をプログラム（ＳＯＫ２４．BAS）を用いて求めてみよう。 ただし、平面座標系は第５系(φo=36º,λo=134º20’)とする。 （解） X=-144654.741 Y= 107365.335 γ =+ 0°40’1”.4312 2.5 平面座標から経緯度への変換 ある地点における平面直角座標(x,y)から緯度経度(φ, λ)および子午線収差(γ)を求めるには、次の式で行える。 _{φ = φ1−}(y m� o)2 2MN t + �y m� o� 4 24MN3 t(5 + 3t2+ η2− 9t2η2) − �y m� o� 6 720MN5t(61 + 90t2+ 45t4) ...(2.5.1) λ = λ_o+ y_m o � Ncosφ− �y m� o� 3 6N3_cosφ(1 + 2t2+ η2) + �y m� o� 5 120N5_cosφ(5 + 28t2+ 24t4) ...(2.5.2) γ = y mo � N t − �y m� o� 3 3N3 (1 + t2− η2) + �y m� o� 5 15N5 t(2 + 5t2+ 3t4) ...(2.5.3) ここで、

φ₁(i)= φ₁(i−1)−S(i−1)_M(i−1)−D

D = a(1 − e2_{) �k}

1φ_o−k_{2 sin2}2 φ_o+k_{4 sin4}3 φ_o−k_{6 sin6}4 φ_o+k_{8 sin8}5 φ_o−k_{10 sin10}6 φ_o+ ⋯ � +_mx′ o S(i)_{= a(1 − e}2_)(k 1φ1− k2sin2φ1(i) 2 + k3sin4φ1(i) 4 − k4sin6φ1(i) 6 + k5sin8φ1(i) 8 − k6sin10φ1(i) 10 + ⋯ ) M(i)₌ a(1−e2) (1−e2_sin2_φ 1 (i)_)3/2 [例題２．５] 点Ａの平面直角座標xyの値は第５系(-144,654.741m,107,365.335m)とする。この点の緯度経度(φ,λ)及び子午 線収差(γ)をプログラム（ＳＯＫ２５．BAS）を用いて求めてみよう。 （解） 緯度Ｂ= 34º41’ 25”.5018 経度Ｌ= 135º30’ 18”.5040

子午線収差（γ）=+ 0°40’1.1573”

2.6 UTM 図法

UTM 投影(Universal Transverse Mercator Projection)による平面座標は世界的に利用されている。わが国で は縮尺1/1 万地形図から 1/25000 基本図，1/5 万編集図における投影に採用されている｡ UTM 座標系は赤道を 6°ずつ東回りに 60 個のゾーンに分割して横メルカトールで投影する図法である。つま り、最初のゾーン1 は東周りに西経(-)180°～174°範囲（中央経線西経 177°）、ゾーン 2 は西経 174º～168º、… とし、ゾーン30 は西経 6°～0°（中央経線西経 3°）、ゾーン 31 は東経 0°～6°、ゾーン 60 は東経 174°～180°（中 央経線東経177°）である。因みに、大阪や沖ノ鳥島はゾーン 53 で、東経 132°～138°、中央経線東経 135°（明 石市）に該当する。（この計算はutm_zone.xlsm で作成している。）日本はゾーン 51～56（中央経線 123º（尖 閣）、129º（沖縄）、135º（大阪）、141º（東京）、147º、153º（北方四島））を採用する。各経度帯内では その中央経線と赤道との交点を座標原点として横メルカトル図法で投影する｡中央経線から東西方向に離れるに 従って距離ひずみが増大するので、平面距離s と球面距離 S の点の投影誤差 ds/ｄS の比率は±4/10000 で収まる ように考える。つまり、中央経線上の縮尺係数は0.9996 に設定するので、中央経線より東西に 180km 離れたと ころの縮率は1.0000、そして同様に 270km 離れた地点の縮率は 1.0004 となる。UTM 投影では中央経線より 270km までは、距離の相対精度が±4/10000 に収まるので、各６°幅の経度帯内で地球を平面とみなして投影する ことができる。（ただし、原点から3°離れると 334ｋｍなので 1.0004 を超える）ＵＴＭ座標においてｘ座標（Ｎ 軸）は中央経線（縦軸）に沿って北半球では0 より増加し、南半球では 10,000km（擬北距＝F.N.）より減少さ せる。ｙ座標（Ｅ軸）は赤道（横軸）を500km（擬東距＝F.E.）としてE 軸に沿い東側へは増加、西側へは減 少する。 Bessel 緯度経度（以下の Ki は GRS80 での値である。）とxy(ｘN,ｙE)座標との関係式は、以下のように表す ことができる。 2.6.1 緯度経度からＵＴＭ座標の変換 _{xN= F. N. +mo[Δφ +}Δλ2 2 Nsinφcosφ + Δλ4 24 Nsinφcos3φ(5 − t2+ 9η2+ 4η4) + Δλ₇₂₀6Nsinφcos5_{φ(61 − 58t}2_{+ t}4_{+ 270η}2_{− 330η}2_t2_{)] …（2.6.1）} yE= F. E. +mo[Δλ ∙ Ncosφ +Δλ 3 6 Ncos3φ(1 − t2+ η2) + Δλ₁₂₀5Ncos5_{φ(5 − 18t}2_{+ t}4_{+ 14η}2_{− 58t}2_η2_{)] …（2.6.2） γ = Δλsinφ +}Δλ3 3 sinφcos2φ�1 + 3η2+ 2η4� + Δλ5 15 sinφcos4φ(2 − t2) …(2.6.3) ここで、 Δφ = a(1 − e2_)[k 1�φ − φ_o� −1₂k2�sin2φ − 2φ_o� +1₄k3�sin4φ − 4φ_o� −

1₆k4�sin6φ − 6φ_o� +1₈k5�sin8φ − 8φ_o� −₁₀1k6�sin10φ − 10φ_o� + ⋯ ] …(2.4.4)

K2 =5.06310862222411D-3 K3 =1.06275902633386D-5 K4 = 2.0820378571650D-8 K5 = 3.9323713711242D-11 K6 = 7.1084534028592D-14 Δλ＝λ-λo、λo＝座標系原点の経度, t=tanφ、η=e’cosφ _{e′ = �}a2−b2 b2 mo=0.9996 a=6,378,137m 1/f=298.257222101 南北半球ではF.E.=500ｋｍ, 北半球では北距 F.N.=0ｍ、南半球 F.N.=10,000ｋｍ [例題２．６．１] ある三角点の GRS80 による世界測地系の緯度・経度は次のとおりとするとき、この点の UTM 座標および子 午線収差（角）をプログラム（ＳＯＫ２６１．ＴＸＴ）により求めてみよう。 φ＝35º39’29”.1572 λ=139º44’28”.8869 また、この点の UTM 座標の原点は次の値である。 φo=0º λo=141º (解) X= 3,946,757.290 Y= 386,070.956 子午線収差γ =- 0o _{44’ 01”.684} 2.6.2 ＵＴＭから緯度経度への変換 φ = φ₁− yE′ 2NMt + yE′4 24MN3t(5 + 3t2+ η2− 9t2η2) − yE′6 720MN5t(61 + 90t2+ 45t4) …（2.6.6） ...（2.2.1） λ = λo+ [_NcosφyE′ − yE′ 3 6N3_cosφ(1 + 2t2+ η2) + yE′5 120N5_cosφ(5 + 28t2+ 24t4)] …（2.6.7） γ =yE′ N t − yE′3 3N3t(1 + t2− η2) + yE ′5 15N5t(2 + t2+ 3t4) ...（2.6.8） ここで、 南北半球では y_E′ _{= (yE− 500,000[m])/mo}

mo = 0 9996.

ただし、γ は子午線収差角である。 φ₁(i)_{= φ1}(i−1)₋S(i−1)−D

M(i−1)

D = a(1 − e2_{) �k}

1φo−k_{2 sin2}2 φo+k_{4 sin4}3 φo−k_{6 sin6}4 φo+k_{8 sin8}5 φo−k_{10 sin10}6 φo+ ⋯ � +_mx′ o S(i)_{= a(1 − e}2_)(k 1φ1− k2sin2φ1(i) 2 + k3sin4φ1(i) 4 − k4sin6φ1(i) 6 + k5sin8φ1(i) 8 − k6sin10φ1(i) 10 + ⋯ ) M(i)₌ a(1−e2) (1−e2_sin2_φ 1 (i)_)3/2 …（2.2.4） ここで、 北半球では _{x′ = xN− F. N. = xN} 南半球では x′ = xS− F. N. = xS− 10,000,000[m] …(2.2.5) である。 [例題２．６．２] ここでは、例題２.６.１の逆計算をプログラム（ＳＯＫ２６２．ＴＸＴ）を用いて行ってみよう。ただし、緯 度経度を求めるために用いるUTM 座標は次のとおりである。 xN=3,946,757.290m yN=386,070.956m (解)

Input data(UTM x,y) Xn= 3,946,757.290ｍ Ye= 386,070.956ｍ LatitudeＢ=35°39’ 29 .1572 ″ LongitudeＬ=139°44′28.8869″ Meridian Convergence γ= - 0°44′01.684″ ３．２ GPS 観測値の平面直角座標への計算 ＧＰＳ観測値であるＷＧＳ８４における緯度・経度・楕円体高 ( , , )ϕ λ h から_{Bessel の(φ, λ, h)に変換する式は} 以下のとおりである。 (1) GPS 観測値 参照点、または求点のDGPS観測から得られたWGS841 _{の経緯度および楕円体高は、次のとおりとする。} GPS = �φλ h�WGS84 ...(3.2.1)

(2) GPS 観測値から WGS84 のXYZへの変換 WGS84 楕円体における３次元直交座標XYZは、次式で算出する。 XWGS84= (N84+ h84)cosφ84cosλ84 YWGS84= (N84+ h84)cosφ84sinλ84 ZWGS84= [N84(1 − e842 ) + h84]sinφ84 � …(3.2.2) なお、以上で用いた要素は、次のとおりである。 卯酉線曲率半径N84=_Wa84₈₄ W84= �1 − e842 sin2φ84 扁平率の逆数1/f84=_a84a_−b84₈₄ 短半径ｂ₈₄= a84(1 − f84) 第一離心率e84= �a84 2 _−b 84 2 a842 ...(3.2.3) また、WGS84 の地球楕円体の値は、次のとおりである。 a84= 6,378,137m _f1 84= 298.257223563 b₈₄= 6,356,752.314m …(3.2.4) (3) WGS84 から ITRF89 系座標への変換 �XY Z�ITRF = RT_�X_Y Z�84 − �−0.0550.541 0.185 � …(3.2.5) (4)ITRF89 系から筑波座標系 92（Bessel）座標への変換 �XY Z�Tsukuba = RT_�X_Y Z�ITRF − �−147.531508.313 680.417 � …(3.2.6) (5) Bessel 経緯度への変換 �XY Z�Bessel = �XY Z�Tsukuba …(3.2.7)

Bessel 経度緯度および正規標高（Orthometric Height）は、これらのXYZ値を用いて、以下のように求める ことができる。

tanλB= Y/X ...(3.2.8)

（近似式）

tanφ_B=_√XZ+(e′)_2+Y22_−ebsin_2acosθ3θ

tanθ =_b√XZ∙a2_+Y2 ...(3.2.9)

tanφ_B=Z+NBeB2sinφB √X2_+Y2 HB=√X 2_+Y2 cosφ_B − NB hB= HB− HG ...(3.2.10) 上に示した Bessel 楕円体上の標高は、次式で求められる。 hB＝（水準標高Ｈ）＋（ジオイド高HG） ...(3.2.11) また、ジオイド高1 _{はジオイド図を参照して求められる。} ここで、Bessel 楕円体は、次のとおりの値である。 aB= 6,377,397.155m bB= 6,356,078.963m _f1 B= 299.152813 ...(3.2.12)

GRS80(geodetic reference system 1980)では

長半径 a= 6 378 137m 及び扁平率 f= 1/298.257 222 101、b＝a(1-f)= b = 6 356 752.314 140 356m

WGS84（GPS 衛星で使用） a=6378137m, f=1/298.257223563

ロシアのGLONASS の使用する楕円体（PZ90）

GLObal’naya NAvigationnaya Sputonikovaya Sistema=英語 Global Navigation Satellite System ロシア宇宙軍の運用するGPS 衛星が GLONASS であり、

GLONASS では座標系として PZ-90 という、北極点の位置を 1900 年から 1905 年までの測定位置の平

均としたものを採用している。これに対して GPS では、WGS84 という 1984 年時点での北極点の測

定位置を元にした測地系を採用している。

“PZ-90”(Parametry Zemli 1990, PE-90; Parameters of the Earth 1990)は、測地系の基準楕円体に

SGS-90 を用いた ECEF 座標系であり、原点を地球重心にとり、Z 軸が IERS の国際基準座標系の X3

軸を慣用極方向として、X 軸はグリニジ子午線と赤道面との交点と地球重心を結ぶ軸、Y 軸は赤道面上

で X 軸に直角な軸で、X-Y-Z は右手系を採用。PE-90 楕円体は長軸半径 a=6378136.0m、扁平率の逆

数 1/f=298.257839303 である。