席新海,朱卫民
(1.河南省地矿局第二地质环境调查院,河南 郑州 450000;2.河南省地矿局测绘地理信息院,河南 郑州 450006)
西藏安多县抱布德铅矿详查测量工作是该区域地质找矿工作流程之一。主要工作内容为测绘覆盖矿区范围1∶10000比例尺地形图;对圈定的详查地段开展1∶2000比例尺地形图测绘和地质工程测量包括地质点、钻孔孔位、探槽起止点、勘探基线1∶1000剖面线测量及基点、端点测量等。
由于矿区位于藏北无人区,平均海拔5200m,工作条件异常艰苦。造成至今还是我国大比例尺测绘成果的盲区。经查询该区域已有测绘资料严重缺乏,基础控制资料缺少,图形资料年代久远。经实地踏勘测区工作条件,很难应用常规的的测量技术手段完成本次测量工作,因此选择合适的测量方案尤为关键,结合工作区实际情况,利用先进的测量技术,在西藏抱布德铅矿详查测量实际工作中进行了一些有益尝试与探讨[1]。
由于西藏安多抱布德铅矿详查矿区缺少必要的已知控制点,无法开展常规控制测量工作。而IGS(International GNSS Service,国际全球卫星定位导航服务组织)站为全球连续运行参考站[1],具有连续观测、精度高、稳定等特点,为无人区控制测量提供了最佳的起算数据,使荒漠地区的测绘工作得以开展。因此西藏抱布德铅矿详查项目就利用5个IGS站,布设矿区C级GPS网,作为本次测量工作的首级控制。
测区内共布设5个C级GPS控制点,编号采用C接自然数顺序编号,如C01、C02、……。实际观测时C级GPS网采用5台仪器同时观测,观测两个时段,采用边连接的方式,每个时段按照静态观测方法连续观测4个小时以上,重复设站数为2。
2.2.1 数据标准化
采用Rinex标准化观测数据文件SITEDAYS.YYO,其中SITE为点位编码,DAY为年积日,S为观测时段号,YY为观测年号,O为观测数据。
2.2.2 收集IGS连续运行站数据
收集IGS数据是为了与外业观测的五个C级GPS点统一计算平差,本次数据处理共收集了KIT3、LHAZ(拉 萨)、SX01(安 乐 河)、ULAB、URUM(乌鲁木齐)五个GPS连续运行站的观测数据。要求是与所布设的GPS点时间一一对应。数据处理采用美国麻省理工学院的GAMIT/GLOBK软件进行基线解算及平差。
2.2.3 数据整理
依据外业观测数据,将同一观测时段的数据进行预处理,数据格式为Rinex格式,并进行以下数据正确性的检验:
点名一致性与正确性;
接收机与天线型号的正确性;
天线高的正确性;
年积日的一致性。
2.2.4 天线高的归算
观测时使用TRIMBLE 5700接收机、TRM39105.00天线及TRM41249.00 天线和TRIMBLE 5800接收机及相应天线。按照天线结构,天线高统一采用观测值归算。在基线解算时由GAMIT软件自动计算天线相位中心位置,归算至标石标志面。
2.3.1 先验坐标的获取
先验坐标采用差分的方法获得,以GPS连续运行站为基准站进行差分,求得GPS观测站的先验坐标,保证坐标可以达到0.1m的精度。
2.3.2 主要参数设置
(1)卫星轨道:IGS快速精密星历,且固定[2];
(2)卫星截至高度角:10s;
(3)数据采集间隔:15s;
(4)坐标约束:GPS连续运行站水平方向约束5cm、垂直方向约束10cm,其他GPS站约束10m;
(5)对流层改正模型:采用Saastamoinen(萨斯塔莫宁)模型进行标准气象改正;
(6)观测值:采用消除电离层后的组合观测值;
(7)数据解算模式:周跳自动修复技术[3]。
2.3.3 GPS基线解算
以GPS Day(年积日)为单位,进行基线解算,见表1。
表1 GPS网环Nrms统计表
由表1可知,GPS网的Nrms值小于0.25周,基线解算时周跳基本剔除干净。
对整网的全部基线结果进行了X2检验,X2检验值小于6,基线全部通过检验,参与平差,见表2。
表2 GPS网X2检验结果统计表
在2000国家大地坐标系下,约束本次数据处理收集的 KIT3、LHAZ(拉萨)、ULAB、URNM(乌鲁木齐)五个GPS连续运行站[4],进行三维约束平差,求出待定GPS网点坐标。
?统计项 最小值/mm 最大值/mm 平均值/mm Xrms 2.7 2.8 2.8 Yrms 4.3 4.8 4.6 Zrms 2.7 3.1 2.9 Nrms 1.2 1.3 1.3 Erms 2.7 2.8 2.8 Urms 4.9 5.6 5.3
由表3可知,GPS点空间直角坐标X分量的中误差平均值为2.8mm,最大值为±2.8mm;Y分量的中误差平均值为±4.6mm,最大值为±4.8mm;Z分量的中误差平均值为±2.9mm,最大值为±3.1mm。GPS点南北分量的中误差平均值为±1.3mm,最大值为±1.3mm;东西分量的中误差平均值为±2.8mm,最大值为±2.8mm;高程分量的中误差平均值为±5.3mm,最大值为±5.6mm。
通过约束平差获得矿区C级GPS点2000坐标下的大地直角坐标,通过坐标转化提供西安80坐标,为后续测绘工作奠定了基础。
TerraSAR-X卫星为德国研制的一颗高分辨率雷达卫星,携带一颗高频率的X波段合成孔径雷达传感器,可以聚束式、条带式和推扫式3种模式成像,并拥有多种极化方式。可全天时、全天候地获取用户要求的任一成像区域的高分辨率影像。Tan-DEM-X于2010年6月21日成功发射,这两颗卫星在3年内将反复扫描整个地球表面,最终绘制出高精度的3D地球数字模型。因此,利用TerraSAR-X卫星数据,基于雷达立体测量技术,快速获得无人区DEM,结合1∶50000国家地形图及外业实测数据等,为开展无人区1∶10000比例尺测图提供了一种新方法,从而摆脱了传统的测量方法,大大减轻了外业测量工作的作业难度,为快速开展无人区地质测量工作提供了新的作业方案[5]。
(1)多分辨率 (1m/3m/18.5m)和覆盖区域:对于特定目标区域采用高分辨率,对大面积覆盖采用中等分辨率;
(2)任何其他的商业星载传感器都无法比拟的几何精度;
(3)极高的辐射精度;
(4)不受天气影响,对地球上的任何地点,重访周期最长2.5天(95%的地区可达到2天重访);
(5)独特的敏捷性(成像和极化模式的快速切换)。
TerraSAR-X卫星数据生成DEM的流程主要包括:SAR复数影像配准、干涉图的生成、去平地相位、干涉相位滤波、相位解缠、基线精确估计及数字高程模型重建[6],见图1。
图1 TerraSAR-X卫星生成DEM流程图
3.2.1 影像的精确配准
用于生成DEM的两幅SAR影像不是同时获取的,相应位置像素点并不是完全对应的,因此必须先对两幅SAR影像进行高精度复图像配准,使得空间各像素点一一对应。影像配准的好坏直接影响干涉条纹的生成,从而影响提取DEM的精度,通常InSAR影像配准精度要求达到1/8个像元精度。
3.2.2 干涉图的生成
在两幅精确配准的SAR复图像中,对应像元基本可以反映同一区域的相干目标特性,将对应的像元复共轭相乘,就可以得到干涉图。
3.2.3 去平地相位
基于卫星轨道参数和影像位置计算平地相位,其根据斜距的分辨率,由影像中心点开始,推算出各个像元在参考椭球面的斜距。利用各像元点间的斜距差得到各像元点间的相位差,进而求得参考椭球面所引起的干涉相位差,将其从原始干涉相位中减去从而去除平地相位的影响。
3.2.4 干涉相位滤波
InSAR干涉图中不可避免的存在各种噪声,造成干涉相位连续性和周期性的破坏,从而影响相位解缠的精度。因此为保证测量结果的准确性必须消除噪声相位影响。
3.2.5 相位解缠
由于干涉相位值在[-π,π]的范围内,只是干涉相位的主值,其整数部分信息丢失,相位解缠操作就是恢复丢失的整数相位信息,以得到真正的相位差值。
3.2.6 精确基线估计
基线是InSAR数据处理中非常重要的一个参数,高精度的基线参数是获取高精度DEM的根本保证。常用的估计方法有:基于轨道信息的基线估计法、基于条纹频率的基线估计方法以及基于地面控制点(GCP)的基线估计方法等。
3.2.7 DEM 生成
通过相位解缠得到相位差真值后,结合高精度的基线数据,就可以计算出每一点的高程值。因为解缠后的相位经过高程估算之后仍然在距离/多普勒坐标系中,需要通过地理编码操作将其转换到地理坐标系下。地理编码不但可以为影像的各像素点确定实际地理位置坐标,而且对影像进行了正射纠正,剔除了干涉影像中的阴影和叠掩部分,得到的是DEM的正射影像图。
(1)利用Global Mapper软件打开DEM数据
(2)在Global Mapper中选取需要的区域,生成所需要的等高线数据。
(3)再将生成的等高线数据输出转成DXF格式。
(4)利用AutoCAD软件打开并编辑生成矿区等高线数据,见图2。
图2 DEM数据获取等高线
工程点测量主要为钻孔位置的测量及基线端点、基点、探槽位置点测量,均采用RTK进行施测。另外由于控制测量仅提供西安80坐标成果,而实际地质工作是在1954年北京坐标系下开展,因此需采用坐标转换方法来提供工程点的1954坐标。
计算模型:“7参数”法(参心平移量:ΔX、ΔY、ΔY;尺度比:M;坐标旋转量:)假设西安80坐标系和北京54坐标系有七个转换参数-3个平移参数、3个旋转参数和1个尺度参数[7],见图3。
图3 “7参数”法模型
假设:
(X80Y80Z80)T为某点在西安80坐标系下的空间直角坐标;
(X54Y54Z54)T为该点在北京54坐标系系下的空间直角坐标;
(ΔX0ΔY0ΔZ0)T为西安80坐标系转换到北京54坐标系的平移参数;
(ωXωYωZ)为西安80坐标系转换到北京54坐标系的旋转参数;
m为西安80坐标系转换到北京54坐标系的尺度参数。
则由西安80坐标系转换到北京54坐标系的转换关系为:
其中:
具体计算步骤如下:
(1)将西安80坐标高斯反算为大地坐标。
(2)将北京54、西安80的大地坐标转换为空间之间坐标,利用7参数法求解西安80坐标系至北京54的7个转换参数,并验证其精度。
(3)利用计算得到的7参数,可将本次测量的各个工程点转换至北京54坐标系下。
(4)将计算得到的北京54空间三维坐标投影为平面坐标。投影时,将空间三维直角坐标转为大地坐标,再进行高斯投影得到平面坐标。
利用5个C级GPS控制点的西安80、北京54坐标,求取西安80到北京54坐标7参数进行数据转换;七个参数包括3个平移参数、3个旋转参数和1个尺度参数(参心平移量:ΔX、ΔY、ΔY;尺度比:M;坐标旋转量:)[8]。
经过转换,转换残差平面最大3mm,最小1mm,精度远远高于地质工程点西安80到北京54坐标的转换要求。
地质工程测量主要包括地质工程点测量(包括钻口、探槽、地质点测量)、地质勘探基线基点测量、勘探线剖面测量等。
主要为钻孔、探槽、地质点的位置测量。以测区内控制点为起算点,输入坐标转换参数,采用美国天宝公司生产的动态GPS对钻孔、探槽、地质点位置直接测定。野外作业时,每次基准站架设好后,均采用流动站检查基准站点坐标及另外已知点坐标,以确保设站无误。每个基准站均进行重测检查,从而保证了GPS测量质量,最后依次测量各工程点,每个工程点测量时都采用固定解,输入正确的棱镜高。
本矿区共计测设钻孔27个,探槽27个,地质点454个,其中钻孔实测最大平面中误差0.024m,限差0.3m,高程最大中误差0.037m,限差0.25m;探槽实测最大平面中误差0.051m,限差1.6m;符合《地质矿产勘查测量规范》(GB/T18341-2001)的规定。
基线测量主要实测基线的起点、终点以及基线与勘探的交叉点,作业过程首先将设计的点采用RTK技术放样到实地,然后再实测设计点的平面位置与高程。
根据地质要求先在相应比例尺的地形图上设计出勘探线两个端点,并在地形图上精确量取勘探线两个端点坐标,野外作业时把勘探线两个端点坐标输入GPS手簿,利用线放样在实地沿着这条线逐一测量地形点和地质分界点。
实测勘探线剖面点最大点位平面中误差0.08m,限差为1.6m,最大高程中误差0.37m,限差为0.67m,均符合《地质矿产勘查测量规范》的规定。
本文从实际生产出发,对西藏安多县抱布德铅矿详查测量技术方案进行了详细的应用探讨,为无人区高海拔开展地质勘探测量项目关键技术提出了解决办法。利用国际IGS站作为起算数据,采用Gamit软件进行数据处理,不仅作业方法先进,而且可以方便快速的布设控制网。另外提出利用TerraSAR-X卫星数据,基于雷达立体测量技术,快速获得无人区DEM,结合1∶50000国家地形图及外业实测数据等,开展无人区1∶10000比例尺测图的技术方法。为荒漠困难地区测绘1∶10000地形图提供了可借鉴的工作经验,是一种有意义的尝试。
[1]包龙海,高佑民,韦书平.无人区域GPS技术作业方案及精度[J].陕西煤炭,2010(1):70-72.
[2]窦和军,武晓忠.应用GAMIT软件解算高精度GPS网[J].西部探矿工程,2007(1):100-102.
[3]赵桂儒,杨天青,徐平.小区域GPS变形监测网GAMIT数据处理结果与IGS站选取的关系探讨[J].地震地磁观测与研究,2006(5):103-106.
[4]赵桂儒.基于GAMIT软件的GPS数据处理框架建设[D].中国地震局地震预测研究所,2007.
[5]侯瑞,谭志祥,黄国满.新型高分辨率星载SAR卫星-TerraSAR-X.中国科技论文在线.
[6]陈艳玲,黄珹,冯天厚.TerraSAR-X卫星及其在地球科学中的应用.中国科学院上海天文台年刊,2007(28):51-57.
[7]柯美忠,张龙其,刘平,邢光成.北京1954坐标系和西安1980坐标系转换的探讨[J].地理空间信息,2005(4):55-57.
[8]束蝉方,李斐,沈飞:空间直角坐标向大地坐标转换的新算法[J].武汉大学学报(信息科学版),2009(5):561-563.