武汉大学姜卫平等:大地测量坐标框架建立的进展与思考
The following article is from 智绘科服 Author 测绘学报
姜卫平1,2
1. 武汉大学卫星导航定位技术研究中心, 湖北 武汉 430079;
2. 湖北珞珈实验室, 湖北 武汉 430079
基金项目:湖北省科技创新人才及服务专项(2022EJD010);国家自然科学基金创新研究群体项目(41721003);湖北珞珈实验室专项基金(220100020)
摘要:坐标框架是描述地球形状及变化、表达地球空间信息的基础, 也是拓展人类活动、促进社会发展的关键地球空间信息基础设施。随着空间大地测量观测技术的发展, 地球科学及相关学科间的交叉渗透融合, 利用其建立全球或区域坐标框架成为当前大地测量的主要任务。研究建立1毫米级坐标框架是国际大地测量学界21世纪的学科目标和重要挑战。本文以当前建立理论最完善、应用最广泛、精度最高的国际地球参考框架为例, 描述了基于空间大地测量观测技术的坐标框架建立方法, 阐述了全球及区域地心坐标框架建设的最新进展及局限性, 最后针对建立1毫米级坐标框架的几个关键问题提出了研究思路。
关键词:坐标系统 坐标框架 空间大地测量技术 非线性变化 多技术组合
引文格式:姜卫平, 李昭, 魏娜, 等. 大地测量坐标框架建立的进展与思考[J]. 测绘学报,2022,51(7):1259-1270. DOI:10.11947/j.AGCS.2022.20220232
JIANG Weiping, LI Zhao, WEI Na, et al. Progress and thoughts on establishment of geodetic coordinate frame[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(7): 1259-1270. DOI: 10.11947/j.AGCS.2022.20220232
引 言
人类活动80%以上的信息与空间位置有关。描述位置信息的前提是建立参照系,也就是坐标系统。大地坐标系统的确定包括选择一个椭球、对椭球进行定位和确定大地起算数据,也就是定义坐标系的原点、轴向和尺度[1]。大地坐标框架(以下简称“坐标框架”)是大地坐标系统的实现,是描述地球形状及其变化和表达地球空间信息的基础,由一组具有坐标(及其随时间变化)的观测站组成。由于科学技术的制约以及其他历史原因,在20世纪末以前,世界上所建立的大地坐标系统及大地坐标框架基本表现为二维、参心特征,采用局域定位和地面网点传递的技术方式提供坐标,未考虑板块运动、地表质量重分布等地球动力学效应对地面点的时变影响,相对精度大约为10-5量级[2-3]。例如,我国的1954年北京坐标系、1980西安坐标系就是如此,其定义包括所采用的地球椭球、大地原点位置及椭球XYZ轴指向,对应的坐标框架则是用经典大地测量技术所测定的全国天文大地网[2-3]。
20世纪末叶,全球导航卫星系统(GNSS)、甚长基线干涉测量(VLBI)、卫星激光测距(SLR)、多普勒无线电定轨定位系统(DORIS)等空间大地测量观测手段成为建立全球或区域坐标框架不可或缺的重要观测技术。随后,地心坐标系及其框架开始逐渐取代传统的参心坐标系统及其坐标框架。例如,GPS采用的世界大地坐标系统(world geodetic system, WGS)及其框架、俄罗斯格洛纳斯(GLONASS)坐标系统及其参考框架PZ-90(Parametry Zelmy 1990)、欧盟伽利略地球参考框架(Galileo terrestrial reference frame, GTRF)、中国北斗导航卫星系统(BeiDou navigation satellite system, BDS)采用的坐标系统(BeiDou coordinate system,BDCS)及其框架(BDS terrestrial reference frame, BTRF)、国际地球参考系统(international terrestrial reference system, ITRS)及其实现国际地球参考框架(international terrestrial reference frame, ITRF),等。其中,ITRF是目前建立理论最完善、应用最广泛、精度最高的全球地心坐标框架,为其他全球性和区域性坐标框架提供统一的空间基准[4]。为了促进可持续发展,2015年联合国通过了采用ITRF作为全球统一大地测量参考框架的决议[5]。
近年来,很多国家也在积极推进区域坐标框架建设的进程,以ITRF为基准,利用GNSS等技术更新了各自的国家/区域地心坐标框架。这些区域坐标框架可以作为ITRF的加密或延伸,不仅支持了ITRF建设,还为全球经济一体化、全球大地测量参考框架综合服务体系的构建奠定了非常坚实的基础。迄今为止,超过80%的国家地心坐标框架与ITRF对准[6],例如,2007年启用的最新韩国大地测量基准2002(Korean Geodetic Datum 2002,KGD2002)与ITRF2000对准[7],2020年启用的最新澳大利亚参考框架2014(Australian Terrestrial Reference Frame 2014, ATRF2014)与ITRF2014一致(https://icsm.gov.au/upgrades-australian-geospatial-reference-system)等。以我国为例,最新的2000国家大地坐标系统(China geodetic coordinate system 2000, CGCS2000)定义与ITRS一致。其实现,也就是CGCS2000坐标框架,代表了我国坐标基准建设的最高水平,精度显著优于我国长期采用的1954北京坐标系、1980西安坐标系[8-9]。然而,严格来说,CGCS2000坐标框架属于区域、静态坐标框架,整体精度相对偏低。目前,一些国家重大工程仍部分依赖于国际数据,难以满足我国测绘基准现代化对高精度坐标框架的需求。
2016年发布的ITRF2014采用4种空间大地观测技术,基于全球第二次数据重新处理计划(repro2)建立。相较于ITRF2008,ITRF2014采用的观测数据及测站数量更多,数据处理模型及策略更先进,并且首次考虑了基准站的非线性运动。因此,其精度优于以往所有ITRF版本,但是长期精度仍为厘米级,无法满足气候变化、地质灾害、地震等大范围或全球尺度毫米级地球系统动态变化监测的需求[4, 10]。尤其是长期海平面变化监测,需要坐标框架的精确度和稳定性水平分别达到1 mm和1 mm/a[11]。2022年4月,最新的ITRF2020正式发布(https://itrf.ign.fr/en/solutions/ITRF2020),不仅提供基准站在参考历元时刻的位置及长期速度,还包括大地震造成的震后形变及周年、半周年参数模型。ITRF2020产品精度较ITRF2014有所提高,但是其量级仍有待进一步评估[6]。研究建立1 mm(中误差1毫米,限差3倍中误差,即1~3 mm)级坐标框架迫在眉睫,是目前大地测量领域面临的一项新任务和新挑战,同时也是全球大地测量观测系统(global geodetic observation system, GGOS)的研究目标,还是国际大地测量学界21世纪的中长期学科目标[12]。
本文首先介绍基于现代空间大地测量技术的坐标框架建立理论与方法,然后详细阐述全球及区域坐标框架的最新进展及其局限性,最后围绕构建毫米级坐标框架的几个关键问题进行了展望,并给出了研究思路。
1 基于空间大地测量技术的坐标框架建立方法
1.1 坐标框架定义及确定方法
基于空间大地测量观测技术建立的坐标框架理论上应该是ITRS的实现。原点、尺度、定向及其随时间的演变是建立坐标框架的必备要素,也就是基准。确定了一套基准,就确定了一个坐标框架。动态坐标框架在数学上表现为14个参数,即7参数(包括3个平移参数、3个旋转参数、1个尺度参数)及其随时间的演变[13]。考虑到ITRF是目前精度最高的大地测量坐标框架,本文以ITRF为例,阐述基于空间大地测量观测技术的坐标框架建立方法。
ITRF的原点理论上位于包括固体地球、海洋和大气的地球质量中心(center of mass, CM),通常利用SLR技术对坐标框架的原点进行约束实现[14]。例如,ITRF2014的原点由SLR长期地球参考框架(terrestrial reference frame, TRF)的原点确定。在进行时间序列堆栈建立SLR长期TRF时,对其平移参数及其速率附加了内部约束,以保证无外部原点基准引入到SLR长期TRF中。值得注意的是,采用的内部约束方法仅顾及一阶项,即仅对ITRF原点相对于CM的长期变化进行约束,并不约束由地表质量负载等引起的季节性变化[15]。严格来讲,目前实现的ITRF原点既不是CM,也不是固体地球的形状中心(center of figure, CF),而是由实际地面观测网维持的地球中心(center of network, CN)。实际布设地面网时都会尽量选择均匀分布的全球网,此时CN近似于CF。由于附加了外部基准约束,ITRF的原点在长时间尺度上近似于CM,近似的程度取决于附加基准约束的质量;在季节性时间尺度上近似于CF,近似的程度取决于观测网空间分布。
坐标框架的尺度通常由光速c、地球总质量、地心引力常数GM及相对论改正模型共同确定。理论上,VLBI、SLR、GNSS和DORIS 4种技术均可以定义ITRF的尺度。GNSS技术除了受GM和c的影响外,还受到如卫星天线相位中心的Z方向偏差、对流层延迟等技术相关系统误差的影响,目前并不适合确定尺度。VLBI数据处理采用几何方法,基本与地球引力场无关,其尺度因子主要取决于光速c,因此VLBI尺度因子的长期稳定度优于SLR和GNSS[13]。鉴于上述原因,从ITRF2000开始,GNSS技术不再参与定义尺度,仅VLBI和SLR用于ITRF尺度定义。
ITRF定向随时间的演变相对于地球表面的水平运动符合相对于地壳整体无残余旋转[13],采用上一版本ITRF提供的实测速度场对定向随时间的演变进行约束[16]。例如,ITRF2014在历元2010.0,相对于ITRF2008的旋转参数和旋转速率为零[4]。基于空间大地测量技术建立的坐标框架定向可以借鉴ITRF的思路,即选择在指定参考历元相对于最新版本ITRF的旋转参数及旋转速率为零来实现。
此外,99瑞典参考框架(Sweden Reference Frame 99, SWEREF99)是ETRS89在瑞典的具体实现[37],由瑞典国家永久GNSS参考站网(Swepos)的30个A级基准站及100个周边国家基准站定义(https://www.lantmateriet.se/en/geodata/gps-geodesi-och-swepos/reference-systems/three-dimensional-systems/SWEREF-99/)。最新的2011日本大地基准(Japanese Geodetic Datum 2011, JGD2011)于2011年10月启用,取代地心三维大地基准JGD2000[38],为日本领土及领海提供坐标基准(https://www.gsi.go.jp/ENGLISH/page_e30030.html)。
3 顾及非线性变化的1毫米级坐标框架构建
从上述全球/区域大地测量坐标框架的现状可以看出,最新的国际地球参考框架仍然无法满足大范围或全球尺度毫米级地球系统动态变化监测的需求。我国现行精度最高的CGCS2000坐标框架是一个区域、静态的框架,整体精度偏低,难以满足我国测绘基准现代化及北斗三号全球导航卫星系统对高精度坐标框架的需求。研究建立顾及非线性变化的1毫米级坐标框架,不仅是大地测量学21世纪的学科挑战,也是我国坐标框架建设的目标,对于构建和维持我国独立自主的北斗全球空间基准具有重要作用。其构建思路主要包括以下几个方面。
3.1 精密空间大地测量数据处理技术
高精度基准站坐标时间序列的获取是建立1毫米级地球坐标框架的必备前提条件。由于数据处理中涉及的地球物理与环境模型及策略不完善,导致空间观测技术类系统误差是限制基准站坐标时间序列精度的最主要因素,会造成基准站“虚假”的非线性位移,与坐标时间序列包括的真实地球物理信号相互作用,从而降低后续坐标框架建立的精度。以GNSS技术为例,忽略高阶电离层延迟会造成GNSS时间序列南北方向的周年、半周年虚假信号[39],实施次分量海洋潮汐改正后多个GPS交点年谐波信号振幅减小[40]等。SLR技术类系统误差包括距离偏差和时间偏差[41]。其中,时间偏差对测站坐标的影响主要体现在水平方向,东西方向最大可达毫米级。距离偏差会造成测站坐标的错误估计,最终影响框架尺度及地心定义。而且,间歇性的距离偏差还会造成测站坐标时间序列“虚假”的跳跃,影响测站的速度。不同分析中心采用的不一致的距离偏差处理方式同样会影响组合产品的精度。针对VLBI技术,系统误差主要包括未模型化或模型化不完善的大气传播误差(对流层、电离层延迟)、测站原子钟误差、天线误差及地球物理模型误差。其中,VLBI天线热膨胀、机械变形及错误的地球物理模型建模导致的测站位置误差均可达毫米级[42]。从多技术组合的角度来看,部分并置站不同观测技术获得的年周期信号之间仍然存在不一致,局部连接差异显著。例如,在ITRF2014的建立过程中,41%的局部连接残差超过5 mm[4],主要是由于技术类系统误差的影响,如未估计SLR距离偏差及VLBI天线重力变形、大规模全球GNSS基准站分子网处理等[26, 42-44]。
研究空间大地测量技术精密数据处理方法,不断精化数学模型,确定最优GNSS/VLBI/SLR/DORIS数据处理模型及策略,并对全球基准站数据进行整网一致性重处理,有助于进一步消除或减弱技术类系统误差的影响,为1毫米级坐标框架的建立提供更准确的输入数据。例如,与实现ITRF2014建立的全球第2次数据重新处理计划(repro2)相比,全球第3次数据重新处理计划(repro3)在数据处理方法及模型方面取得了很大的改进,包括引入了IERS建议的平均线性极潮新模型、太阳辐射压模型、海洋潮汐模型FES2014b、高频EOP新模型及未指北跟踪天线相位中心改正旋转,更新了GLONASS-M/GLONASS-K卫星的地球反照率模型及Galileo卫星天线相位中心偏差,并建议各分析中心将所有行星纳入三体运动模型等(http://acc.igs.org/repro3/repro3.html)。研究确定数据处理模型变更导致的GNSS技术类系统误差量级及其对坐标时间序列的贡献,有助于削弱“虚假”位移对真实地球物理信号的淹没效应,为后续坐标框架的建立提供可靠“干净”的数据源。
国际激光测距服务(international laser ranging service, ILRS)分析常务委员会(Analysis Standing Committee, ASC)通过同时估计测站坐标及距离偏差,并引入新的质心目标特征(center of mass target signature)模型来消除系统误差的影响,有效降低了ITRF2020中单VLBI/SLR技术获得的系统尺度差异(0.22×10-9),精度较ITRF2014(1.37×10-9)显著提高[6],但是系统误差仍然存在,而且并未考虑时间偏差的影响。2019年起,国际VLBI服务(international VLBI service, IVS)引入新的平均极潮模型、高频EOP模型,银河系光行差模型及天线重力变形模型针对1979—2020年的全球VLBI观测数据进行重新处理,发现4种新模型造成的VLBI台站高程变化超过15 mm[45]。尽管如此,VLBI数据处理将观测到的河外射电源作为点源处理。最新的研究结果表明,射电源结构,也就是发射角分布导致的误差占全部VLBI误差方差的40%,与测站相关误差几乎为同一量级[46]。完善SLR距离偏差模型,并深入研究忽略时间偏差、射电源结构造成的全球SLR、VLBI台站“虚假”运动特征,有望大幅度削弱SLR/VLBI技术类系统误差,进一步提高坐标框架建立的精度。
3.2 基准站非线性运动建模
由于仅模型化了周年、半周年季节变化及大地震造成的测站震后形变,无法精确描述基准站呈现的复杂非线性运动,这是导致ITRF2014精度难以达到毫米级的一个主要原因。地球物理信号不仅表现为周年、半周年运动,同时呈现年间及短周期变化(如日变化)特征[47-48]。学者评估DTRF2014发现,水文负载改正能显著减小基于VLBI技术得到的尺度周年信号,同时改正非潮汐大气压及水文负载效应则能使地球参考框架原点的年周期信号消失,且大幅度减小基于SLR及VLBI技术获得的尺度年周期信号[23]。然而,新一代地球参考框架ITRF2020的建立仍然采用周年、半周年季节变化模型描述基准站非线性运动,并未将环境负载改正模型纳入最终产品[6]。而且,除环境负载外的其他地球物理因素同样可能造成基准站的非线性变化。例如,非季节温度变化造成的测站垂直位移某些区域可达3 mm[49],地表覆盖的土壤层同样可能导致基准站热膨胀效应呈现相位滞后[50],水位变化造成的孔隙弹性形变等。
针对地球物理效应造成的基准站非线性位移对于坐标框架建立的贡献当前仍然处于起步阶段,还有许多科学问题值得研究,主要体现在4个方面:①水文负载对于GNSS基准站垂直非线性运动的贡献差异较大,且与地表水平位移弱相关;②环境负载位移模型的建立忽略了地球的滞弹性特征及区域地理环境差异的影响;③非季节温度变化造成的水平地表位移及其他地球物理因素对基准站非线性运动的贡献目前尚无研究涉及;④融合地球物理模型建立坐标框架在地球形状中心参考框架下实施,造成了坐标框架定义参数及产品的扭曲。评估不同时空分辨率全球及区域水文模型的精度,建立精密水文模型,并与地表气压、洋底压力模型联合,建立全球GNSS基准站的精密环境负载位移时间序列;发展基于滞弹性地球模型及测站格林函数建立环境负载模型的新方法,构建顾及地表土壤厚度的全周期GNSS基准站三维热弹性应变新模型,并探索可能造成基准站非线性运动的其他地球物理因素;最后,融合环境负载、热膨胀及其他模型,建立毫米级基准站地球物理运动模型,并评估其对地球参考框架的贡献具有重要的理论价值,有望进一步提高产品精度,为实现1毫米级坐标框架建立提供方法支持。
3.3 空间大地测量技术组合
利用单一技术难以同时解决框架的原点、尺度精度问题及站点分布密度问题,综合GNSS、VLBI、SLR、DORIS等多源空间大地测量观测技术是实现1毫米级地球参考框架的重要基础。其挑战和难点主要表现为3个方面:首先,由于不同测量技术观测精度参差不齐、系统误差显著,而且组合时未顾及季节信号等问题,导致基于多源观测技术建立的地球参考框架精度受限;其次,多源观测技术的组合主要依赖地面并置站,并置站的数量、精度、全球分布、局部连接的质量等因素对技术组合解的精度也起到决定性的作用[51];第三,局部连接一般通过地面常规测量方法获得,与空间大地测量手段获得的结果存在难以精确解释的差异,同样制约着多技术组合的精度和可靠性,例如,参与建立ITRF2014的并置站中,仍有41%的局部连接残差超过5 mm,21%的局部连接残差大于10 mm[4]。
目前ITRF的建立还是基于坐标层面[4],而DTRF的建立也是基于法方程层面[22]。为了保证各技术之间的一致性,更严密的组合方法是基于不同空间技术原始观测值层面的组合,也就是说,采用同一数据处理软件、相同常数及误差改正模型,统一处理GNSS/VLBI/SLR/DORIS原始观测数据[13, 9]。而且,我国北斗三号全球导航系统目前采用的坐标框架仅用GPS数据建立[32],精度仅为厘米级,与ITRF框架还存在较大的差距。
识别并减弱GNSS/VLBI/SLR/DORIS等空间大地测量技术之间的系统误差, 确定不同技术及并置站的权重,实现在观测值层面的多源空间观测技术精密融合,可以为1毫米级地球参考框架建立提供技术支撑。具体可考虑从以下3个方面进行:①利用紧组合方式,实现同种技术内多分析中心基准站坐标时间序列组合解的精确获取;②深入研究技术间组合的输入/输出标准、模型、组合策略及不同技术系统误差的处理、定权方式,提高并置站和局部连接精度及可靠性,从算法上改进技术间组合方法;③对输入的技术内组合坐标时间序列进行分析,加入选取的符合要求的并置站和局部连接信息,消除系统误差,并对技术间组合的法方程进行基准定义,以获取最终的组合解。
考虑到SLR站不仅观测LAGEOS(laser geodynamics satellite)等激光测地卫星,还观测GNSS卫星和低轨卫星、GNSS卫星和低轨卫星还能同时作为GNSS和SLR的空间并置站,可以增加有效的共同观测量,为解决地面并置站空间分布不合理,且与空间大地测量手段结果存在差异这一问题提供可能的有效途径。
3.4 地心运动
基于空间大地测量观测技术确定的坐标框架原点实际上为分布全球的基准站观测网中心CN,在观测网均匀分布的情况下,CN近似于CF。然而,根据IERS协议,ITRS的原点应该位于CM,通常将CF相对于CM的变化定义为地心运动[15]。季节性地心运动也是目前建立与维持地心坐标框架的主要误差源。国际上主要采用网平移法或一阶形变法,利用单一或多源空间大地测量数据估计季节性地心运动[52-53],但是目前为止尚未提供一套公认的地心运动模型。
从造成地心运动的根本原因入手,联合多种空间大地测量技术及地球物理模型反演地心运动是今后的发展趋势,其结果可最大限度地综合利用不同技术的优势,以解决坐标框架原点定义与实现间的不一致问题,为地球动力学研究提供更高精度的地心参考框架。特别是随着精度的不断提高,GNSS技术对于反演地心运动的贡献将越来越大[10, 13, 53-54]。
3.5 其他问题
除上述关键问题外,坐标框架建立的精度还与基准站网中观测站点的数量、质量、选取、分布、均匀性及密度、数据处理方法等有关。在观测站点方面,对现有空间大地测量地面监测网进行升级改造,是建立1毫米级坐标框架的一项重要任务。一方面,现有IGS测站普遍存在近场相位多路径效应[55],而且部分SLR测站的仪器存在老化,在后续基准站网建立过程中应加强基础设施的升级。另一方面,应用于新一代坐标框架ITRF2020建立的GNSS基准站全球均匀覆盖程度较ITRF2014显著改善,测站数达到约1100个,但是在非洲北部仍然存在空白,而且南半球的VLBI/SLR站点仍然过于稀疏[6],从而限制坐标框架原点、尺度的精度,后续还需要继续加密南半球的地面监测网,尤其是扩大VLBI/SLR站点的全球均匀覆盖率。而且,作为未来导航、通信、遥感融合发展的平台,低轨卫星(low earth orbit, LEO)星座发挥着至关重要的作用,是当前和未来的发展趋势。发展以北斗为代表的GNSS、LEO等高中低轨异构星座与地面增强的空天地一体化数据处理新技术和新方法,并建立GNSS/LEO/SLR空间连接(space tie),有望弥补地面并置站全球分布不均匀及地面/空间观测技术精度不匹配的不足。
在数据处理方面,最新ITRF建立采用的数学模型仅涉及函数模型,随机模型假设为白噪声[6]。然而,国际公认的GNSS坐标时间序列的随机特征表现为白噪声加有色噪声,纯白噪声假设会导致基准站速度不确定度的过低估计[56]。不断优化与完善坐标框架数学模型,并将有色噪声引入技术内组合及技术间组合模型,是将来关注的一个研究方向。
当前全球用户依赖的GNSS卫星导航定位系统信号具备天然的脆弱性,易受环境阻碍或被干扰、劫持。随着量子物理的快速发展,无须卫星的量子定位导航技术正在稳步推进,未来将可以实现新型定位导航。量子系统运行的方式最接近宇宙,会比其他系统以10~1000的倍数变得更快、更精准、更灵敏、更强大。此外,GNSS采用铷原子钟和铯原子钟(微波钟)提供高精度时间频率基准(原子秒),其天稳定度在10-14~10-15,限制了卫星定位的精度在米级。随着以更高稳定度的光钟(基于光学频率的原子钟)为代表的下一代原子钟技术不断突破,例如我国光钟最高稳定度为10-18,国际水平已达10-19,科学界将于2026年重新讨论秒的定义[57]。光频比微波频率高4~5个量级,由此可增加频率的精度和带宽(携带信息的能力)。如果未来用光钟替代现在使用的铷钟和铯钟,卫星定位导航及坐标框架的精度可能将显著提升。
4 结束语
第一作者简介:姜卫平(1972-), 男, 博士, 教授, 博士生导师, 研究方向为卫星大地测量学理论与方法及工程应用研究。
E-mail: wpjiang@whu.edu.cn
点击文末“阅读原文”即可查看原文章。
荐读
测绘大讲堂分享 | 姜卫平主任:北斗赋能万物互联 助力数字经济发展
姜卫平教授:北斗让万物更互联
《慧天地》敬告
《慧天地》公众号聚焦国内外时空信息科技前沿、行业发展动态、跨界融合趋势,发现企业核心竞争力,传播测绘地理信息文化,为时空信息类相关专业学子提供日常学习、考研就业一站式服务,打造政产学研金服用精准对接的平台。
《慧天地》借鉴《读者》办刊理念,把时空信息领域的精华内容汇聚到平台上。我们高度重视版权,对于精选的每一篇推文,都会在文章开头显著注明出处,以表达对作者和推文引用平台版权的充分尊重和感谢;对于来源于网络作者不明的优质作品,转载时如出现侵权,请后台留言,我们会及时删除。感谢大家一直以来对《慧天地》的关注和支持!
——《慧天地》运营团队
投稿、转载、商务等合作请联系
微信号:huitiandi321
邮箱:geomaticshtd@163.com