子午線長さの計算


 
楕円孤の長さ
 

 長円半径 a 短円半径 b の楕円定義は次の2次関数となる。
x2/a2+y2/b2=1

 各点の接線勾配は、次の式。
dy/dx=-b2/a2 × x/y

 緯度の微小な増分による楕円孤長の増分d l は

d l 2 = d x 2 + d y 2

 dy を dx で表現して
d l 2 = d x 2 + (-b2/a2 × x/y )2× d x 2

 y を x で表現して
d l 2 = ( 1 + b2 / a2 × ( x/a )2 / ( 1 - ( x/a )2 ) ) × d x 2

 x/a を X として
d l 2 = ( 1 - b2 / a2 × X2 / ( 1 - X2 ) ) × a 2  d X 2


 1 - b2/a2 を e 2 として整理すると

d l =[ ( 1 - e2 2 ) / ( 1 - X2 ) ]1/2 × a  d X



 水平投影距離 x が 0 からx までの楕円孤長は下記となる。

L = a ∫[( 1 - k2 2 ) / ( 1 - X2 ) ]1/2 d X  ( X = 0 〜x/a )


 ベッセルの第2種楕円積分公式を利用する。


 
ベッセルの第2種楕円積分
 


L = a ∫[( 1 - e2 2 ) / ( 1 - X2 ) ]1/2 d X  ( X = 0 〜x/a )



 上式は極座標表現で以下となる。

L = a ∫( 1 - e2sin 2 θ ) 1/2 d θ  ( θ = 0 〜φ )


 この場合θは北極から赤道に向けての角度であることに注意。

楕円孤の長さは初等函数では表せず、以下の級数展開となる。

2p = ∫( 1 - e2 sin 2 θ ) 1/2 d θ  ( θ = 0 〜φ )として


E( φ, e ) = φ - Σ [(2(p-1)-1) ! ! / (2p) ! ! ] × e 2p ×S 2p (φ)
 ( p = 1 〜∞ )


 (2p) ! ! は偶数の2重階乗。例えば8 × 6 × 4 × 2

 (2(p-1)-1) ! ! は奇数の2重階乗。例えば7 × 5 × 3 × 1


 S 2p (φ) は∫ ( sin2p θ ) d θ  ( θ = 0 〜φ ) で


 2P C P / 2 2P)φ +Σ { (-1) (p+r) 2pr / [2 2p(p-r) ]×sin[2(p-r)φ] }
 ( r = 0 〜 p-1 )

2pr は二項係数。2p ! / (r ! (2p-r) ! )
例えば63 = 6 × 5 × 4 × 3 × 2 × 1 / (3 × 2 × 1 ) /(3 × 2 × 1 )


気の遠くなる計算式だが、整数の係数とサイン関数で展開できる。 これはルジャンドルとかヤコブとかベッセルらの19世紀ヨ-ロッパ数学の叡智。 級数展開は無限で正解なのだが、工学的には5、6項目で有効数字以下となる。 と言い聞かせて体力勝負の展開をする。ただし、EXCELの力を借りる。


 
子午線の長さ
 

 地球子午線全周の1/4を求める。積分範囲は0〜π/2。
この場合の積分定数がベッセル第2種完全楕円積分E( e )。
L = a E( π/2, e )

E( π/2, e ) = E( e ) のsin成分はすべてゼロとなり、
第1項のみで求まる

地球の寸法(GRS80楕円体)
長円半径 a = 6378137m 短円半径b=6356752.314m
扁平率 f= 1/298.257222101 
離心率 e 2 = 0.00669438

子午線長(緯度0〜90度)
 L =10、001、965.7292m


 1792年より7年の歳月をかけて、ほぼ同経度にある2つの都市、ダンケルク(ベルギー)とバルセロナ(スペイン)間の距離をフランス人が測量した。 フランス革命中に、地球の北極から赤道間までの距離の1千万分の1を1mと定義するためであった。その歴史が、いまでも子午線長に残っている。

 
ベッセル第2種楕円積分の級数展開
 

 ベッセル第2種楕円積分の第一項は上の様に子午線全長で検証した。次に中間値の検証をしよう。sin2φ、sin4φ、sin6φ....のどこまでが必要なのだろうか。
 sin16φまでの展開式を整理した。


 今、地球の子午線距離を問題にしている。sinや角度ラジアンはせいぜい少数以上一桁。1万kmの計算で0.1mmまで求められれば工学的には理論値。GRS80回転楕円体の緒元を代入してφとsinの係数を小数点以下4桁になるまでの展開項数は?

 
結論はsin6φまで

 地球の1/4子午線距離をさらに30°、60°、90°の1/3区間で計算してみよう。それを合計して1/4子午線距離になれば計算式の照査ができる。

北極より0〜30°
3338617.5226m
33.3796%
30〜60°
3333989.5489m
33.3333%
60〜90°
3329358.6578m
33.2871%
合計
10001965.7293m
100.00%

 学生時代の「数学公式T」岩波全書221を久しぶりに開いた。曲面アーチの計算では よく利用した懐かしい本だ。あの当時、メモ用紙を用意して腕力計算をしたが、もうその気力と体力はない。しかし、使いなれたマイクロソフトEXCELで華麗に計算結果を整理した。いい時代になったものだ。EXCELには二項係数、2重階乗、は関数として完備されている。大した物だ。


 
赤道からの子午線長
 

 ベッセル第2種楕円積分の場合、X=sinθの変数置換を利用している。この場合のθは北極からの角度である。緯度は赤道面からの角度であるから都合が悪い。指定された緯度φの赤道面からの子午線長は次ぎのように求める。

 L = ∫ (θ = 0 〜π/2 )- ∫ (θ = 0 〜π/2 - φ )

 つまり、1/4子午線長から、その地点の北極からの子午線長を引くことになる。級数展開式を 眺めてみると、第1項目は定数でありsinの項は含まれない。第2項目のsin2θ等の項は
sin2(π/2 - φ)=sin( π - 2φ)= sin( 2φ)
sin4(π/2 - φ)=sin(2π - 4φ)=-sin( 4φ)
であるから整理できる。sin( 10φ)までの級数展開として、直接緯度φの赤道面からの子午線長の計算式を誘導しておく。


 この式を使い、赤道から緯度9度ごとの子午線長を計算してみる。
赤道面より0〜9°
998543.856m
9.983%
9〜18°
998705.823m
9.985%
18〜27°
999013.831m
9.988%
27〜36°
999437.614m
9.992%
36〜45°
999935.573m
9.997%
45〜54°
1000458.892m
10.003%
54〜63°
1000956.347m
10.008%
63〜72°
1001379.315m
10.012%
72〜81°
1001686.508m
10.015%
81〜90°
1001847.971m
10.017%
合計
10001965.729m
100.00%

 まだ地球の形が定かでなかった時代、このような緯度による子午線の長さを各所で測量し、最小二乗法で長半径、短半径、扁平度を推定したのだろう。

 
地心緯度と測地緯度
 

 2000年を節目として、日本は測量の原点である地球の形を、ベッセル回転楕円体からGRS80楕円体に変更した。回転楕円体では、地球の中心軸を通る鉛直面は、どこでも同じ楕円形である。水平面で輪切りにすれば、どこでも円であることになる。            
準拠楕円体 GRS80楕円体ベッセル楕円体
長半径 a 6,378,137km6,377,397km
短半径 b 6,356,752km6,356,078km
扁平率 f=(a-b)/a 1/298.257222101 1/299.152813
離心率e=(a 2-b 2) -/2/a 1/12.222 1/12.240

 1840年頃から、測量データにより色々な楕円体が提唱されてきた。日本の近代測量術は、ベッセル楕円体を採用して出発している。楕円体といっても、限りなく円に近い。もし、正確な地球を黒板に書いてみろと言われたらコンパスで円を描く。 中学の教師の言葉が今でも真実だろう。
 楕円は2つの焦点を持つ円とも言える。楕円中心から焦点までの距離は長半径 a の離心率倍の位置にある。a /12 程の所にある。
 地球が円ではなく楕円だということは、近代測量術の基本であった。だが、日本国という限定地域で考える場合、それに近似する円を想定して測量計算式が展開されていた。僕もそのような教科書で講義を受けた。グローバルな全地球的時代には、「測量学再入門」が必要なのだろう。
 もう一度、測地緯度とは何か、図を見て確認しよう。


   測量術の準拠楕円体とは何か。地球重力圏に捕捉された水圏の理想表面ともいえる。これをジオイドと呼ぶ。水平面を定義して標高の拠り所である。重力はこの面に直角に働く。測量機械に吊るされた「下げ振り」の方向が「地心」である。「その先に地球の重心位置がある」とするのだが、それは長年の習慣でそう思うだけの嘘。楕円体と決めた時点で重心方向の角度は地心緯度ではなくなる。測地緯度、あるいは地理学的緯度と区別される。測量と地図の緯度は測地緯度。GPSやGISでも測地緯度を用いる。地球から外の天体に目を向ける天文学の緯度とは区別される。

上の図では、θ=測地緯度 φ=地心緯度。

tan φ = (1 -e2) ×tan(θ) の関係がある。

測地緯度
地心緯度
00°00′00″
00°00′00″
09°00′00″
08°56′27″
18°00′00″
17°53′14″
27°00′00″
26°51′41″
36°00′00″
35°49′02″
45°00′00″
44°48′47″
54°00′00″
53°49′01″
63°00′00″
62°50′38″
72°00′00″
71°53′12″
81°00′00″
80°56′25″
90°00′00″
90°00′00″

 日本の地図の緯度は国土地理院から測地緯度として提供される。楕円積分の場合、地心緯度に一旦変換して子午線の長さを計算することになる。

 
測地緯度からの子午線長
 

 測量計算において緯度から子午線長を求めるということは、測地緯度から直接子午線長を計算することである。国土地理院の「精密測地基準点作業規定」に計算式が掲載されている。

S0 = a(1-e2) [Aφ0-B/ 2×sin2φ0+C/ 4×sin4φ0-D/ 6×sin6φ0
+E/ 8×sin8 φ0 -F/10×sin10 φ0]

a 6378137 m A 1.005052501813080
e2 0.00669438   B 0.005063108622224
        C 0.000010627590263
        D 0.000000020820379
        E 0.000000000039324
        F 0.000000000000071

 XMLデータとして公開されている数値地図25000の経緯度座標値を公共座標XYに変換するには、この式を利用する。このページは、この計算式の物理的意味と精度、適用限界を考察することを目的としている。
 XMLデータをXYに変換するには膨大な点数の計算をする。その数値はCADに取込まれるのであるが、現在のパソコン、CADソフトは3次元を扱える。GPSが主流になる 21世紀を見据えれば、地球3次元XYZデータに変換して蓄積したほうが賢明なのかもしれない。 3次元CADはグローバル座標と共にユーザー座標を持っている。指定したユーザー座標平面に瞬時に3次元データを投影する。地図投影法はパソコンの回路とソフトがするのである。地図投影法により変形したデータは極力排し、ピユアなXYZデータをCADに送り込むべきだろう。地図投影法は曲面を平面に投影する知恵と工夫の叡智なのだが、新しい時代の新しい器には、発想転換したデータが必要ではないだろうか。

 さて、「精密測地基準点作業規定」のBLXY座標変換式に用いる子午線長計算式は、 GRS80の場合次ぎのような式になる。

So(φ) = 6367449.1459 φ -16038.5087 sin(2φ) +16.8326 sin(4φ) -0.0220 sin(6φ) +0.0000 sin(8φ) +0.0000 sin(10φ)

一方、楕円積分の式は下記。

Lo(φ) = 6367449.1458 φ -5346.1696 sin(2φ) -1.1222 sin(4φ) -0.0006 sin(6φ) +0.0000 sin(8φ) +0.0000 sin(10φ)


測地緯度
地心緯度
Soの子午線長
Loの子午線長
00°00′00″
00°00′00″
0m
0m
10°00′00″
09°56′04″
1105854.8
1102220.9
36°00′00″
35°49′02″
3985542.7
3975394.9
40°00′00″
39°48′38″
4429529.0
4419011.6
45°00′00″
44°28′27″
4984944.4
4974252.2
90°00′00″
90°00′00″
10001965.7
10001965.7

 両式には誤差が生じている。以下考察をしてみる。

子午線長計算式で緯度90度の場合1/4地球周長を与える。両者は一致する。 この数字が一致することは、a、e2の数値とその取扱いに 間違いがないことになる。なを、 子午線長計算式の緯度90度の場合、第1項のみが有効で以下sinを含む項は0となる。

 子午線長計算式の緯度45度の場合、第2項以下はsinπ/2、sinπ、sin2πとなり 第3項以下は0と考えて差し支えない。緯度45度の場合第2項の係数を評価することになる。両者で大きく異なる。これは、測地緯度と地心緯度の関係を調整しても解消しない大きさである。
 So(φ)式はLo(φ)の級数展開項を元に,測地緯度から地心緯度への変換を全項に(1 - e2)倍することで処理している。だが、φの項とsin2φ等の項への重みは一律にすべきではない。おそらく旧日本測地系にフィットする係数を最小二乗法で決定して運用されていたのであろうが、その場合、緯度90度での子午線長 を犠牲にして緯度40度近辺の精度を確保する係数を決定すべきではなかったか。

 So(φ)式は日本の公共測量大座標XYの計算に用いられるが、系の基準点と求点の子午線長差として運用される。赤道面からの子午線長として単独では用いられない。この場合、SoとLoの差は1/1200程度である。

 So(φ)式は旧公共測量大座標値との整合、地図展開のための係数等を考慮しての計算式と推定される。部分的にこの計算式を引用して地球規模の計算に用いることは結果の保証をするものだはないと推測する。

 以上は私個人の見解である事を付記しておきます。

 
海峡を渡った大三角網
 

 本来、測地成果は国家主権の範疇である。国境で経緯度が異なることはあり得る。海峡で囲まれた日本国においては、国境の経緯度はそれ程逼迫した問題にはなり得ないが、一度だけ大三角網が海峡を越えたことがある。

 日本の近代測量術は、1871年(明治4年)にイギリス人技師の指導を受け内務省で始まった。その後、1878年に陸軍の参謀本部に地図課と測量課が誕生した。陸地測量部の前身である。1982年にドイツに3年間留学した田坂大尉が帰国したことにより、測量はドイツ式に変わった。BL計算という用語にドイツ語Breite(緯度)が今でも残る。
 殖産興業、富国強兵の国家事業の一環として、辺長約45kmの一等三角点(本点)から辺長約1.5kmの四等三角点まであって、全国に整備されて行く。 この三角網は海峡を渡り朝鮮半島にまで広がった。1932年(昭和7年)に満州国長春 (当時の新京)に基準点を設けて測量を開始し、その三角網が朝鮮との国境の鴨緑江に達したとき、2つの三角網に大きな差があることが判明した。経度が17.0秒、緯度が9.4秒と無視できない大きさであった。この原因は当時不明であった。天体観測の結果、日本の基準点に誤差があると判明したのは終戦前後の時期である。

 日本列島の南東には深度8000mに及ぶ海溝がある。大陸の縁辺部に位置する日本においては重力方向は地下の大陸塊に引き付けられ傾いている。全地球規模の理想球体に準拠すれば特殊事情が生じるのである。これは日本国の宿命である。

 昭和7年の頃、「不可解」として処理された誤差は、GPSの測量により明白のものとなった。この歴史的誤差を解消したことがJGD2000の持つ本質である。理想楕円体からの重力線の局所ひずみを「鉛直線偏差」と呼ぶ。17秒の角度を、明治から始まった日本国測量プロジェクトは見事に検知していたのである。
 世界標準を採用して21世紀に乗り出した日本は、「鉛直線偏差」を国内事情として抱え、他国より状況は格段に厳しい。だが、先達の 足跡の上に測量界の叡智と努力が脈々と継承されて行くと信じるのである。




2003.11.30
by Kon