1. 引言
由于非线性系统对其系统参数和初始值的极端敏感性，任何观测噪声和动力噪声都会在在非线性系统的作用下以指数级放大。因此传统的滤波或去噪算法非但难以奏效，往往还会改变非线性序列的固有特性，使其不再反映原系统的任何属性。当系统独立状态变量和相关变量都有误差时，处理不当还将带来所谓的变量误差问题 [1] 。非双曲型非线性系统同宿切面点和同宿横截点的存在，使得这一问题变得更具复杂性和挑战性。非线性系统的流形分析对于非线性复杂大系统分析设计 [2] 具有非常重要的理论和现实意义。
在非线性时间序列的去噪或轨迹重影算法中，经常要考虑系统的双曲特性。系统稳定流形与不稳定流形的交点称为同宿横截点，稳定流形与不稳定流形平行的点称为同宿切面点。处处双曲型非线性系统是指，系统的稳定流形和不稳定流形永远不会平行，因此没有同宿切面点，也就是说对系统的每一个点，系统稳定流形和不稳定流形间的夹角要远大于零。非双曲型非线性系统是指有同宿切面点的非线性系统。然而真正意义上的同宿切面点是很少出现的，通常所说的同宿切面点也多指准同宿切面点。若在某点处，其局部不稳定流形中的单位向量
与其稳定流形中的单位向量
间的夹角
不大于某指定的较小角度
，则称该点为
准同宿切面点。
2. 轨迹重影基本问题
考虑由下式确定的自治非线性动力系统
(1)
其中
为由系统状态变量组成的向量，
为系统的动力学函数。该系统的受扰观测向量如下
(2)
其中
和
分别表示实际序列向量和观测噪声向量。设
为来自于系统内部的热噪声，如下式所示
(3)
式中
通常也被称为系统动力噪声。
在观测噪声下，系统的单步动力误差如下
(4)
其中
、
分别是系统在
时刻和
时刻的观测向量。对式(4)而言，若满足
(5)
则称序列
为动力方程
所决定的一个确定性序列，其中
为观测序列的长度。
轨迹重影就是要在机器精度将所有
降至0，这是一个典型的高维非线性寻根问题。
3. 同宿切面点对轨迹重影算法的影响分析
3.1. Newton类方法的不稳定性分析
求解高维寻根问题的一个实用方法就是Newton-Raphson算法。由于Newton-Raphson算法用一阶泰勒展开进行非线性系统的局部线性化，当线性化不够准确时，算法的不稳定性便在所难免。而且该方法对解的初始猜测值有着较为苛刻的要求 [3] [4] 。当初始猜测值在实际解的邻近区域之内时，解的迭代过程为
(6)
其中
为
维动力误差向量，
和
分别为迭代前后的观测序列向量，
为如下
阶矩阵
(7)
其中
为系统
在第
个采样点的Jacobian矩阵。由于
非方阵，为使式(6)可解，通常辅以约束条件
(8)
对非双曲型非线性系统而言，Schreiber给出了一种在上述约束条件下求解式(6)简单有效的流形分解 [5] [6] 方法。该方法所依据的条件是能将系统在每一点处的切空间分解为稳定子空间和不稳定子空间。该方法求解式(6)的迭代公式如下
(9)
其中上标
和
分别代表系统在每一点处的稳定方向和不稳定方向。
设
点为
准同宿切面点，若
，则矩阵
的条件数有下界，且该下界为点
局部稳定流形方向
与局部不稳定流形方向
间夹角的函数。不难证明 [7]
(10)
结合式(6)可得
(11)
然而大多数非线性系统并不具备处处双曲特性，这意味着观测序列中某些点处的稳定流形和不稳定流形将是平行的，可见因系统同宿切面点或
准同宿切面点的出现，使得矩阵
成为一个坏条件或降秩矩阵。这将使得Newton类算法在这些同宿切面点处失去稳定性，这是造成Newton类算法不稳定的根本原因。
3.2. Gradient Descent算法对同宿切面点的免疫性
轨迹重影问题的另一种求解办法是Gradient Descent算法 [8] [9] 。Gradient Descent算法把轨迹重影问题看成是一个最小化问题。由式(4)可得该最小化问题的代价函数为
(12)
Gradient Descent算法的迭代公式如下
(13)
其中
是算法的迭代步长，
为代价函数
的梯度向量，不难证明 [6]
(14)
由上式可以看出，Gradient Descent算法中代价函数的梯度下降方向总和
的奇异方向垂直，即梯度下降方向在
的奇异方向上的映射分量为
，这是Gradient Descent算法不受系统同宿切面点影响的关键因素与根本原因。因此当代价函数
梯度方向变得近乎平坦，即平行于同宿切面方向时，算法在这些点处将不能以线性率减少误差，但不会不稳定。
其实，还可以从另一角度，即从Gradient Descent算法本身来解释此现象。最小化策略往往是以找到一个噪声水平更低的序列为目标，算法本身并不要求受扰序列与系统的一个确定序列在同一个线性邻域内，因而算法不会因非线性系统同宿切面点的出现而被迫变得不稳定，这是对Gradeient Descent算法免疫于系统同宿切面点影响的最直观解释。
3.3. 仿真分析
本节通过仿真实验来验证上述分析的正确性。被广泛研究的非双曲型非线性系统Hénon映射的动力学可由下式确定
(15)
上式中当
、
时系统将处于混沌态。由于此处仅进行Newton方法和Gradient Descent方法的性能比较，故假设系统的数学模型已知，无需进行模型选取与参数估计。以
，可得变量
长度
的序列向量
，取观测向量
的分量为
，
(16)
其中
为向量
的第
个分量，
为高斯白噪声向量
的第
个分量，
的均值为0，标准差取序列
标准差的1%，长度与
相同。以时间延迟1、潜入维2重构序列
和
的相空间 [10] [11] 分别如图1(a)、图1(b)所示。
为评估算法的收敛性能，定义如下两个判断准则
(17)
(18)
其中式(17)中
表示系统的阶数，
为单变量观测序列在第
点处的单步动力误差，
和
分别为该单变量在第
时刻的观测值和实际值。式(17)用来描述系统的均方动力误差和，式(18)用来描述观测序列与实际序列的均方误差和。
实验1：Newton方法的轨迹重影仿真实验
将Newton方法用于上述受扰序列
，迭代10次后的结果如图2所示，其中(a)、(b)两图的纵轴为对数坐标。由图2可见，算法在迭代10次后快速收敛，最大动力误差约为
。图(a)中估计轨迹的波峰处为序列的同宿切面点或准同宿切面点，在这些点处，矩阵
是坏条件的，算法不仅在这些点处无效，甚至是不稳定的，使得误差不仅不降低，反而因序列的非线性而被放大。而在其它点处，算法在机器精度内逼近了原未受扰序列。从观测误差的角度分析，区间
和
上的去噪序列即为原未受扰序列的部分重影轨迹。
实验2：Gradient Descent方法的轨迹重影仿真实验
由Gradient Descent方法对上述观测序列
进行轨迹重影的结果如图3所示，其中(a)、(b)两图的纵轴为对数坐标。由于Gradient Descent方法只有线性阶的收敛速度，故为达到图2的估计精度，算法需运行400,000次（还可进一步迭代下去），可见收敛速度是相当慢的。但由图3(b)可以看出，动力误差几乎全部收敛为零，该算法几乎不受系统同宿切面点的影响，因此所得重影序列更符合系统动力学特性，因而也更接近实际序列
。或者说仅从序列的确定性，即动力误差的角度分析，上述去噪序列就是原未受扰序列
。观察图3(a)和图3(b)，可以看到一个有趣的现象，即动力误差和观测误差是不完全同步的，也就是说并不是在所有同宿切面点处都存在动力误差，而且从图3(c)还可以看出，均方误差和达到最小时算法并没有收敛，这个有趣问题目前还没有理论解释。
旨在弥补上述两类算法的不足，文献 [12] [13] 提出一种既拥有Gradient Descent方法的良好稳定性又不失Newton方法快速收敛性的新方法，并改良了同宿切面点和准同宿切面点处矩阵的条件数，进而可用线性阶的Cholesky分解技术对其进行求逆运算，快速稳定地在机器精度内实现了非双曲型非线性系统观测序列的轨迹重影。
4. 同宿横截点对轨迹重影算法的影响分析
试验发现对上述
的观测序列，仍使用均值为零、均方差为0.01的高斯白噪声，使用文献 [12] [13] 描述的去噪算法，对于此噪声型的20条不同实现序列，唯有序列第628点左右的估计结果出现了较大的差别，其余点的估计结果如图4所示。图5给出了第628点左右11个点在上述20条不同噪声实现序列下的去噪结果。有趣的是，图5中第628点处的估计结果清晰地被分为上下两部分，且这两部分间的平均距离与噪声的均方差十分接近。注意，图5的纵轴为对数坐标，同样的误差被放缩的程度不同，实际上两部分有着大致相同的形状。第628点的估计结果及其全域稳定流形与不稳定流形示于图6之中。
(a) (b)
(c) (d)
Figure 2. Estimation results of Newton methods after 10 iterations. (a) Measurements error, dotted: original, solid: after iterations; (b) dynamic error, dotted: original, solid: after iterations; (c) Emeas against the iterate number k; (d) Edyna against the iterate number k
图2. Newton方法类方法迭代10次的结果。(a) 观测误差，点线：原始观测噪声，实线：迭代后；(b)系统动力误差，点线：原受扰序列，实线：迭代后；(c) 均方观测误差收敛效果；(d) 均方动力误差收敛效果
(a) (b)
(c) (d)
Figure 3. Estimation results of Gradient Descent methods after 400,000 iterations. (a) Measurements error, dotted: original, solid: after iterations; (b) dynamic error, dotted: original, solid: after iterations; (c) Emeas against the iterate number k; (d) Edyna of the preceding 5000 points against the iterate number k
图3. Gradient Descent类方法迭代400,000次的结果。(a) 观测误差，点线：原始观测噪声，实线：迭代后；(b) 系统动力误差，点线：原受扰序列，实线：迭代后；(c) 均方观测误差迭代效果；(d) 前5000点的均方动力误差迭代效果
Figure 4. Results of measurement error after 12 iterations
图4. 序列的去噪结果
Figure 5. Estimation results of eleven points around the 628th points
图5. 第628点左右11点的估计结果
Figure 6. Spatial plot of the 628th points of the noise-reduced trajectories of Figure 6 (+) and the 628th point of the correct trajectory (□) with its stable (dotted) and unstable (solid) manifolds
图6. 第628点的估计结果及其稳定流形和不稳定流形。“□”为其实际位置，“+”为其估计位置，点线：不稳定流形，实线：稳定流形
从图6可以看出，此时第628点为同宿横截点。这两个同宿横截点间的距离为0.0095，和噪声的均方差0.01十分接近。这样算法选择通过哪一点完全由在该点的噪声幅度所决定，或者说仅仅依据噪声的幅度，算法无法判断其中哪一点为系统的真实轨迹点。因噪声序列来自于同一种分布的噪声模型，且通过该点的噪声是随机的，故算法各以50%的概率通过系统真实轨迹点和同宿横截点。仿真中20条噪声序列，其中11条通过系统真实点，另外9条则选择了另一个同宿横截点。试验表明当噪声的均方差降至0.005时，上述现象不复存在，算法选择错误同宿横截点的概率已变得非常小，几乎可以忽略不计。由此可知，由于同宿切面点和同宿横截点的出现，仅当噪声水平小于系统同宿横截点距离的最小值时算法才是真正有效的，因此即便是不受系统同宿切面影响的算法，当噪声水平相当于或大于系统同宿横截点间距离的最小值时，算法至少在这些点处是无效的，这或许就是Grebogi和Walker方法 [14] [15] 在某些点处失败的原因。故为得到非双曲型非线性系统的重影轨迹，须要求噪声的均方差小于同宿横截点间距离的最小值。因此若事先比较噪声均方差与系统同宿横截点间距离最小值两者间的关系，即可确定算法可成功应用的范围。由于同宿横截点间距离的最小值一般都非常小，所以由受较大幅度噪声序列估计系统重影轨迹一般是不可实现的，现有文献只是单纯地将此现象归因于非线性系统对初值的敏感性，本文量化地给出了新的解释。
下面再讨论一下上述结论与噪声分布之间的关系。由于系统的真实点和其同宿横截点距离很近，算法选择哪一点完全由噪声的具体实现而定。对本文实验所用正态分布噪声而言，当噪声的均方差小于系统同宿横截点间的距离时，算法已能高概率判断出其真实点，而不通过其同宿横截点。然而试考虑如下情形的随机噪声
(19)
其均方差为3.1623e+003，远远大于系统同宿横截点间的距离0.0095，它显然不符合“噪声均方差小于同宿横截点间最小距离”这样的结论，然而算法对这样的噪声是成功的。因为去掉点“10000”后，该噪声的分布为正态的，噪声的均方差为“0.0043”，小于0.0095，也就是说这样的噪声只是在某一时刻使得算法无所适从，或者说噪声的分布只会影响到算法通过系统真实点的概率。显然本文所得“噪声均方差小于同宿横截点间距离”结论最适合于均值为0、均方差小于同宿横截点间距离的均匀分布噪声，因为对上述噪声，算法可以概率1通过其真实轨迹点，而完全不受系统同宿横截点的影响。因而可认为上述结论的成立和噪声的分布是有关的。然而测量误差通常可被看作高斯分布，故该结论仍具有重要的现实意义，为研究一些原无法解释的现象提供了新的理论依据。因为该结论对所有拥有同宿横截点的非线性系统都成立，故又具有普遍的理论意义。
5. 结论
定性的分析了同宿切面点可能对算法的稳定性和快速性所造成的影响，研究了同宿横截点间的距离和干扰噪声均方差二者间的定量关系，得出了仅当噪声水平小于系统同宿横截点距离的最小值时，非双曲型非线性系统的轨迹重影或去噪算法才是真正有效的这一重要结论，为同类算法提供了新的启发点。