高分辨率格网模型在地质灾害风险普查数据坐标转换中的应用研究

胡志瑞, 范宏涛, 马 欣, 黄 玮*, 赵大江,李 鹏, 李小琼

(1.宁夏国土资源调查监测院,宁夏 银川 750002; 2.宁夏地质灾害应急中心,宁夏 银川 750002;3.宁夏矿产资源储量评审中心,宁夏 银川 750002; 4.自然资源部大地测量数据处理中心,陕西 西安 710054;5.宁夏国土整治修复中心,宁夏 银川 750002)

第一次全国自然灾害综合风险普查是一次重大的国情国力调查[1],地质灾害风险普查作为其中一项,主要目标是掌握我国地质灾害风险隐患底数,为深入开展地质灾害防灾减灾工作奠定基础.地质灾害数据是进行地质灾害风险普查的基础.长期以来,地质灾害数据采用1954年北京坐标系(简称BJS54)或1980年西安坐标系(简称XAS80).BJS54,XAS80的测绘基准是参心、二维、低精度、静态的,已无法满足智慧交通、现代物流、防灾减灾等领域对高精度测绘地理信息服务的需求,制约高新技术的应用.CGCS2000具有三维、高精度、地心、动态、统一和实用等诸多优势,采用统一坐标系对经济建设、国防建设、社会发展和科学研究十分必要[2].我国从2008年开始启用2000国家大地坐标系(简称CGCS2000),并经10 a的过渡期进行推广应用[3].过渡期间,原测绘部门完成了基础测绘成果和基础地理信息数据库的坐标转换,三、四等天文大地网的平差,建立了独立坐标系与2000国家大地坐标系的联系,基本完成了现行大地坐标系向CGCS2000的过渡.

近年来,学者针对BJS54或XAS80向CGCS2000转换进行了大量研究.从转换方法看,主要采用参数法[4—5];从数据格式看,主要有CAD,ArcGIS平台格式[6—7].但基于国产MapGIS平台数据、针对地质灾害领域的坐标转换研究相对较少.空间位置是否准确统一、属性信息是否完整翔实,影响地质灾害数据质量[8],而高质量数据也是构建防灾减灾数据体系、发挥地质灾害数据价值的基础.在全国各省地质灾害风险普查工作中,收集到的各类数据需预处理,而数据坐标系统转换是模型重要的环节,关系到各类数据的套合,不仅影响数据的精度,还决定后期空间分析和评价结果的准确性.因此有必要建立一套适用于省级范围的坐标转换模型.笔者研究地质灾害数据由BJS54或XAS80向CGCS2000的高精度转换方法,以期为其他地区数据的坐标转换提供参考.

1 坐标转换方法

1.1 参数法

参数法是基于方程的变换方法,根据转换区域大小选取适当的转换模型,通过选取重合点并经过粗差分析后计算转换参数,然后进行坐标转换的一种转换方法.对于面积较大区域,一般采用分区转换.二维七参数模型、三维七参数模型以及椭球面的多项式拟合算法适用于全国或椭球面的经纬度差在3°及以上的省级大地坐标的转换;布尔莎模型适用于空间直角坐标的转换;莫洛金斯基模型适用于省级及以下范围大地坐标的转换;三维四参数模型适用于2°×2°以内局部区域空间直角坐标的转换;局部区域平面坐标的转换可用二维四参数模型[9].

1.2 格网法

格网法的基本思路是,将待转换区域以一定网格划分,根据转换前后坐标系统间公共点的坐标差,通过一定的算法计算各格网结点的坐标差.然后用坐标差来内插网格内任意点的值求得改正量.这样就实现了坐标系的变换.算法主要采用协方差推估、最小曲率、加权平均等数学方法.国外应用格网法转换坐标较早,美国采用最小曲率法[10]、日本采用克里格插值法[11]、澳大利亚采用最小二乘法[12],分别建立了适用于本国的大地坐标转换模型.我国学者近年对格网模型进行研究[13],开发了一些坐标转换软件[14].

目前,BJS54,XAS80与CGCS2000坐标数据的相互转换,针对测绘数据的研究较多,方法也相对成熟[15],但针对地质灾害数据的研究较少.

采用格网模型的坐标转换优势在于不受制图比例尺和范围限制;相同地理位置的转换改正量保持不变,而且转换改正量是连续的,弥补了分片区转换后相邻分区之间改正量不连续的缺陷.格网模型适合连续分幅的地图转换且便于后期图幅接边.笔者采用文献[13]中的方法,对宁夏地质灾害数据的坐标转换进行研究.

2 数据与流程

2.1 实验数据

地质灾害数据具有专业性强、多源异构、属性信息复杂、数据量大等特点:从比例尺看,各种比例尺并存;从坐标系看,有BJS54,XAS80且以BJS54为主;从数据类型看,有矢量数据、栅格数据和文本数据,矢量数据是在MapGIS平台产生的,以*.wt(点),*.wl(线),*.wp(面)格式为主.笔者以宁夏地质资料馆馆藏地质灾害详查数据、区域地质图等地质灾害风险普查所需基础数据为研究对象(表1).

表1 实验数据

2.2 转换流程

研究表明,BJS54转换到CGCS2000,其图幅X方向平移-29～-62 m、Y方向平移-56～84 m;XAS80转换到CGCS2000,其图幅X方向平移-9～43 m、Y方向平移76～119 m[9].各类地形图因坐标系转换造成内部任意两点的长度和方位变动均在制图精度内,可忽略不计,但由坐标系更换引起图廓点坐标的变化超出制图精度范围.以1∶5万图幅为例,BJS54转换到CGCS2000,图廓X方向平移约1.2 mm、Y方向平移约1.7 mm;XAS80转换到CGCS2000,图廓X方向平移约0.8 mm、Y方向平移约2.4 mm.同时,不同坐标系下转换后数据接边与重合不容忽略.

数据转换点位的平均精度应小于图上0.1 mm.1∶50万图的方里网和图廓线平均位移0.1 mm,而1∶25万图的方里网和图廓线平均位移达0.2 mm,因此大于等于1∶25万地形图中点(图廓点)的位置变动已超过制图精度,需重新标记.对于小于1∶25万的地质图,由坐标系更换引起图廓点坐标的变动在制图精度以内,可忽略其影响.所以对于大于或等于1∶25万的地质图、各类成果数据,必须进行坐标转换.

通过模型计算宁夏XAS80向CGCS2000转换的高精度、高分辨率格网改正量(dB,dL),可得出格网点CGCS2000下大地坐标.BJS54转换到XAS80坐标系已有成熟的方法,因此BJS54转CGCS2000可先转为XAS80,再通过模型计算.宁夏XAS80向CGCS2000转换的流程见图1.

图1 XAS80到CGCS2000的转换流程

3 高分辨率格网模型的建立

3.1 控制点选择

范围为宁夏境内2000国家GPS大地控制网点,GPS C级网点,国家天文大地网与高精度GPS2000网联合平差点,全国三、四等三角网在CGCS2000下平差的三、四等三角网点.为保证宁夏省界处的改正量精度,还搜集了相邻省份(甘肃、内蒙古、陕西)与宁夏行政界线相距10 km以内的控制点.

3.2 控制点分析

控制点的精度决定着坐标转换改正量的精度,所以需要对参与计算的控制点的来源、精度、获取年代进行分析.计算时优先选择获取年代相对较近的点位.

控制点按精度选取先后顺序,即与国家三角点重合的2000国家GPS大地控制网点,与国家三角点重合的GPS C级网点,国家天文大地网与高精度GPS2000网联合平差点,全国三、四等三角网在CGCS2000下平差的三、四等三角网点.国家三角点布测于半个世纪前,由于人类活动等各种因素的影响,部分点位发生了变化,因而存在粗差点,如果这些粗差点参与改正量计算,会影响结果的精度.

因此,在选取控制点后,需对选取的各控制点进行试算,直至选取的控制点不含有粗差点.首先求取与控制点距离最近的至少3个控制点的改正量.然后与该点的改正量比较,若该点的改正量绝对值介于周边几个控制点改正量绝对值的平均值的1/3与3倍,则认为该控制点正常,否则认为是异常点.对异常点根据分布图形进一步判定,确定后剔除.对每个控制点遍历以上步骤后,最终确定用于计算宁夏全区高分辨率格网坐标转换改正量的控制点为2 429个.

3.3 改正量计算模型

对2 429个控制点空间上的密度进行分析,通过反距离加权模型获得连续高精度的趋势[16].反距离权重法假定每个测量点都受一种局部影响,而这种影响随着距离的增大而减小.当为任何未测量的位置预测时,反距离权重法预测位置周围的测量值.与距离预测位置较远的测量值相比,距离预测位置最近的测量值对预测值的影响较大.这种方法为距离预测位置最近的点分配的权重较大,且权重与反距离(数据点与预测位置之间)的P次幂成正比.具体计算公式:

(1)

(2)

B2000=B80+dB,

L2000=L80+dL,

(3)

式中:ΔBi=Bi,2000-Bi,80,ΔLi=Li,2000-Li,80;Bi,Li为控制点大地坐标,B0,L0为格网点大地坐标;B2000,L2000为目标CGCS2000大地坐标,B80,L80为XAS80大地坐标.

3.4 格网分辨率的确定

分辨率设定思路是使每个格网内有1个控制点,格网的分辨率越高,内插值越接近于真值,坐标转换精度也就越高.由于该过程的运算量大,影响转换速度,但当格网分辨率达到一定值时,转换精度基本趋于稳定.根据测绘一、二、三、四等控制点的实际间隔[13],将格网分辨率设为3″×3″可满足计算要求,即网格大小为92.7 m.

3.5 改正量的计算

3.5.1XAS80向CGCS2000转换改正量的计算 在椭球面上,以XAS80下各个格网点为中心,在一定的搜索区域选择已剔除局部粗差点的一、二、三、四等控制点,依次求取各控制点的大地坐标改正量.然后采用反距离加权平均法计算宁夏境内XAS80坐标系到CGCS2000的高分辨率、3″×3″间隔格网改正量dB80T2000,dL80T2000.搜索区域一般为矩形、圆形和椭圆形.研究中,需考虑在各个方向上均等的数据点,因此将搜索区域定义为圆形.

3.5.2BJS54向CGCS2000转换改正量的计算 BJS54是通过逐级控制、分级分区平差的,没有经过整体平差,因而精度较低.因此在不同的区域,坐标精度不同,有的地方相差甚至达到米级.

XAS80建立后,基本比例尺地形图更新时使用XAS80.BJS54转换到XAS80已有较成熟、高精度的方法.因此,BJS54向CGCS2000转换高精度改正量的计算步骤:第一步,计算BJS54向XAS80转换的改正量dB54T80,dL54T80;第二步,计算XAS80向CGCS2000转换的改正量dB80T2000,dL80T2000;第三步,将BJS54向XAS80转换的改正量与高分辨率格网XAS80向CGCS2000转换的改正量叠加,得到高精度、高分辨率格网BJS54向CGCS2000转换的改正量dB54T2000,dL54T2000,即dB54T2000=dB54T80+dB80T2000,dL54T2000=dL54T80+dL80T2000.

通过以上方法得到XAS80向CGCS2000转换的改正量和BJS54向CGCS2000转换的改正量.在实际转换中,先通过MapGIS获得工程文件中的点、线、面文件.点文件通过改正量换算为转换后的坐标,线、面文件也是由若干点组成的.因此先读取构成线或面文件中各个点的坐标,然后将每个点坐标转换更新,从而实现线、面文件的转换.

4 精度评价

坐标转换精度评价有内符合和外符合2种评价方法.依据计算转换改正量的重合点残差中误差评估坐标转换精度,若点位残差小于3倍点位中误差,满足精度要求.笔者采用内符合评价方法,计算转换前后坐标残差中误差,对坐标转换精度进行评估[17].对于n个评估点,坐标转换模型精度计算公式:

对于空间直角坐标X残差中误差,

(4)

对于空间直角坐标Y残差中误差,

(5)

式中:n为点位个数;v为重合点残差;v′为重合点转换坐标值与重合点已知坐标值的差.

平面的点位中误差:

(6)

4.1 XAS80向CGCS2000转换的精度评估

对宁夏具有XAS80和CGCS2000的三角点,选取257个重合点进行精度评估,点位坐标残差频率分布见图2～图3,横坐标是评估点的残差范围,纵坐标是对应范围点的个数.

图2 X坐标残差直方图

图3 Y坐标残差直方图

由图2～图3可知,X坐标残差基本呈正态分布,残差最大值为0.408,最小值为-0.638,11.3%的点在(-0.002,-0.001],9.7%的点在(-0.001,0],8.6%的点在(0,0.001],9.3%的点在(0.001,0.002].Y坐标残差最大值为0.289,最小值为-0.501,10.1%的点在(-0.002,-0.001],14.4%的点在(0.001,0],10.9%的点在(0,0.001],8.9%的点在(0.001,0.002].经过计算,X坐标的残差中误差MX=0.003 m,Y坐标的残差中误差MY=0.003 m,平面的点位的中误差Mp=0.004 m.

4.2 BJS54向CGCS2000转换的精度评估

从宁夏具有BJS54和CGCS2000的三角点中选取257个重合点进行精度评估,点位坐标残差频率分布见图4～图5.

图4 X坐标残差直方图

图5 Y坐标残差直方图

由图4～图5可知,X坐标残差频率呈正态分布,25.3%的点在(0,0.1].残差最大值为0.408,最小值为-0.638;Y坐标残差最大值为0.289,最小值为-0.501,50.6%的点在(0,0.1].经过计算,X坐标的残差中误差MX=0.181 m,Y坐标的残差中误差MY=0.142 m,平面的点位中误差Mp=0.230 m.

4.3 BJS54和XAS80向CGCS2000转换的精度对比

将BJS54和XAS80下的257个重合点分别转换到CGCS2000,进行残差对比分析,点位精度分布见图6～图7.对比发现,XAS80转换CGCS2000的精度高于BJS54转换的.地质风险普查是以县为单元开展的,比例尺为1∶5万.按照精度要求[9],1∶5万空间数据库坐标转换精度小于等于5.0 m.根据控制点内符合精度,文中所建BJS54和XAS80转换CGCS2000的高分辨率格网模型完全满足县级地质灾害风险普查的精度要求.

图6 XAS80,BJS54下X坐标的转换残差

图7 XAS80,BJS54下Y坐标的转换残差

5 结论

1)探讨地质灾害风险普查相关数据由XAS80及BJS54分别转换到CGCS2000的流程方法,建立省级高精度、高分辨率格网模型,并进行精度评估.结果表明,应用该方法,BJS54,XAS80转换到CGCS2000的数据,均满足精度要求,其中,XAS80转换CGCS2000的精度优于BJS54转换CGCS2000的.

2)在宁夏完成了县级MapGIS格式地质灾害数据坐标转换,并很好地解决区域不同图幅因采用参数法转换坐标造成的接边问题.

3)部分地质图中的高程采用1956年确定的黄海高程系,并在MapGIS点文件中标注,需与1985年国家高程基准加以区分.地质灾害工作中经常使用WGS84坐标系.WGS84坐标系与CGCS2000的参考椭球体相近,而历元和参考框架不一致,二者存在分米级(同一历元下在厘米级)差别[18],若测量精度在此范围内,可视为WGS84坐标系与CGCS2000一致.

5)在宁夏实践地质灾害数据坐标转换方法,可为其他区域数据转换提供参考.