mikeo_410
  1. 2.地球
    1. 1.地球楕円体上の点
    2. 2.地球儀と楕円体の描画
    3. 3.3Dの直線と平面
    4. 4.卯酉線
    5. 5.楕円の標準形
    6. 6.5点で決まる楕円
    7. 7.曲率半径と平均曲率
    8. A1.地図、地球儀(Rスクリプト、データ)
    9. A2. 付録6.計算式集の公式

曲率半径と平均曲率

1. 楕円の曲率半径

楕円を描いている曲線は、短い区間を円で近似するようです。この円の半径が曲率半径で、その逆数が曲率のようです。

二階微分可能な曲線  の点  における曲率半径 R は、

        

で、 のときは、曲率半径 ∞ です。

楕円を、

        

と表すと、

        

maximaで微分して、

        

したがって、曲率半径 R の2乗は、

                        1.1

2. 平均曲率半径

準拠楕円体上の1点における曲面の曲率を考えます。地表のある地点で、その点を通る鉛直方向の直線を含む平面で楕円体を切断します。断面は楕円で、切断する平面の向きによって曲率が変わります。その平均値が平均曲率だと考えられます。その逆数は平均曲率半径です。

方位盤は水平に置かれ、その地点の接面です。錘を下げて決まる鉛直方向は法線です。

法線を含む平面は法面で、法線と、方位盤の方向を示す直線で決まる面です。

この法面と楕円体の交線は、垂直截線と呼ばれるようです。方位盤の方位を示す目盛りは、垂直截線の一部です。

方位は様々な表現方法がありますが、真北から時計回りに360度まで測る「方位角」が共通の基準のようです。方位角は、その地点を通る経線と垂直截線の成す角度であり、経線が作る平面と、法面の成す角角度です。

測量法に基づいた「作業規程の準則」の付則「付録6計算式集」には、平均曲率半径として、 6,370,000 m がでています。また、公式が示されいて、緯度36度の平均曲率半径は、6371488.6206 m となります。

GRS80の赤道半径は、6,378,137 m です。極半径は、6,356,752.314140356 m となります。赤道は真円ですが、経線は楕円を描きます。

経線上では、赤道との交点で最も曲率半径が短く、極で最長となることは明らかです。式1.1によって計算すると、赤道で 6,335,439.327083874 m、極で 6,399,593.625864022 m となります。

極を通る経線はすべて同じ楕円ですから、極における平均曲率半径は、6,399,593.6259 m です。極以外の地点では方位角によって曲率が変わることになります。

「付録6計算式集」の平均曲率半径は、子午線曲率半径と卯酉線曲率半径の幾何平均です。

つまり、方位角0° と 90° の垂直截線の曲率半径の幾何平均です。

準拠楕円体は極軸に回転した回転楕円体なので赤道は真円です。赤道の曲率半径は赤道半径そのものです。子午線曲率半径と卯酉線曲率半径の幾何平均は経度に無関係に、緯度だけで決まります。どの地点でも、その地点を通る経線が最小の曲率半径となり、卯酉線が最大の曲率半径となることが前提のようです。

緯度±90°の極では、6,399,593.6259 m です。緯度0°の赤道上の点では、経線の曲率半径と赤道半径から、6,356,752.3141 m となります。この値は、極半径と同じ値のようです。

3. 平面直角座標9系の原点における平均曲率半径

平面直角座標9系の原点は、緯度北緯36度、東経139度50分と定められています。

この地点に方位盤を置くことを考えます。真北を定めて。方位盤を水平に置きます。方位盤の中心を通る鉛直方向の直線は法線で、赤道面と成す角が 36° です。方位盤の南北の線は、東経139度50分の経線です。

したがって、東西の直線は、東経139度50分の経線に直交する卯酉線です。他の方位を示す直線も全て、法線を含む法面と、方位盤の交線です。

方位角は、その地点を通る経線が作る平面を基準として、法線の周りに回転する方面との成す角度で表されます。真北を0°として、時計回りに360°まで測られます。

全ての方位に対して、準拠楕円体の法面による断面が考えられ、それぞれ曲率半径があります。

準拠楕円体の中心から、東経139度50分の経線と赤道の交点へ向かう直線をX軸とします。準拠楕円体の中心から、自転軸の北の極へ向かう直線をZ軸とします。準拠楕円体の中心から、西経130度40分の経線と赤道の交点へ向かう直線をY軸とします。

卯酉線のイメージは、下図のようになります。地球は、ほとんど球ですから、経線の描く楕円を紙に書いても円にしか見えません。図は大きく強調したものです。

「卯酉線」に従って、平面直角座標9系の原点は、経度0°方向をX軸としたXYZ座標で、

(-3947708.577205533,  3332137.506016432,  3728191.675729482)

と計算されます。

  1. r$> func <- function(φ, λ){
  2.         f = 1/298.257222101
  3.         b <- 1-f
  4.         β <- atan(b*tan(φ))
  5.         list(x=cos(β)*cos(λ),y=cos(β)*sin(λ),z=b*sin(β))
  6.     }
  7.     φ <- 36/180*pi; λ <- (139+50/60)/180*pi
  8.     u <- func(φ, λ)
  9.     a <- 6378137
  10.     c(u$x,u$y,u$z)*a
  11. [1] -3947708.577205533  3332137.506016432  3728191.675729482

平面直角座標9系の原点の経度を0°としたXは、以下のようになります。

  1. r$> sqrt(u$x^2+u$y^2)*a
  2. [1] 5165998.777539881

あるいは、

  1. r$> u <- func(φ, 0)
  2.     u$x*a
  3. [1] 5165998.77753988
3.1. 経線の曲率半径

式1.1に、この値を適用すると、曲率半径 6357482.4375 を得ます。平面直角座標9系の原点における、経線の曲率半径です。

3.2. 卯酉線の曲率半径

準拠楕円体は回転楕円体なので、経線は全て同じ楕円です。平面直角座標9系の原点を通る経線を0°としても不都合はありません。これで、平面直角座標9系の原点を通る経線が作る平面は、XZ平面と考えることができます。

平面直角座標9系の原点Aは、

となります。

計算を簡単にしたいので、赤道半径を1として計算します。GRS80の赤道半径 aで割って、

        a = 6378137

点Aにおける法線の切片(地軸Zとの交点)s は、

と表せ、法線は (0, 0, s)を通ります。

卯酉線が作る平面は法線を含み、、(0, 0, s)を含みます。また、XZ平面にある経線に直交するので、Y軸に平行です。したがって、任意のyに対して(0, y, s)を通ります。

、(0, 0, s)、(0, 1, s)の3点で決まる卯酉線が作る平面を、

とします。

maxima で、L、M、N、O を求めると左図のようになります。

回転楕円体は、

        

です。平面の式と、回転楕円体の式の連立方程式を解いて、卯酉線上の点の座標を4つ求めます。

まず、卯酉線と赤道面の交点2つを求めます。z = 0 として、x、 y を求めます。

次に、x = 0、z = s として、直線ベクトル(0,1,s)との交点を求めます。

maxima の表示はコピーしてRのスクリプトに貼り付けられます。

a=1、b=0.6 と強調して図示すれば下図のように描けます。

卯酉線を卯酉線の作る平面上で考えると、準拠楕円体が地軸を軸とした回転楕円体であることから、法線に対して左右対称です。卯酉線は、法線を短軸とした楕円です。

地球楕円体の比率で描くと、左図のようになります。地球楕円体は、ほとんど真球なので、法線が回転楕円体の中心を通るように見え、卯酉線と赤道面との交点と、卯酉線と法線とZ軸の交点を通るY軸と同じ方向ベクトルの直線との交点は区別できません。

オレンジ色の線は、卯酉線を、赤道面との交線を軸にXY平面まで回転したものですが、赤道と重なって図では区別ができません。

法線と赤道面のなす角が緯度なので、法線が作る平面を、赤道面との交線を軸に、緯度だけ回転してXY平面に重ねます。

  1. r$> x<-cos(β)-(b*sin(β))/tan(φ)
  2.     y<-sqrt((1-cos(β)^2)*tan(φ)^2+2*b*cos(β)*sin(β)*tan(φ)-b^2*sin(β)^2)/tan(φ)  
  3.     y2 <- sqrt(-cos(β)^2*tan(φ)^2+2*b*cos(β)*sin(β)*tan(φ)-b^2*sin(β)^2+b^2)/b
  4.     x2 <- 0
  5.     x <- c(Ax, x,  x, x2,   0)
  6.     y <- c(0,  y, -y, y2, -y2)
  7.     z <- c(Az, 0,  0,  s,   s)
  8.     u <- Rotate(L, φ, list(x=x,y=y,z=z))
  9.     u
  10. $x
  11. [1]  0.999878441669715401  0.005422141138500591  0.005422141138500591
  12. [4] -0.001279993892284230 -0.001279993892284230
  13. $y
  14. [1]  0.0000000000000000  0.9999853000846934 -0.9999853000846934
  15. [4]  0.9999921881746419 -0.9999921881746419
  16. $z
  17. [1] -6.808789643208968e-17  0.000000000000000e+00  0.000000000000000e+00
  18. [4]  1.604619215278547e-17  1.604619215278547e-17

卯酉線上の5点は回転によって全て z がゼロになります。点AはX軸上に移動します。

この5点で決まる楕円の式の係数から、楕円の標準形を求めます。

  1. r$> v <- EllipseCoeff(u$x,u$y)
  2.     S <- StandardEllipse(v$A,v$B,v$C,v$D,v$E,v$F)
  3.     S
  4. $a
  5. [1] 0.9999948989782719
  6. $b
  7. [1] 0.9988327146454221
  8. $rotation
  9. [1] 1.570796326794897
  10. $cx
  11. [1] 0.001045727024292981
  12. $cy
  13. [1] -3.420291656950161e-17

長軸 u$a をx軸半径、短軸 u$b をy軸半径として、横長の楕円を描けば、点Aの曲率半径は、x=0のときの値です。これに、赤道半径を乗じれば、地球規模の点Aにおける卯酉線の曲率半径になります。

  1. r$> R <- CurvatureRadius(S$a, S$b, 0)
  2.     R
  3.     a <- 6378137
  4.     R*a
  5. [1] 1.001158435562008
  6. [1] 6385525.660720157

緯度36°における卯酉線の曲率半径は 6385525.6607 m と計算されます。

3.3. 平均曲率半径

緯度36°における方位角0°、90°、つまり経線方向と卯酉線方向の2つの曲率の幾何平均は、

6371488.620 m

となりました。単純平均は、

6371504.049 m

です。

3.4. 公式

        

        

緯度φ=36°における、経線方向の曲率半径は、

緯度φ=36°における、卯酉線方向の曲率半径は、

        

緯度φ=36°における、方位角α の曲率半径は、

        

緯度φ=36°における、方位角45° の曲率半径は、

        

また、「付録 6 計算式集」に従って、

        

        

4. 全方位の曲率半径

ある緯度φ の地点をAとすると、地点Aにおける方位角α の方位線の曲率半径を計算します。

準拠楕円体は地軸の周りに回転した回転楕円体なので、等緯度の地点の状態は全て同じです。地軸をZ軸、回転楕円体の中心を原点、Aを通る経線の作る面をXZ平面、赤道面をXY平面として、普通のXYZ座標で考えます。

計算を簡単にするために赤道半径を1として計算します。

点Aの座標は、

方位線が作る、回転楕円体の切断面は、その地点を通る経線の作る面を、法線を軸に、時計回りにα回転したものです。

点Aを通る経線がXZ平面になるように座標軸を決めたので、経線による断面は左図のようになります。直線ASは法線で、点Sの座標は (0, 0, s) です。

        

        

点Aにおける法線は、(0, 0, s) を通り、方向ベクトル

        

        

です。

 は、 なので、

        

です。したがって、

        

また、曲線 y=f(x) の点 Q(Xq,Yq) における法線の公式は、

 

のようです。

yをzと読み替えて当て嵌めてみます。

楕円、

        

の式をzについて解けば、

        

となります。点Qは第1象限にあると考えているので、正の式を採って微分すると、

        

A(Ax,Az) における法線は、

となります。切片は、

        

となります。

点 A を通る経線が作る平面はXZ平面です。この平面は、点A を含む、法線ベクトル (0, 1, 0) の面と表されます。式では、y=0 と表されます。平面の式の係数L=0、M=1、N=0、O=0 で、式に現れない x、z は任意の値だと言うことのようです。

y=0 で表される平面を、点Aにおける法線を軸に、時計回りにα回転した平面は、点A、Sを含み、その法線の向きはα回転しています。

任意の直線の周りの回転は、各軸の周りの回転を組み合わせて行います。

直線や平面を表す方向ベクトルは、逆向きでも良いですし、始点や含む点にも制約はありません。これらの回転を整然と説明することができないので、かつて作成した回転関数に合わせて1例としたいと思います。

手順を説明するために、関数にしたものを掲げます。

  1. AzimuthCurvatureRadius <- function(φ, az, f=1/298.257222101) {
  2.     b <- 1 - f          # 極半径
  3.     β <- atan(b*tan(φ)) # 更成緯度
  4.     A0 <- c(cos(β), 0, b*sin(β))  # 緯度36度の地点Aの座標
  5.     S <- c(0, 0, (b-1/b)*sin(β))  # 法線の切片
  6.     if(abs(az)<24*.Machine$double.eps){
  7.         # 方位角がゼロの場合、XY平面では直線となるので別に計算
  8.         R <- CurvatureRadius(1, b, A0[1]) # 経線上の曲率半径
  9.         k <- list(a=1, b=b, rotation=0, cx=0, cy=0)
  10.         A <- c(A0[1],A0[3]) # 標準形の点Aの座標
  11.         O <- c(-1,0)
  12.         x <- A0[1]          # 曲率を求めたx
  13.     } else {
  14.         Normal <- Line2(A0, S)           # 法線
  15.         P <- list(x=c(0),y=c(0),z=c(1))  # XZ平面上の点
  16.         u <- Rotate(Normal, az, P)       # 平面を方位角だけ回転
  17.         PS <- Plane(A0,S,c(u$x,u$y,u$z)) # 方位角azの法面
  18.         # 垂直截線のXY平面への投影
  19.         qc <- ProjectEllipsoidSectionXY(PS$L,PS$M,PS$N,PS$O,b)
  20.         xy <- StandardEllipseXY(qc, 6) # 投影楕円の周の点列
  21.         z <- -(PS$L*xy$x+PS$M*xy$y+PS$O)/PS$N # x,yからz算出
  22.         # 切断面と赤道面の交線
  23.         IL <- Intersection2Planes(
  24.             list(L=PS$L,M=PS$M,N=PS$N,O=PS$O), list(L=0,M=0,N=1,O=0))
  25.         w <- AngleXY(PS) # 切断面と赤道面のなす角
  26.         xy$x[6] <- A0[1]; xy$y[6] <- A0[2]; z[6] <- A0[3]; xy$z <- z
  27.         u <- Rotate(IL, -w, xy) # 垂直截線をXY平面まで回転
  28.         A <- c(u$x[6], u$y[6])  # 回転後のA0のXY平面での座標
  29.         v <- EllipseCoeff(u$x,u$y) # 5点で決まる楕円の式の係数
  30.         k <- StandardEllipse(v$A,v$B,v$C,v$D,v$E,v$F) # その標準形
  31.         # 楕円の長軸(標準形のX軸)のXY平面での位置
  32.         O <- OrignCoeff(k, -1, 0) # 座標(-1,0)の行先
  33.         # 2直線の交点と、それぞれの直線状の点から成す角を求める
  34.         θ <- OriginAngle(c(k$cx,k$cy), O, A)
  35.         # 点Aから楕円の長軸に垂線を下して曲率を求めるxとする
  36.         x <- sqrt((A[1]-k$cx)^2+(A[2]-k$cy)^2)*cos(θ)
  37.         R <- CurvatureRadius(k$a, k$b, x) # 曲率半径
  38.     }
  39.     list(R=R, ellipse=k, A=A, x=x, S=S)
  40. }

緯度0° から 90° を考えます。方位角も0° から 90° を考えます。方位角による曲率の変化は、この範囲が分かれば、後は東西、南北に折り返して同じです。

方位角を15°刻みに計算すると以下のようになります。

  1. r$> a <- 6378137
  2.     t <- seq(0,pi/2,length.out=7)
  3.     φ <- 36/180*pi
  4.     for(az in t){
  5.         R <- AzimuthCurvatureRadius(φ, az)$R
  6.         cat(az/pi*180,R*a,"\n")
  7.     }
  8. 0 6357482.437549671 
  9. 15 6359353.277720708
  10. 30 6364470.126018431 
  11. 45 6371473.192120609 
  12. 60 6378491.686673882 
  13. 75 6383639.391985905
  14. 90 6385525.660720106

下図の2つの平面は、XZ平面と、点Aにおける方位角45°の法面です。XZ平面は点Aを通る経線の作る平面と同じです。

点Aにおける方位角45°の法面には、垂直截線を青で描きました。オレンジ色の楕円は垂直截線をXY平面まで回転したものです。

スクリプトを説明すると、更成緯度βを求めて、緯度φに相当する地点AのXYZ座標A、法線とZ軸の交点Sを求めます。

Normal は、直線ASです。

点Aを通る経線はXZ平面上にあります。XZ平面上の点(0,0,1)を、方位角azだけ、直線ASの周りに回転します。

その点と、点A、Sの3点で決まる平面PSは、方位角azの法面で、この面による断面が曲率を与える垂直截線です。

垂直截線は、回転楕円体の式と、法面の式からzを消去して、XY平面に投影した楕円の式を得ることで描きます。したがって、法面が赤道面に垂直の場合は計算できません。XY平面に投影した楕円の式を満たす点列を生成して、法面の式からzを計算することで、垂直截線を法面に描くことができます。

qcは、XY平面に投影した楕円の式の係数で、xyは、その楕円の標準形を求めて発生させた楕円の周上の点列です。5点あれば良いのですが、0から360°を等分して、最初は0°、最後は360°となるので、最後を除いて使います。

ILは、法面と赤道面の交線です。wは、法面と赤道面の成す角です。

xyを-w回転して、XY平面に重ねます。この際、点Aの座標も回転します。正しければ、回転結果のzが全てゼロになっているはずです。

XY平面まで回転した方位角45°の垂直截線を描くと左図のようになります。楕円の曲率半径は、式1.1によるので、標準形の楕円の長軸上の長さxから計算します。

準拠楕円体は赤道半径と扁平率で与えられた回転楕円体です。赤道は真円です。経線は全て同じ楕円を描きます。ある地点の曲率は緯度だけで決まります。

経線が描く楕円の長軸半径は赤道半径です。短軸半径は極半径で、21384.6859 m だけ赤道より短くなっています。

曲率半径が最も長いのは極で、6399593.6259 m と算出されます。曲率半径が最も短いのは赤道上の南北方向で、6335439.3271 m となります。差は、64154.2987 m です。

平面直角座標9系の原点の緯度36°における方位角と曲率半径は左図のようになります。

東西の卯酉線曲率半径は、6385525.6607 m で最も長く、南北方向の子午線曲率半径が最も短い、6357482.4375 m となります。差は、28043.2232 m です。

平均曲率半径は、2つの値の幾何平均で、6371488.6206 m です。

この2つの値の単純平均は、6371504.0491 m で、15°刻みの7方向の単純平均は、6371490.8247 m となります。1° 刻みの単純平均は、6371488.7901 m でした。0.1°、0.01° 刻みを計算すると、それぞれ、6371488.6377 m、6371488.6223 m となって、卯酉線曲率半径と子午線曲率半径の幾何平均に近づきます。


題目一覧へmikeo_410@hotmail.com(updated: 2022/09/08)