【技术】基于CGCS2000新框架的地物坐标更新策略研究与分析

《测绘通报》2024年第3期，审图号：GS京（2024）0499号

基于CGCS2000新框架的地物坐标更新策略研究与分析

成英燕1, 常春涛2, 王虎1, 徐彦田1, 许长辉1

1. 中国测绘科学研究院，北京 100036;

2. 武汉大学，湖北 武汉 430079

基金项目：国家自然科学基金(42274044); 崂山实验室科技创新项目(LSKJ202205100; LSKJ202205105); 国家重点研发计划(2020YFB0505802)

关键词：CGCS2000 坐标产品动态更新 天文大地网精度提升 最佳分布密度

引文格式：成英燕, 常春涛, 王虎, 等. 基于CGCS2000新框架的地物坐标更新策略研究与分析[J]. 测绘通报，2024(3)：54-62. DOI: 10.13474/j.cnki.11-2246.2024.0310.

摘要

摘要 ：本文针对未来CGCS2000框架更新后，由此衍生的全国各类实体坐标成果的更新策略和方法进行分析研究，为未来CGCS2000框架维持在GPS布点、数据处理方法及达到的精度方面提供思路和应对策略。采用地形及面积与我国大致相同，且可以公开下载的国外密集的CORS站数据，设计、处理形成了研究所需的不同时长的GPS观测数据集及厘米级、分米级CGCS2000坐标数据仿真集。分别讨论了天文大地网CGCS2000坐标从分米级至厘米级的提升，以及基于CORS的CGCS2000新框架衍生的实体坐标转换的实施方法和策略。对转换到新框架下影响坐标转换精度的因素：格网内插方法、控制点密度、搜索范围、格网间距等进行充分试验分析，确定了：①天文大地点CGCS2000精度从分米级到厘米级提升的最大控制点间距为2.2°，相应为230 km，即全国需要布363个点；但顾及地方框架的自身维护，各省(市)至少有4个站，则全国需在2178个天文大地点进行GPS观测。②要实现新的CGCS2000框架下实体坐标间隔10年的动态更新，若布设985个CORS站，在我国东部和中部控制点间距分别为482 km、西部地区为107 km时可达到5 cm的精度；而全国布置2170个CORS站，对应我国东部、中部、西部地区控制点密度分别为118、374、75 km时可达到3 cm的精度。因试验区为地理近似而非完全一致，本文所得结论和实际会有差异，仅为我国框架衍生的各类实体坐标成果更新机制提供参考和借鉴。

正文

国际地球参考框架(ITRF)是目前国际上公认的影响最大、精度最高的全球地心坐标参考框架，目前已发布了14个版本[1-4]，最新发布的为ITRF2020[5]。CGCS2000地心坐标框架是基于国际ITRF97框架建立的，表达在2000历元的框架点坐标，对准的精度在3 cm以内[6]。CGCS2000框架包括2600多个20世纪90年代建立的站点[7-9]，除25个连续运行基准站(平均间隔600 km)外，其他均为定期复测的点，且分布在地震活动活跃的断层区域。由于板块运动，CGCS2000与随后发布的ITRF框架差异越来越大[10]，主要体现在与ITRF框架精细实现的坐标系的原点、尺度、方向的系统性差异，以及表达地物坐标的现势性和精度的差异[11]。CGCS2000启用后，我国又新增和建成了覆盖中国全域的3000个左右连续运行基准站(CORS)。由基于少数CORS站控制的CGCS2000过渡到完全由CORS站维护的CGCS2000新基准框架势在必行。为区别于现CGCS2000，本文以CGCS2000CORS命名新框架，主要讨论基于CGCS2000CORS框架衍生的各类实体坐标的更新和维护问题。

首先对新框架CGCS2000CORS与现CGCS2000框架基础设施的异同进行分析。现CGCS2000框架包括CGCS2000骨干网(即2000国家GPS大地控制网，本文以CGCS2000Old命名)和纳入CGCS2000的天文大地网(本文以CGCS2000Astro命名，即天文大地点采用经典大地测量与CGCS2000联合平差后得到的坐标)两部分组成。

(1) 更新框架CGCS2000CORS与CGCS2000Old之间只有少数几个CORS站重合(不到25个，因部分CORS站已不再维护)。

(2) CGCS2000CORS完全由CORS站组成，CGCS2000Old框架大部分站点为定期复测站(其中1000站不再复测)，两框架大部分站不重叠。

其次对基于框架衍生的实体坐标更新进行分析。众所周知，CGCS2000启用后，原历史测绘坐标成果到CGCS2000的转换是基于同时具有3套坐标(1954北京坐标系、1980西安坐标系、CGCS2000Astro坐标)的天文大地点生成格网改正数完成的[12-13]，实现由1980西安坐标系到CGCS2000坐标转换的平均精度为0.11 m[14]。这期间CGCS2000Old大部分站并未直接参与测绘坐标成果的转换。CGCS2000框架各网间相互关系如图 1所示。

图 1 CGCS2000框架各网间相互关系

图选项

CGCS2000Old框架点中只有336个点与CGCS2000Astro重合，两个网联合平差时，CGCS2000Old的336个点以全方差形式参与平差，使得336个重合点同时具有CGCS2000下的两套坐标：CGCS2000Old框架下的坐标(三网平差)及CGCS2000Astro框架下的坐标(两网平差)，这两套坐标的差异在厘米级。而只有CGCS2000Astro在后续测绘坐标成果转换中发挥了作用。因此在启用CGCS2000CORS后，如何将基于CGCS2000Astro转换的坐标成果升级变更为基于CGCS2000CORS基础设施的坐标成果并持续更新，是本文讨论的问题。

1 方案

基于CGCS2000CORS框架网实现CGCS2000坐标框架衍生的实体坐标成果的持续更新和维护过程中，天文大地网CGCS2000Astro在转换中起承前启后的作用。原历史成果基于天文大地网坐标成果完成了到CGCS2000的转换，当将来CGCS2000Old框架更新到CGCS2000CORS后，需要将基于CGCS2000Astro的实体坐标成果转换到CGCS2000CORS框架下。显而易见，要建立联系需要在选取的部分天文大地网CGCS2000Astro站上与CGCS2000CORS站进行同期GPS观测(如图 1(b)中的第2层)，以便得到两个网在同框架同历元下的坐标。但是需要观测多少点(即密度)、多长时间确定的站坐标才能起到很好的传递作用？众所周知，5万天文大地点损毁严重，尚存的可观测的站点不多，况且从成本考虑也不可能都观测，因此需设计一定密度和分布合理的天文大地点，并通过一些数学方法完成实体坐标从CGCS2000Astro到CGCS2000CORS的转换，并达到期望的精度，解决框架断层衔接问题。

要想实现过渡和转换，两网之间要有共同点或确定的转换关系。由图 1(a)可以看出，CGCS2000Astro与CGCS2000Old，以及CGCS2000Old和CGCS2000CORS基本没有共同点，同时两两之间不仅有技术手段的差异还有精度的差异。因此本文提出基于CGCS2000Astro框架的测绘成果整体精度的提升和过渡的思路和方法，如图 1(b)所示。基于下载的密集站点集，分别按不同格网大小设计观测站密度，并按以下两步进行坐标成果转换及升级：①将CGCS2000Astro框架精度提升到与CGCS2000Old同等精度。具体做法为，选择一定数量和分布密度的天文大地点(按GPS C级点的观测要求)与全国CORS站联测，获得天文大地点和CORS站在观测历元ti、ITRF瞬时框架下的坐标(如ITRF2020、历元2025框架)；经板块运动改正[11]分别得到天文大地点及CORS站在CGCS2000(ITRF97、2000历元)下的坐标，分别标记为CGCS2000ti→2000Astro＿GPS(即天文大地点上进行GPS观测经改正得到的CGCS2000下的坐标)及CGCS2000CORS(CORS站的CGCS2000坐标)，坐标精度可达厘米级；计算CGCS2000ti→2000Astro＿GPS与CGCS2000Astro坐标的差值，作为控制数据内插得到格网改正数模型；基于坐标差值格网模型可将已转换到2000历元下实体地物的CGCS2000的坐标精度提升到厘米级。②基于CORS网速度构建格网速度场，将步骤①获得的2000历元下厘米级精度的实体地物坐标，基于格网速度场内插获得实体地物速度，由此速度将地物坐标推算得到在新CGCS2000CORS框架下坐标。具体为，将步骤①得到的CGCS2000ti→2000CORS坐标与如10、20或25年后CGCS2000CORS在2010、2020、2025历元下的坐标CGCS20002010CORS、CGCS20002020CORS、CGCS20002025CORS分别形成不同的站点坐标差数集，除去各站随时间以线性速度运动引起的坐标差，即线性运动坐标差集，得到各站非线性变化坐标差值，即非线性运动坐标差集，以这两套坐标差数据集作为控制数据，分别形成两套格网改正数，对待定点进行线性运动及非线性运动改正，并分析坐标精度。

由于目前国内无法收集可公开下载的高密度GNSS跟踪站数据供本文进行框架精化研究，因此，选择可以免费获取、地形大致与国内相同的NOAA连续运行参考站网(The NOAA Continuously Operating Reference Stations Network，NCN)，下载2019年1月的NCN数据。经处理形成3套数据：

(1) CGCS2000(Astro＿GPS)坐标：NCN数据按C级网观测时长切分模拟，在天文大地点上进行GPS C级网观测数据，经处理得到站点在ITRF2014、2019历元下坐标，再经板块模型改正[15-18]及框架转换，得到ITRF97、2000历元坐标，精度为厘米级[19]。

(2) CGCS2000Astro坐标：在上述结果基础上加入分米级误差(来自我国天文大地网真实数据)，即模拟CGCS2000下天文大地网坐标。

(3) CGCS2000CORS坐标：对一个月的NCN数据进行处理，得到仿真CGCS2000CORS新框架坐标成果[20-21]。

以上述坐标数据集为基础，分析讨论我国现有CGCS2000实体坐标成果更新到新框架CGCS2000CORS所需站分布密度、更新周期、更新方法及可能达到的精度。

2 天文大地点CGCS2000坐标精度精化策略与分析2.1 控制点密度设计与算法实现

控制点分布密度会直接影响内插的格网改正数模型。为快速挑选出分析所需的控制点密度，本文设计了基于格网的控制点快速抽稀及站点优选算法。基本流程为：将区域按照经纬度划分成一定间距的格网，并根据站点坐标集获取所有点位的格网范围；通过格网遍历搜索抽稀点，以格网范围确定格网内坐标集，计算格网内所有点与格网中心点间的距离；采用冒泡法确定最靠近格网中心点的测站，即为抽稀点，此方法可按设计控制点密度实现控制点位分布均匀化。图 2为分别以0.8°、1.6°两个密度为例的抽稀效果。

图 2 基于格网的最优抽稀效果

图选项

抽稀点作为构建坐标差值格网模型的控制点，采用分析确定的内插方法计算各个格网点的坐标差值，在此基础上以板块为划分单元，以速度和坐标差为判断指标，采用监督聚类法对站点进行进一步优选，从而分析不同点密度达到的插值精度。最后利用双线性内插计算待定点坐标差值，分析建模精度。

2.2 影响格网模型构建因素试验分析

在控制点数一定的情况下，影响天文大地点精细化格网模型构建精度的因素主要有：①空间插值数学方法选择；②初始格网的间距；③控制点搜索范围。

采用1502个天文大地点在2000历元下的坐标差CGCS2000Astro＿GPS-CGCS2000Astro为初始控制网集，通过对这些点逐步进行均匀抽稀，得到不同密度下的控制点集，其他点为未知点作为精度检核，研究分析天文大地点精细化格网模型构建精度。

2.2.1 插值方法选择对模型精度的影响

本文主要比较分析反距离加权法、克里金法及基于Delaunay三角剖分的坐标插值方法[22-23]，并确定最优插值策略，以构建最优精细化格网模型。

首先按0.5°分布密度选择控制点，分布如图 3所示，其中固定点826个，待插点676个。格网模型为0.5°，采用3种插值方法进行插值试验。反距离加权和克里金格网点在内插时以1.0°范围搜索控制点，其内插结果如图 4所示。可以看出，3种插值方法中反距离加权法精度最高，达到厘米级精度的未知点点数占比为92.16%(两分量均达到厘米级的站点统计)，其中N方向达到厘米级精度的点数占比95.12%，E方向占比96.30%；受克里金插值算法限制，搜索范围内固定点较少时导致部分待插点周围4个格网点有部分插值失败，以至于无法进行双线性内插，使得达到厘米级精度的未知点点数占比略低，为88.91%，其中N方向达到厘米级精度的点数占比91.86%，E方向占比92.01%；基于三角剖分的线性插值法在内插时，当格网点未包含在三角网内或外部三角网时，导致部分未知点周围4个格网点有部分未被内插，以至于无法进行双线性内插，经精度检核达到厘米级精度的未知点点数占比89.64%；其中N方向达到厘米级精度的点数占比91.86%，E方向占比92.75%。后文只针对反距离加权法进行讨论。

图 3 控制点间距0.5°格网1.0°搜索范围站分布

图选项

图 4 控制点0.5°间距插值精度对比

图选项

2.2.2 搜索范围对模型精度的影响

在进行反距离加权插值时，设置搜索范围会直接影响插值的精度，设计0.5°~2.5°的不同控制点密度，分析与最优搜索范围的关系，其结果如图 5所示。可以看出，最优搜索范围与控制点位间距之间存在一定的线性关系。

图 5 搜索范围对插值精度的影响

图选项

针对最优搜索范围与控制点位间距线性关系，采用最小二乘原理对其进行线性拟合，可得出关系为

式中，Tspacing为控制点位间距；Soptimal为最优搜索范围。

2.2.3 格网间距对模型精度的影响

格网间距也是影响插值精度的重要因素，试验分析分别设计站点间隔为0.5°~2.5°，按0.1°递增，分析最优精度时的格网大小，对应关系见表 1。可以得出，不同控制点间距内插的最优格网大小都在0.5°~0.7°。

表选项

2.2.4 控制点密度对精度影响分析

为获得最佳的插值精度，控制点密度越高越好。但控制点越密意味着观测量大且成本很高，因此需分析达到一定的精度条件下最优级的控制点密度。对控制点由0.5°至2.5°逐步抽稀点位间距，分析计算达到厘米级精度时，控制点的最大间距，结果如图 6所示。可以看出，控制点最大的点位间隔为2.2°，此时最优格网大小为0.7°。如果按站点的分布间距为实地约230 km，全国需要363个点，考虑各省至少有4个站，则全国需要2178个天文大地点。

图 6 控制点间距与精度关系

图选项

3 CGCS2000CORS框架产品动态更新策略与分析3.1 速度场数据源

在将天文大地点的CGCS2000Astro坐标成果的精度提升至厘米级的基础上，还需考虑历经20多年CGCS2000框架更新后，依赖于框架的衍生实体坐标成果的转换与更新。首先需考虑框架点的动态运动，重点考量的是站点的速度信息。由我国已有的框架点线性速度可以看出，我国板块运动的量级平均为3~4 cm/a[15]，按ITRF更新的频率5、10、15年，板块运动的影响量级分别为20、40、60 cm。分析不同站点分布密度下生成的格网速度内插待定点速度，由线性速度推算得到待定点坐标所能达到的精度，在此基础上进一步顾及非线性因素影响，以10年为例，分析更新框架坐标所能达到的精度，并确定相应精度国家CORS站应布点的站数及分布密度。

框架站点运动包括线性运动(板块的运动)及非线性运动，如潮汐负荷、大气负载、雪和土壤湿度、地下水开采等具有周期性变化及地震构造运动阶跃式跳变引起的运动[24]。由于本文格网建模点主要是基于框架中稳定点进行的，因此分析序列中不包括地震等因素引起的站点坐标跳变。站点在任意时刻的位置可表示为

式中，X0、 分别为参考历元测站的位置和线性速度；XR(t)是测站在任意时刻的位置；Xnon-linear为非线性因素引起的站点位移。

下载2009.0—2019.0的坐标时间序列及站点速度数据，首先根据站点的位置和板块信息进行站板块投影(如图 7所示)，对框架点速度插值进行分析。利用监督聚类法剔除粗差点，见表 2。

图 7 板块分区

图选项

3.2 框架线性维持精度分析

分别对5年(2014.0—2019.0)和10年(2009.0—2019.0)两个不同时间跨度的坐标时序进行线性框架维持精度分析。选取一部分站点的速度值作为控制数据，另一部分作为待定点，由已知控制点速度进行空间插值得到待定点速度值[17]，分别外推5年和10年得到坐标差值，与已知坐标差值进行比较，从而分析插值的精度。站点位移计算公式为

式中，XR(·)为相应年待定站点位置；Vmodel为空间插值所得速度值。

对于各板块分区，控制点密度在站点数据集中试验，设计抽稀间隔从0.5°至4.5°，以0.1°为步长均匀抽稀。采用反距离加权插值方法对待定点进行速度插值，分析由位于不同地形复杂度或不同板块区域控制点插值待定点速度，推算不同时间跨度坐标，与已知坐标比较精度分别达到5、3、1 cm时，控制点最佳密度。分析结果和达到相应内插精度时待定点的占比情况见表 3。

SW＿CONUS整个区域地形高程跨越0~2000 m，地形起伏较大，E＿CONUS区域板块除西部小部分区域为山区、高程在2000~4000 m外，其余地势均较平坦，这两个区域基本代表两个极端地形。

计算分析可知，对于大部分地形较平坦的E＿CONUS板块，当最大控制站点间距为4.5°，由速度场外推10年，99%以上的待定点坐标精度仍可达5 cm，即可满足1∶500大比例尺实体坐标成果转换的精度要求。若要达到优于1 cm的精度或缩小外推年限(5年)，88%以上的待定点可达到优于1 cm的精度；或加大控制点密度，控制点间距为1.1°，80%以上的待定点优于1 cm，而未达标的点基本均位于山区。此区域可对应我国东部地区的地形，即对于我国东部大部分区域，采用本文方法10年框架维持可以达到毫米级的精度。

SW＿CONUS区域对应于我国西部板块活动复杂区域，在仅考虑线性因素的前提下，当控制点间距为1.8°时，由速度外推5年，90%以上的待定点坐标精度可达5 cm；控制点间距为1.1°时，80%以上的待定点可达3 cm。外推10年、控制点间距为1.0°时，80%以上的待定点可达5 cm；控制点密度进一步加大至0.6°时，70%以上的待定点可达3 cm。

NW＿CONUS对应于我国中部地区，当控制点间距为4.5°时，外推10年，98%以上的待定点精度可达5 cm；控制点间距为3.5°时，91%以上的待定点精度可达3 cm；控制点间距为0.8°时，80%以上的待定点精度可达毫米级。即对于我国中部大部分区域，采用0.8°控制点间距、10年框架维持可以达到毫米级的精度。

由以上分析可知，对于山区地形起伏大的地区，若只考虑线性速度场进行框架维持，即使外推5年也都无法使站坐标精度优于1 cm，同时若考虑10年框架维持且和其他区达到相同精度时，待定点的占比仍较低。

3.3 顾及非线性因素的坐标插值精度分析

在上述利用外推速度得到坐标差值的基础上，将站点真实的坐标差与外推速度计算的坐标差进行插值，即为考虑到非线性因素的影响下的坐标差值。以板块运动最为活跃的SW＿CONUS区域为例，分析插值的精度。以10年框架更新利用外推速度得到坐标差值的基础上，对非线性维持的站点位移进行分析，站点位移计算公式为

式中，XR(2019)为计算所得未知站点位置；X(2009)为2009年站点位置；Vmodel为空间插值所得速度值；Xnon-linear为非线性因素所引起的位移量。控制点上顾及非线性因素，以坐标二次差为控制，对待定点进行内插求得坐标的改正数，其结果与线性维持精度分别达到5和3 cm结果进行比较，如图 8—图 9所示。

图 8 SW＿CONUS区顾及非线性时5 cm精度占比

图选项

图 9 SW＿CONUS区顾及非线性时3 cm精度占比

图选项

可以看出，对于10年框架维持，当顾及非线性因素且控制站点间距为0.7°时，SW＿CONUS区达到5 cm精度的测站占比由84.80%提升至87.14%，提升2.34%；达到3 cm精度的测站占比由66.96%提升至75.02%，提升8.06%。

综合以上分析，考虑实体坐标产品更新周期为10年，则不同板块能达到的最优精度及达到此精度的最大站分布间距见表 4。达到最大精度时控制点与待定点分布情况如图 10所示。

图 10 10年框架维持各区点插值精度达到最优时控制点(大三角)与检核点(小三角)的分布情况

图选项

3.4 插值精度与高程变化及板块边界的相关性分析

不同插值模型对控制点所在空间特征要求不同。当控制点密度较大、空间特征较简单时，插值模型普遍能取得较高的插值精度(如E＿CONUS)。当控制点集不足时，反距离加权插值模型要求控制点集尽量分布于地形特征点上。为探究站点位移空间插值与高程变化和板块边界的相关性，绘制PB2002板块划分与试验地区地势如图 11所示。

图 11 PB2002板块划分与试验地区地势

图选项

对时间间隔5年和10年统一采用控制点密度4.5°间隔内插得到的站点速度外推坐标差进行精度分析，待插点精度如图 12所示。可以看出，对于地形起伏平缓、海拔较低区域的大多数站采用线性维持坐标精度均能达到优于1 cm的精度，坐标维持精度较差的点大都分布在西部海拔较高的区域，即站点坐标插值精度与站点高程相关。同时通过比较分析板块边界位置发现，站点坐标变化与其所在的板块位置息息相关。以全球板块详细边界模型(PB2002模型)[25]对全球板块边界进行分析，由图 11、图 12可知，误差较大点主要分布在PB2002板块划分的板块交界处，处于本试验区域划分的SW＿CONUS区，可以得出坐标插值精度与板块边界活动不稳定有关。

图 12 10年跨度控制点间距4.5°待插点精度

图选项

综合地区分析可知，对于地形起伏平缓、海拔较低的区域，考虑线性因素并采用合理的控制点位间距进行插值时精度可达毫米级；区域位于板块内部且地形起伏较大时，考虑线性因素并采用合理的控制点位间距进行插值时精度可达3 cm；区域位于板块边界区域且地形起伏较大时，需同时考虑线性与非线性因素，并采用合理的控制点位间距进行插值时精度可达3 cm。

4 结语

本文选择地理空间特征与我国相似且有上千CORS站数据可供公开下载的国外区域作为研究对象，设计了基于不同测站分布密度、不同GPS观测时长获得坐标基准，对采用不同插值方法、不同格网间距、不同控制点探索范围对实体坐标精度升级和更新可能达到的精度进行分析，为我国坐标框架升级继而进行实体坐标动态更新制定应对策略提供依据。通过分析可知，反距离加权为天文大地点精细化格网模型构建的最优插值算法，插值的最优策略为，设置格网大小为0.7°，且控制点最优搜索范围与控制点位间距有线性关系。同时在天文大地点厘米级精度实现时，最大控制点间距为2.2°，对应实地230 km，全国需要363个点，考虑便于各省自身框架维护，至少需要4个站，全国需设2178个天文大地点。而对于基于框架点实体坐标产品的动态维持，将上述区域相应空间地理特征匹配映射我国相应区域：E＿CONUS对应东部地区，NW＿CONUS对应中部地区，SW＿CONUS对应西部地区。以时间间隔10年更新分析，根据不同区域地理特征不同，精度分别达到最优和5 cm时，CORS站布站策略见表 5。

对应表中最优精度，全国需布置CORS站2170个。具体到我国东部地区，控制点间距为118 km时，80%以上的点可达毫米级的精度；中部地区控制点间距为374 km时，90%以上的点可达3 cm的精度；西部地区控制点间距为75 km时，75%以上的点可达3 cm的精度。

框架点更新时间间隔10年，若需达到5 cm的精度，全国需布设CORS站985个。东部、中部地区控制点间距为482 km时，90%以上的点可达5 cm的精度，西部地区控制点间距为107 km时，85%以上的点可达5 cm的精度。值得说明的是：本文虽然讨论和分析的最长时长为10年，但30年坐标变化趋势仍可分为线性和非线性运动，根据我们对国内IGS站的时序分析结论，外推到目前仍然成立；同时考虑上述试验是针对国外类似地形区域进行的，与我国实际地形差异大的主要为西部山区，我国CORS站有在海拔5000多米的地方，而美国最高的地方只有4000多米。因此西部地区的布站相同时，精度会比本文预期的精度差。但分析方法及上述结论仍可对我国实体坐标产品更新、坐标框架改正模型构建所需控制点分布密度和更新策略有较好的借鉴作用。

作者简介

作者简介：成英燕(1964—)，女，硕士，研究员，主要研究方向为CGCS2000框架构建与维持、GNSS数据处理。E-mail：yycheng@casm.ac.cn

通信作者：常春涛。E-mail：chuntaochang@whu.edu.cn

社群交流/原创投稿/ 商务合作

来源：智绘科服