７．曲率半径と平均曲率

曲率半径と平均曲率

1. 楕円の曲率半径

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

二階微分可能な曲線 $y=f\left({x}\right)$ の点 $A\left({q,f\left({q}\right)}\right)$ における曲率半径 R は、

$R=\frac{{\left({1+{f'\left({q}\right)}^{2}}\right)}^{\frac{3}{2}}}{\left|{f''\left({q}\right)}\right|}$

で、 $f''\left({a}\right)=0$ のときは、曲率半径 ∞ です。

楕円を、

$\frac{{x}^{2}}{{a}^{2}}+\frac{{y}^{2}}{{b}^{2}}=1$

と表すと、

$y=\pm{}\frac{b\sqrt{{a}^{2}-{x}^{2}}}{a}$

maximaで微分して、

$y'=\pm{}\frac{{bx}^{\ }}{a\sqrt{{a}^{2}-{x}^{2}}}$

$y''=\pm{}\frac{ab}{\left({{a}^{2}-{x}^{2}}\right)\sqrt{{a}^{2}-{x}^{2}}}$

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

${R}^{2}=\frac{{\left({1-{\left({y'}\right)}^{2}}\right)}^{3}}{{\left({y''}\right)}^{2}}=\frac{{\left({\left({{b}^{2}-{a}^{2}}\right){x}^{2}+{a}^{4}}\right)}^{3}}{{a}^{8}{b}^{2}}$ 1.1

2. 平均曲率半径

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

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

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

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

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

測量法に基づいた「作業規程の準則」の付則「付録６計算式集」には、平均曲率半径として、 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 です。極以外の地点では方位角によって曲率が変わることになります。

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

つまり、方位角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)

と計算されます。

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

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

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

あるいは、

r$> u <- func(φ, 0)
u$x*a
[1] 5165998.77753988

3.1. 経線の曲率半径

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

3.2. 卯酉線の曲率半径

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

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

$\left({{A}_{x}{\ A}_{y}\ {A}_{z}}\right)=\left({5165998.777539880,\ 0,\ 3728191.675729482}\right)$

となります。

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

a = 6378137

$\left({{A}_{x}{\ A}_{y}\ {A}_{z}}\right)=\left({5165998.777539880/a,\ 0/a,\ 3728191.675729482/a}\right)$

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

$s=b\cdot{}sin\beta{}-cos\beta{}\cdot{}tan\varphi{}$

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

卯酉線が作る平面は法線を含み、 $\left({{A}_{x}{\ A}_{y}\ {A}_{z}}\right)\$ 、（0, 0, s）を含みます。また、XZ平面にある経線に直交するので、Y軸に平行です。したがって、任意のyに対して（0, y, s）を通ります。

$\left({{A}_{x}{\ A}_{y}\ {A}_{z}}\right)\$ 、（0, 0, s）、（0, 1, s）の3点で決まる卯酉線が作る平面を、

とします。

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

回転楕円体は、

${x}^{2}+{y}^{2}+\frac{{z}^{2}}{{b}^{2}}=1$

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

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

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

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

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

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

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

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

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

r$> x<-cos(β)-(b*sin(β))/tan(φ)
y<-sqrt((1-cos(β)^2)*tan(φ)^2+2*b*cos(β)*sin(β)*tan(φ)-b^2*sin(β)^2)/tan(φ)
y2 <- sqrt(-cos(β)^2*tan(φ)^2+2*b*cos(β)*sin(β)*tan(φ)-b^2*sin(β)^2+b^2)/b
x2 <- 0
x <- c(Ax, x, x, x2, 0)
y <- c(0, y, -y, y2, -y2)
z <- c(Az, 0, 0, s, s)
u <- Rotate(L, φ, list(x=x,y=y,z=z))
u
$x
[1] 0.999878441669715401 0.005422141138500591 0.005422141138500591
[4] -0.001279993892284230 -0.001279993892284230
$y
[1] 0.0000000000000000 0.9999853000846934 -0.9999853000846934
[4] 0.9999921881746419 -0.9999921881746419
$z
[1] -6.808789643208968e-17 0.000000000000000e+00 0.000000000000000e+00
[4] 1.604619215278547e-17 1.604619215278547e-17

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

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

r$> v <- EllipseCoeff(u$x,u$y)
S <- StandardEllipse(v$A,v$B,v$C,v$D,v$E,v$F)
S
$a
[1] 0.9999948989782719
$b
[1] 0.9988327146454221
$rotation
[1] 1.570796326794897
$cx
[1] 0.001045727024292981
$cy
[1] -3.420291656950161e-17

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

r$> R <- CurvatureRadius(S$a, S$b, 0)
R
a <- 6378137
R*a
[1] 1.001158435562008
[1] 6385525.660720157

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

3.3. 平均曲率半径

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

6371488.620 m

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

6371504.049 m

です。

3.4. 公式

$f=\frac{1}{298.257222101}$

$a=6378137\ ,\ \ b=a\left({1-f}\right)$

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

${R}_{N}=\frac{{\left({ab}\right)}^{2}}{{\left({{\left({a\cdot{}cos\varphi{}}\right)}^{2}+{\left({b\cdot{}sin\varphi{}}\right)}^{2}}\right)}^{\frac{3}{2}}}=6357482.4375$

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

${R}_{E}=\frac{{{a}^{2}}^{\ }}{\sqrt{{\left({a\cdot{}cos\varphi{}}\right)}^{2}+{\left({b\cdot{}sin\varphi{}}\right)}^{2}}}=6385525.6607$

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

${R}_{\alpha{}}=\frac{{R}_{N}\cdot{}{R}_{E}}{{R}_{N}{sin}^{2}\alpha{}+{R}_{E}{cos}^{2}\alpha{}}$

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

${R}_{45}=\frac{{R}_{N}\cdot{}{R}_{E}}{\frac{{R}_{N}}{2}+\frac{{R}_{E}}{2}}=6371473.1921$

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

$f=\frac{1}{298.257222101}$

$e=\sqrt{2f-{f}^{2}}$

$W=\sqrt{1-{e}^{2}\cdot{}{sin}^{2}\varphi{}}$

$M=\frac{a\left({1-{e}^{2}}\right)}{{W}^{3}}=6357482.4375$

$N=\frac{a}{W}=6385525.6607$

$R=\frac{a(1-f)}{{W}^{2}}=6371488.6206$

4. 全方位の曲率半径

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

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

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

点Aの座標は、

$\left({{A}_{x}{\ A}_{y}\ {A}_{z}}\right)=\left({cos\beta{},\ 0,\ b\cdot{}sin\beta{}}\right)$

$tan\beta{}=b\cdot{}tan\varphi{}$

$f=\frac{1}{298.257222101}$

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

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

$tan\varphi{}=\frac{\overline{AT}}{\overline{ST}}=\frac{\overline{AR}+\overline{RT}}{\overline{ST}}=\frac{b\cdot{}sin\beta{}-s}{cos\beta{}}$

$s=b\cdot{}sin\beta{}-cos\beta{}\cdot{}tan\varphi{}$

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

$\left({{Ax,\ 0,\ Az-s}_{\ }}\right)=\left({cos\beta{},\ 0,\ cos\beta{}\cdot{}tan\varphi{}}\right)$

$=\left({cos\beta{},\ 0,\ \frac{sin\beta{}}{b}}\right)$

です。

$cos\beta{}\cdot{}tan\varphi{}$ は、 $tan\beta{}=b\cdot{}tan\varphi{}$ なので、

$cos\beta{}\cdot{}tan\varphi{}=cos\beta{}\frac{tan\beta{}}{b}=cos\beta{}\frac{\frac{sin\beta{}}{cos\beta{}}}{b}=\frac{sin\beta{}}{b}$

です。したがって、

$s=b\cdot{}sin\beta{}-cos\beta{}\cdot{}tan\varphi{}=b\cdot{}sin\beta{}-\frac{sin\beta{}}{b}=\left({b-\frac{1}{b}}\right)sin\beta{}$

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

$y-{Y}_{q}=-\frac{x-{X}_{q}}{{f}^{'}\left({{X}_{q}}\right)}$

のようです。

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

楕円、

${x}^{2}+\frac{{z}^{2}}{{b}^{2}}=1$

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

$z=\pm{}b\sqrt{1-{x}^{2}}$

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

$z'=-\frac{{bx}^{\ }}{\sqrt{1-{x}^{2}}}$

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

$z-{A}_{z}=-\frac{x-{A}_{x}}{-\frac{{b{A}_{x}}^{\ }}{\sqrt{1-{{A}_{x}}^{2}}}}$

$z=\frac{\sqrt{1-{{A}_{x}}^{2}}}{b{A}_{x}}x+{A}_{z}-\frac{\sqrt{1-{{A}_{x}}^{2}}}{b}$

となります。切片は、

${A}_{z}-\frac{\sqrt{1-{{A}_{x}}^{2}}}{b}=b\cdot{}sin\beta{}-\frac{\sqrt{1-{cos}^{2}\beta{}}}{b}=b\cdot{}sin\beta{}-\frac{sin\beta{}}{b}=\left({b-\frac{1}{b}}\right)sin\beta{}$

となります。

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

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

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

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

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

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

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

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

r$> a <- 6378137
t <- seq(0,pi/2,length.out=7)
φ <- 36/180*pi
for(az in t){
R <- AzimuthCurvatureRadius(φ, az)$R
cat(az/pi*180,R*a,"\n")
}
0 6357482.437549671
15 6359353.277720708
30 6364470.126018431
45 6371473.192120609
60 6378491.686673882
75 6383639.391985905
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 となって、卯酉線曲率半径と子午線曲率半径の幾何平均に近づきます。