投稿问答最小化  关闭

万维书刊APP下载

基于重力场产品的水下匹配导航发展现状及影响因素分析

2023/7/10 14:19:07  阅读:35 发布者:

以下文章来源于经纬石旁话遥测 ,作者万晓云

本文改编自学术论文

《基于重力场产品的水下匹配导航发展现状及影响因素分析》

已刊载于《武汉大学学报(信息科学版)》2023年第6

万晓云1  吴云龙2   郭恒洋1 李明1

1. 中国地质大学(北京)土地科学技术学院,北京,100089

2. 中国地质大学(武汉)地理与信息工程学院,湖北 武汉,430078

万晓云

博士,副教授,主要从事重力场数据处理及应用研究。wxy191954@126.com

因高精度、无源性的特点,水下匹配导航技术已成为水下运载体导航的主要发展方向之一。分析水下匹配导航存在的问题,详细介绍了基于重力场产品的匹配导航发展现状,并对制约该技术定位精度的因素进行了分析。通过对比南海和太平洋研究区的重力异常、水深以及重力梯度特征参数,验证了重力辅助匹配导航技术更适合在重力特征变化明显的区域使用;当重力场产品格网大小为 l时,匹配导航的最高精度为根号6/6l

万晓云,吴云龙,郭恒洋,等 . 基于重力场产品的水下匹配导航发展现状及影响因素分析[J]. 武汉大学学报(信息科学版), 2023, 48(6):879-890. DOI:10.13203/j.whugis20220568

惯性导航是一种重要的水下导航方式,但是惯性导航系统误差会随着时间的积累不断增大,在水下运载体航行的过程中需要对惯性导航系统产生的误差进行不断地修正务的水下运载体不宜经常浮出水面校正误差。为了克服惯性导航的系统性偏差,水下运载体多采用组合导航方式,开展基于重力场信息辅助的水下匹配导航研究具有重要应用价值。

重力匹配辅助导航的基本原理是通过实时测量运动载体所在位置的重力场参数,并与事先观测和存储的高精度背景场图进行比对,从而准确确定当前载体所在的位置。目前,系统中用到的重力场特征参数主要为重力异常值或重力梯度值。重力梯度数据分辨率更高,不会受到水下运载体的加速度影响,能够反映更为精细的局部地球物理特征,因此重力梯度更适合作为匹配导航的背景场数据。然而重力梯度背景场数据的构建以及实时观测均比重力异常复杂。除了背景图的构建,数据之间的匹配算法也对匹配导航成败发挥了关键作用。

另一方面,受制于背景场特征的强弱以及仪器的噪声水平,基于重力场产品的匹配导航也并不总是有效,甚至在许多区域难以发挥作用,这也是重力场匹配导航仅是作为辅助导航手段的原因。然而随着量子重力仪、冷原子重力仪等高精度重力仪器的研制,构建超高精度的背景场数据以及实时测量高精度的重力场产品变得更为容易。在此种形势下,有必要对重力场匹配导航的发展现状及制约因素进行梳理与分析,这是本文工作的根本出发点。

本文主要对基于重力场产品的匹配导航发展现状如重力场测量仪器、匹配算法等进行总结分析,并对影响匹配导航的因素,如匹配量、适配区、重力场产品格网大小和精度等进行讨论,以期为该项技术未来的发展提供参考。

1 国内外研究现状

1.1

匹配导航系统研究现状

20 世纪 70 年代起,美国就已经开始了对水下无源导航技术的研究;80 年代初,Lockheed Martin 公 司 成 功 研 发 了 重 力 敏 感 系 统(gravity sensitive systemGSS),通过计算实时垂线偏差来修正惯性导航系统中因扰动重力而引起的误差;90 年代初,该公司研发了无源重力辅助匹配惯性导航系统,由 GSS、陀螺导航器、重力基准图和深度探测器一起配合工作;90 年代末,该公司 又 研 发 了 重 力 模 块(universal gravity moduleUGM),通过已有的重力异常数据、重力梯度数据以及地形数据与惯性导航系统的相关数据进行匹配,可实现水下重力匹配惯性导航。文献[10]提出通过多种插值方法得到高精度的海洋重力异常图,并在波罗的海区域进行测试;文献[11]总结了水下自主航行器(autonomous underwater vehicleAUV)匹配辅助导航的历史;文献[12]认为 AUV 的通信、定位和导航在将来面临挑战,并正在为新的、更复杂的应用开辟可能性。

随着精密测量科学的飞速发展以及处理测高卫星数据的理论和方法的成熟,中国在重力辅助匹配导航研究领域进步明显。国家地震局研制的 DZY-2 型重力仪、中国航天科技集团公司第九研究院第十三研究所研制的 SAG-2M 重力仪、中国科学院测量与地球物理研究所研制的 CHZ型重力仪等都已通过了实验阶段。通过对高精度地球重力场模型的研究,能够计算出高分辨率的重力异常数据和重力梯度数据,为建立高精度海洋重力异常基准图和重力梯度基准图奠定基础,促进重力辅助匹配导航技术的研究。国内多所高校与科研单位先后推出了一系列高阶地球重力场模型,为制备高精度海洋重力异常基准图和重力梯度基准图提供了条件。近年来,国内部分单位已开展了重力匹配导航的海上试验,如在国家自然科学基金的支持下,中国科学院测量与地球物理研究所的王虎彪团队开展了两次海上匹配导航试验,定位精度达到 2~3 海里;中国船舶集团有限公司李晓平团队在南海开展海洋重力匹配导航试验。以上研究表明,中国在重力场辅助惯性导航方面已取得可喜的进步。

1.2

重力场测量仪器研究现状

重力辅助匹配导航的精度与背景场精度有关,也与导航时的重力场数据实时测量精度有关,因此重力场测量仪器的测量精度对匹配导航至为关键。

近年来,中国的重力仪器研制发展迅速。2017 年第十届全球绝对重力仪国际比对中,国内多家单位如华中科技大学、清华大学、中国科学院测量与地球物理研究所等研发的仪器表现不俗 ,其 精 度 达 到 或 接 近 国 际 先 进 水 平。2022年,中国科学院精密测量科学与技术创新研究院研制的 CHZ-Ⅱ型海洋重力仪在南海测试,结果中重力仪零漂为 1.05 mGal/月,交叉点内符合中误差为 0.68 mGal,平差后优于 0.2 mGal。上述进展表明,中国的地面及海洋重力仪器研制已经达到国外同类仪器水平。值得注意的是,在国产仪器的工程化应用以及新体制的重力仪器研制方面,中国与国际领先国家还有不小差距,很多仪器还处于样机阶段,对于运动平台的大范围高精度观测等也仅处于研发阶段,这些是相关仪器用于匹配导航必须要解决的问题。

此外,重力梯度测量也是一项困难的大地测量任务。从第一台扭秤式重力梯度仪的诞生到量子梯度仪,重力梯度仪经历了多次重大改进与创新。国外部分具有代表性的重力仪和重力梯度仪信息统计见表 1

1 重力仪和重力梯度仪信息统计

注:1 E = 10-9/s2

GOCE 搭载的静电式重力梯度仪的极限测量精度可达到 0.2 mE Hz ,但受制于生产精度与测量噪声,其实际在轨精度仅为 10 mE Hz。寻求新型的测量载荷是未来提升重力场探测空间分辨率的重要方向之一。目前先进的重力梯度仪还有超导重力梯度仪和原子干涉重力梯度仪,超导重力梯度仪的精度可达 1 mE Hz ,但其依赖于低温环境,且技术成熟度尚未达到工程化;量子重力梯度仪利用超低温原子团代替惯性质量来感应重力作用,通过探测原子的内态布居数来获得重力梯度信息。华中科技大学研制的重力梯度仪在采样率为 0.25 Hz 时,灵敏度达到 670 E/根号Hz;中国科学院武汉物理与数学研究所在“十二五”期间研制出垂向原子重力梯度仪和水平原子干涉重力梯度仪,测量精度分别达到 7.5 E7.4 E1 E = 10-9/s2);吉林大学开展了多种类型的重力梯度仪研制,并取得显著进展,其研制的基于冷原子干涉原理的重力梯度仪 分 辨 率 达 到 32 E@200 s(移 动 平 台 环 境),MEMSmicro-electro-mechanical system)重 力 梯度仪样机分辨率优于 30 E@120 s,超导重力梯度仪噪声水平达到 0.87 mE/ 根号Hz

需要说明的是,不同文献所提供的重力场测量仪器的精度测试方法、测试环境等均有所差异,难以直接比较。对于不同的观测环境和观测平台,重力仪和重力梯度的精度水平表现也会不同。整体而言,中国重力梯度仪器的研制水平与国外相比差距要大于重力仪器的研制。未来,随着各种高精度仪器的广泛应用,重力场相关数据的信息更加精确,导航背景场信息进一步丰富,匹配导航的适用性有望增强。

此外,除了仪器本身的测量精度外,观测环境给重力场高精度观测带来的干扰也是实际应用需要重点关注的问题。这是因为水下匹配导航载体往往处于运动状态,这给重力数据的高精度获取带来了不利影响。因此需要研究和攻克水下移动平台的高精度重力测量及数据处理技术。文献[44]建立了水下移动重力测量误差模型,给出了 1 mGal 精度重力测量对水下定位设备的性能要求,相关研究表明,水下高精度移动测量对水下观测平台、数据处理算法等均提出了较高要求。尽管如此,国防科技大学、海军工程大学等国内多家单位均进行了水下测量试验,表明中国在该领域取得较大进步。

1.3

匹配导航算法研究现状

除了高精度重力观测仪器以及高精度匹配导航背景场数据以外,重力辅助匹配导航技术的另一个关键核心就是匹配算法。国外对匹配导航算法的研究比较早。20 世纪 70 年代,美国学者提出地 形 轮 廓 匹 配 算 法(terrain contour matchingTERCOM)和圣地亚惯性地形辅助导航算法(Sandia inertial terrain aided navigationSITAN),并在飞行试验中得到应用。虽然这两种算法针对地形匹配提出,但同样适用于重力匹配导航。20世纪80年代,克拉索夫斯基提出了重力场数据的相关极值和匹配算法。1999 年,Kamgar提出了最近点迭代等值线算法(iterative closest contour pointICCP)。近 年 来 ,Allotta 等提 出 了 一 种 专 为AUV 水下导航设计的卡尔曼滤波器(unscented Kalman filterUKF),并实现了海上测试;MehdiMohammad针对 AUV 提高水下导航系统的性能,提出一种基于 H-无穷大滤波器的集成导航系统;Parth通过研究重力异常参数与颗粒过滤器SLAM 之间的关系,为 AUV 提供了最佳路径,并实现目标位置的最小定位误差;Elmokadem等为AUV 的动态定位和航点跟踪开发了一种鲁棒控制方案;Salavasidis等等提出了地形辅助导航定位方法,使用测深测量来限制惯性导航误差的增长。

国内在重力辅助匹配导航算法领域的研究取得了丰富的成果。许大欣将 SITAN 算法应用于重力辅助匹配导航中,利用重力异常数据作为匹配导航的背景场数据。闫利等通过仿真试验分析了 TERCOM 算法对不同重力场分布区域的适应性。柴华等基于四元数和卡尔曼滤波,设计了一个 12 状态滤波器,并在地固系下完成了惯性导航系统的初始精对准。魏二虎等针对重力辅助匹配导航提出了一种改进的 TERCOM 匹配算法。武凛等通过引入平均重力差(average gravity differenceAGD)参 数 ,表 示 某海域重力变化水平,研究发现在南中国海,北印度 洋 等 地 区 重 力 匹 配 导 航 更 适 用 。欧 阳 明 达等提出了基于改进 A*算法的重力辅助匹配导航算法,利用海域重力异常和粗糙度等指标对该海域适配性进行评价,有效提高了匹配导航路径规划的可行性。邹嘉盛等提出一种基于约束条 件 的 重 力 匹 配 导 航 算 法 ,该 算 法 基 于 PNN probabilistic neural network)匹 配 算 法 ,利 用 了惯导提供的轨迹方位和航距信息,提高了算法的匹配速度和匹配精度。刘慧等提出了一种新的基于多尺度的重力匹配导航算法,测量精度和匹配要求较低的同时表现出较高的导航精度、匹配成功率和匹配效率。王成龙等提出了一种基于 Delaunay 三角测量的匹配区选择算法。该算法建立了三维模型,通过特征提取和聚类分析来选择匹配区域。与传统算法相比,该算法选择匹配区域更有效。Wan Yu指出,相比于重力梯度张量数据,重力梯度不变量数据作为匹配导航背景场数据时,也具有较高的匹配精度,特别是当重力梯度张量数据存在姿态误差时,采用重力梯度张量不变量作为匹配导航背景场数据更具优势。更进一步地,于锦海等提出在水下导航定位中引入不变量坐标系,利用观测得到的重力梯度和重力梯度不变量可以直接构建关于观测点的方程组,该算法稳定,能较好的对点位坐标进行恢复。欧阳明达等针对 SITAN 算法在粗差探测方面的不足,引入抗差估计方法,实验表明,抗差估计方法用于水下重力匹配导航,改善了导航精度,增强了可靠性。欧阳明达等提出了将粒子滤波与粒子群优化算法相结合,提出了改进粒子滤波的水下重力匹配导航算法,有效提高水下重力匹配导航精度。

1.4

存在的问题及发展趋势

从国内外重力匹配辅助导航技术的研究现状来看,国外在仪器设备、匹配导航算法等领域的研究均远早于国内,在重力辅助匹配导航研究领域的技术手段和理论方法相对成熟,有些研究通过了水下实验阶段,已经开始为潜艇等水下航行器服务。相较于国外,国内对重力辅助匹配导航的研究起步较晚,特别是高精度的实时重力测量仪器的研发,与国外存在着一定的差距,但是随着近年来中国学者的不断努力,在高精度重力辅助匹配导航的高精度背景场数据制备、重力辅助匹配导航算法以及高精度的重力场测量仪器等方面都取得了丰硕的研究成果。

从当前的研究现状来看,研发出高精度的实时重力场数据测量仪器,特别是成本低、精度高的仪器,将会推动重力辅助匹配导航技术的发展。随着卫星测高技术和重力卫星的发展,计算出高精度全球重力场模型成为当前的研究热点。一方面,通过单一的测量手段获取海洋重力场数据已经无法满足匹配导航背景场数据精度要求,融合多种高精度重力场数据才可制备出高精度重力场辅助匹配导航背景场数据。另一方面,寻找新的方法和手段对已存储的重力数据进行插值,提高重力数据的近似精度也可以制备出满足精度要求的匹配导航背景场数据。

关于重力匹配辅助导航算法的研究,前人已经做过了大量的工作。总结前人的工作,不难发现重力匹配辅助导航算法的工作效率、可靠性以及其适用性是评价一种匹配导航算法是否优异的重要标准。而提高重力辅助匹配导航算法的效率,一方面需要深刻理解和掌握所运用的重力辅助匹配导航背景场数据的特征。前人的研究中所运用的背景场数据大多是重力异常数据,基于重力梯度数据的研究相对较少,而基于重力梯度不变量数据的研究则更少。重力梯度数据具有 5 个独立分量,受地形因素影响较大,理论上更适用于重力场辅助匹配导航,但需要研究与之相适应的匹配算法。另一方面,各种传统的重力辅助匹配导航的算法有其自身的优点和缺陷,如ICCP 算法不适用于初始误差较大的情况、TERCOM 算法实时性比较差等。因此,进一步探索和研究出新的适用于不同环境的重力辅助匹配导航算法,如人工智能算法等,也将会成为未来研究的重点。

综上所述 ,基于重力场产品的水下匹配导航还存在如下问题:(1)精度高、成本低的重力场测量仪器的研制和工程化应用有待加强,特别是高精度全张量重力梯度仪 ;(2)对 于 匹 配量的挖掘和合理利用有待进一步加强,例如不同 量 的 组 合 量 等 ;(3)匹配算法 也需要进一步研究,如随着背景场精度和实时测量精度的不断提升,是否能够实现基于引力场信息的空间位置坐标直接解算以及人工智能算法的应用等。

2 讨论与分析

2.1

匹配量的比较

选择研究区为 112°E~114°E17°N~19°N 的中国南海海域,使用 XGM2019e_2159 模型球谐展开到 2159 阶,计算出中国南海研究区域 1′×1′格网的重力异常及扰动重力梯度的 6 个分量TxxTyyTzzTxyTyzTxzTxy = TyxTyz = TzyTxz =Tzx),结果如图 1、图 2 所示。

1 中国南海研究区域重力异常

2 中国南海研究区域扰动重力梯度

在重力梯度张量基础上,计算扰动重力梯度 张量不变量,计算式为:

式中,I1I2I3为重力梯度不变量,其中 I1作为评定球谐展开计算扰动重力梯度的精度指标,其理论值应 0。按式(1)计算出I1I2I3,结果如图 3所示。

3 中国南海研究区重力梯度不变量

由于I1的理论值为 0,本次计算得到的值量级为 1024,主要反映了计算误差,因此在图 3 中并未展示不变量的I1分量。对比图 1~3 可知,不同的重力场产品在同一区域呈现的特征是显著不同的,如图 2 和图 3 中极值点的位置不同,表明不同分量的重力场产品对于载体位置的敏感度不同;同类型的重力场产品的特征有一定的相似性,如梯度分量 Tyy Tyz 的纹理具有一定的相似性,梯度不变量 I2I3的极值点的分布也较为相似。以上分析表明,不同的重力场产品分量对于载体位置均具有一定的敏感度,理论上均可用于匹配导航;从特征分析表明,不同分量的匹配效果会有所差异。除了上述讨论的匹配量外,也有学者提出采用两两采样点间的重力异常之差作为新的匹配量,仿真试验表明,重力异常之差匹配量的精度优于绝对重力异常匹配量。事实上,重力异常之差可反映重力异常梯度信息,因此该匹配的优越性间接证明了重力梯度用于匹配导航的巨大潜力。

2.2

适配区的选择

为了分析区域特征对匹配导航的适应性影响,选择研究区域为 95°W~97°W5°S~7°S 的太平洋海域进行对比分析。为了反映区域特性,表 2 统计对比了中国南海研究区与太平洋研究区的水深数据。

2 不同研究区域水深统计/m

由表 2 可知,南海研究区的水深变化剧烈,表明该区域的海底地形复杂;太平洋研究区的水深变化量较小,表明该区域的海底平缓。使用XGM2019e_2159 模型球谐展开到 2159 阶,计算太平洋研究区 1′×1′格网的重力异常、重力梯度及其不变量,结果分别如图 4~6 所示。为了便于对比,已将图 4~6 的色标范围与图 1~3调整一致。

4 太平洋研究区域重力异常

重力场的特征参数是影响匹配精度的重要因素,常用的重力梯度特征参数包括重力梯度分量和不变量的均值、标准差以及粗糙度等。重力梯度标准差反映了局部重力梯度数据的离散程度;粗糙度反映的是局部重力梯度张量的平均光滑程度。特征参数的值越大,表征重力梯度信息越丰富,越有利于匹配。统计南海和太平洋研究区域扰动重力梯度分量和梯度不变量,特征参数见表 3

对比表 3 中重力梯度 6 个分量及不变量的均值和标准差,发现南海研究区的重力梯度数据分布更加离散;对比绝对粗糙度,发现南海研究区的背景场起伏程度更大;对比经纬度方向的粗糙度,发现重力梯度数据的特征在经纬度方向上的较大差异。通过以上特征参数分析,结合表 2 两个研究区的水深变化,发现南海研究区的地形比太平洋研究区更复杂,南海研究区的重力数据信号变化更加剧烈,更适合重力辅助匹配导航。

不同的背景场特征对匹配算法提出了不同的要求,因此在背景场特征明确的前提下,需要研究匹配区的特征性差异对不同匹配算法的要求,如在背景场特征突出的情况下,普通算法就可能实现较高的精度;而在背景场特征不太明显的情况下,需要研究更优的算法;在特定方向特征明显的情况下,需要研究相适应的匹配算法。

5 太平洋研究区域扰动重力梯度

6 太平洋研究区域重力梯度不变量

3 不同研究区域重力梯度分量特征参数值

2.3

重力场产品格网大小的影响

单从匹配导航本身来讲,背景场格网的大小意味着最终的匹配点位于格网任意点的概率一样,因此定位误差可认为是个二维均匀分布。图 7为均匀分布示意图。

7 均匀分布示意图

假定在 x 方向,格网的边界点坐标是 ab,则匹配点位于 a b 之间的概率 f ( x) 1/(b - a),据此可推导其落点 P 的均值 E 和方差 V,计算式分别为:

事实上,以上结果对应的是一维均匀分布的均值和标准差。对于匹配导航而言,可以将 xy的定位误差看作两个一维分布。根据以上结论,定位结果(xy)的均值将位于格网的中心点,其方差的计算式为:

l= b - a为格网边长,则匹配导航定位的标准差为 根号6/6*l。格网大小为 1~10 km 对应的匹配导航定位误差标准差见表 4。由表 4 可知,在匹配导航时,格网越小(即重力场分辨率越高),匹配导航定位误差越小,标准差略小于格网边长的一半。

4 不同格网边长对应的误差标准差/m

2.4

重力场产品精度的影响

在分辨率提高的同时,也需要保证重力场产品的精度。本节以重力异常为例,分析重力场产品精度对匹配导航的影响 。由于

则重力异常误差 δΔg对导航的定位误差可表示为:

由式(5)可知,重力异常精度越高,则匹配导航的精度也越高,但对于不同区域,相同的重力异常产品精度,其匹配导航的精度也将不同,这主要取决于匹配区域的重力场变化程度。事实上,式(5)中 TzxTzyTzz即量化表征了重力异常变化的剧烈程度。根据图 2,在南海研究区内,上述值的绝对值均值分别为 10.34 E6.64 E12.47 E;而在太平洋区内(根据图 5),上述值的绝对值均值分别为 2.18 E2.42 E3.26 E。因此,若重力异常的误差为 1 mGal,则其在南海研究区 3 个方向产 生 的 定 位 误 差 分 别 为 966.81 m1 505.85 m802.15 m;在太平洋研究区产生的定位误差分别为 4 597.16 m4 134.12 m3 069.49 m。反之若 3个方向的目标定位精度均为 500 m,则重力异常的精度在南海研究区需要达到 0.332 mGal,在太平洋研究区需要达到 0.109 mGal。因此,在同等定位精度要求下,重力场变化剧烈的区域对重力场产品的精度要求要更低;在同等重力场产品精度下,重力场变化剧烈的区域定位精度更高。

需要说明的是,本文的研究区仅是南海的一个局部区域,因此以上结论并不适合整个南海区域。此外,以上结论也仅适用于基于重力异常的匹配导航。若匹配量为绝对重力,结论也将不同,原因是地球表面绝对重力的梯度值超过 3 000 E,若按式(5),则 1 mGal 的重力变化可敏感约 3 m的位置变化,显然这是目前的匹配导航技术难以达到的。一方面是因为绝对重力的背景图难以获得,而重力异常图可以通过测高卫星间接获得;另一方面,绝对重力测量比相对重力测量更加困难,但随着仪器研制技术的进步,如高精度冷原子绝对重力仪等,未来的重力场测量有望在导航领域发挥巨大作用。

3 结语

凭借重力场数据无源性和稳定性等特点,重力辅助匹配导航逐渐成为水下运载体导航定位的主要手段之一。本文总结了国内外在匹配导航领域的发展现状,分析了存在问题,并进一步以南海和太平洋部分区域为例,从不同角度分析了影响匹配导航的因素,包括匹配量、适配区、重力场产品格网大小、重力场产品精度等。通过对比两个研究区的重力场特征参数,以及同等定位精度要求下对重力场产品的精度要求,发现复杂海域的重力信号变化更剧烈,匹配导航更适用于重力异常变化明显的海域。

虽然匹配导航精度受制于背景场格网大小、重力场产品精度以及匹配区域的重力场特征等诸多因素,但随着高精度重力仪器的成功研制,特别是基于量子技术的高精度重力及重力梯度仪器的研制,将有望为匹配导航提供更准确的数据,从而极大促进匹配导航技术的发展。另一方面,随着航天技术的发展,基于卫星观测数据反演得到全球海洋重力场产品的精度越来越高,中国在该领域已经取得快速发展;此外,美国于2022 年底发射了 SWOTsurface water and ocean topography)卫星,也将对海洋重力场产品精度和分辨率的提升发挥巨大作用,这些进展将极大促进匹配导航背景场数据库的建设。在重力仪器和背景场数据显著提升的情况下,高效的匹配算法也将是未来研究的重点,人工智能技术有望在该领域发挥重要作用。匹配导航最终受多种因素的联合影响,除了以上单一因素的进展之外,未来还需要对各因素之间的耦合效应开展研究,才能最终使得匹配导航达到理想的效果。

转自:“测绘学术资讯”微信公众号

如有侵权,请联系本站删除!


  • 万维QQ投稿交流群    招募志愿者

    版权所有 Copyright@2009-2015豫ICP证合字09037080号

     纯自助论文投稿平台    E-mail:eshukan@163.com