An Approach for Real-Time GPS/GLONASS Satellite Clock Estimation with GLONASS Code Inter-Frequency Biases Compensation
-
摘要: 在进行GPS/GLONASS联合卫星钟差估计时，GLONASS码频间偏差（inter-frequency bias，IFB）因卫星频率间的差异而无法被测站接收机钟差参数吸收，其一部分将进入GLONASS卫星钟差估值中。通过引入多个"时频偏差"参数（inter-system and inter-frequency bias，ISFB）及附加基准约束对测站GLONASS码IFB进行函数模型补偿，实现其与待估卫星钟差参数的有效分离，并对所估计实时卫星钟差和实时精度单点定位（real-time precise point positioning，RT-PPP）进行精度评估。结果表明，在卫星钟差估计观测方程中忽略码IFB，会明显降低GLONASS卫星钟差估值精度；新方法能有效避免码IFB对卫星钟差估值的影响，所获得GPS、GLONASS卫星钟差与ESA（European Space Agency）事后精密钟差产品偏差平均均方根值分别小于0.2 ns、0.3 ns。利用实时估计卫星钟差进行静态RT-PPP，当观测时段长为2 h时，GPS单系统、GPS/GLONASS组合系统的3D定位精度优于10 cm，GLONASS单系统3D定位精度约为15 cm；三种模式24 h单天解的3D定位精度均优于5 cm。Abstract: In GPS/GLONASS satellite clock estimation, GLONASS code inter-frequency biases (IFB) can not be assimilated into the receiver clock parameters due to the discrepancy of GLONASS frequencies, which is partly involved in the GLONASS satellite clock estimates. In order to separate the GLONASS IFBs and the satellite clock estimates effectually, an approach for real-time GPS/GLONASS satellite clock estimation with multiple inter-system and inter-frequency bias (ISFB) parameters and additional restriction is proposed. The estimated satellite clocks and real-time precise point positioning results are then evaluated. The results show that ignoring GLONASS code IFBs could dramatically degrade the GLONASS satellite clock accuracy. The new approach can efficiently avoid the adverse impact of GLONASS code IFBs on satellite clock estimates. Compared with the ESA final satellite clock products, the average RMS of clock biases are less than 0.2 ns for GPS and 0.3 ns for GLONASS respectively. RT-PPP experiments were conducted with different session lengths. The 3D positioning accuracy of RT-PPP using the estimated satellite clocks was better than 10 cm for GPS-only and combined GPS/GLONASS modes, 15 cm for GLONASS-only mode with 2-h observations, and better than 5 cm for all the three modes with 24-h observations.
-
高分影像多分辨率的本质是卫星遥感传感器放大了对地观测的尺度，提高了影像的分辨能力。高分辨率意味着影像可以获得更多地物细节，因此，地物之间的尺度关系更加丰富，尺度差异性更加显著。而地物分布的复杂性主要在于它是多尺度的统一，这就导致了地物信息提取中多尺度分析的问题。如果没有对高分影像的多尺度分析和多尺度特征的结构化建模，那么影像特征计算与信息提取中的尺度确定只具有单一的选择性，其自动化程度必然很低。
已有的高分影像多尺度分析方法大致包括两类。一是基于特征向量的分类器分类方法。特征计算方法主要有灰度共生矩阵[1]、像元形状指数(pixel shape index, PSI)[2]以及数学形态学方法[3]等；分类器主要有人工免疫分类器(artificial immune classifier, AIC)[4]和支持向量机[5-6]等。这些方法已经取得了一定的效果。但特征计算只能表达影像地物的多尺度结构，本身并没有尺度识别和选择的能力；同时，分类器本身也不具有尺度识别和选择的能力。因此，通过向量累加构成多尺度高维特征空间(向量)，再用分类器进行分类的特征融合方法，缺乏精准的尺度之间的关联或传递机制，是一种对特征的“平均融合”处理，所以那些个性突出的尺度特征依然会偏离融合中心，造成尺度偏移和椒盐噪声等现象。从这个角度看，通过特征计算与分类器进行多尺度融合具有一定的粗放性，而且其有效程度依赖于高维特征空间(向量)的可分程度[7]。二是基于多层结构的多尺度对象融合方法。多层结构主要指金字塔结构[8-9]。利用面向对象方法[10-11]进行多尺度处理存在两个困境：①处理精度受制于对象的产生方法(即影像分割方法)；②多尺度对象之间如何融合。文献[12]讨论用均值漂移(mean shift, MS)方法进行像元的高维特征聚类获得对象，再以对象的平均特征值作为融合后的对象特征与单个像元进行混合分类，这个方法将面向像元和面向对象(object-based image analysis, OBIA)的混合分析方法成功应用在城市高光谱多尺度分类中。文献[13]认为，基于MS方法的多尺度分割精度明显高于分形网络进化方法(fractal net evolution approach, FNEA)方法[11]，并指出在多层OBIA方法中，MS是一种非常有潜力的技术。对于多层结构的多尺度对象特征而言，MS方法的性能要优于FNEA[7]。但是，这些方法与基于特征向量方法本质上基本相似。为此，文献[8]用小波金字塔建立多尺度结构，从粗尺度到细尺度，逐级使用MS进行特征空间聚类，在多层结构内建立了比较精准的多尺度对象的融合机制，其多尺度特征传递策略具有重要意义。文献[14]分析了尺度传递中特征空间稀疏性和复杂性的变化规律以及带宽尺度选择的敏感性问题；文献[15]分析了小波特征的可分性与多尺度结构的表达合理性问题；在此基础上，文献[14-15]提出了新的约束MS方法，讨论了带宽尺度敏感性问题。文献[10]通过细节亚对象(detailed sub-objects, DSB)与粗略超对象(coarser super-objects, CSB)的聚集效应，文献[16]通过尺度选择的类别比例(class proportion, CP)效应，分别研究了尺度的不确定性问题。
从以上分析可以看出，目前多尺度分析方法在特征计算与融合策略上都取得了一些进展，为后续研究奠定了坚实基础。但是这些方法还存在以下共同缺陷:①在利用多尺度特征互补性的同时，无法兼顾尺度的个性特征；②没有充分利用多尺度特征的层次性，从而无法在结构化特征内部建立多尺度特征的关联与传递途径；③没有利用多尺度特征的冗余性，从而削弱了多尺度信息的语义一致性表达能力。
基于以上分析，本文提出基于特征结构化多尺度分析建模方法研究，主要包括两个方面。一是面向高分影像地物结构的各向异性特征和多尺度特征，针对性地研究特征变换方法和特征结构化表达形式。二是用结构化特征表达形式引导多尺度分析建模，包括3个结构相关但形式不同的模型：①下采样金字塔模型结构；②非下采样直塔模型结构；③金字塔与直塔的混合模型结构。并针对非下采样直塔模型结构进行了具体分析与实验，验证了建模方法的可行性。
1 多尺度分析建模方法概述
多尺度分析建模方法包括3个主要部分，其技术总路线如图 1所示。图 1中第一部分为影像几何、结构和光谱特征变换与提取。采用变换域或光谱分析等特征描述算法，对高分影像进行特征提取，获得多尺度、多方向、多层级等多通道的特征表达形式，为构成多层结构的多尺度特征空间奠定基础。第二部分为构造多尺度特征空间。主要运用小波与多尺度几何分析(multi-scale geometric analysis, MGA)等变换方法，构建多尺度、多层级特征的内在联系，将影像特征在尺度维最大化衍生，产生内在联系紧密的结构化影像特征序列，为建立多尺度分析模型奠定基础，其中红色“×”号表示不建议使用这个技术路线构建直塔结构。第三部分为建立多尺度分析模型。在尺度维进行空间多尺度变换，包括下采样、上采样以及非下采样等技术，将特征尺度和空间尺度的内在关联性并行纳入塔式结构之中，建立多尺度分析模型。
MGA方法不仅提供了优良的多尺度变换域特征，而且多尺度变换域的统计模型能够有效地描述变换域系数在尺度间、尺度内和方向间的统计相关性。多尺度分析模型的研究目的在于同时利用高分影像在不同尺度下的变换特征，为不同分布特征的目标自适应找到合适的尺度，主要通过3种多尺度分析模型来实现，包括金字塔、直塔以及混合塔模型结构。
1.1 金字塔模型结构
既有研究成果运用小波金字塔实现影像多尺度聚类分析[14-15]，综合考虑小波基的正交性、紧致性和消失矩的阶数，其算法选用Daubechies4小波基，算法原理与流程说明如下。
1) 小波变换包括两步。①对原始影像逐层进行小波分解，从第1至第n层，由于下采样使影像尺度逐层以1/4递减变粗，第1层尺寸是原始影像的1/4，这里n=4；②0层是原始影像进行小波变换后非下采样，与原始影像同尺寸。
2) 由0~n层构成金字塔结构。
3) 在金字塔的每一层，分别对高、低频提取标准偏差和均值特征，构成特征空间。
4) 多尺度间聚类信息传递与均值漂移：①对金字塔第ni层(上层)特征空间，进行约束均值漂移聚类(当ni=n时即为顶层，此时无约束条件)，对第ni层聚类影像作边缘标记，均质区聚类值即为第ni-1层(下层)均值漂移聚类的约束路径值；②将第ni层聚类影像标记后上采样(扩展4倍)，其像元与第ni-1层相对应，在第ni-1层中对与第ni层均质像元相应的像元特征取平均值，为第ni-1层的新特征值，边缘像元特征值保持不变，它们共同构成第ni-1层新的特征空间，实现第ni层聚类尺度信息对第ni-1层的传递，对ni-1层新的特征空间进行约束均值漂移聚类，再对第ni-1层聚类影像作边缘标记，进入下一轮迭代；③依次向下，重复步骤②的操作过程，直至0层，获得最终分割结果。
上述算法针对纹理影像聚类问题，对纹理内部特征的多尺度逐级综合能力很强，因此具有很好的区域分割能力。但其局限性也很明显，它假设多尺度对象之间的不确定性只存在于边缘附近，这种假设显然缺乏一般性。高分影像地物内部的空间结构特征、地物间的尺寸形状及其近邻的细节关系特征等都非常丰富和复杂，多尺度对象之间的不确定性不仅存在于边缘，更多的是相互间存在重叠、交错与切分等现象，这些都是多尺度分析中必须解决的互补与融合问题。
一维小波分析所具有的优异特性并不能简单推广到二维或更高维，虽然高分影像的细节特征更为清晰，为高维稀疏特征分析提供了数据基础，但是小波分析的局限性却更加明显。与小波变换相比较，MGA则是对线奇异或面奇异的高维函数最优或最稀疏的表示方法。
基于这个事实，为了充分利用高分影像的细节特征，融合多尺度空间结构信息，进一步改进和发展既有的研究成果，本文提出下采样多尺度几何分析(sampled multiscale geometric analysis, SMGA)的金字塔模型结构，如图 2所示。其基本思路是用SMGA变换替代小波变换，并将多尺度对象之间的边缘标记方法替换为对象整合和特征融合等方法。首先针对不同尺度特点，用下采样变换压缩空间几何结构，突出像元邻域关系特征，获得空间多尺度、特征多通道的金字塔模型；再使用特征互补与融合等方法，对每个尺度层进行独立的聚类操作，获得聚类对象，并使这些对象逐级接受上层粗尺度对象的整合，逐级获得尺度传递和特征融合，为目标选择最佳类别。这样一方面可以避免粗尺度对象的不纯净性对细尺度的影响，另一方面可以使高分影像的丰富细节特征在尺度间充分继承。算法的一般性原理与流程说明如下。
![]()
图 2 多尺度几何分析特征变换的金字塔模型结构递归算法流程图Figure 2. Flowchart of Recursive Algorithm of Feature Transformation for Pyramid Structure Model Using Multi-scale Geometric Analysis1) 采用下采样方式进行多尺度特征变换，构成多尺度金字塔模型。
2) 进行MS聚类(或者其他聚类模型，下同)，生成多尺度空间压缩对象层，再对多尺度空间压缩对象层进行上采样，获得空间扩展后的多尺度对象层的CSB或DSB。
3) 从顶层开始，逐层动态进行CSB对DSB间的对象整合，获得新的具有特征融合的对象层。
4) 对新对象层的融合特征空间进行MS聚类，动态获得新的CSB层，如图 2中箭头流向与说明所示。
5) 新的CSB层再对其下层的DSB层进行整合，重复步骤1)至4)，依次向下直到0层，结束。
1.2 直塔模型结构
高分影像的最大优势即精细的对地观测能力，所以在多尺度分析中，应该具有保持影像细节的能力。然而，下采样过程会损失一些细节特征，SMGA金字塔模型存在微小结构丢失的风险，限制了影像高分辨率的优势。本节提出非下采样多尺度几何分析(non-sampled multiscale geometric analysis, NSMGA)模型，即直塔模型结构，可以保持细节特征，避免微小结构丢失。直塔模型结构的自下而上多通道特征空间逐级稀疏，背景数据的冗余性逐级加大，虽然增加了统计模型的运用难度，但是冗余性有利于多通道间的特征互补，只要针对稀疏性和冗余性特点适当改进聚类方式或过程，就可以在细节保持与宏观结构间得到最优效果。本文在§2详细讨论直塔模型的一种具体实现过程，它可以实现高分辨影像线结构特征的稀疏表示。
1.3 混合塔模型结构
金字塔结构的重要意义在于，通过下采样压缩不同尺度层的空间冗余信息，增强近邻像元的空间相关性，将尺度信息隐含于压缩的近邻像元关系之中，可以更好地适应于各种统计模型。直塔模型结构采用非下采样的尺度变换，各尺度层间的像元尺度关系变化不依赖空间压缩操作，产生的冗余性有利于特征互补，但是削弱了尺度的逐级空间关联，使统计模型处理产生困难。
基于这个事实，本文进一步提出金字塔和直塔结构相结合的混合塔模型，如图 3所示。混合塔模型既可以保持大目标细节特征以及避免小目标丢失，也可以同时在时(空)频域中进行尺度选择操作，将多尺度的影像空间和多通道的特征空间有机结合。图 3算法的一般性原理与流程说明如下。
![]()
图 3 多尺度几何分析特征变换的混合塔模型结构递归算法流程图Figure 3. Flowchart of Recursive Algorithm of Feature Transformation for Hybrid Tower Structure Model Using Multi-scale Geometric Analysis1) 对金字塔结构的各层聚类结果分别进行上采样(空间扩展)，得到空间尺度相同的对象层。
2) 将上采样金字塔各尺度层的对象分别与直塔各对应尺度层的对象进行整合，生成多尺度对象层。
该算法的后续处理步骤同金字塔模型。
上述3种模型的特征通道都是多维变换，除MGA变换方法之外，也可以采用其他方法，其模型框架具有一般性，具体模块的作用与内涵都应该与应用问题相结合，进行具体改变，而不必局限于本文流程图。
2 直塔模型结构实验研究
金字塔模型研究已经取得了一些成果[8, 15, 17]，本节对直塔模型进行具体实验研究，以验证模型方法的正确性和有效性。
2.1 影像线结构特征变换与提取
MGA包括多种具体变换域方法，例如，Ridgelet变换、单尺度Ridgelet变换、Curvelet变换、Bandelet变换以及Contourlet变换等。在建模中，需要结合高分影像的特点与提取信息的任务目标，有针对性地进行具体选择和实践。Ridgelet变换技术的应用目前主要集中于影像压缩与去噪方面[18-20]，针对边缘检测或稀疏表示问题的研究还很不充分，本节基于改进Ridgelet变换的直塔模型，实现影像线结构特征的稀疏表示。
Radon变换是Ridgelet变换的基础，为此，简单回顾Radon变换原理。若快速衰减的光滑函数为f(x)=f(x1, x2)，x∈R，设Rf是函数f的Radon变换，用狄拉克δ函数表示，沿直线L(θ, t)积分，则Radon变换如下[21]：
$$ ({R_f})\left( {t, \theta } \right) = \int {f\left( x \right)\delta ({x_{\rm{1}}}{\rm{cos}}\theta + {x_2}{\rm{sin}}\theta - t){\rm{d}}\mathit{x}} $$ (1) 式中，θ∈[0, 2π)，t∈R。
在与影像对偶的Radon变换域，设脊方向的规范正交基为Uμ(t)，t∈L2(R)，角方向的基为Vν(θ)，θ∈L2[0, 2π)，则数字Ridgelet变换的基函数Wλ可以写成如下形式[21]：
$$ \begin{array}{l} {W_\lambda } = {P_A}\left[ {{U_\mu }\left( t \right) \otimes V{{\left( \theta \right)}_{\nu |\mu }}} \right], \\ \;\;\;\;\;\;\forall \lambda = \left( {\mu , \nu } \right), \mu , \nu \in \mathit{\boldsymbol{Z}} \end{array} $$ (2) 式中，角方向的基(Vν|μ:ν)依赖于μ；A是服从反对称性L2(dtdθ)空间上的闭合子空间；PA是从L2(dtdθ)到A的正交投影；⊗表示张量级。文献[21]给出数字Ridgelet变换包的基函数定义如下。
定义1 基函数Uμ和Vν(或Vν|μ，如果μ和ν存在函数关系)分别是从小波包中选择的基函数，其中Vν也可以从余弦包中选择，因此，称式(2)是脊波包的基函数。
由定义1可以推论出Radon变换域中对偶正交基框架的有关性质。
推论1 根据定义1，基函数Uμ和Vν可以分别从小波包中选择父小波和母小波基函数，式(2)依然是脊波包的基函数。
推论2 根据定义1，如果μ, ν=1，式(2)依然成立。
文献[21]对式(2)的Wλ构造原理给出了严格条件以及弱化条件，虽然弱化条件可能导致基函数冗余性增大以及空间衰减性较差，但其基函数结构更具普适性。因此，本节基于Wλ构造原理的弱化条件，提出基于Ridgelet变换特征的多尺度综合框架，即基于改进Ridgelet变换特征的多尺度直塔模型。
根据Radon变换原理，从式(1)可以看出，如果频率域径向方向θ=θ0的正交方向θ0+π/2存在一条直线x1=x2tan (θ0+π/2)+t′0，则在θ0径向方向上的脊函数值x1cosθ0+x2sinθ0=t0是一个常数，这个常数的积分形式为脊函数Rf(t0, θ0)= ∫f(x)dx，与这些脊函数相正交的小波束共同构成Ridgelet变换基。
显然，在Radon变换域中，直线目标在点(t0(θ0), θ0)的邻域，同时表现出脊方向与角方向奇异性(迅速衰减)，即点奇异。然而，一般非直线目标在Radon变换域中，由于Rf(t, θ)随着角度变化具有连续不间断的平滑性[22]，点(t0(θ), θ)的邻域表现为一条随径向角度θ而变化的光滑曲线t0(θ)，即曲线奇异。
由于式(2)构成了具有反对称性的小波张量基，设R*表示Radon等距逆变换，则Ridgelet变换单元ρλ可以写成：
$$ {\rho _\lambda } = {{\bar R}^*}[{W_\lambda }] $$ (3) 可见，Radon变换域对线结构的特征具有可测性，根据推论1、推论2以及式(2)的结论，对Radon变换域的对偶正交基可以进行如下改进：
1) 根据推论1，为检测脊方向线结构奇异性间断点，基函数Uμ采用db2母小波，可以获得脊方向高频奇异点特征；为逼近角方向线结构连续性变化，基函数Vν采用Meyer父小波[23]，可以获得角方向低频近似连续特征。
2) 根据推论2，式(2)依然成立，其瓦片实际上已经退化为一个点，即在式(2)中Radon投影变换不是以瓦片区域为单元进行，而是以径向脊的函数U1(t)和角的函数V1(θ)的局部点(U1(t), V1(θ))特征进行(t∈L2(R), θ∈L2[0, 2π))。
3) 根据推论1、推论2，式(2)可以写为
$$ \begin{array}{l} {W_{\lambda = 1}} = {P_A}\left[ {{U_1}\left( t \right) \otimes {V_1}\left( \theta \right)} \right], \\ \;\;\;\;\;\forall \lambda = \left( {\mu = 1, \nu = 1} \right) \end{array} $$ (4) 进一步写为：
$$ \begin{array}{l} {W_{\lambda = 1}}\left( {t, \theta } \right) = (\mathit{\boldsymbol{\psi }}\left( t \right) \otimes \mathit{\boldsymbol{w}}\left( \theta \right) + \\ \;\;\;\;\;\mathit{\boldsymbol{\psi }}\left( { - t} \right) \otimes \mathit{\boldsymbol{w}}(\theta + \pi ))/2 \end{array} $$ (5) 式中，ψ(·)是对脊方向逐点U1(t)进行db2母小波卷积的向量；w(·)是对角方向逐点V1(θ)进行Meyer父小波卷积的向量。
从以上分析可以看出，线结构特征存在于Wλ=1变换域中，设PRWλ是对变换域Wλ=1的模式识别算子，BλLine为线结构子带空间，可得出如下关系：
$$ B_\lambda ^{{\rm{Line}}} = {\rm{P}}{{\rm{R}}_{{W_\lambda }}}({W_{\lambda = 1}}) $$ (6) 式中，λ的含义等同式(2)，区别仅在于BλLine中的λ是为进行线结构模式识别所做的瓦片划分。
因为式(6)中BλLine是线结构特征模式在Wλ=1变换域中的子带空间，可以使用式(3)的Radon等距逆变换R*进行空间域线结构特征模式的投影变换，得到线结构变换的空间域系数ρλLine：
$$ \rho _\lambda ^{{\rm{Line}}} = {{\bar R}^*}\left[ {B_\lambda ^{{\rm{Line}}}} \right] $$ (7) 2.2 多尺度特征空间构成
基于式(7)，线结构系数变换的多尺度特征空间数学模型如下。
设IN×N∈L2(R)为输入影像，其中N是影像尺寸，必须满足N=2n, n=1, 2…nlit，nlit是有限正整数，则多尺度分解子带的直塔结构Pdec的数学表示为：
$$ {P_{{\rm{dec}}}} = \bigcup\limits_{H = 1}^4 {(P_ \bot ^H(\bigcup\limits_{r = 1, c = 1}^{{N_H}, {N_H}} {I_{r, c}^{{A_H}}} ))} $$ (8) 式中，P⊥H表示垂直结构单尺度子带，H是尺度层序号；Ir, cAH表示单尺度子带内的影像块，r、c分别是影像块位置坐标；NH=N/SIr, cAH，SIr, cAH是Ir, cAH的边长；AH是与尺度层有关并服从反对称性L2(dtdθ)空间上的闭合子空间。本实验中SIr, cAH=2n, n=4, 5, 6, 7，即分别是从16×16到128×128二进制影像块划分构成的多尺度分解子带集，即直塔结构。
在式(8)基础上，通过对单尺度子带及其旋转后的子带进行影像块划分格网的位置平移，构成多尺度空间位移子带集，其数学表示为：
$$ P_{{\rm{dec}}}^{{\Delta _H}} = \bigcup\limits_{H = 1}^4 {(P_ \bot ^H(\bigcup\limits_{\begin{array}{*{20}{c}} {{\Delta _H} = 0, 1, 2, 3}\\ {{\theta _H} = 0, 1, 2, 3} \end{array}} {(\bigcup\limits_{r = 1, c = 1}^{{N_H}, {N_H}} {I_{r, c}^{{A}}} )} ))} $$ (9) 式中，ΔH表示单尺度子带直接平移的位移方式；θH表示单尺度子带先旋转再平移的位移方式，本实验采用正交位移量为SIr, cAH边长的一半，对角位移动量为其$ \sqrt 2 $倍，位移方向分别是位置0(左上方，即原位置)、位置1(正右方)、位置2(正下方)、位置3(右下方)，按位移后的格网重新划分影像块；再对原始影像90度旋转，并重复4种平移的位移方式，重新划分影像块，共同构成8个一组的多尺度空间位移子带集。
用式(7)对每一个影像块进行线结构系数变换，并对其系数进行二值化处理，二值化阈值为系数均值的1.5倍，对各类影像均适用，表示为：
$$ I_{1.5}^{{\rm{Line}}} = \rho _\lambda ^{{\rm{Line}}}{(\bigcup\limits_{r = 1, c = 1}^{{N_H}, {N_H}} {I_{r, c}^{{A}}} )_{1.5}} $$ (10) 综合式(8)至(10)，得到多尺度合成子带集的数学表达式：
$$ P_{{\rm{dec}}}^{\sum\nolimits_H {} } = \bigcup\limits_{H = 1}^4 {(P_ \bot ^H(\sum\nolimits_H {(\bigcup\limits_{\begin{array}{*{20}{c}} {{{\bar \Delta }_H} = 0, 1, 2, 3}\\ {{{\bar \theta }_H} = 0, 1, 2, 3} \end{array}} {(I_{1.5}^{{\rm{Line}}})} } )))} $$ (11) 式中，ΔH和θH表示对变换后的影像进行相应的位移恢复；∑H表示按照尺度进行归一化叠加计算。
则输出影像IOut是对多尺度合成子带再进行归一化叠加计算：
$$ {I_{{\rm{Out}}}} = \sum\limits_{H = 1}^4 {P_{{\rm{dec}}}^{\sum\nolimits_H {} }} $$ (12) 2.3 直塔模型结构与实验分析
为保持建模方法研究的系统性和一般性，本节对直塔模型结构依然采用与金字塔以及混合塔模型结构类似的模块与符号表示，以便于对比说明建模方法的可行性，如图 4所示。
![]()
图 4 多尺度几何分析特征变换的直塔模型结构递归算法流程图Figure 4. Flowchart of Recursive Algorithm for Straight Tower Structural Model of Feature Transformation Based on Multi-scale Geometric Analysis图 4中特征变换1~n是采用第一种非下采样方式得到的影像多尺度、多通道特征空间，构成直塔模型结构，结合线结构稀疏表示的具体算法原理与流程，说明如下。
1) 输入影像，按照式(8)进行多尺度分解。如图 4中横轴的第1、2项处理方式，纵轴是与尺度相关的特征1~n的变换，采用非下采样方式的多尺度、多通道特征空间，构成多尺度分解子带集Pdec，即直塔模型。本实验最小影像块为16×16，最大影像为128×128。
2) 按照式(9)对Pdec进行位移操作，构成多尺度空间位移子带集PdecΔH，即对直塔模型的多尺度、多通道空间做进一步的特征变换，也是横轴第2项处理方式的继续。
3) 在多尺度空间位移子带集PdecΔH(直塔模型)中，用式(10)对PdecΔH中的单尺度层影像块逐个进行I1.5Line计算，即进行横轴第3项处理，这里采用最简单的阈值二值化处理(二类分割)，得到单尺度影像的8个空间位移子带，然后进入横轴第4项处理方式，即用式(11)对8个单尺度二值化空间位移子带归一化叠加计算，得到多尺度合成子带集Pdec∑H。
4) 进入横轴第5至7项处理，用式(12)对多尺度合成子带集Pdec∑H进行多尺度归一化叠加计算，通过多尺度、多空间位移信息的互补与融合，输出稀疏影像IOut。
基于式(7)，线结构变换系数的非下采样直塔模型中，系数的冗余性主要存在于各尺度线结构特征的空间位置附近，并在各尺度子带通道内呈现出稀疏分布，而且越往塔顶的较粗尺度层冗余性越强。由于算法采用各通道内二值分割以及各通道间归一化叠加的计算策略，不对特征空间做阶段性整合操作，计算方法简单，因此效率较好。
采用一组合成影像说明模型性能，如图 5所示；采用一组自然影像验证模型应用效果，如图 6所示。
![]()
图 5 基于脊波变换特征的多尺度直塔模型的合成影像稀疏表示Figure 5. Synthetic Image Sparse Representation for Multi-scale Straight Tower Structural Model Based on Ridgelet Transform Feature![]()
图 6 基于脊波变换特征的多尺度直塔模型的影像稀疏表示Figure 6. Image Sparse Representation for Multi-scale Straight Tower Structural Model Based on Ridgelet Transform Feature图 5为合成影像(512×512像素)及其对应的实验结果。图 5(a)、5(b)是纹理背景与两个不同尺度圆线的合成影像，图 5(c)、5(d)是圆环纹理的合成影像，图 5(e)、5(f)是两个不同尺度圆线添加20%噪声的合成影像。§2.1中，式(7)是基于Radon变换与小波分析的多尺度特征分析，它继承了Radon变换对影像形状结构的旋转不变性与平移不变性。同时，在影像多尺度分解子带构成的直塔模型中，子带内非线结构的微小结构具有方向、大小的随机性特征，而多尺度子带间的线结构则具有方向、大小的相关性、连续性以及稳定性。所以对尺度维进行多尺度子带之间的归一化系数叠加合成运算后，在增强线结构的同时，有效抵消了随机微小结构之间的累加性，即抑制了随机的各向异性特征叠加，而有效地互补增强了各子带间稳定的线结构信息，实现影像整体线结构的稀疏表示。从图 5实验结果可以看出，3幅影像都可以很好地稀疏表示出圆环的线结构，具有抗纹理或噪声干扰的能力，验证了上述分析。同时，还存在一些细节效果上的差异，由于噪声影像的随机性更强，所以它的背景抑制效果更好，如图 5(f)所示；如果影像纹理的随机性较强，则其效果也较好，如图 5(b)所示；如果影像纹理具有一定的规律性，则其纹理的微结构也会被识别出来，但是其总体强度依然远远小于线结构特征，如图 5(d)所示。
图 6是景物影像(512×512像素)及其对应的实验结果。图 6(a)、6(b)是脸谱挂件影像，图 6(c)、6(d)是航空影像，分辨率为0.5 m，图 6(e)、6(f)是卫星影像，分辨率为1.5 m左右。从图 6的实验结果可以看出，3幅影像都可以很好地稀疏表示出影像的线结构，对非线结构区域都具有较好的抑制能力。由于脸谱中间凸起，造成影像平面受光不均匀，通常在高亮区域会对线结构目标识别产生影响，但是在多尺度直塔模型结构中，依然可以获得较好的稀疏表示效果，如图 6(b)所示；一般算法很难同时适用于直线与曲线稀疏表示问题[21]，但是在多尺度直塔模型结构中，航空影像中同时包含直线与曲线结构，其效果依然很好，如图 6(d)所示；对卫星影像中包含多种尺度线结构的情况，其效果基本一致，都可以较好地稀疏表示，如图 6(f)所示。对比图 6(d)、6(f)可以看出，影像分辨率越高，其稀疏表示的线结构越细，因此，多尺度直塔模型结构更适应于高分影像线结构的稀疏表示与提取。
3 结语
3种形式多尺度分析的塔式模型结构可以有效解决高分影像的多尺度信息稀疏表示与提取问题。其中金字塔模型的实践已经比较成熟[8, 15, 17]，并主要解决区域分割问题。本文对直塔模型进行了比较详细的研究，并主要解决线结构特征提取问题，实验验证了其有效性。混合塔模型是在这两者基础上融合提出的，因此其理论与实践都具有可靠的基础，可以进一步开展技术细节的研究。在多尺度间的信息整合与特征融合问题中，将单尺度(像素级)和多尺度(超像素级)分析技术混合集成在塔式模型结构中，可以增强多尺度分析的灵活性，并为更深入的技术研究提供较好的模型基础与参考。
本文研究存在的问题主要为，在直塔模型实验中，对影像双线结构的提取与表示需要继续研究，双线的间距较近时只能表示其中较强的一条，如图 6(b)、6(d)所示；在两条线的交叉处只能表示其中较强的一条，如图 6(d)所示，其原因是由于式(6)只设计了识别单线结构的模式，下一步将对此进行改进。
-
![]()
图 1 卫星钟差实时估计采用IGS测站分布图
Figure 1. Distribution of the IGS Stations Selected for Real-Time Satellite Clock Estimation
![]()
图 2 实时估计初始钟偏差与ESA事后产品值互差
Figure 2. Differences of the Initial Clock Biases Between the Estimates and the ESA Final Products
![]()
图 4 忽略码IFB的GLONASS实时钟差偏差
Figure 4. Biases of the GLONASS Real-Time Satellite Clock Estimates Ignoring IFBs
![]()
图 5 顾及码IFB的GLONASS实时钟差偏差
Figure 5. Biases of the GLONASS Real-Time Satellite Clock Estimates with IFBs Compensation
![]()
图 6 GPS/GLONASS实时钟差精度统计
Figure 6. Comparison of RMS of the GPS/GLONASS Real-Time Satellite Clock Estimate Biases with/without IFB Compensation
![]()
图 7 不同模式下RT-PPP定位结果
Figure 7. RT-PPP N/E/U Biases with Different Processing Modes Using the Estimated Satellite Clocks
![]()
图 8 各测站3D偏差RMS值
Figure 8. RMS of RT-PPP 3D Biases with Different Observation Modes Using the Estimated Satellite Clocks
表 1 用于钟差估计结果分析的卫星信息
Table 1 Satellite Information Selected for Satellite Clock Estimate Analysis
系统 卫星 卫星及频标类型[频率号] GPS G03 Block IIA Cs G26 Block IIA Rb G11 Block IIR-A Rb G19 Block IIR-B Rb G07 Block IIR-M Rb G24 Block IIF Cs G01 Block IIF Rb GLONASS R14 Cs[-7] R06 Cs[-4] R22 Cs[-3] R13 Cs[-2] R16 Cs[-1] R15 Cs[0] R05 Cs[1] R24 Cs[2] R23 Cs[3] R21 Cs[4] R07 Cs[5] R08 Cs[6] 
下载: 导出CSV
-
[1] Hauschild A, Montenbruck O. Kalman-Filter-Based GPS Clock Estimation for Near Real-Time Positioning[J]. GPS Solutions, 2009, 13(3):173-182 doi: 10.1007/s10291-008-0110-3
[2] Laurichesse D, Mercier F, Berthias J P. Real Time Precise GPS Constellation Orbits and Clocks Estimation Using Zero-Difference Integer Ambiguity Fixing[C]. The 22nd International Technical Meeting of the Satellite Division of the Institute of Navigation, Savannah, USA, 2009
[3] Zhang X H, Li X X, Guo F. Satellite Clock Estimation at 1 Hz for Realtime Kinematic PPP Applications[J]. GPS Solutions, 2011, 15(4):315-324 doi: 10.1007/s10291-010-0191-7
[4] 张小红, 李星星, 郭斐, 等.基于服务系统的实时精密单点定位技术及应用研究[J].地球物理学报, 2010, 53(6):1308-1314 http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201006010.htm Zhang Xiaohong, Li Xingxing, Guo Fei, et al. Service-Based Real-Time Precise Point Positioning and Its Application[J]. Chinese Journal of Geophysics, 2010, 53(6):1308-1314 http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201006010.htm
[5] 李浩军, 王解先, 陈俊平, 等.基于GNSS网络的实时精密单点定位及精度分析[J].地球物理学报, 2010, 53(6):1302-1307 http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201006009.htm Li Haojun, Wang Jiexian, Chen Junping, et al.The Realization and Analysis of GNSS Network Based Real-Time Precise Point Positioning[J]. Chinese Journal of Geophysics, 2010, 53(6):1302-1307 http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201006009.htm
[6] Hadas T, Bosy J. IGS RTS Precise Orbits and Clocks Verification and Quality Degradation Over Time[J]. GPS Solutions, 2014, 19(1):93-105 doi: 10.1007/s10291-014-0369-5
[7] Rossbach U. Positioning and Navigation Using Russian Satellite System GLONASS[D]. Germany:University of the Federal Armed Forces Munich, 2001
[8] Yamada H, Takasu T, Kubo N, et al. Evaluation and Calibration of Receiver Inter-channel Biases for RTK-GPS/GLONASS[C]. The 23rd International Technical Meeting of the Satellite Division of the Institute of Navigation, Portland, USA, 2010
[9] Meindl M, Dach R, Jean Y. International GNSS Service Technical Report 2011[R]. Switzerland:University of Bern, 2012
[10] Shi C, Yi W T, Song W W, et al. GLONASS Pseudorange Inter-Channel Biases and Their Effects on Combined GPS/GLONASS Precise Point Positioning[J]. GPS Solutions, 2013, 17(4):439-451 doi: 10.1007/s10291-013-0332-x
[11] 刘志强, 王解先, 段兵兵.单站多参数GLONASS码频间偏差估计及其对组合精密单点定位的影响[J].测绘学报, 2015, 44(2):150-159 doi: 10.11947/j.AGCS.2015.20130800 Liu Zhiqiang, Wang Jiexian, Duan Bingbing. Estimation of GLONASS Code Inter-Frequency Biases with Multiple Parameters Based on a Single Station and Its Impact on Combined Precise Point Positioning[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(2):150-159 doi: 10.11947/j.AGCS.2015.20130800
[12] Springer T. IGSMAIL-6107:ESA/ESOC Ultra-Rapid and Rapid GNSS Solutions[EB/OL]. http://igscb.jpl.nasa.gov/pipermail/igsmail/2010/006178.html, 2010-03-07
[13] Song W W, Yi W T, Lou Y D, et al. Impact of GLONASS Pseudorange Inter-Channel Biases on Satellite Clock Corrections[J]. GPS Solutions, 2014, 18(3):323-333 doi: 10.1007/s10291-014-0371-y
[14] Cai C S, Gao Y. Modeling and Assessment of Combined GPS/GLONASS Precise Point Positioning[J]. GPS Solutions, 2013, 17(2):223-236 doi: 10.1007/s10291-012-0273-9
[15] 孟祥广, 郭际明. GPS/GLONASS及其组合精密单点定位研究[J].武汉大学学报·信息科学版, 2010, 35(12):1409-1413 http://ch.whu.edu.cn/CN/abstract/abstract823.shtml Meng Xiangguang, Guo Jiming.GPS-GLONASS and Their Combined Precise Point Positioning[J]. Geomatics and Information Science of Wuhan University, 2010, 35(12):1409-1413 http://ch.whu.edu.cn/CN/abstract/abstract823.shtml
[16] Ge M, Chen J, Douša J, et al. A Computationally Efficient Approach for Estimating High-Rate Satellite Clock Corrections in Realtime[J]. GPS Solutions, 2012, 16(1):9-17 doi: 10.1007/s10291-011-0206-z
-
期刊类型引用(3)
1. 周会超，郑晋军，王海红，毕少筠. 地月空间天基时空基准的构建与传递. 中国空间科学技术(中英文). 2024(03): 51-59 . 
百度学术
2. 杨阳，稂龙亚，吕玉祥，张孜豪. 一种用于电力终端的BDS定位漂移校正算法. 信息技术. 2019(02): 125-128 . 百度学术
3. 张继海，董绍武，卢鋆，袁海波，广伟. 北斗卫星天线相位中心改正在PPP中的应用研究. 仪器仪表学报. 2019(07): 90-96 . 百度学术
其他类型引用(8)
计量
- 文章访问数: 2290
- HTML全文浏览量: 172
- PDF下载量: 432
- 被引次数: 11
