利用测地坐标计算椭球面上凸多边形面积的算法研究

2 利用测地坐标计算椭球面三角形的面积
2
中国科技论文在线 ___________________________________________________________________________ http://www.paper.edu.cn
P3
(s x
3
,sy3
)
F
P2
(s x
2
,sy2
F=
1 m 1 s yi + s yi +1 s xi − s xi +1 + ∑ 2 i =1 24 R02
(
)(
)
∑ s 3yi s xi −2 − s xi −1 +
i =1
m
(
)
1 8R02
∑s
i =1
m
yi
s yi −1 s yi + s yi −1 s xi − s xi −1
（8）
(
)(
F = s y3 s x3 − s x1 −
(
)
s3 y3 6 R02
(s
x3
− s x1
)
⎡ + ⎢ s y1 − s y3 s x3 − s x1 ⎢ ⎣
(
)(
)
2 ⎛ ⎛ ⎞⎞ s y ⎜ ⎜1 − 1 ⎟ ⎟ 3 ⎤ ⎜ ⎜ 2 R02 ⎟ ⎟ s3 y − s y3 s x3 − s x1 ⎠⎟ 1 ⎥ ÷ ⎜1 + ⎝ − 2 2 6 R0 ⎛ ⎞ ⎥ ⎦ ⎜ ⎜1 − s y3 ⎟ ⎟ ⎜ ⎜ ⎜ 2 R02 ⎟ ⎟ ⎟ ⎠⎠ ⎝ ⎝
(
)
(
)⎤
⎥ds x = s y2 − s y1 s x2 − s x1 − ⎥ ⎦
(
)(
)
(s
3 y2
− s3 y1
2 0
)
6R
(s
x2
− s x1
) (5)
上式中, n 及 R 只取至二次项,这在一定范围内已能满足精度要求。由式（5）不难看出，利用测地坐标 依测地格网来计算椭球面上的面积要比利用大地经纬度依大地格网进行计算要简单方便得多。
(
)(
)
)(
)
( (
) ( ) (
)(
)
)
(
)
( − (s
1 2 2 2 3 2 [ s3 y1 + s y1 s y3 − 2 s y1 s y3 + s y1 s y3 + s y3 − 2 s y1 s y3 s x3 − s x1 8 R02
3 y2 2 2 2 2 s y3 + s y s y3 + s 3 + s y2 s y − 2s y y3 − 2 s y 2 s y 3 3 2 2
∂r ∂r ds x ds y , PL = PQ = ∂s x ∂s y
→ ⎛ → ⎞ ∂r ⎜ ∂r ⎟ dσ = ⎜ ds y × ds x ⎟ ∂s x ⎜ ∂s y ⎟ ⎝ ⎠
→
→
(1)
微分四边形 PQTL 的面积可用来代表该测地格网在椭球面上的微分面积 dσ ，即有
(2)
椭球面上由坐标曲线 u1 = s y1 , u 2 = s y2 v1 = s x1 , v 2 = s x2 所构成的测地格网的四边形 D 的面积为
目（03-04-02） 作者简介：施一民(1942), 男 ,浙江宁波人, 教授，博士生导师。.E-mail：yimshi@citiz.net 1
中国科技论文在线 ___________________________________________________________________________ http://www.paper.edu.cn
∂r ∂r × ∂s y ∂s x
→ →
→
→
2
∂r ∂r ∂r ∂r =( × )⋅( × ) ∂s y ∂s x ∂s y ∂s x
→ → → → → →
→
→
→
→
∂r ∂r ∂r ∂r ∂r ∂r ∂r ∂r )⋅( )−( )⋅( ) = EG − F 2 = n 2 ⋅ ⋅ ⋅ ⋅ =( ∂s x ∂s y ∂s y ∂s x ∂s x ∂s x ∂s y ∂s y
___________________________________________________________________________ http://www.paper.edu.cn
利用测地坐标计算椭球面上凸多边形面积的算法研究
施一民 1 , 朱紫阳 2
(同济大学 测量与国土信息工程系,上海,200092；2.广东省国土资源厅测绘院, 广州,510500)
由各顶点的平面直角坐标计算凸多边形面积已有简易可行的熟知公式，然而算得的高斯平面上 的面积与所需的椭球面上的面积相比，存在较可观的高斯投影变形，尤其是当所求面积与中央子午线
相距较远时。若能在与测区平均高程面相接近的区域性椭球面 上计算图形面积，则可因杜绝高斯投 影变形而提高面积量算的实际精度，然而按现有方法计算甚是不易：对于由大地经纬线所包围的曲边 梯形尚可由积分求得；对于椭球面三角形还有公式可解；但对于由大地线所构成的任意凸多边形就需 分割成若干个椭球面三角形，分别求面积再取总和。如果已知各顶点的大地经纬度，为了计算图形面 积还需涉及一系列的大地主题反解。实际工作中是难以进行这样繁复计算的。本文正是要探求如何用 文献［2］ 、 ［3］ 、 ［4］所提出的测地坐标来方便地计算椭球面上凸多边形的面积。
(
)(
)
(
)(
)
[ (
)
(
)
(
)]
+
1 s y s y s y + s y3 s x1 − s x3 + s y2 s y1 s y2 + s y1 s x2 − s x1 + s y3 s y2 s y3 + s y2 s x3 − s x2 8R02 1 3 1
[
(
)(
)
(
)(
)
(
[7] 由式（7）可知,利用三顶点的测地坐标所算得的椭球面三角形面积可分为主项 F0 以及含 小项 δF 两部分。主项 F0 的表示式完全等同于用平面坐标计算面积的常用公式。
L
r ′ d sx
T
dσ
r ( s x , s y ) d sy P
Q
(B0 , L 0 )
P0
图 1：椭球面上基于测地格网的微分面积
Fig.1 Differential area based on geodesic Grid on ellipsoidal surface
由向量运算中的 Lagrange 恒等式得到：
摘 要：本文推导出用三顶点的测地坐标计算地球椭球面上三角形面积的公式，公式表明，其 主项的表示式与按平面坐标求面积的计算式完全一致，而附加项的表示式亦有规律可循，因此 该公式的适用范围可由椭球面三角形推广至椭球面上任意凸多边形。与高斯平面上计算的面积 相比，由于它不再蒙受投影变形的影响，更接近于实际面积,这就为至今难以实施的椭球面上面 积计算开拓了一个新的途径, 有利于作出更客观的 GIS 空间量度与分析。实际数据的验算充分 证实了该算法的正确性和有效性。 关键词：测地坐标 椭球面三角形 凸多边形 面积公式
3 利用测地坐标计算椭球面上凸多边形的面积
式（7）完全可以推广到椭球面上凸多边形的计算，设为任意凸 m（m≥3）边形 P1P2 P3 LPm ， 若各顶点由 1 至 m 按顺时针排列，则面积计算公式为:
4
中国科技论文在线 ___________________________________________________________________________ http://www.paper.edu.cn
(
)
(
)
(
)
( )( ) )( )]
1 的 R02
24 R
2 0
1 1 1 s y1 + s y2 s x1 − s x2 + s y2 + s y3 s x2 − s x3 + s y3 + s y1 s x3 − s x1 2 2 2 1 3 3 s3 + y s x − s x2 + s y 2 s x1 − s x3 + s y3 s x2 − s x1 24 R02 1 3
− s y2 s x2 − s x1 +
(
)
s3 y2 6R
2 0
(s
x2
⎡ − s x1 − ⎢ s y1 − s y2 s x2 − s x1 ⎢ ⎣
) (
)(
)
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
− s y2 s x3 − s x2 +
(
)
s
3 y2
6 R02
(s
x3
⎡ − s x2 − ⎢ s y3 − s y2 s x3 − s x2 ⎢ ⎣
)
P1 s x 1 , s y1
(
)
(B0 , L 0 )
P0
图 2：测地坐标系中的椭球面三角形 Fig.2 Ellipsoidal triangle in geodesic coordinate system
设椭球面上由三条大地线所围成的三角形 P 1 P2 P 3 的 三 顶 点 的 测 地 坐 标 分 别 为
) (
)(
)
(6) 对式(6)进一步整理化简，得出
F = F0 + δF
= s y3 s x3 − s x1 + − s y2 s x3 − s x2 +
(
)
(
) (
1 1 s y1 − s y3 s x3 − s x1 − s y2 s x2 − s x1 − s y1 − s y2 s x2 − s x1 2 2 3 3 s y1 s x2 − s x3 + s 3 1 y 2 s x3 − s x1 + s y3 s x1 − s x2 − s y3 − s y2 s x3 − s x2 + 2 12 R02
→ ⎛ → ⎜ ∂r ∂r × A = ∫∫ dσ = ∫∫ ⎜ ⎜ ∂s y ∂s x D ⎝
⎞ ⎟ ds y ds x ⎟ ⎟ ⎠
（3）
---------------------------------------------------------------------------------------------------------------基金项目：国家自然科学基金资助项目(40471114) ,地球空间环境与大地测量教育部重点实验室开放基金资助项
(s
x1
, s y1 , s x2 , s y2 , s x3 , s y3 。过 P1 , P2 , P3 三点分别作测地格网线（图 2）, 若把两内角均为直角的
) (
) (
)
三个曲边“梯形” 各分成曲边“矩形” 和“直角三角形” 两部分，而在“直角三角形” 面积的计 算中，将其看成所相应的曲边“矩形” 的一半左右。因这两个曲边“直角三角形” 仅有一条边不等 长，故其面积之比应为长度归化因子之比。则由式（5）可得出该椭球面三角形的面积为
由（3） 、 （4）两式，得出
(4)
A = ∫ EG ds y ds x = ∫ nds y ds x = ∫
D D
s x2
s x1
∫
s y2
s y1
⎛ s2 ⎜1 − y ⎜ 2R 2 0 ⎝
⎞ ⎟ds y ds x ⎟ ⎠
=∫
s x2
s x1
3 ⎡ s3 y 2 − s y1 ⎢ s y2 − s y1 − 6 R02 ⎢ ⎣
[1]
1 利用测地坐标计算椭球面上测地格网曲边“矩形”的面积
按一定间隔的一族大地线及一族测地平行线构成了地球椭球面上的正交参数曲线格网（测地格 网）
[2 ]
,借此可度量椭球面上由两对正交的坐标曲线所围成的四内角均为直角的四边形面积。在图 1
所示的测地格网的一个微分网格 PQTL 中，当 Q，L 点无限接近于 P 点，则有
(
)(
)
2 2 2 3 2 − s3 y1 + s y1 s y 2 − 2 s y1 s y 2 + s y1 s y 2 + s y 2 − 2 s y1 s y 2 s x2 − s x1 x3
)( )(s
− s x2
) )]
+{
=
3 3 s3 y1 s x3 − s x2 + s y 2 s x1 − s x3 + s y3 s x2 − s x1
(
)(
Baidu Nhomakorabea
)
3
中国科技论文在线 ___________________________________________________________________________ http://www.paper.edu.cn
2 ⎡ sy 1 ⎢ 1− 3 2 ⎤ − s3 s 2 R y y 0 − 1 2 2 ⎥ ÷ ⎢1 + 2 6 R0 ⎦ sy ⎥ ⎢ ⎢ 1 − 22 2 R0 ⎢ ⎣ 2 ⎡ sy 3 − 1 ⎢ 3 ⎤ − s3 s 2 R02 y y − 3 2 2 ⎥ ÷ ⎢1 + 2 ⎢ 6 R0 ⎥ ⎦ ⎢ 1 − s y2 2 R02 ⎢ ⎣
)
若将上式分成主项及次要项两部分，则有
F0 =
1 m ∑ s y + s yi +1 s xi − s xi +1 2 i =1 i