面向多源卫星测高数据的动态窗口约束交叉点算法

毛瀚辰 ,  郑伟 ,  张文松 ,  沈祎凡 ,  吴沐府 ,  田万辰

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (30) :12902-12914

PDF (14531KB)
科学技术与工程 ›› 2026, Vol. 26 ›› Issue (30) :12902-12914 DOI: 10.12404/j.issn.1671-1815.2601072
论文：天文学、地球科学

面向多源卫星测高数据的动态窗口约束交叉点算法

作者信息 +

Dynamic Window Constrained Crossover Algorithm for Multi-mission Satellite Altimetry

Author information +
文章历史 +

摘要

为解决海量多源卫星测高数据交叉点定位难以兼顾精度与效率的问题,通过多项式拟合与轨道对称性原理,构建了一种新型动态窗口约束交叉点算法。该方法结合空间预筛选与计算机图形学逻辑构建多层级约束机制以消除冗余计算,通过动态限制计算窗口边界策略将解算约束在极小固定范围,提升运算效能并规避了长弧段拟合导致的几何偏移,且利用多种测高卫星数据对该方法的适配性与定位精度进行了验证。结果表明:在计算效率方面,新型动态窗口约束交叉点算法较传统分段拟合法提升至少50%;在定位精度方面,有效克服了几何偏移与漏算现象,实现了点位精准计算;在平均海面高度建模实验中,动态窗口约束交叉点算法计算结果与SDUST2020及CLS15模型具有高度一致性。其中,Jason-1 GM数据精度最高,其均方根误差优于0.08 m;而CryoSat-2与Sentinel-3A的RMSE经证实由几何基准差异引起,而非算法误差。可见,动态窗口约束交叉点算法在保证高精度的同时显著提升了计算效率,为海面高度的精密获取提供了高效手段。

Abstract

To address the problem that traditional crossover positioning methods for massive multi-source satellite altimetry data struggle to balance accuracy and efficiency, a novel DWCC(dynamic window constrained crossover) algorithm was developed by integrating polynomial fitting with the principle of orbital symmetry. A multi-level constraint mechanism, combining spatial pre-screening and computer graphics logic, was constructed to eliminate redundant computations. The solution was constrained to an extremely small fixed range through a dynamic window boundary strategy, thereby improving computational efficiency and avoiding geometric offsets caused by long-arc fitting. The adaptability and positioning accuracy of the proposed method were rigorously validated using multiple satellite altimetry datasets. The results indicate that the computational efficiency of the DWCC algorithm is improved by at least 50% compared with the traditional piecewise fitting method. In terms of positioning accuracy, the geometric offset and missing point issues are effectively overcome, and precise crossover point positioning is achieved. In MSS(mean sea surface) modeling experiment, the crossover results computed by the DWCC algorithm are found to be highly consistent with the SDUST2020 and CLS15 MSS models. Among all datasets, the Jason-1 GM data are shown to achieve the highest accuracy, with RMSE(root mean square error) better than 0.08 m. The relatively larger RMSE values of CryoSat-2 and Sentinel-3A are confirmed to be caused by systematic differences in geometric datums rather than by algorithm errors. Therefore, the DWCC algorithm is demonstrated to significantly improve computational efficiency while maintaining high accuracy, providing an efficient and reliable technical means for the precise acquisition of sea surface height.

Graphical abstract

关键词

测高卫星 / 卫星轨迹交叉点 / 动态窗口约束交叉点算法 / 平均海面高度

Key words

satellite altimetry / satellite track crossovers / dynamic window constrained crossover algorithm / mean sea surface

引用本文

引用格式 ▾
毛瀚辰, 郑伟, 张文松, 沈祎凡, 吴沐府, 田万辰. 面向多源卫星测高数据的动态窗口约束交叉点算法[J]. 科学技术与工程, 2026, 26(30): 12902-12914 DOI:10.12404/j.issn.1671-1815.2601072
MAO Han-chen, ZHENG Wei, ZHANG Wen-song, SHEN Yi-fan, WU Mu-fu, TIAN Wan-chen. Dynamic Window Constrained Crossover Algorithm for Multi-mission Satellite Altimetry[J]. Science Technology and Engineering, 2026, 26(30): 12902-12914 DOI:10.12404/j.issn.1671-1815.2601072

登录浏览全文

4963

注册一个新账户 忘记密码

海平面是指在无潮、无波、无风的理想状态下,海洋所能保持的水平面。它是联系海洋、大气、陆地和冰川系统的重要参照面[1]。平均海面高度(mean sea surface, MSS)是通过长期观测得到的海面高度平均值,广泛应用于海洋学、地球物理学等领域的研究[2],也为海洋重力场变化率模型[3]、海洋重力异常反演[4]、气候变化与动态海面地形研究等领域[5]提供了重要的研究基础,也是开展国际海底区域资源勘查、提升深海高新科技创新能力及维护国家深海权益的关键保障[6]。
平均海面高度的计算主要依赖于卫星测高技术。自SEASAT(sea satellite)和GEOSAT(geodetic satellite)卫星首次实现全球海洋中尺度海平面变化的连续观测以来,已经连续发射了数十颗测高卫星[7]。对于海面测高研究,这些不同类型、不同任务的卫星数据为平均海面建模[8]提供了重要的数据基础。在进行海面测高过程中,测高卫星通过雷达技术实现海面点的高度测量,例如,Jason系列则通过多代卫星长期观测建立了全球海平面变化的高精度基准数据集[9],且在其寿命后期的大地测量阶段(geodetic mission, GM)通过改变轨道周期,获取了极高密度的全球空间覆盖数据;Sentinel系列卫星则凭借其高频次、高分辨率的观测能力,为现代海面高度监测提供了关键数据。再者,CryoSat-2卫星还利用合成孔径雷达进一步提升了更广阔区域的海洋监测能力[10];而以导航卫星反射信号GNSS-R为代表的微波遥感技术,凭借其信号源丰富、高时空分辨率等优势,也成为获取全球高精度海面高度及海洋重力场信息的重要补充手段[11]。与此同时,中国近年来在卫星测高领域取得了长足进步,构建了以海洋系列(HaiYang-2, HY-2)为核心的自主测高星座。HY-2系列卫星(如HY-2B、HY-2C、HY-2D)的成功运行,标志着中国已具备全球海洋动力环境的长效监测能力[12]。与国际主流的Jason系列卫星相比,中国海洋卫星在轨道设计上兼顾了重复周期与覆盖密度,为高精度平均海面高度建模提供了关键的国产数据源[13]。此外,高分七号(GaoFen-7, GF-7)等搭载的激光测高载荷,进一步提升了中国在复杂海域及极地冰盖高度测量中的自主可控水平[14]。
然而,利用这些不同任务、不同轨道密度的卫星数据获取高精度海平面的关键在于获取精确的地面卫星轨迹交叉点。交叉点是卫星轨迹上升弧段与下降弧段轨道的交会位置,在轨道重复周期内,上升弧段和下降弧段在交叉点处测得的海面高度参数应完全一致,但由于观测系统误差、仪器偏差、环境干扰等因素[15],实际数据在交叉点处常出现显著差异,这种不符值会对海平面变化监测的可靠性构成挑战,因而影响到海平面高度模型构建。因此,研究高精度的交叉点定位方法对平均海面高度建模精度和可靠性具有重要意义。
为了提升交叉点的计算精度,科学家已开展了广泛研究。在交叉点定位方面,Tai[16]于1989年提出一种FI算法(fixed iteration),该算法通过轨道拟合求解交叉点位置,但拟合精度不足导致结果偏差较大;周晓光[17]设计了11种分段类型的分段拟合法(piecewise fitting method)提升了计算精度与适用范围,但需要人工判断最优分段类型,存在自动化程度低的缺陷。Greene等[18]提出基于计算机图形学原理的RST算法(rapid rejection and straddle test)有效解决了自动化问题,但算法复杂度导致计算效率偏低;另一方面,为了提升计算效率,格网法将研究区划分为规则格网,并将同时存在上升轨道和下降轨道测点的格网中心点视为交叉点,因其计算耗时低而颇具吸引力。Zhu等[19]在其研究中发现该方法精度较差,影响海平面的高精度建模,纬度差法虽然也能减少部分计算量,但在处理大规模数据时仍然面临计算耗时长的问题等[20]。赵子龙[21]研究发现,当数据量较大或交叉点数量过大时,传统算法计算耗时会随数据量显著增加,进一步影响整体计算效率,限制了其在高频次、大体量卫星数据中的应用。
针对目前传统交叉点算法普遍存在的精度与效率难以兼顾的问题,现提出一种兼顾计算精度与效率的新型动态窗口约束交叉点算法。在多项式拟合与交叉点计算原理基础上实现一系列成果:首先,构建多层级约束机制,通过原始数据空间预筛选并结合计算机图形学中的快速排斥与跨立测试,从物理空间到几何拓扑逻辑层层剔除冗余计算;其次,创新性地提出动态定长计算窗口策略,动态调整计算区域并将解算过程严格约束在极小固定范围内,在显著提升矩阵运算效能的同时,有效规避长弧段拟合带来的几何偏移;最后,通过结合CryoSat-2、Jason-1 GM及Sentinel-3A多源卫星数据进行验证,证实该方法在不同轨道特性与海量数据体下均具有极高的适配性与定位精度。实验结果表明,DWCC算法在不同轨道特性和数据体量下均表现出极高的适配性,在确保交叉点定位高精度的同时,提升计算效率,为高效获取海平面高度提供可行的手段。

1 动态窗口约束交叉点算法的构建

卫星轨迹交叉点是指卫星在同一海域内,上升轨道与下降轨道在空间上的几何交点。对于CryoSat-2、Jason-1和 Sentinel-3A 等传统测高卫星,其观测数据以沿轨点列的形式分布。交叉点定位的核心任务是寻找两组离散时间序列在经纬度空间中的交会位置。在实际计算中,由于采样频率限制,观测点很难精确落在几何交叉点上。因此,通常需要对交点坐标进行求解。

在以往研究中,各种交叉点算法被广泛应用于确定交叉点的存在,但面对数据体量日益增长的卫星测高数据,这些方法处理能力已显露出明显的局限性。基于多项式拟合与交叉点计算原理,引入多重影响因子进行约束,在确保交叉点定位精度前提下大幅提高计算效率。图1为DWCC的构建流程图。

如图1所示,该方法主要由原始数据空间预筛选、快速排斥与跨立测试、交叉点的计算3个基本模块组成,各模块的实现步骤详细阐述如下。

1.1 原始数据空间预筛选

卫星测高数据通常以离散的沿轨点列形式分布。在进行交叉点解算之前,首先需要对输入的原始上升轨道和下降轨道数据进行空间预筛选,以剔除无效数据并缩小计算范围。若直接使用研究区域的严格几何边界对轨道数据进行裁剪,由于卫星离散采样的特性,可能会出现某些轨道在空间上确实穿越了研究区域,但其相邻的实测点恰好落在了边界之外的现象。这种情况会导致有效轨道被错误剔除,进而造成交叉点的漏算。为了解决这一问题,DWCC算法引入了筛选缓冲带机制。具体而言,算法首先统计所处理卫星数据的相邻测点最大间距参数k。随后,将研究区域的经纬度边界分别向外扩展k°,构建出一个外扩的计算缓冲区域。在此基础上进行点分布判断,遍历所有原始轨道,仅保留至少有一个观测点落入该扩展缓冲区域内的上升与下降轨道数据,其余未落入该区域的轨道则直接予以剔除。这一预处理模块从物理空间上杜绝了边缘处漏判现象,同时最大程度地减少了冗余轨道带来的无效计算,为后续的高效定位奠定了数据基础。

1.2 快速排斥与跨立测试

快速排斥测试(rapid rejection test)是判断由轨道端点构成的对角矩形包围盒是否重合,可以划分为3种情况:为矩形重合且轨道相交、矩形重合但轨道不相交、矩形不重合且轨道不相交,如图2所示。

图2(a)中,矩形R1与R2相交,线段M1M2和N1N2相交;图2(b)中,矩形R1与R2相交,线段M1M2和N1N2不相交;图2(c)中,矩形R1与R2不相交,线段M1M2和N1N2不相交。

图2(b)展示了矩形包围盒重合但两条轨道并不相交的情况,而如图3所示的跨立测试(straddle test)可以进一步排除此情况,图3(a)与图3(b)分别展示了M1M2跨立于N1N2以及M1M2不跨立于N1N2的两种情况。如图3所示,若线段M1M2的两端点位于N1N2的两边,则M1M2跨立N1N2。跨立测试通过计算两线段的矢量叉积来判断线段的跨立情况,若两线段互相跨立,则两线段相交。

其中,判断线段M1M2跨立N1N2且N1N2跨立M1M2的条件为

$\left\{\begin{array}{l}({\mathit{M}}_{2}-{\mathit{N}}_{2})\times ({\mathit{N}}_{1}-{\mathit{N}}_{2})·({\mathit{N}}_{1}-{\mathit{N}}_{2})\times \\ \mathit{　}\mathit{　}({\mathit{M}}_{1}-{\mathit{N}}_{2})\ge 0\\ ({\mathit{N}}_{2}-{\mathit{M}}_{2})\times ({\mathit{M}}_{1}-{\mathit{M}}_{2})·({\mathit{M}}_{1}-{\mathit{M}}_{2})\times \\ \mathit{　}\mathit{　}({\mathit{N}}_{1}-{\mathit{M}}_{2})\ge 0\end{array}\right.$

式(1)中:M2-N2表示N2为起点、M2为终点的向量,其他向量同理。

通过快速排斥跨立实验对交叉点的存在进行判断,能够避免大量冗余计算。

1.3 交叉点的计算

首先需快速确定概略交叉点的经度。卫星绕地球一周的轨道地面轨迹可分为两个弧段,分别为上升轨道和下降轨道。卫星地面轨迹在赤道两侧呈现高度的对称性。利用这一几何特征,可快速锁定交叉点的概略位置,为后续精确计算提供初值约束。假设上升轨道和下降轨道的分布如图4所示,同一卫星的上升轨道与下降轨道呈现对称性。

(1)完整轨道情形:如图4(a)所示,当测点数据连续时,利用上升轨道端点Pa1、Pa2与下降轨道端点Pb1、Pb2构建线性方程,其连线交点的经度即为概略经度λ0。

${\mathit{\lambda }}_{0}=\left[\right({\mathit{\lambda }}_{\mathrm{b}1}{\mathit{\varphi }}_{\mathrm{b}2}-{\mathit{\lambda }}_{\mathrm{b}2}{\mathit{\varphi }}_{\mathrm{b}1}\left)\right({\mathit{\lambda }}_{\mathrm{a}2}-{\mathit{\lambda }}_{\mathrm{a}1})-({\mathit{\lambda }}_{\mathrm{a}1}{\mathit{\varphi }}_{\mathrm{a}2}-{\mathit{\lambda }}_{\mathrm{a}2}{\mathit{\varphi }}_{\mathrm{a}1}\left)\right({\mathit{\lambda }}_{\mathrm{b}2}-{\mathit{\lambda }}_{\mathrm{b}1}\left)\right]·\left[\right({\mathit{\lambda }}_{\mathrm{a}2}-{\mathit{\lambda }}_{\mathrm{a}1}\left)\right({\mathit{\varphi }}_{\mathrm{b}2}-{\mathit{\varphi }}_{\mathrm{b}1})-({\mathit{\lambda }}_{\mathrm{b}2}-{\mathit{\lambda }}_{\mathrm{b}1}\left)\right({\mathit{\varphi }}_{\mathrm{a}2}-{\mathit{\varphi }}_{\mathrm{a}1}{\left)\right]}^{-1}$

式(2)中:λa1、λa2、λb1、λb2分别为点Pa1、Pa2、Pb1、Pb2的经度;${\mathit{\varphi }}_{\mathrm{a}1}$、${\mathit{\varphi }}_{\mathrm{a}2}$、${\mathit{\varphi }}_{\mathrm{b}1}$、${\mathit{\varphi }}_{\mathrm{b}2}$分别为点Pa1、Pa2、Pb1、Pb2的纬度。

(2)非完整轨道情形:针对图4(b)所示的因数据质量控制剔除或卫星轨迹跨越陆地等情况导致的数据中断情况,上述端点法失效。此时,基于轨道对称性原理,选取两条轨道上纬度相近的观测点,取其经度均值作为概略经度λ0。

${\mathit{\lambda }}_{0}=\frac{1}{2}({\mathit{\lambda }}_{\mathrm{a}}+{\mathit{\lambda }}_{\mathrm{b}})$

该方法能有效规避数据缺失带来的计算失效问题,将概略经度的估算误差控制在0.5°以内。

通过轨道对称性确定概略交叉点经度λ0后,以λ0为中心,从上升轨道和下降轨道上沿着经度方向选择相邻的观测点,从而限制所选观测点的经度范围,形成动态计算窗口,对选取测点的经度范围形成约束。随后,选取概略交叉点经度左右大约10个观测点数据形成动态定长计算窗口,这样既保证了精确交叉点的经度落于动态计算窗口内,也能够尽可能地减少参与计算的数据。对动态定长窗口内的上升轨道进行二项式拟合,即

$\mathit{\varphi }=\mathit{a}{\mathit{\lambda }}^{2}+\mathit{b}\mathit{\lambda }+\mathit{c}$

式(4)中:λ、$\mathit{\varphi }$为经纬度坐标;a、b、c分别为用于拟合上升轨道迹数据的二次多项式的系数。

此过程通过概略交叉点经度动态确定计算窗口的位置,并以固定大小确保高计算效率,为后续通过动态计算窗口内的数据确定交叉点的准确坐标奠定基础。在此采用了二次多项式,因为它能够充分捕捉卫星轨道在小的受限窗口内的局部曲率,同时避免过度拟合并确保数值稳定性,这与之前研究中的做法一致[22]。

得到上升轨道的拟合结果后,将下降轨道上的观测点的经度代入与上升轨道相匹配的多项式中,以获得相应的纬度值$\mathit{\varphi }\mathit{\text{'}},$这些纬度值代表了在相同经度位置中上升轨道上的纬度值。在$\mathit{\varphi }\mathit{\text{'}}$与下降轨道上同经度点的纬度值$\mathit{\varphi }$之间做差从而形成一系列偏差。通过差值正负符号变化的两组点,假设为第i组和第i+1组点,可以计算出精确的交叉点位置为

$\left\{\begin{array}{l}{\mathit{\lambda }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{e}}={\mathit{\lambda }}_{\mathit{i}}+\frac{({\mathit{\varphi }}_{\mathit{i}}-\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}})({\mathit{\lambda }}_{\mathit{i}+1}-{\mathit{\lambda }}_{\mathit{i}})}{(\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}+1}-\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}})-({\mathit{\varphi }}_{\mathit{i}+1}-{\mathit{\varphi }}_{\mathit{i}})}\\ {\mathit{\varphi }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{e}}=\frac{{\mathit{\varphi }}_{\mathit{i}}\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}+1}-{\mathit{\varphi }}_{\mathit{i}+1}\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}}}{(\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}+1}-\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}})-({\mathit{\varphi }}_{\mathit{i}+1}-{\mathit{\varphi }}_{\mathit{i}})}\end{array}\right.$

式(5)中:${\mathit{\lambda }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{e}}\mathit{、}{\mathit{\varphi }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{e}}$分别为交叉点的经度和纬度;${\mathit{\varphi }}_{\mathit{i}}$、$\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}}$分别为第 i组点中,轨迹1和轨迹2的纬度;${\mathit{\varphi }}_{\mathit{i}+1}$、$\mathit{\varphi }\mathit{\text{'}}{\mathrm{ }}_{\mathit{i}+1}$分别为第 i+1组点中,轨迹1和轨迹2的纬度;λi、λi+1 分别为第i组和第i+1组点的经度。

图5展示了该过程的计算步骤。在前面的步骤中,基于已构建的上升轨道多项式拟合模型,将下降轨道实测点的经度代入解算得到上升轨道在同经度处的理论纬度值。该步骤实现了升、降轨数据在经度维度的空间对齐,将交叉点搜索转化为同经度下的纬度差值计算,显著降低了算法复杂度。

2 动态窗口约束交叉点算法的验证

2.1 研究数据及区域

为了验证所提出的DWCC算法在不同轨道特性与数据体量下的有效性,选取了CryoSat-2、Jason-1 GM(大地测量任务阶段)以及Sentinel-3A三颗测高卫星的沿轨数据作为实验对象。具体选用的数据时间段涵盖了不同的观测任务周期:CryoSat-2数据选取时段为2010年7月16日—2011年9月16日;Jason-1 GM数据选取时段为2012年5月7日—2013年1月1日,该阶段提供了极高密度的空间覆盖;Sentinel-3A数据选取时段为2017年8月7日—2018年3月21日。图6展示了上述3个数据集在全球范围内的卫星地面轨迹分布情况,这些数据涵盖了从稀疏到密集的多种轨道分布模式,为评估算法性能提供了全面的数据基础。

为验证计算结果的可靠性,选取研究区域范围为 100°E~150°E, 0°N~20°N。该海域位于西太平洋低纬度带,地理环境较为复杂,不仅包含开阔大洋,还密布着众多的岛屿与群岛,在此区域进行交叉点定位实验,能够全面验证DWCC算法在不同地理环境下的稳健性与适应性。

2.2 动态窗口约束交叉点算法的实现

首先,对输入的轨道数据进行空间预筛选。考虑到卫星测高数据的离散采样特性,若直接使用研究区域的几何边界进行裁剪,可能会导致部分实际上穿越该区域但测点恰好落在边界外的轨道被错误剔除(即“边缘效应”)。为此,对CryoSat-2、Jason-1 GM及Sentinel-3A数据的相邻测点间距进行了统计分析,结果表明全部相邻测点的经纬度差值绝对值均收敛于0.3°以内。基于此,将研究区域的经纬度边界分别向外扩展0.3°建立筛选缓冲带,仅保留落入该扩展区域内的上升与下降轨道数据,从而在杜绝漏判的同时最大程度减少无效计算。

其次,快速确定交叉点的概略位置并构建动态计算窗口。利用前文所述的轨道对称性原理,快速解算概略交叉点经度,如图7所示。

获取概略交叉点经度后,在此基础上,以λ0为中心构建精确动态计算窗口,如图8所示。

假设lat1和lat2为分段拟合法中某一分段的界限,且交叉点刚好落于该分段中,传统分段拟合法通常在较大的分段区间K1内进行拟合,且窗口大小随数据量增加而扩大;而DWCC算法仅需根据概略交叉点经度确定的动态计算窗口K2内进行多项式拟合与差值计算。这种定长窗口策略不仅显著降低了矩阵运算的维度,还有效避免了长弧段拟合带来的误差,实现了计算效率与稳定性的提升。

最后,在约束窗口内完成交叉点的精确解算。在确定的K2窗口内,对上升轨道数据进行二次多项式拟合,并将下降轨道各测点的经度代入拟合方程,反演得到上升轨道在同经度处的理论纬度值。随后,计算同经度下上升与下降轨道的纬度差值序列,通过检测差值序列的正负符号变化来锁定交叉点,并利用线性插值法最终求解出交叉点的精确经纬度坐标。

2.3 交叉点位置计算结果

为了验证DWCC算法的性能,将其与传统的分段拟合法及格网法进行了对比实验,且3种方法均基于完全一致的数据预处理流程。其中,分段拟合法所使用的分段方法为拟合精度较好的12段分段(分段方式为:-90°~-65°~-60°~-50°~-40°~-20°~0°~20°~40°~50°~60°~65°~90°),而格网法构建了分辨率为1'×1'的网格进行计算。如图9所示,以CryoSat-2卫星数据为例,展示了3种方法解算的交叉点位置与原始测点及卫星轨迹的局部效果。

图9(a)为DWCC的结果;图9(b)为分段拟合法的结果;图9(c)为格网法的结果。通过观察黑色交叉点与红蓝轨迹线的相对位置,可以直观地对比不同方法的计算效果。在交叉点位置上,从图9(a)中可以清楚地看到DWCC计算出的交叉点紧密贴合在上升与下降轨道的实际交汇处,说明该方法能很好地还原交叉点的真实位置,相对来说分段拟合法则存在一定的偏移;而对于格网法,交叉点位置完全偏离了真实位置,因此该方法只能作为近似交叉点的计算方式。

为了更全面地评估算法性能,除了比较位置精度,还统计了3种方法在计算耗时和交叉点获取数量上的表现,如表1所示。

从交叉点数量来看,DWCC算法表现最佳。以CryoSat-2数据为例,DWCC算法成功提取出42 298个交叉点,而分段拟合法识别出了38 204个。这种数量上的差距主要源于分段拟合法的固有缺陷:当交叉点位于分段区间的边界处时,拟合的不连续性容易导致部分交叉点被漏算。格网法的局限性则更为突出,其获取的交叉点数量远低于实际情况,如在Sentinel-3A数据中仅识别出超过1 000个点,存在严重的遗漏现象。这也证实了格网法虽然运算速度快,但受限于网格分辨率,无法完整捕捉细微的轨道交汇情况,难以满足高精度海面建模的需求。综上所述,DWCC算法不仅定位精度高,且能确保获取交叉点时基本没有遗漏,综合优势显著。

根据表1的耗时统计结果可知,DWCC算法在计算效率上展现出显著优势。相较于传统的分段拟合法,DWCC算法通过剔除冗余计算和优化窗口搜索策略,在保证高精度定位且基本无交叉点遗漏的前提下,大幅缩减了运算时间。尽管格网法的计算耗时最少,但正如前文所述,其以牺牲巨大的定位精度和点位完整性为代价,因此不具备实际应用价值。

为了量化评估DWCC算法相对于分段拟合法的效率提升幅度,通过时间提升率来直观反映效率的提升,计算公式为

$\mathit{R}=\frac{{\mathit{T}}_{\mathrm{p}\mathrm{i}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}-\mathrm{f}\mathrm{i}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}}-{\mathit{T}}_{\mathrm{F}\mathrm{C}\mathrm{M}\mathrm{C}}}{{\mathit{T}}_{\mathrm{p}\mathrm{i}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}-\mathrm{f}\mathrm{i}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}}}\times 100\mathrm{\%}$

式(6)中:Tpiecewise-fitting为分段拟合法用时;TFCMC为DWCC用时;R为时间提升率。

基于式(6)计算得出的统计结果如表2所示。

在定位精度与交叉点完整获取方面,格网法由于受制于离散化的算法原理,其计算出的交叉点位置完全偏离了真实位置,仅能作为一种粗略的近似估计,且伴随严重的交叉点漏算现象,无法满足精密应用需求;分段拟合法虽然改善了定位精度,但相对而言仍存在一定的偏移,且因分段问题在边界处存在点位遗漏。相比之下,DWCC算法展现出最优的综合性能,其计算结果紧密贴合在上升与下降轨道的实际交汇处,能很好地还原交叉点的真实位置,更实现了对交叉点基本完整的获取,有效解决了传统方法中精度与完备性难以两全的问题。

在计算效率方面,表2的统计数据表明,DWCC算法相对于传统的分段拟合法具有显著的性能优势。通过引入动态窗口约束机制,DWCC算法在降低拟合维度的同时剔除了大量冗余运算,从而实现了效率的跨越式提升。实验数据显示,在3颗卫星的任务处理中,DWCC算法的时间提升率均超过了50%,分别为54.20%、51.03%和55.97%。综上所述,DWCC算法在确保高精度定位和无遗漏搜寻的基础上,大幅降低了计算成本,为处理海量卫星测高数据提供了一种高效、可靠的技术手段。

3 基于动态窗口约束交叉点算法的平均海面高度计算3.1 海面高计算中参数选取

在确定交叉点的精确位置后,需利用周边测点的海面高度数据推算交叉点处的瞬时海面高,可使用反距离加权法对交叉点位置附近测点数据插值计算得到交叉点处海面高度(sea surface height,SSH)数据为

$\mathit{S}\mathit{S}\mathit{H}=\frac{\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}\frac{1}{{\mathit{s}}_{\mathit{i}}}{\mathit{h}}_{\mathit{i}}}{\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}\frac{1}{{\mathit{s}}_{\mathit{i}}}}$

式(7)中:hi为第i点的海面高度;si为第i点与待算点之间的球面距离。

传统的反距离加权法对所有测点一视同仁,但在实际海洋环境中,距离交叉点越近的测点,其物理相关性越强。为了避免较远距离点对交叉点处海面高度数值产生影响,采用指数衰减型反距离加权模型,通过引入衰减因子来调节权重的空间分布,可以使得远距离测点的影响迅速减弱。

$\left\{\begin{array}{l}\mathrm{S}\mathrm{S}\mathrm{H}=\frac{{\mathit{h}}_{1}{\mathit{p}}_{1}+{\mathit{h}}_{2}{\mathit{p}}_{2}+\dots +{\mathit{h}}_{\mathit{n}}{\mathit{p}}_{\mathit{n}}}{{\mathit{p}}_{1}+{\mathit{p}}_{2}+\dots +{\mathit{p}}_{\mathit{n}}}\\ {\mathit{p}}_{\mathit{i}}={\mathrm{e}}^{-\mathit{\alpha }{\mathit{s}}_{\mathit{i}}}\end{array}\right.$

式(8)中:pi为第i点的权重系数;α为衰减因子。

图10描绘了单个点权重受不同衰减因子的影响。α取值越大,远距离测点对交叉点处的影响衰减越快。

衰减因子α的选取直接决定了插值结果的平滑度与准确性。为了确定最优的α,选取地形变化较为显著的局部海域(113.5°E~115.5°E, 15°N~17°N)进行测试。SDUST2020MSS在该区域中的MSS等值线如图11所示。

以点(114°E,16°N)为例,该处的SDUST2020MSS模型参考值为14.640 7 m。在α取不同整数值(1、3、5、7、9)的情况下进行试算,结果显示:当α较小(α<3)时,由于过度引入了远距离测点,计算结果显著偏离参考值(如α=1时的15.345 3 m),表现为平滑过渡;当逐渐增大时,计算结果迅速收敛至参考真值;在α=5时,计算结果为14.649 3 m,与参考值的偏差最小,精度达到最优;当α过大(α>7)时,权重过度集中于最近的单个测点,失去了加权平均的抗差能力,导致结果反而出现微小偏差(如α=9时的14.615 9 m)。基于上述分析,为了在压制远场噪声与保留局部特征之间取得平衡,最终选取α=5作为后续计算的参数设置。

3.2 结果评定指标

为了全面评估DWCC算法的计算精度与稳定性,选取标准差(standard deviation, STD)和均方根误差(root mean square error, RMSE)作为核心评价指标。

标准差用于衡量观测值或计算差值的离散程度,反映了数据的内部稳定性。在交叉点不符值分析中,STD越小,表明算法计算得到的交叉点海面高度差异波动越小。计算公式为

$\mathit{S}\mathit{T}\mathit{D}=\sqrt{\frac{1}{\mathit{n}-1}\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}({\mathit{\Delta }}_{\mathit{i}}{-\stackrel{-}{\mathit{\Delta }})}^{2}}$

式(9)中:n为参与统计的样本数量;Δi为第i个交叉点处的上升轨道与下降轨道MSS差异值;$\stackrel{-}{\mathit{\Delta }}$为所有差异值的平均值。

均方根误差用于衡量计算结果与参考模型(SDUST2020MSS、CLS15MSS)之间的绝对偏差,RMSE值越接近0,说明计算结果与权威模型越吻合,计算公式为

$\mathit{R}\mathit{M}\mathit{S}\mathit{E}=\sqrt{\frac{1}{\mathit{n}}\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}(\mathrm{M}\mathrm{S}{\mathrm{S}}_{\mathrm{c}\mathrm{a}\mathrm{l}}-\mathit{M}\mathit{S}{\mathrm{S}}_{\mathrm{r}\mathrm{e}\mathrm{f}}{)}^{2}}$

式(10)中:n为参与统计的样本数量;MSScal为基于DWCC算法计算得到的平均海面高度;MSSref为对应位置参考模型的平均海面高度值。

3.3 平均海面高计算

基于指数衰减型反距离加权方法,分别计算了由DWCC算法定位的交叉点处上升与下降轨道的平均海面高度MSS。理论上,同一地理位置的上升轨道、下降轨道之间MSS数值应当保持严格一致,但由于海洋环境时刻处于动态变化之中,且卫星观测系统不可避免地存在微量误差,二者之间往往存在微小的差异。为此,统计了CryoSat-2、Jason-1 GM及Sentinel-3A三颗卫星在研究区域100°E~150°E, 0°N~20°N内的交叉点MSS差异信息,具体的统计结果如表3所示。

表3展示了CryoSat-2、Jason-1 GM及Sentinel-3A卫星在研究区域(100°E~150°E, 0°N~20°N)内的升降轨交叉点MSS差异统计信息。从系统性偏差Mean来看,3个任务的升降轨平均差异均表现出极高的一致性。其中,Jason-1 GM的平均差异最小,仅为-0.000 6 m;Sentinel-3A次之,为-0.002 3 m;即便是覆盖密度较高的CryoSat-2数据,其平均偏差也控制在-0.056 7 m的极低水平。这一结果表明,在该研究区域内,上升轨道与下降轨道数据之间不存在显著的系统性偏移,这种内符合精度的稳定性意味着原始沿轨数据在几何基准上高度统一。

在数据离散度方面,Jason-1 GM 表现出最佳的稳定性,其标准差与均方根误差均为0.337 m,反映了其大地测量任务(GM)的高精度特性。相比之下,CryoSat-2和Sentinel-3A的STD和RMSE数值相对较大。这种较大的离散度并非完全源于轨道误差,更多地反映了该海域复杂的中海洋环境以及近岸或岛礁区域不可避免的波形重跟踪噪声。

基于上述交叉点分析,进行了MSS模型的构建。通过升降轨数据的平均偏差对系统性误差进行改正,并利用指数衰减型反距离加权法对研究区域中MSS进行计算,结果如图12所示。

图12直观展示了利用DWCC算法计算得到的平均海面高度(MSS)。为了进一步定量评估计算结果的精度,将DWCC算法与传统的分段拟合法计算得到的MSS结果,分别与国际通用的CLS15 MSS模型[23]、SDUST2020 MSS模型[24]参考基准进行了比对。需要特别说明的是,由第3节的交叉点位置验证结果可知,格网法受限于其离散化的算法原理,计算得出的交叉点位置存在严重偏差且伴随大量点位漏算,仅能作为一种初步的粗略计算方法,无法满足当今高精度海面高度模型构建的需求,因此本节未将其加入最终的平均海面高度计算与结果对比中。各方法在不同卫星数据下的差异统计结果如表4与表5所示。

从表4与表5的统计数据可以看出,在与多组参考模型的对比中,DWCC算法的整体计算精度均优于分段拟合法。在CryoSat-2、Jason-1 GM以及Sentinel-3A共3组数据测试中,DWCC计算结果的均方根误差(RMSE)与标准差(STD)均普遍小于分段拟合法。这表明分段拟合法在交叉点定位阶段存在几何偏移与点位遗漏,导致其最终的海面高度计算精度略差;而DWCC算法则能够更加精准地还原真实海面起伏特征。

此外,对比不同卫星的计算结果发现,Jason-1 GM数据的计算结果精度最高,其与CLS15和SDUST2020模型的均方根误差(RMSE)分别仅为0.067 2 m和0.082 9 m,且平均偏差极小,说明其结果能准确反映海面起伏特征,与参考模型具有高度的一致性。相比之下,CryoSat-2和Sentinel-3A的计算结果与模型之间存在约-0.7 m的系统性偏差。经分析,该偏差主要源于数据产品所采用的参考椭球基准不一致,SDUST2020和CLS15模型均基于Topex/Poseidon(T/P)椭球构建,而CryoSat-2和Sentinel-3A数据则基于WGS84椭球,两者之间存在约0.7 m的理论几何差异[25]。在去除该系统偏差后,其标准差(STD)仍保持在合理范围内,说明DWCC算法能够有效恢复海面高度的变化特征。

通过对比两种方法的精度,证实了DWCC算法优于分段拟合法。为了直观地剖析DWCC算法的误差分布特征,将DWCC算法的计算偏差进行了可视化展示,如图13所示。图13展示了DWCC计算结果与参考模型差异的概率分布图。图13中彩色柱状图展示了差异值的概率密度直方图,反映了误差的实际分布频率;而平滑实线则为拟合的概率密度曲线,用于直观揭示误差分布的整体趋势与形态特征。

从分布形态来看,3组卫星数据的偏差分布均呈现出显著的正态分布特征,曲线光滑且对称。这表明DWCC算法的计算残差主要由随机噪声构成,未出现明显的异常值干扰或非线性畸变,验证了算法在处理复杂海面信号时的稳健性;从峰值位置来看,图13(b)中Jason-1 GM的分布峰值几乎重合于零值线,这直观地验证了其与参考模型在参考椭球基准上的高度一致性。相比之下,图13(a)中CryoSat-2与图13(c)中Sentinel-3A的分布曲线整体相对于零值线有一定偏移,峰值中心稳定落在-0.7~-0.75 m,这主要源于数据产品所采用的参考椭球基准不一致。

综上所述,通过表4与图13的验证,不仅量化了DWCC算法的高精度特性,更清晰地分离了系统性基准偏差与随机观测误差。证明了在统一参考椭球基准后,该算法具备构建高精度、高分辨率平均海面高模型的能力。

4 结论

高效而精准的交叉点定位方法是快速获取卫星轨迹交叉点的关键手段,对于平均海面高度模型构建也具有重要影响。当前常用的交叉点算法在分析中显示,它们均在计算效率或定位精度方面存在一定不足。针对卫星测高数据中交叉点定位精度与计算效率难以兼顾的挑战,提出了一种新型动态窗口约束交叉点算法(DWCC),旨在高效处理卫星数据并获取高精度交叉点结果,能够有效保障交叉点位置的高精度。

(1)动态窗口约束交叉点算法构建。 基于多项式拟合与交叉点计算原理,引入快速排斥与跨立测试机制进行交叉点的存在性检测,并通过约束精确交叉点动态计算窗口策略降低冗余计算,提出了一种新型动态窗口约束交叉点算法。该方法通过动态窗口自适应调整计算区域,并使用计算窗口定长策略限制计算量,在无交叉点遗漏的同时,高效获取高精度交叉点位置,为平均海面高度模型构建奠定基础。未来研究可进一步优化多项式拟合阶数与窗口大小的自适应调整策略,以增强算法在极地轨道或大幅度机动轨道任务中的普适性。

(2)动态窗口约束交叉点算法验证。在不同研究区域中选择不同方法进行交叉点位置验证,结果表明:DWCC算法在实现交叉点基本无遗漏获取的前提下,能够精准定位交叉点,并显著提升计算效率。其中,DWCC算法较分段拟合法计算效率提升超过50%;而在定位精度方面,相比于分段拟合法存在的几何偏移与边界漏算现象,DWCC算法实现了交叉点位置与真实轨迹交汇处的紧密贴合,从而验证了该算法在兼顾高精度、高完备性与高效率方面的显著优势。这些优势为未来处理下一代宽刈幅测高卫星产生的海量数据提供了技术储备,有望实现全球尺度下更短时间分辨率的交叉点动态监控。

(3)平均海面高度模型构建与精度验证。结果显示,在与参考模型的对比中,DWCC算法的整体计算精度显著优于分段拟合法,能够更加精准地还原真实海面起伏特征。DWCC算法计算得到的平均海面高度与SDUST2020及CLS15 MSS模型呈现出高度的一致性。其中,Jason-1 GM数据的均方根误差(RMSE)优于0.08 m,验证了算法的高精度特性;同时,研究明确了部分卫星数据存在的系统性偏差源于参考椭球定义的几何基准差异,而非计算误差。该方法为高精度海面动态监测与建模提供了可靠的数据基础,具有广阔的应用潜力。在未来工作中,将进一步探究该算法在近岸复杂海域及多源数据融合中的适用性,并通过统一基准变换模型,进一步精化近岸及复杂海域的MSS建模,为海平面变化监测提供更精准的几何参考基准。

(4)算法的多领域应用与智能化展望。本算法获取的高精度交叉点不仅可用于海面高度建模,还将为海洋重力异常反演、海底地形建模以及海洋重力场变化率研究提供关键的数据支持。未来将探索融合深度学习技术,进一步提升多源卫星数据融合处理的智能化水平,为构建高空间分辨率、高精度的全球数字海洋模型提供高效的数据处理手段。

参考文献

原文顺序 | 出版日期 | 本文引用
[1]

Cazenave A, Remy F. Sea level and climate: measurements and causes of changes[J]. WIREs Climate Change, 2011, 2(5): 647-662.
[2]

Andersen O B, Stenseng L, Piccioni G, et al. The DTU15 MSS (mean sea surface) and DTU15LAT (lowest astronomical tide) reference surface[C]// Proceedings of ESA Living Planet Symposium. Prague: ESA, 2016: 1579.
[3]

Zhu F, Guo J, Zhang H, et al. SDUST2020MGCR: a global marine gravity change rate model determined from multi-satellite alti-meter data[J]. Earth System Science Data, 2024, 16: 2281-2296.
[4]

Sandwell D T, Smith W H F. Marine gravity anomaly from Geosat and ERS-1 satellite altimetry[J]. Journal of Geophysical Research: Solid Earth, 1997, 102(B5): 10039-10054.
[5]

Nerem R S, Beckley B D, Fasullo J T, et al. Climate-change-driven accelerated sea-level rise detected in the altimeter era[J]. Proceedings of the National Academy of Sciences of the United States of America, 2018, 115(9): 2022-2025.
[6]

李志军, 宁靓, 陈任霖. 中国国际海底区域资源勘查开发: 现状分析与战略研究[J]. 科学技术与工程, 2024, 24(27): 11558-11564.
[7]

Li Zhijun, Ning Jing, Chen Renlin. Exploration and development of international seabed area resources in China: current situation analysis and strategic research[J]. Science Technology and Engineering, 2024, 24(27): 11558-11564.
[8]

Morrow R, Le Traon P Y. Recent advances in observing mesoscale ocean dynamics with satellite altimetry[J]. Advances in Space Research, 2012, 50(8): 1062-1076.
[9]

Andersen O B, Zhang S, Sandwell D T, et al. The unique role of the Jason geodetic missions for high resolution gravity field and mean sea surface modelling[J]. Remote Sensing, 2021, 13(4): 646.
[10]

Beckley B D, Callahan P S, Hancock D W, et al. On the ‘cal-mode’ correction to TOPEX satellite altimetry and its effect on the global mean sea level time series[J]. Journal of Geophysical Research: Oceans, 2017, 122: 8371-8384.
[11]

Wingham D J, Francis C R, Baker S, et al. CryoSat: a mission to determine the fluctuations in Earth’s land and marine ice fields[J]. Advances in Space Research, 2006, 37(4): 841-871.
[12]

郑伟, 李钊伟, 吴凡. 基于新一代GNSS-R星座海面测高原理提高水下惯性/重力组合导航精度研究进展[J]. 科学技术与工程, 2019, 19(36): 21-36.
[13]

Zheng Wei, Li Zhaowei, Wu Fan. Research progress in improving the accuracy of underwater inertial/gravity integrated navigation based on the new generation of GNSS-R constellation sea surface altimetry principle[J]. Science Technology and Engineering, 2019, 19(36): 21-36.
[14]

Yang J Z, Jia Y J, Fan C Q, et al. Preliminary results of the global ocean tide derived from HY-2A radar altimeter data[J]. Acta Oceanologica Sinica, 2023, 42(2): 65-73.
[15]

Wang Y, Guo J Y, Ya S X, et al. Impact of pseudo-stochastic pulse and phase center variation on precision orbit determination of Haiyang-2A from experimental HY2 receiver GPS data[J]. Remote Sensing, 2024, 16(8): 1336.
[16]

李国元, 唐新明, 陈继溢, 等. 高分七号卫星激光测高数据处理与精度初步验证[J]. 测绘学报, 2021, 50(10): 1338-1348.
[17]

Li Guoyuan, Tang Xinming, Chen Jiyi, et al. Data processing and preliminary accuracy verification of Gaofen-7 satellite laser altimetry[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(10): 1338-1348.
[18]

Mashburn J R. Analysis of GNSS-R observations for altimetry and characterization of earth surfaces[D]. Boulder: University of Colorado, 2018.
[19]

Tai C. An efficient algorithm for computing the crossovers in satellite altimetry[J]. Journal of Atmospheric and Oceanic Technology, 1989, 6: 13-18.
[20]

周晓光, 苗洪利, 王云海, 等. 卫星地面轨迹分段拟合确定交叉点的方法研究[J]. 测绘学报, 2012, 41(6): 811-815.
[21]

Zhou Xiaoguang, Miao Hongli, Wang Yunhai, et al. Research on crossover method based on piecewise fitting of satellite ground tracks[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(6): 811-815.
[22]

Greene C A, Gwyther D E, Blankenship D D. Antarctic mapping tools for MATLAB[J]. Computers & Geosciences, 2017, 104: 151-157.
[23]

Zhu C C, Li Z, Guo J Y, et al. The calculation and analysis of along-track and cross-track crossover SSHs discrepancies from wide-swath altimetry data[J]. SSRN, 2024. DOI: 10.2139/ssrn.4882946.
[24]

Yuan J J, Guo J Y, Niu Y P, et al. Mean sea surface model over the Sea of Japan determinedfrom multi-satellite altimeter data and tide gauge records[J]. Remote Sensing, 2020, 12(24): 4168.
[25]

赵子龙. 卫星交叉点快速确定方法研究[D]. 西安: 长安大学, 2021.
[26]

Zhao Zilong.Research on fast determination method of satellite crossover[D]. Xi’an: Chang’an University, 2021.
[27]

Li X, Zhang S, Geng T, et al. An improved algorithm for extracting crossovers of satellite ground tracks[J]. Computers & Geosciences, 2022, 166: 105179.
[28]

Pujol M I, Schaeffer P, Faugère Y, et al. Gauging the improvement of recent mean sea surface models: a new approach for identifying and quantifying their errors[J]. Journal of Geophysical Research: Oceans, 2018, 123: 5889-5911.
[29]

Yuan J, Guo J, Zhu C, et al. SDUST2020 MSS: a global 1'×1' mean sea surface model determined from multi-satellite altimetry data[J]. Earth System Science Data, 2023, 15(1): 155-169.
[30]

Feng Y, Yang L, Zhan P, et al. Synthesis of the ICESat/ICESat-2 and CryoSat-2 observations to reconstruct time series of lake level[J]. International Journal of Digital Earth, 2023, 16(1): 183-209.

基金资助

国家自然科学基金重点项目(42274119)

国家自然科学基金重点项目(U25A20783)

国家重点研发计划(2022YFF1400500)

国家高层次领军人才计划(X2023081)

“十四五”综合实验区重大专项(WZC25110122)

“十四五”综合实验区重大专项(CECOM-2025-BDJC02-SY-0737/06)

国家领军人才山东省配套项目(IDFA10202024)

哈尔滨工业大学(威海)青年教师发展基金(理学工学类)(IDGA10002224)

河南省测绘科学与技术学科“双一流”建设科研项目(BZCG202303)

辽宁省科技厅博士科研启动基金计划(2025080690-JH3/101)

辽宁省青年科技人才托举工程(2025-Z0233)

国家自然科学基金(42508023)

测绘遥感信息工程国家重点实验室(LIESMARS)开放基金(25R04)

河南省科学院基本科研业务费项目(20260625001)

PDF (14531KB)

7

访问

0

被引

详细
导航
相关文章

/

〈 〉

京ICP备05035734号-4
版权所有 © 《科学技术与工程》编辑部
地址： 北京市海淀区学院南路86号 《科学技术与工程》编辑部　邮编：100081
电话：(010)62118920 　E-mail：stae@vip.163.com
技术支持：北京玛格泰克科技发展有限公司 　
网站管理员登录