1 第 章位置の計算 測量において地球上の水平位置は経度と緯度 平面座標 又は極座標の距離と方向角 ( 球面距離と方位角 ) で表 される 以下では 経緯度と平面座標の関係等について計算方法について例題 ( プログラム ) を用いて示すことにしよう 以下の計算での数値は Bessel 楕円体から GRS80 楕円体に変更してあるので注意のこと.1 経緯度および方位角図.1.1 に示すように与点 Aの緯度と経度 1 (φ, λ) は既知とする いま 点 Aと新点 Bとの球面距離 (S) と点 Aにおける点 Bの方位角 (α) を測定したとすれば 点 Bの緯度と経度 (φ, λ ) と点 BにおけるAの方位角 ( 反方位角 ) は 次に示す式で求められる φ = φ + W3 1 e u 3W4 e (1 e ) u sin φ cos φ W4 tan φ (1 e ) v W5 e (1 e ) 3 [1 sin φ e sin φ(5 6 sin φ)]u 3 W5 [1 + 3 6(1 e ) tan φ e tan φ(13 10 sin φ)]uv (.1.1) λ λ cosφ = Wv + W uv + W3 3 1 + 3 tan φ + e 1 e cos φ u v 1 3 (W3 tan φ)v 3 (.1.) Np dλ 90-φ 90-φ B α S A α φ =φ+δφ-δφ φ λ 赤道 λ 図.1.1 緯度 経度と方位角 α = α + (Wtanφ)v + W (1 + tan φ + e 1 e cos φ)uv + 1 6 W3 tanφ 5 + 6 tan φ + e 1 e cos φ (1 e ) cos4 φ u v 4e 1 6 W3 tanφ 1 + tan φ + 1 e cos φ v 3 + 180 o e (.1.3) 1latitude( 緯度 ) longitude( 経度 ): 楕円体上の地点の位置 azimuthal angle( 方位角 ): 楕円体上 ある位置の子午線の北を基準に他の位置まで右回りに測った角
u = S a cosα, v = S a sinα W = 1 e sin φ e = a b a [ 例題.1] 新点 Bの緯度経度 (φ B, λ B ) および反方位角 (α ) をプログラム (SOK1.BAS) を用いて求めてみよう ただし 与点 Aの緯度経度 (φ A, λ A ) AB 間の球面距離 (S) および点 Aにおける点 Bの方位角 (αは次の値とする φa =34 o 41 5.0000, λa =135 o 30 19.0000 S=14,999.930m, α=134 o 9 9.7 ( 解 ) 緯度 B = 34º35 43.6660 経度 L = 135º37 18.914 方位角 α = 314º33 8.404. 点の経緯度より球面距離および方位角の計算 点の緯度経度が与えられている場合 前節の式.1.1~ 式.1.3 を用い繰り返し反復計算によって (u, v) を求めると これらの 点間の球面距離 (S) および方位角 ( α) は次式から解ける S = a u + v α = tan 1 v u (..1) ( 第 1 近似 ) u f 1 v g 1 ( 第 近似 ) u f 1 + f u 1 + f 3 v 1 v g 1 g u v 1 ( 第 3 近似 ) u f 1 + f u + f 3 v + f 4 u 3 + f 5 u v v g 1 g uv g 3 u v + g 4 v 3...(..) ここに f 1 = (φ φ) 1 e W 3 f = 3We (1 e ) sinφcosφ f 3 = Wtanφ f 4 = W e (1 e ) [1 sin φ e sin φ(5 6sin φ)]
3 f 5 = W 6(1 e ) [1 + 3tan φ e tan φ(13 10sin φ)] g 1 = (λ λ) cosφ W g = Wtanφ g 3 = W 3 (1 + 3tan φ + e 1 e cos φ) g 4 = W 3 tan φ (..3) [ 例題.] 点 A,Bの緯度経度が与えられているものとする これらの値からプログラム (SOK.BAS) によって球面距離 (S) と点 Aにおける点 Bの方位角 (α) を求めてみよう ただし 点 Aの緯度経度 (φ A, λ A ) および点 Bの緯度経度 (φ B, λ B ) は 次のとおりとする φa=34º41 5. 0000 φb=34º35 43. 6660 λa=135º30 19. 0000 λb=135º37 18. 914 ( 解 ) 方位 α =134º9 9.6880 球面距離 S= 14,999.931m.3 点の平面座標による球面距離および方位角の計算点 A,Bの平面直角座標 1 はそれぞれA(x A, y A ), B(x B, y A ) とすれば AB 間の平面距離 (s) および点 Aにおける点 Bの方向角 (t) は次式から計算できる ( 平面距離 ) s = (x B x A ) + (y B y A ) (.3.1) ( 方向角 ) t = tan 1 y B y A x B x A ( 球面距離 ) s ( s S ) (.3.) (.3.3) ( 縮尺係数 ) 縮尺係数は その点の縮率であり 日本の場合楕円体に横円筒をかぶせたときの ds/ds は次式で求められる y m = m o 1 + M 1 N 1 m + y 4 o 4M 1 N 1 m 4 o...(.3.4) 平面直角座標系の場合 原点の縮尺係数 m 0 = 0.9999 UTM の場合 m 0 = 0.9996とする 1plane rectangular coordinates: 地球楕円体に横円筒をかぶせて楕円体の中心を投影原点とし 地点を円筒面 に投影し この円筒面を切り開いて平面にした座標 direction angle: この円筒面の真上である座標の北を基準として他の点まで右回りに測った角
4 ( 距離の補正 )s/s は平面距離 s を求めたり 球面距離 S を求めたりする係数である S s = 1 1 1 m o 6r m m (y 1 + y 1 y + y ) + o η t 6r 3 m o 3 (x x 1 )(y y 1 ) = 1 1 1 m o 8r m m (y 1 + y ) 1 o 4r m m (y y 1 ) + o 6r 3 3 m (x x 1 )(y y 1 ) o η t ( 方位角 ) (.3.5) α = t + γ + T 1 t 1 γ:p における子午線収差 T 1 t 1 は方向補正 α = t + γ + T t (.3.6) ( 方向補正 ) T 1 t 1 = x x 1 (y 6r 1 + y ) η t 1 9r 3 (x x 1 ) (y y 1 ) + η t 6r 3 (y y 1 )(3y 1 + y 1 y + y ) T t = x x 1 (y 6r 1 + y ) + η t 1 9r 3 (x x 1 ) (y y 1 ) η t 6r 3 (y y 1 )(y 1 + y 1 y + 3y ) R = MN (.3.7) M = a(1 e ) W 3 N = a W W = 1 e sin φ e = a b a [ 例題.3] 平面直角座標系における 点 A,Bの値が与えられているものとする これらの値を用いてプログラム (SO K3.BAS) により 平面距離 (s) 球面距離(S) 方向角(t) および点 A における点 B の方位角 ( α) を求めてみよう ただし 点 A(-144,654.741m,107,365.335m), 緯度 φa=34º41 5. 0000 点 B(-155,04.18m,118,187.713m) および点 Aにおける子午線収差は γ=+0º43 54.09 である ( 解 ) 平面距離 s= 15,000.785 方向角 t= 133 49 31.016 球面距離 S= 14,999.931m 方位角 α= 133 49 8.7
5.4 経緯度から平面座標の計算 地球の基準楕円体 1 に横円筒をかぶせて各地点をその円筒上に投影しておき その円筒を切り開けば 楕円体 から平面への写像が直接的に得ることができる これは横メルカトール図法 と呼ばれている 測量で多く用い られる横メルカトール図法としての平面直角座標系では 投影原点の縮率を mo=0.9999 として距離の精度が 1/10,000 以内に収まるように配慮されている 地点の緯度 経度 (φ, λ) から平面座標 (x,y) および子午線収差 (γ) を求める式は 次のとおりである x=m o [ φ + λ 4 Δλ Nsinφcosφ + Nsinφ 4 cos3 φ(5 t + 9η + 4η 4 ) + λ6 70 Nsinφcos5 φ(61 58t + t 4 + 70η 330t η )] (.4.1) y = m o ΔλNcosφ + Δλ3 N 6 cos3 φ(1 t + η + λ 5 10 Ncos5 φ(5 18t + t 4 + 14η 58t η ] (.4.) γ = Δλsinφ + Δλ3 3 sinφcos φ(1 + 3η + η 4 ) + Δλ 5 15 sinφcos4 φ( t ) (.4.3) Np x Δλ=λ-λo 原子午線 子午線 y γ P φo φ Δφ δφ x y φ λ 赤道 λo 平行圏 図.4.1 緯度経度と平面座標 Δφ = a(1 e )[k 1 φ φ o 1 k sinφ φ o + 1 k 4 3 sin4φ 4φ o 1 6 k 4 sin6φ 6φ o + 1 8 k 5 sin8φ 8φ o 1 10 k 6 sin10φ 10φ o + ] (.4.4) k1=1.00505501813087 k=0.00506310864 k3=0.0000106759063 k4=0.0000000080379 1 reference ellipsoid: 準拠楕円体 transverse mercator projection: 横メルカトル投影
6 k5=0.00000000003934 k6=0.000000000000071 Δλ=λ-λ o λ o = 座標系原点の経度 t=tanφ η=e cosφ e = a b b [ 例題.4] 点 A の緯度経度は (φ=34º41 5. 0000, λ=135º30 19. 0000) とする この点の平面直角座標 (x, y) の値および 子午線収差 (γ) をプログラム (SOK4.BAS) を用いて求めてみよう ただし 平面座標系は第 5 系 (φo=36º,λo=134º0 ) とする ( 解 ) =-144654.741 Y= 107365.335 γ =+ 0 40 1.431.5 平面座標から経緯度への変換ある地点における平面直角座標 (x,y) から緯度経度 (φ, λ) および子午線収差 (γ) を求めるには 次の式で行える ) MN φ = φ 1 (y m o λ = λ o + γ = y mo y mo 3 Ncosφ y mo N t y mo 3 3N 3 t + y mo 4 t(5 + 4MN 3 3t + η 9t η ) y mo 6 t(61 + 70MN 5 90t + 45t 4 )...(.5.1) y mo 5 (1 + 6N 3 cosφ t + η ) + (5 + 10N 5 cosφ 8t + 4t 4 )...(.5.) φ 1 = φ1 (i 1) S (i 1) D M (i 1) (1 + t η ) + y mo 5 t( + 5t + 3t 4 )...(.5.3) 15N 5 D = a(1 e ) k 1 φ o k sinφ o + k 3 4 sin4φ o k 4 6 sin6φ o + k 5 8 sin8φ o k 6 10 sin10φ x o + + m o S = a(1 e )(k 1 φ 1 k sinφ 1 M = a(1 e ) (1 e sin φ 1 ) 3/ + k 3sin4φ 1 4 k 4sin6φ 1 6 + k 5sin8φ 1 8 k 6sin10φ 1 10 + ) [ 例題.5] 点 A の平面直角座標 xy の値は第 5 系 (-144,654.741m,107,365.335m) とする この点の緯度経度 (φ,λ) 及び子午 線収差 (γ) をプログラム (SOK5.BAS) を用いて求めてみよう ( 解 ) 緯度 B= 34º41 5.5018 経度 L= 135º30 18.5040
7 子午線収差 (γ)=+ 0 40 1.1573.6 UTM 図法 UTM 投影 (Universal Transverse Mercator Projection) による平面座標は世界的に利用されている わが国では縮尺 1/1 万地形図から 1/5000 基本図,1/5 万編集図における投影に採用されている UTM 座標系は赤道を 6 ずつ東回りに 60 個のゾーンに分割して横メルカトールで投影する図法である つまり 最初のゾーン 1 は東周りに西経 (-)180 ~174 範囲 ( 中央経線西経 177 ) ゾーン は西経 174º~168º とし ゾーン 30 は西経 6 ~0 ( 中央経線西経 3 ) ゾーン 31 は東経 0 ~6 ゾーン 60 は東経 174 ~180 ( 中央経線東経 177 ) である 因みに 大阪や沖ノ鳥島はゾーン 53 で 東経 13 ~138 中央経線東経 135 ( 明石市 ) に該当する ( この計算は utm_zone.xlsm で作成している ) 日本はゾーン 51~56( 中央経線 13º( 尖閣 ) 19º( 沖縄 ) 135º( 大阪 ) 141º( 東京 ) 147º 153º( 北方四島 )) を採用する 各経度帯内ではその中央経線と赤道との交点を座標原点として横メルカトル図法で投影する 中央経線から東西方向に離れるに従って距離ひずみが増大するので 平面距離 s と球面距離 S の点の投影誤差 ds/ds の比率は ±4/10000 で収まるように考える つまり 中央経線上の縮尺係数は 0.9996 に設定するので 中央経線より東西に 180km 離れたところの縮率は 1.0000 そして同様に 70km 離れた地点の縮率は 1.0004 となる UTM 投影では中央経線より 70km までは 距離の相対精度が ±4/10000 に収まるので 各 6 幅の経度帯内で地球を平面とみなして投影することができる ( ただし 原点から 3 離れると 334kmなので 1.0004 を超える )UTM 座標においてx 座標 (N 軸 ) は中央経線 ( 縦軸 ) に沿って北半球では 0 より増加し 南半球では 10,000km( 擬北距 =F.N.) より減少させる y 座標 (E 軸 ) は赤道 ( 横軸 ) を 500km( 擬東距 =F.E.) として E 軸に沿い東側へは増加 西側へは減少する Bessel 緯度経度 ( 以下の Ki は GRS80 での値である ) と xy(xn,ye) 座標との関係式は 以下のように表すことができる.6.1 緯度経度から UTM 座標の変換 x N = F. N. +m o [Δφ + Δλ Δλ4 Nsinφcosφ + 4 Nsinφcos3 φ(5 t + 9η + 4η 4 ) + Δλ 6 70 Nsinφcos5 φ(61 58t + t 4 + 70η 330η t )] (.6.1) y E = F. E. +m o [Δλ Ncosφ + Δλ3 6 Ncos3 φ(1 t + η ) + Δλ 5 10 Ncos5 φ(5 18t + t 4 + 14η 58t η )] (.6.) γ = Δλsinφ + Δλ3 3 sinφcos φ 1 + 3η + η 4 + Δλ 5 15 sinφcos4 φ( t ) (.6.3) Δφ = a(1 e )[k 1 φ φ o 1 k sinφ φ o + 1 4 k 3 sin4φ 4φ o 1 k 6 4 sin6φ 6φ o + 1 k 8 5 sin8φ 8φ o 1 k 10 6 sin10φ 10φ o + ] (.4.4) K1=1.005 05 501 813 086 71
8 K =5.0631086411D -3 K3 =1.067590633386D -5 K4 =.080378571650D -8 K5 = 3.9337137114D -11 K6 = 7.10845340859D -14 Δλ=λ-λo λo= 座標系原点の経度, t=tanφ η=e cosφ e = a b b mo=0.9996 a=6,378,137m 1/f=98.57101 南北半球では F.E.=500km, 北半球では北距 F.N.=0m 南半球 F.N.=10,000km [ 例題.6.1] ある三角点の GRS80 による世界測地系の緯度 経度は次のとおりとするとき この点の UTM 座標および子午線収差 ( 角 ) をプログラム (SOK61.TT) により求めてみよう φ=35º39 9.157 λ=139º44 8.8869 また この点の UTM 座標の原点は次の値である φo=0º λo=141º ( 解 ) = 3,946,757.90 Y= 386,070.956 子午線収差 γ =- 0 o 44 01.684.6. UTM から緯度経度への変換 φ = φ 1 y E NM t + y E 4 4MN 3 t(5 + 3t + η 9t η ) y E 6 70MN 5 t(61 + 90t + 45t 4 ) (.6.6)...(..1) λ = λ o + [ y E Ncosφ y E 3 6N 3 cosφ (1 + t + η ) + y E 5 10N 5 cosφ (5 + 8t + 4t 4 )] (.6.7) γ = y E N t y E 3 3N 3 t(1 + t η ) + y E 5 15N 5 t( + t + 3t 4 )...(.6.8) 南北半球では y E = (y E 500,000[m])/m o
9 m o = 0. 9996 ただし γ は子午線収差角である φ 1 = φ1 (i 1) S (i 1) D M (i 1) D = a(1 e ) k 1 φ o k sinφ o + k 3 4 sin4φ o k 4 6 sin6φ o + k 5 8 sin8φ o k 6 10 sin10φ x o + + m o S = a(1 e )(k 1 φ 1 k sinφ 1 + k 3sin4φ 1 4 k 4sin6φ 1 6 + k 5sin8φ 1 8 k 6sin10φ 1 10 + ) M = a(1 e ) (1 e sin φ 1 ) 3/ (..4) 北半球では x = x N F. N. = x N 南半球では x = x S F. N. = x S 10,000,000[m] (..5) である [ 例題.6.] ここでは 例題.6.1の逆計算をプログラム (SOK6.TT) を用いて行ってみよう ただし 緯度経度を求めるために用いる UTM 座標は次のとおりである xn=3,946,757.90m yn=386,070.956m ( 解 ) Input data(utm x,y) n= 3,946,757.90m Ye= 386,070.956m LatitudeB=35 39 9.157 LongitudeL=139 44 8.8869 Meridian Convergence γ= - 0 44 01.684 3. GPS 観測値の平面直角座標への計算 GPS 観測値であるWGS84における緯度 経度 楕円体高 ( ϕ, λ, h) から Bessel の (φ, λ, h) に変換する式は以下のとおりである (1) GPS 観測値参照点 または求点のDGPS 観測から得られたWGS84 1 の経緯度および楕円体高は 次のとおりとする φ GPS = λ...(3..1) h WGS84 1 World Geodetic System(1984):1984 年に人工衛星を利用して決められた回転楕円体
10 () GPS 観測値から WGS84 の YZ への変換 WGS84 楕円体における 3 次元直交座標 YZ は 次式で算出する WGS84 = (N 84 + h 84 )cosφ 84 cosλ 84 Y WGS84 = (N 84 + h 84 )cosφ 84 sinλ 84 (3..) Z WGS84 = [N 84 (1 e 84 ) + h 84 ]sinφ 84 なお 以上で用いた要素は 次のとおりである 卯酉線曲率半径 N 84 = a 84 W 84 W 84 = 1 e 84 sin φ 84 扁平率の逆数 1/f 84 = a 84 a 84 b 84 短半径 b 84 = a 84 (1 f 84 ) 第一離心率 e 84 = a 84 b 84 a...(3..3) 84 また WGS84 の地球楕円体の値は 次のとおりである a 84 = 6,378,137m 1 f 84 = 98.573563 b 84 = 6,356,75.314m (3..4) (3) WGS84 から ITRF89 系座標への変換 Y = R T Y 0.055 0.541 (3..5) Z ITRF Z 84 0.185 (4)ITRF89 系から筑波座標系 9(Bessel) 座標への変換 Y = R T Y 147.531 508.313 (3..6) Z Tsukuba Z ITRF 680.417 (5) Bessel 経緯度への変換 Y = Y (3..7) Z Bessel Z Tsukuba Bessel 経度緯度および正規標高 (Orthometric Height) は これらの YZ 値を用いて 以下のように求める ことができる tanλ B = Y/...(3..8) ( 近似式 ) tanφ B = tanθ = ( 厳密式 ) Z+(e ) bsin 3 θ +Y e acosθ Z a b +Y...(3..9)
11 tanφ B = Z+N Be B sinφ B +Y H B = +Y N cosφ B B h B = H B H G...(3..10) 上に示した Bessel 楕円体上の標高は 次式で求められる hb=( 水準標高 H)+( ジオイド高 HG)...(3..11) また ジオイド高 1 はジオイド図を参照して求められる Bessel 楕円体は 次のとおりの値である a B = 6,377,397.155m b B = 6,356,078.963m 1 f B = 99.15813...(3..1) GRS80(geodetic reference system 1980) では 長半径 a= 6 378 137m 及び扁平率 f= 1/98.57 101 b=a(1-f)= b = 6 356 75.314 140 356m WGS84(GPS 衛星で使用 ) a=6378137m, f=1/98.573563 ロシアの 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 の国際基準座標系の 3 軸を慣用極方向として 軸はグリニジ子午線と赤道面との交点と地球重心を結ぶ軸 Y 軸は赤道面上で 軸に直角な軸で -Y-Z は右手系を採用 PE-90 楕円体は長軸半径 a=6378136.0m 扁平率の逆数 1/f=98.57839303 である 1 基準楕円体表面の標高をゼロとしたときの平均海面高 つまり楕円体面からのジオイド面の高さ