多波束数据处理中的UTM直角坐标与经纬度坐标转换技术.doc

第 24 卷第 2 期2006 年 4 月海 洋 科 学 进 展ADV A N C ES IN MA R IN E SC I EN C EVol . 24 No . 2Ap ril ,2006多波束数据处理中的 UTM直角坐标3与经纬度坐标转换技术熊应乾1 ,2,杜德文1 ,姜莉莉3 ,杨刚1,吕海龙1 ,闫仕娟1(1 . 国家海洋局 第一海洋研究所 ,山东 青岛 266061 ;2 . 中国海洋大学 ,山东 青岛 266071 ;3 . 中国科学院 地理科学与资源研究所地图室 , 北京 100101)摘 要 :简要介绍了地图投影的一般概念 ,重点描述了 U TM 的优点及其与大地坐标之间的相互转换 ,该转换方法软件化后形成的技术可以为多波束数据后处理商业软件提供坐标转换的补丁技术 。最后附上部分经实际检验的Fo rt ra n 程序供参考 。关键词 :多波束测量 ;通用横轴墨卡托投影 ( U TM) ;大地坐标 ;坐标转换中图分类号 : P716文献标识码 : A文章编号 :167126647 (2006) 0220246208多波束测深在当代各国的科研 、航海 、河道港口的规划整治及军事等方面发挥了越来越大的作用 。然而由于多波束勘测仪器型号及后处理软件千差万别 ,各种数据格式及椭球/ 投影方式的采用使数据的兼容和整 合存在困难 ,从而影响了对多波束数据的综合应用 ,这个问题在国内外已经引起了广泛的关注 , GM T ( ge2 ne ric mappi ng tool s) 及 MB Cha r t 1 出现的部分原因即在于此 。在我国使用较多的多波束数据采集及后处理软件包括 Ko ng sber g Si mra d 公司的 Nep t une 软件包 ( 采 集 : Me r maid/ Me rli n ,后处理 : Nep t u ne 和 Cf loo r Ro xa r 公司) 和 Trito n Elic s 公司的 Dat a Fu sio n ( 包括 I2 si s , Bat hyp ro , Delp h map , Delp h Sei smic , S GIS) 等 。Nep t une 主要基于工作站使用 ,整个系统比较成熟稳 定 ,处理能力强 ,但软硬件费用高昂 ,U M X 语言操作难于为其他科学工作者学习和使用 ,在功能上也有所欠 缺 ,如不能提供网格经纬度数据 (其最新版本已出现微机版但其性能有待改进 ,如数据量大时由于内存使用 效率过低无法计算 ,在 Po sei do n 中作后散射图批处理时同一点在不同批次的数据块中灰度值不完全统一) ; Dat a Fu sio n 为微机版 ( Wi ndo w s 操作系统) ,适用于多波束 、浅地层剖面仪 、旁侧声纳的数据采集及后处理 , 其室内工作时所需的软硬件性价比大为提高 ,且操作平台和软件本身相对更易于使用 ,适用范围更广 ,但也 存在缺点 ,目前还不能提供原始经纬度数据 ,也不能提供从直角坐标转换到经纬度的功能 ,给实际工作带来 极大不便 ,这也是促成本文的原因之一 。本文提出一个简单快捷 、切实易行的方法实现多波束测深数据的兼容和整合 。由于多波束数据采集时 需运用 GPS 导航并提供中央波束的经纬度数据 ,而 GPS 使用 W GS84 椭球 ; U TM 是一种精度高 、变形幅度小而分布均匀 、适用范围广的投影方式 ,因此本文主要基于 W GS84/ U TM 考虑多波束数据的兼容和整合 。先由各种后处理软件直接获得原始离散点经纬度数据 ,若不能输出经纬度数据则可通过输出墨卡托投影数据后转换成经纬度数据 (我国多采用墨卡托投影 ,与经纬度数据转换的方法及转换公式在诸多文献书刊中都 有出现 2 ) ,也可输出 U TM 直角坐标数据后转换成经纬度数据 ( 与经纬度数据的转换在本文详述) ,然后再 将需要转换的经纬度数据转换成 U TM 直角坐标数据而实现不同航次 、不同仪器和软件采集处理的多波束 数据的兼容和整合 。3 收稿日期 : 2005205220基金项目 : 国家高技术发展研究规划项目 海底底质分类的多参数识别技术 ( 2003 A A613040)作者简介 : 熊应乾 ( 19732) ,男 ,湖南沅江人 ,助理研究员 ,从事海洋地质及海底信息技术研究 。(杜素兰编辑)1简介地球是一个表面很复杂的球体 ,人们以假想的平均静止的海水面形成的“大地体”为参照 ,求出近似的椭球体 ,理论和实践证明 ,该椭球体近似一个以地球短轴为轴的椭园而旋转的椭球面 ,这个椭球面可用数学公 式表达 ,通过将自然表面上的点归化到这个椭球面上即可 。对应于各椭球建立的大地坐标系按其原点相对地球质心的位置可分为局部坐标系和地心坐标系 。局部 坐标系的原点偏离地心可能达几十到几百米 ,而地心坐标系的原点理论上与地心重合 (实际上与地心难免有 些偏离) 。局部坐标系一般是以经典测量技术为基础建立的 ,而地心坐标系则是以卫星大地测量为基础建立 的 。由于局部坐标系存在一定的局限 ,随着科学技术的发展尤其是 GPS 技术的广泛应用 ,大地坐标系应用 的深度和广度已今非昔比 3 ,而在多波束测深过程中离不开 GP S 定位技术 ,因此选用 W GS 84 这样一个成熟的高精度的地心坐标系 3 ,4 就成为理所当然的事情了 。表 1 列出了几种常用的椭球体参数 。表 1 常用的一些椭球体及参数Ta ble 1 So me f requently used ref erence ellip soids a nd t heir pa ra meter s椭球名称长半轴/ m1/ 扁率适用范围海福特椭球 ( 1910)Kra sso vsky ( 1940)6 378 3882971952 年以前北京 54 坐标系6 378 245298 . 3中国国际大地测量协会推荐椭球 I.U . G. G. ( 1975)6 378 140298 . 257西安 80 坐标系W GS 84GRS 67Kra sso vsky ( 1940) Int er natio nal ( 1924) Cla r ke ( 1880)Cla r ke ( 1866) Be ssel ( 1841)Ai r y ( 1830)6 378 1376 378 1606 378 2456 378 3886 378 2496 378 2066 377 3976 376 542298 . 257298 . 25298 . 3297293294 . 98299 . 15299全球 ( GPS)澳大利亚 1966 ,南美 1969俄罗斯 欧洲法国 , 非洲北美 德国 英国 Eve re st ( 1830) 6 377 276 300 印度注 :此表部分内容引自电子文献 投影 : 地图投影是把地球椭球面上的经纬线网和地理要素表示到平面上的数学法则 。按投影面的形 状 ,地图投影分为方位投影 、圆柱投影 、圆锥投影 ;按投影变形的性质 ,地图投影又可分为等角投影 、等面积投 影 、等距离投影 。表 2 为几种常用等角投影方式一览表 。 ht t p :/ / www . sp envi s. o ma . be/ sp envi s/ help / backgro und/ coo rt ra n . ht ml248海 洋 科 学 进 展24 卷图 1 常用投影系统Frequently used p rojectio n systemsFig. 1(a) Mercato r p rojectio n system , ( b) U niver sal Tra nsver se Mercato r p rojectio n system2U T M 投影许多世界地图所使用的投影是墨卡托投影 一种圆柱投影 ,墨卡托投影的圆柱在赤道与地球相切 ,但实际上并不是一定要求在赤道上与地球相切 ,它也可以在某一经度上与地球相切 ,此时一般称之为横向墨 卡托投影 。由于横向墨卡托投影在较窄的区域内精度很高 , 因而成为了一种称之为通用横轴墨卡托投影 ( U TM) 的全球性坐标体系的基础 。墨卡托投影是一种等角投影 ,地球上的角度及一些小的地形以相同的角度和形状投影到地图上 ,适合制 作航海用图 ,其代价就是远离地图中央的部分变形大 ,如格林兰在采用墨卡托投影的地图上与南美的的面积 相当 ,但实际上只有南美的 1/ 8 。直接用经纬度坐标生成的“真实”地形也常常与实际地形有差异 ,尤其是当 纬度偏高时 (见图 2) 。U TM 将全球分为 60 个投影带 ,每 6 一带 ,每带则使用横向墨卡托投影 ,其中央经线为 6 带中 部的经线 ; 投影带编号 1 60 , 与百万分之一地形 图的有关规定一致 ,适用介于 80 S 到 84 N 的区 域 。南北半球从赤道处起以 8 分带 ( 图 1 b) ,除此 之外的两 极地 区 则 采 用 通 用 极 球 面 投 影 ( U P S) 。 U TM 的椭圆柱面 与地 球的 交 线在 对称 于 距中 央 子午线约 1 40 处 ,投影后两条割线无变形 ,中央 子午线的长度比 m0 小于 1 ,设定 m0 = 0 . 999 6 。U TM 属于横轴等角割圆柱投影 ,因其直角网 格和十进制特性而比经纬度更便于使用 ,但无法独图 2 塔斯玛尼亚地形图Fig. 2 Topo grap hic map of Ta sma nia立地确定 网格 的位 置 。U TM 在同 一纬 线 上 长 度变形随远离中央子午线而增大 ; 在同一经线上长度 ht t p :/ / www . vi sualizatio n sof t wa re . co m/ 2de m/ ut mi nfo . ht ml变形则随纬度的增高而减小 ,在两条交线以内为负变形 ,以外为正变形 ,总体而言改善了变形的分布 ,低纬度的长度变形得到进一步的改善 ,特别是高纬地区的变形比高斯 - 克吕格投影小 2 ,3 ,5 。图 2 中左图运用原始大地坐标由 GTO PO30 软件生成 ,右图则由 3D EM 软件以 U TM 对同一个 D EM (数字地形模型) 重新采样得到 ,它更接近真实的地形 。3U T M 与经纬度的互换由于 W GS84 椭球有精度高和适用范围广等优点 ,而采集多波束数据时需运用 GPS 导航并提供中央波束的经纬度数据 ,故 GPS 使用 W GS84 椭球 ; U TM 是一种精度高 、变形幅度小而分布均匀 、适用范围广的投 影方式 ,因此本文建议使用 W GS84/ U TM 整合多波束数据 ,并在文献 6 的基础上提供了 U TM 与大地坐 标的之间相互转换的公式 。本文对公式有所改进 ,经商业软件 Coo r dt ra n s (提供单点各种形式坐标之间的转换) 检验 ,坐标转换精度提高到分米级 。3 . 1从经纬度坐标转换到 UTM 直角坐标如图 3 所示 , P = 考察点 ; F = 考察点 P 在中央经线上的投影 , F 点的纬度称之为考察点在中央经线上的投影纬度 f p ; O = 原点 ; O Z = 中央经线 ; L P = 过 P 点的纬度线 ; Z P = 过 P 点的子午线 ; OL = 始自赤道k0 S的弧长 ; L F = 曲率的纵坐标 ; O F = N = 网格北向距离 ; F P = E =央经线的网格距离 ; G N = 网格北 ; C = 真北和网格北之间的夹角 。(1) 从经纬度转换到 U TM这些公式由美国军方 (1973 年) 公式略作改动 :lat = 考察点纬度lo ng = 考察点经度lo ng0 = 分带中央经度与中图 3 考察点位置示意图Fig. 3 Schematic diagra m of t he survey point po sitio nk0 =e2 =e2 =n =r hoscale alo ng lo ng0 =0 . 999 6(1 - b2 / a2 ) = 0 . 006 694 379 990 为第一偏心率(ea/ b) 2 = e2 / (1 - e2 ) = 0 . 006 739 496 742 为第二偏心率(a - b) / (a + b)= a (1 - e2 ) / 1 - e2 si n2 (lat ) 3/ 2 为子午圈平面曲率半径nu = a/ 1 - e2 si n2 (lat ) 1/ 2 为与子午圈垂直的平面的曲率半径 ,同时也是考察点到极轴的垂直距离p =(lo ng - lo ng0 )si n1 = pi/ (180 60 60) = 4 . 8 481 368 x 10 - 6 为 1 弧度的正弦值(2) 计算子午线弧长S 过考察点子午圈弧度 ,地球表面自赤道至考察点距离 ,所有角度以弧度表示 。S = Alat - Bsi n (2lat ) + Csi n (4lat ) - Dsi n (6lat )+ Esi n (8lat )A = a 1 - n + (5/ 4) ( n2 - n3 )+ (81/ 64) ( n4 -n5 ). . . - n5 ) . . . = (3a n/ 2) 1 - n + (7/ 8) ( n2- n3 ) + (55/ 64) ( n4BCDE= (15a n2 / 16) 1= (35a n3 / 48) 1+ (3/ 4) ( n2 - n3 ) . . . + (11/ 16) ( n2 - n3 ) . . . -nn= (315a n4 / 51) 1 - n . . . 3 . 2由 UTM 直角坐标转换到经纬度y = 网格纵坐标 ; x = 网格东向横坐标 。相对中央经线而言 ,应从传统 U TM 坐标值减去 500 ,000 。250海洋 科 学 进展24 卷(1) 计算子午线弧长M = y/ k0 .计算考察点的投影纬度 fpmu = M/ a (1e1 = 1 - (1- e2 / 4 -3e4 / 645e6 / 256 . . . )e2 ) 1/ 2 + J 3 si n (6 mu)- e2 ) 1/ 2 / 1+ (1 -fp = mu + J 1 si n (2 mu) + J 2 si n (4 mu)+ J 4 si n (8 mu)J 1J 2J 3J 4=(3e1 / 2-27e13 / 32 . . )= (21e12 / 16 - 55e14 / 32 . . )= (151e13 / 96 . . )=(1097e14 / 512 . . )(2) 计算经纬度e2C1T1R1N1= (ea/ b) 2 = e2 / (1 - e2 )= e2 co s2 (fp )= t a n2 (fp )= a (1 - e2 ) / (1 - e2 si n2 (fp ) ) 3/ 2= a/ (1 - e2 si n2 (lat ) ) 1/ 2D =x/ ( N1 k0 )lat = fp - Q1 ( Q2 - Q3 + Q4)= N1 t a n (fp ) / R1= (D2 / 2)= (5 + 3 T1 + 10C1 - 4 C12 - 9e2 ) D4 / 24Q1Q2Q3Q4= (61 + 90 T1298C1 + 45 T12 - 3C12 - 252e2 ) D6 / 720Q6 + Q7) / co s (fp )+-lo ng = lo ng0 + ( Q5Q5Q6Q7= D= (1= (5C1) D3 / 6+-2 T12 C1+ 28 T1 - 3 C12 + 8e2 + 24 T12 ) D5 / 1204结语本文所述坐标变换技术原理简单 ,相应软件运行可靠 ,速度快 ,适合做大批量的数据转换 。但在理论上的精度仍略有欠缺 ,只能达到厘米级 ,有待在此基础上进一步完善提高 。 附录 :实用简化程序Co nver ti ng U TM to L atit ude a nd Lo ngit udecy a nd x su bt ract 500 ,000 f ro m co nve ntio nal U TM coo r di nat edi me n sio n xyz (3000000 ,3)c ha ract e r 3 15 ol dfile ,newfileop e n (5 ,file =p a ra . t xt , st at u s =ol d)rea d (5 , (a) ol dfile rea d (5 , (a) newfile rea d (5 , 3 ) nnop e n (3 ,file = ol dfile , st at u s =ol d)op e n (4 ,file = newfile , st at u s =un k no w n)rea d (3 , 3 ) ( ( xyz (i ,j) ,j = 1 ,3) ,i = 1 ,nn)do 120 i = 1 , nnxyz (i ,1) = xyz (i ,1) - 500000 . 0 co nti nuet lo ng0 = 123 . 0e = 0 . 00669437999 t k0 = 0 . 9996e12 = 0 . 0067394967421332 a = 6378137 . 0b = 6356752 . 3142 p e = sqr t (1 . 0 - e)e1 = (1 . 0 - p e) / (1 + p e)120= e1 3 1 . 5 - 27 . 0 3 e1 3 3 3/ 32 . 0= 21 . 0 3 e1 3 3 2/ 16 . 0 - 55 . 0 3 e1 3= 151 . 0 3 e1 3 3 3/ 96 . 0= 1097 . 0 3 e1 3 3 4/ 512 . 0tj1tj2 tj3 tj43 4/ 32 . 0gM = a 3 (1 . 0 - e/ 4 . 0 - 3 . 0 3 e 3 3 2/ 64 . 0 - 5 . 0 3 e 3 3 3/ 256 . 0)do 100 i = 1 ,n nt M = xyz (i ,2) / t k0 t M u = t M/ gMtfp = t M u + tj1 3 si n (2 3 t M u) + tj2 3 si n (4 3 t M u) + tj3 3 si n (6 3 t M u)+ tj4 3 si n (8 3 t M u)= e12 3 (co s ( tfp ) ) 3 3 2= ( t a n ( tfp ) ) 3 3 2= a 3 (1 - e) / sqr t ( (1 - e 3 si n ( tfp ) 3 si n ( tfp ) ) 3 3 3)#c1 t 1 r1t n1 = a/ sqr t (1 - e 3 si n ( tfp ) 3 si n ( tfp ) )d = xyz (i ,1) / ( t n1 3 t k0)q1q2 q3= t n1 3 t a n ( tfp ) / r1= d 3 3 2/ 2= (5 + 3 3 t 1 + 10 3 c1 - 4 3 c1 3 3 2 - 9 3 e12) 3 d 3 3 4/ 24 . 0q4 = (61 . 0 + 90 . 0 3 t 1 + 298 . 0 3 c1 + 45 . 0 3 t 1 3 3 2 - 3 3 c1 3 3 2 - 252 . 0 3 e12) 3 d 33 6/ 720 . 0q5q6 q7= d= (1 + 2 . 0 3 t 1 + c1) 3 d 3 3 3/ 6 . 0= (5 . 0 - 2 . 0 3 c1 + 28 . 0 3 t 1 - 3 . 0 3 c1 3 3 2 + 8 . 0 3 e12 + 24 . 0 3 t 1 3 3 2) 3 d 3 3 5/ 120 . 0xyz (i ,2) = ( tfp - q1 3 ( q2 - q3 + q4) ) 3 180 . 0/ 3 . 1415926xyz (i ,1) = tlo ng0 + ( ( q5 - q6 + q7) / co s ( tfp ) ) 3 180 . 0/ 3 . 1415926 co nti nuedo 150 i = 1 ,n nw rit e (4 ,1400) ( xyzl (i ,j) ,j = 1 ,3)co nti nuefo r mat (2f 15 . 7 ,f 10 . 2)e nd1001501400252海 洋 科 学 进 展24 卷运用本程序已一次性成功转换上百万个数据 ,以下为所附程序应用实例 (右边为输出数据) :668296 . 21668288 . 46668282 . 33668275 . 83668270 . 05668263 . 73668258 . 04668252 . 52668246 . 56668240 . 16124 . 7448807124 . 7447968124 . 7447357124 . 7446671124 . 7446136124 . 7445450124 . 7444839124 . 7444305124 . 7443695124 . 74430083320488 . 943320492 . 523320495 . 353320498 . 363320501 . 033320503 . 953320506 . 583320509 . 123320511 . 883320514 . 8458 . 2757 . 6157 . 3357 . 2957 . 1457 . 0757 . 1457 . 2457 . 2057 . 0030 . 003810930 . 003841430 . 003866230 . 003894830 . 003919630 . 003948230 . 003973030 . 003994030 . 004024530 . 004053158 . 2757 . 6157 . 3357 . 2957 . 1457 . 0757 . 1457 . 2457 . 2057 . 00参考文献 ( Ref e re nc e s) : 1 GAO J Y , FA N G Y X , XU S Y.Tra n sfo r matio n a nd unificatio n of vario us fo r mat s of raw multi bea m dat a . Ma ri ne Science Bulleti n J ,2002 , 21 ( 6) : 68274 . 高金耀 ,方银霞 ,徐赛英. 原始多波束数据的格式转换与统一 J . 海洋通报 ,2002 ,21 ( 6) :68274 .L I R C , WA N G Z Y. Mappi ng p rojectio n M . Wu ha n : Chi na U niver sit y of Geo science s Pre ss ,1992 . 李汝昌 ,王祖英. 地图投影 M . 武 汉 :中国地质大学出版社 ,1992 .W EI Z Q . Generatio n cha nge of Chi ne se geo detic coo r di nat e syst e m J . J o ur nal of Wu ha n U niver sit y : Info r matio n Science s , 2003 , 28 ( 2) :1382143 . 魏子卿. 我国大地坐标系的换代问题J . 武汉大学学报 :信息科学版 , 2003 , 28 :1382143 .N IN G J S. Mo der n geo detic ref erence syst e msJ . Sur veyi ng a nd Mappi ng Bulleti n , 2002 , 6 : 125 . 宁津生. 现代大地测量参考系统J .测绘通报 ,2002 ,6 :125 .ZHAO Y H . Experi ment al applicatio n of mo dified U TM i n ur ba n sur veyJ . Cit y Reco nnai ssa nce , 1994 , 18 : 43245 . 赵已鐄. U TM 投影 的改良用于城市测量的探索J . 城市勘测 ,1994 ,18 : 43245 .DU TC H S. Co nverti ng U TM to latit ude a nd lo ngit ude (o r vice ver sa) EB/ OL . 2005205219 . ht t p : www . uwgb . edu/ dut chs/ U sef ul2Dat a/ U TM Fo r mula s. ht m 2 3 4 5 6 Conversion of Geodetic Coordinates to UTM Rectangular Coordinatesin Processing the Multi2bea m Bathymetric DataXION G Yi ng2qia n1 ,2 , D U De2