在现代科学计算中，求非线性方程或极值问题是非常复杂的工作，传统的求解非线性问题通常是先将其转换成线性问题，然后解析求解，或者采用迭代计算，如Newton-Raphson算法、最陡梯度算法^[1-2]，但在计算中需要进行微分计算，或者说，需要计算函数的导数，然而，实际中有很多非线性方程无法进行微分计算，因此无法采用传统的微分算法求解。幸运的是，采用无微分算法可以克服这种困难。

目前，有很多无需进行微分就可以对非线性方程或极值求解的算法，如Nelder-Mead Simplex算法 (或称为NMS算法)，Taxi Cab算法，Powell算法等^[3-4]。这些算法的共同点是，通过迭代计算求解，很容易在计算机上实现，为解非线性问题提供了极大方便。这需要考虑两个问题：一是计算速度，二是计算精度。如果这两个方面都达到最优，则计算方案最佳。但很多无微分算法都存在某方面的缺陷，仍然需要做进一步改进。

Nelder-Mead单纯算法作为一种简单的无微分极值算法，在数值计算中得到了广泛应用，如在Matlab中使用了这种算法^[5]，但这种算法存在一定不足，如收敛速度慢，计算精度不高，也需要进行改进。

在气象科学研究中，需要进行各种计算，如通过非线性曲线拟合寻找各个气象要素变量之间的关系^[6-7]；根据观测资料建立经验公式，并确定其中的参数^[8-9]；解非线性方程组^[10]等，这些计算非常复杂，很多问题是非线性的，很难使用传统的计算方法求解。传统上，很多气象物理方程的求解，需要对方程进行微分计算，这限制了一些问题的求解，而现代无微分算法则可以克服这方面的不足。为此，本文基于对NMS算法的改进，作为一个此类算法的应用实例，为气象科学计算提出一种新的思路。

Nelder-Mead单纯形算法是Nelder等^[11]提出的一种犹如下山一样查找一个函数的局部最小值点的方法，其中函数变量可以是多元的。关于NMS算法的详细介绍可参见文献[4, 11]。

在几何上，单纯形或者N-单纯形是和三角形类似的N维几何体。精确地说，单纯形是某个N维以上的欧几里空间中的N+1个仿射无关的点的集合的凸包。例如，0-单纯形就是点，1-单纯形就是线段，2-单纯形就是三角形，3-单纯形就是四面体，…，如此类推，N-单纯形就是N+1面体。

NMS算法是使用单纯形的思想，对于N维空间可以构成N+1顶点的多面体。亦即对于拥有N个变量的函数，可以构成N+1个点。

在计算开始时，先计算出每个点的函数值，查找出最大值的点 (G) 和最小值的点 (S)，再计算除G点外的其他点的中点 (M，一般取平均)，在GM连线方向上，计算G的反射点、扩展点、压缩点。若没有找到这些点，则在G至S的线段上计算一个收缩点若在这些点中，存在一个最佳点，则使用这个点作为一个新的有效点，替代G点，重现构成一个新的N+1多面体。如此进行，直到找到函数的最小值点。

从NMS算法的计算过程可知，NMS算法计算较繁琐，需要进行各种判别，计算精度不高，收敛较慢。在计算这些点时，因为线段长度固定，很容易陷于局部收敛，出现计算早熟，这很像遗传算法中出现的早熟现象。

目前，有很多人提出了对NMS算法的改进方案^[12-14]，但没有从根本上解决NMS算法的收敛问题。

为了避免原有NMS算法中存在的问题，可以将此算法中的反射点、扩展点、压缩点以及收缩点的计算转化成使用适应参数计算最佳点的方法，这与最陡梯度算法和Powell算法中计算下一个迭代点的方法相同，但此时查找最小值点的方向为NMS算法的下山方向，即最大值的点G与中点M的连线方向 (图 1)。

图 1

图 1. 改进的NMS下山方向设置示意图 Fig 1. Searching direction setting for the improved NMS

这种设置方向类似于梯度方向，但梯度算法中需要分别对每个变量计算偏导数。

对于一个拥有N个变量的函数，可以构成一个有N+1个顶点的单纯形，其中点坐标可以表示为

(1)

式 (1) 中, i为顶点索引号，k为函数中变量的索引号，m为中点索引号，g为最大函数值所在点的索引号。从最大函数值的坐标点G至中点M的方向向量可以表示为

(2)

式 (2) 中, f_M为中点M的函数值，f_G为顶点G的函数值。

计算下山方向向量后, 设置一个适应参数α，结合G和向量E_k，可以是用线性插值方法计算出下一个迭代点, 用它来替代原来NMS算法中的反射点、扩展点、压缩点或收缩点。

(3)

将式 (3) 代入到极值函数中，则可以得到一个以α为自变量的函数f (α)。计算此函数取极值时的α，然后，再一次将α代入式 (3)，则可以得到迭代点坐标。重复此计算，直到查找到极值点为止。判别极值终止的一般方法是，先计算这N+1个顶点的函数值的均方差，然后与一个事先给定的误差值进行比较，若小于此误差值，则迭代终止，否则继续进行迭代。

可以使用单元函数极值的无微分算法来计算α。这种算法有很多，如黄金分割搜寻法、斐波那契搜寻法等，因此，改进后的NMS算法仍然为一种无微分算法。

这种引用适应参数计算下一个迭代点的方法，与Powell算法中的相同，但搜索极值点的方向不相同。本改进算法的搜索方向是准梯度方向，而Powell算法的搜索方向是先依次沿着已知的N个共轭方向搜索，得一个最佳点，然后沿着本轮迭代的初始点与该最佳点连线方向进行搜索，求得这一阶段的最佳点。再用最后的搜索方向取代前N个方向之一，开始下一阶段的迭代。由此可见，与Powell算法相比，本改进算法搜索过程较简便。

NMS算法作为一种无微分算法，可以应用于气象上的非线性问题的计算。在传统的数值计算中，通常是将非线性方程转化为线性方程，然后进行求解。然而，实际问题中的非线性方程式是多样的，不易求解，此时可以采用无微分极值算法来完成此项工作。

曲线拟合是气象上经常要做的工作。例如，对气象观测资料的时间或空间变化曲线进行拟合，得到一个经验公式^[8]，用于数值模式中。通常这些曲线不能使用线性公式拟合，否则所得到的经验公式不能准确地描述实际物理量变化。同样，如果采用有限的、可以转换成线性公式或非线性公式，如指数公式等，也可能达不到实际要求。因此，可以根据需要设定拟合公式模型，然后确定其中的系数。

使用NMS算法进行非线性曲线拟合的具体方法：① 根据实际观测资料的变化规律，给定一个经验拟合公式模型 (可以根据需要给定)，例如, z=f(a, b, c, x, y)，其中f是一个具体的表达式。② 设计最小二乘法计算式：

(4)

式 (4) 中, x_i, y_i和z_i分别为观测数据，N为观测资料的组数。③ 使用NMS算法计算函数式 (4) 最小值，确定其中的系数，即可以得到所需要经验拟合公式。

可以看出，使用类似于NMS算法的无微分算法进行非线性曲线拟合，不需要像传统拟合方法那样，求解线性方程组。不仅设定拟合公式有一定的自由性，而且算法非常简单。另外，也可以使用Matlab进行曲线拟合计算。

无论是线性还是非线性方程组，其中不包括微分方程组，都可以使用无微分求极值的算法进行求解。例如，对于气象学上的非线性动力、化学方程组^{[10, 15]}，可以使用NMS算法，或其他求解无微分极值的方法计算。

对于下面的一个方程组：

(5)

令δ₁=f₁²(…)，δ₂=f₂²(…)，δ₃=f₃²(…)，δ₄=f₄²(…)，并构造以下函数：

(6)

当δ₁=δ₂=δ₃=δ₄=0时，以上函数达到最小值零，此时的解也是式 (5) 的解。也就是说，求解方程组 (5) 可以转换成计算函数式 (6) 的极值。这样，就可以使用NMS算法求解非线性方程组。

在现有的气象数值模式中，涉及到很多方程组求解计算。基于以上方法，可以直接使用无微分极值算法，求解气象动力学方程组。例如，用NMS算法代替牛顿迭代法。

为了获得气象物理量变化的经验公式，需要对气象观测资料进行曲线拟合，在数值模式计算中，也需要进行曲线拟合。在数值计算中，如果某个计算公式非常复杂，需要花费大量计算时间，此时可以对它的计算结果进行拟合，以提高计算效率。当然，也可以使用查表的形式来提高计算效率，如概率统计中的t-函数表。

单纯形初始化对优化计算结果有一定影响，有很多关于如何确定这种初始值的研究^[16]。合理地给出初始值，可以提高函数收敛速度，并避免计算结果局地收敛的发生。

在使用改进的NMS算法或原NMS算法计算时，有两种类型的初始值需要设定，一是变量的初值，二是单纯形顶点的初值。对于变量的初值，也是迭代求解的初值。这种初值一般是猜想值，根据实际情况，应合理地给出估计值，然后根据这个初值进行迭代求解。

对于单纯形顶点，由于顶点不能重合，因此初始值不能相同。一般是将初始单纯形的顶点在空间上看成是均匀分布的，即各个顶点之间的距离是相等的。

如何构建在空间上具有均匀分布顶点的多面体，对于N个变量，假定所对应的单纯形N+1个顶点的单位向量为

(7)

式 (7) 中，a为待定常数。从式 (7) 可以看出，当j, k≠0时，；而当j≠0时，|v_j-v₀|=(1-2a+N·a²)^1/2。也就是说，除了v₀外，其余的点之间的距离都是相等的。如果令 (1-2a+N·a²)^1/2=，则所有这些点之间的距离均相等。

通过解方程，可以得到

(8)

如果通过空间平移，将v₀变成v₀=(0, 0, …, 0)，并对顶点单位标准化，则可以得到初始单纯形的N+1个顶点的标准单位向量：

(9)

式 (9) 中，，。

设由初始变量值构成的一个空间点为X₀=(x₁, x₂, …, x_N)，每维空间上的步长都为s，则可以构成单纯形N+1个顶点的初始坐标，即

(10)

对陆面过程中的一些参数，如零平面位移 (D₀)、地面与冠层空隙空气之间的空气动力阻 (R_d) 等所表示的是植被类型、植被高度和叶面积指数 (I_LA) 的函数，计算过程非常复杂^[17-18]。在早期的SSiB陆面过程模式中，由于I_LA取值不是连续的，D₀与R_d的值是通过一个查找表获得。如果I_LA的值是连续变化的，D₀与R_d的值可以使用一个经验拟合公式来表示。

为了验证改进的NMS算法计算能力，这里先根据原始的、计算近地面参数D₀和R_d的函数，对于某种特定植被类型，通过改变不同的植被高度、叶面积指数，生成用于曲线拟合的数据，其中Bayesian函数的计算可以借用一些软件中 (如Microsoft Visual Fortran等) 提供的IMSL函数库。

图 2显示的是对于热带阔叶林关于地面零平面位移D₀和空气阻力原始计算数据点及其拟合曲线。在本示例中，分别取植被高度Z₂=4 m，Z₂=20 m和Z₂=35 m，并按线性等间隔取I_LA的值，由原始的计算D₀和R_d计算公式，得到原始数据点，由图中的点表示。在进行曲线拟合之前，绘制原始数据点，通过原始数据的变化情况，合理地设计拟合D₀和R_d的数学模型，分别设置为

(11)

(12)

图 2

图 2. 阔叶林近地面参数曲线拟合结果 (a) 零平面位移D0, (b) 空气动力阻力Rd Fig 2. Simulations on the variations of ground parameters for the broad leaf tree (a) zero plane displacement D0, (b) aerodynamic resistance Rd

式 (11) 和 (12) 中, c₁, …, c₇和b₁, …, b₇分别代表拟合D₀和R_d变化曲线的经验公式中的系数。利用式 (4) 和改进的NMS算法，可以确定这些系数 (具体数值略)，从而可以得到D₀和R_d计算的经验公式。再通过拟合的经验公式，可以绘制如图 2中所示的拟合曲线，由图 2中的光滑曲线表示。

在实际应用中，拟合方程中参数一般是通过物理试验来确定的，所要做的是根据这些参数查找一个经验公式，确定经验公式的系数。系数的个数根据曲线拟合情况而定。所选系数个数的原则是拟合曲线与原始点越吻合越好。

从设计的计算D₀和R_d的数学模型可以看出，这两个公式的结构是非线性的，很难转化为线性的计算公式来求解线性方程组。而使用极值算法，如NMS算法和改进的NMS算法，可容易确定其中的系数。从图 2可以看出，D₀和R_d对植被高度和I_LA变化的拟合曲线与原始的数据点变化十分接近，说明使用改进的NMS算法能够很好地对本示例中D₀和R_d的变化进行拟合。

本示例说明了如何使用改进的NMS算法对非线性变化曲线进行拟合计算。也可以采用类似计算步骤，使用其他的极值算法，如现代遗传算法、模拟退火算法和蚁群算法等，进行非线性变化曲线拟合。

NMS算法是一种不需求解微分、可以查找多元函数最小值点坐标的技术，在数值计算中，得到了广泛应用，也可以在气象科学研究领域中得到应用。但它存在一些缺点，如收敛速度较慢、容易产生局地收敛等。为此，本文提出了一种对NMS算法改进的方法。在改进的计算方案中，最小值点的搜索方向仍然为原NMS算法的搜索方向，但通过引入适应参数α，来计算一个新顶点，用于替换单纯形多面体中具有最大函数值的顶点。这种计算方法与最陡梯度算法和Powell算法中确定下一迭代点的方法相同。

改进的NMS算法与Powell算法有一定相似性，但它们搜索最小值点的方向不同。通过比较可以看出，改进的NMS算法搜索过程简单，因而搜索路径快捷。

无微分极值算法有很多种，NMS算法只是其中的一种。本文以NMS算法为例，简要介绍了在气象上如何使用无微分极值算法进行非线性曲线拟合，以及如何使用它来求解非线性方程组。本文还介绍了如何对单纯形计算进行初始化。

为了检验改进后NMS算法的计算能力，以及说明具体使用无微分极值算法进行曲线拟合，本文以陆面过程中的参数为例，对热带阔叶林地面零平面位移和空气阻力参数进行了经验曲线拟合，结果表明，改进的NMS算法具有较高的非线性曲线拟合能力。