从水电工程图中精确量算面积的方法（何永明,朱贤博）
|idpi
关 键 词: 三次样条函数法;加权平差方法;投影交换;工程图纸
中图分类号: TU20 文献标识码: A
1 精确量算地图的面积
在水利水电工程施工中，经常需要在地图上量取各种地理信息，如面积、坐标等，但任何地图都不可避免地存在因地图投影而产生的变形。除等积投影外，在图上量算出的图面面积与主比例尺分母平方的乘积并不等于相应的椭球体表面上的面积，而椭球体表面上规则区域的面积是可积的，如某种梯形。基于这种情况，在地图上量算面积，可通过如下方法实现。
（1）运用地图投影变换的理论，把通过数字化仪采集的图上所求区域的边界曲线各点的平面直角坐标x，y反算为大地坐标φ，λ。
（2）运用三次样条函数法，可以根据实地边界曲线的变化来选取支距点，采用符合边界曲线更的数学模型，应用于内外业面积量算中，能提高面积的计算精度。
（3）应用测量中的平差方法，当所量各子区域满布一个由经纬线所构成的大区域时，则对面积量算误差进行以子区域面积为权的加权平差。这样即可使总区域大小不变，又可使子区域量算精度提高一个数量级。
2 几个基本原理的应用过程
2.1 地图投影反解变换
已知地图投影函数的一般公式是:
x=f1 （φ，λ）
y=f2 （φ，λ） （1）
则其反函数式可写作:
φ=f1-1 （x，y）
λ=f2-1 （x，y） （2）
当式（1）确定时，式（2）也随之确定。其中x，y为图上坐标;φ，λ为大地纬度、经度;f 1 、f 2 为单值、连续函数。
如对默卡托投影，则正算式为:
x=r0/c0（q-q 1 ）
y=r0/c0（λ-λ 1 ） （3）
q=th -1 sinφ-e th -1 （esinφ） （4）
r=acosφ/（1-e 2 sin 2 φ）
q=（xc 0 /r 0 ）+q 1 （5）
λ=（yc 0 /r 0 ）+λ 1
ψ=sin -1 thq （6）
φ=ψ+0.192288sin 2 ψ+0.000376sin 4 ψ+…
式（3）～（6）中，q为等量纬度;r 0 为基准纬圈半径;a为椭球体长半轴;e为椭球体第一偏心率。
一般地图投影的解算中，对已知投影的反算式多有介绍。但对于未知投影的地图，或虽知投影，但由式（2）不易求得的投影，只要图上有精确的纬线网，可通过数值变换的方法，将其上各点反算出来。如变换多项式:
x=f 1-1 （φ，λ）=a 0 +a 1 φ+a 2 λ+a 3 φ 2 +a 4 φλ+a 5 λ 2 +a 6 φ 2 λ+a 7 φλ 2 +a 8 φ+a 9 λ 3 （7）
y=f 2-1 （φ，λ）=b 0 +b 1 φ+b 2 λ+b 3 φ 2 +b 4 φλ+b 5 λ 2 +b 6 φ 3 +b 7 φ 2 λ+b 8 φλ 2 +b 9 λ 3
可先在图上均匀采集10个经纬线交点（φ，λ）及相应的（x，y）值，组成线性方程组，用主元消去法或最小二乘逼近来求解各系数ai，bi。
待求面积的图上边界曲线各点，用手工方法量取其图上坐标，并反解经纬度，工作量是非常大的，采用以计算机控制的数字化仪时，由于数字化仪坐标系与地图平面直角坐标系的不一致，加之纸张的变形，必须通过一种称之为仿射变换的方法，对所取点进行旋转、平移和伸缩的变换处理，才能得到图上正确的坐标。
2.2 三次样条函数法求面积
用三次样条函数对不规则曲线进行拟合，可使模拟曲线光滑至有连续二阶导数，即任意两相邻多项式的峰值点得到了圆滑处理，因而其数学模型能更好地与地面曲线相符。
对每一曲线段采用三次多项式Pi（X）使之在区间［ Xi，Xi+1］内，而每一个三次多项式Pi（X）在其两端点Xi及Xi+1处必须和已知函数在该点的值相等，则:
Pi（X）=（fi/hi-hici/6）（Xi+1-X）+（fi/hi-hici+1/6）（X-Xi）+ci/6hi （Xi+1-X）3+ci+1/6hi（X-Xi）3
（i=1，2，…，n-1）
根据积分原理可得其面积为
Ai =∫Pi（X）dx=∫［（fi/hi-hici/6）（Xi+1-X）+（fi/hi-hici+1/6）（X-Xi ）+ci/6hi（Xi+1-X）3+ci+1/6hi（X-Xi）3 ］dx=0.5（fi/hi-hici/6）（Xi+1-X）2+0.5（fi+1/hi-hici+1/6）（Xi+1-Xi）2+ci/24hi（Xi+1-Xi）4+ci+1/24hi（Xi+1-Xi）4
将hi=Xi+1-Xi代入并整理得:
Ai=［（fi+fi+1）/2-h2i（ci+ci+1）/12］ki+k3i（ci+ci+1）/24
（i=1，2，…，n-1） （8）
式（7）即为曲线与X轴在区间［ Xi，Xi+1］上所围面积的计算公式。显然，曲线与X轴在区间［X1，Xn］上所围面积的计算公式为:
A=∫Pi（X）dx=∫P1（X）dx+∫P2（X）dx+…+ ∫Pn-1（X）dx=A1+A2+…+An-1 =∑Ai （9）
以上即为应用三次样条函数法计算不规则图形的面积计算公式。下面给出应用三次样条函数法在计算机上计算不规则图形面积的处理步骤:
（1）输入支距点数n，数组说明，输入支距点坐标Xi，Fi，计算Hi=Xi+1-Xi;
（2）按hi-1ci-1+2（hi+hi-1）+hici+1=6［（fi+1-fi）/hi-（fi-fi-1）hi-1］（i=2，3，…，n-1）且c1=cn=0，组成求算三次样条多项式系数ci的线性方程组，用矩阵形式表示为M×C=W;
（3）解算线性方程组求三次样条多项式系数C=M -1 W;
（4）按（8）式和（9）式计算图形的面积A=∑｛［0.5（fi+fi+1）-h2i（ci+ci+1）/12］hi+h3i（ci+ci+1）/24｝。
2.3 椭球面上梯形面积的计算
椭球面上的梯形面积是指由两条纬线和两条经线所围成区域的面积，可用积分方法导出，推导过程如下:
已知沿纬线微分线段dsn=rdλ，沿经线微分线段dsm=Mdφ，则椭球面上微分梯形面积可得:
dT=Mrdφdλ=b2cosφ/（1-e2sin2φ）2dφdλ （10）
对上式取定积分有:
T= ∫∫b2cosφ/（1-e2sin2φ）2dφdλ （11）
积分后得:
T=b2（λ2-λ1）［sinφ2/（1-e2sin2φ2）+th-1（esinφ2）-sinφ1/（1-e2sin2φ1）-th-1（esinφ1）］/2 （12）
式中b为椭球体短半轴。
2.4 面积加权平差的方法
由于在地图上绘制，数字化仪采点，以等高椭球面梯形代替不等高椭球面梯形的过程中，不同程度地存在着误差，因而使量算出的面积结果也必然含有一定误差。对这些误差不进行处理，必然降低量算精度。因此对量算结果进行平差处理是提高精度的一个重要途径。椭球面上梯形，在地图上是由经纬线所围成的区域，其真实面积可由式（12）求得。若所量各子区域满布于这块椭球面梯形区域内，且量算无误差，则各子区域面积和应等于这块椭球面梯形面积。设总区域范围为######φ1-φ2，λ1-λ2，其真实面积为T，其内各子区域Di的面积量得为Si，其和为 ∑Si，显然由于量算误差的存在，T≠∑Si，其差为:
ΔS= ∑Si-T
为此，可以采用测量中的加权平差方法，以面积为权，把这些误差分配到各子区域中，从而修正面积量算结果为较真实的结果。加权平差后各子区域面积可用正式表达:
S′i =Si×T/∑Si
其满足 ∑S′i=T。
2.5 精度分析
精度分析分两个方面:①经纬度转换误差;②面积量算误差。经测试，经纬度转换误差在图上达到不大于0.4mm。其引起的原因是被量地图的绘制印刷误差和数字化取点误差，设其分别为0.3mm和0.2mm，按误差传播定律可得:
δd=sqrt（0.32+0.22）=0.36（mm）
因而证明本程序算法正确，误差符合要求。
3 结语
本文中所介绍的面积量算算法正确，精度较高，求得的面积是可靠的。但本文所讨论的误差，只是就图而言。至于图上所量区域的边界线，描绘是否精确，概括是否合理，是不顾及的。因而若想进一步提高地表某一区域的面积量算精度，还需从测绘边界线入手，量算时尽量选用较大比例尺且误差较小的地图。
作者简介: 何永明，男，武警水电第二总队第七支队，助理工程师。