基于贪心算法的气象修正表自动分区拟合算法设计

赵明阳，杨揆，王璟旻，陈捷

(1.中国船舶集团有限公司系统工程研究院，北京 100094；2.海军装备部，陕西 西安 710030)

火炮实际射击时的弹道、气象条件和标准条件往往不一致，即非标准射击条件。这些偏差使得按基本射表计算的射击诸元击发的弹丸可能不能命中目标，为保证实际射击条件的计算，通常需要进行弹道、气象等非标条件的修正[1-2]。为避免计算机存储大量气象修正表数据以便快速准确计算气象修正量,目前常用连续解析函数尽可能准确地逼近气象修正表，然后把这些函数存入计算机代替气象修正表[3-5]。实际逼近过程中，有时无论怎样提高解析函数阶数都无法通过一组多项式系数在整个定义域内计算出满足精度要求的修正量，此时需要分区处理[6]。不同类型的气象参数，分区方式也不尽相同，区域的划分是一个动态过程，对数据量庞大的修正表进行人工分区需要进行反复的试探，效率低下，不易控制分区数量和残差，目前业内尚未对自动分区方法进行过深入研究。笔者设计了基于贪心算法的自动分区拟合算法，实现了气象修正表逼近函数系数的自动分区计算，并且残差可控，确保了逼近准确性。

1 算法设计

逼近气象修正表的主要工作是对修正表的修正量离散数据进行拟合，求解逼近函数各个高阶多项式的系数[7-8]，需要根据水平距离(或斜距离)或高度中的一个自变量分段拟合，或者选择水平距离(或斜距离)和高度两个自变量的方式分区拟合，分段拟合其实是分区拟合的一个特例。

1.1 算法选择

算法设计的原则是在保证解析函数精度的情况下，尽量减少分区，使解析函数用单组拟合系数尽可能地逼近气象修正表定义域中更大区域上的修正量数据，符合残差要求的最大区域作为一个分区。在搜索最大区域时，采用自顶向下的设计，只需在每个决策点做出在当前问题中的最优选择，直接加入局部最优解集合，而不必考虑子问题的解。如果有剩余子问题，再将子问题作为当前问题开始下一次的最大区域搜索。基于以上设计思想，适合选用贪心算法进行自动分区。贪心算法的显著特点是每次对未知区域搜索都采取当前状态下的局部最优解，通过选择局部最优产生出全部定义域上的全局最优解。

贪心算法通用步骤如下[9]：

1)从全部定义域上的离散数据出发，根据贪心策略先得出当前状态的最优解，缩小问题规模。

2)采用迭代方法，向剩余部分的定义域前进，根据贪心策略得到剩余子问题的局部最优解，逐步缩小问题规模。

3)将所有解综合，验证结果。

1.2 分区原则

贪心算法所得的局部解是否最优，可用拟合残差和拟合区域面积大小来判断。拟合结果的验收依据是残差极值和残差均方差，可以针对每种气象参数设置一个残差极限门限和均方差门限。当修正表当前选择区域的数据拟合残差低于残差门限和均方差门限时，可以扩大选择区域再次计算拟合系数和拟合残差，直到残差超过门限为止，满足精度要求的最大区域为气象修正表的一个分区，这就是贪心算法的贪心策略。

残差门限应合理设置，气象修正表是通过实测和理论计算综合处理得出，本身存在误差，因此残差门限设置过低不会带来更多收益。一般选择均方差不大于气象修正表数据精度的5倍，最大残差控制在不大于均方差的5倍[6]。

1.3 拟合方法

设Y=(y1,y2,y3,…,yp)T是修正表中的一组离散修正量数据，是自变量X=(X1,X2,X3,…,Xm;m≪p)的连续函数，自变量X中的元素均由DaHb(a≥0,b≥0,a+b≤n)形式的二元基函数组成的向量[6]，n为逼近函数的阶数，D为某个修正量数据所对应的距离或斜距离，视计算结果的分区数量或拟合残差的情况选择距离或斜距离，H为相同修正量数据所对应的高度。

根据维尔斯拉斯(Weierstrass)逼近定理，逼近函数的多项式项数越多，其逼近效果越好[10]，但实际计算中发现，各自变量多项式对于逼近函数的贡献并不相同，部分项对结果的贡献极其微弱，影响计算速度，而且项数过多时，虽然在逼近离散修正量时可降低残差，但在计算修正表未定义的离散点修正量时，有时会出现修正量明显正常范围的不合理波动[6]，这种不合理波动需要通过验算发现。为避免过多的高阶项，n一般不大于7，逼近函数多项式的项数不超过36。当选择4阶逼近函数时，m=15，函数表达式为

y=f(D,H)=β1+β2D+β3H+β4D2+β5DH+
β6H2+β7D3+β8D2H+β9DH2+β10H3+β11D4+
β12D3H+β13D2H2+β14DH3+β15H4+ε，

(1)

式中，ε为逼近残差。

B=(β1,β2,β3,…,β15)T是一组待求系数，写为矩阵形式为Y=XB+ε。

(2)

采用最小二乘法拟合数据。最小二乘法是使残差平方和Q(B)达到最小。根据极值定理B应满足：

(3)

可得非齐次线性方程组系数矩阵：

(4)

1.4 自动分区

分区拟合过程中，每次迭代开始前，要在气象修正表内选定一块尚未拟合的连续区域M，为便于规定该区域的定义域范围，连续区域最好是规则形状，本文选择矩形区域。再根据区域的距离和高度范围，筛选区域M内的修正量数据，作为本次迭代的Y向量，然后对筛选的数据进行拟合并计算残差。如果在区域M内，存在宽度和长度足够大的连续子区域Pi(Pi⊆M)，使Pi内的全部数据的拟合精度满足要求，则将精度满足的区域标记为已拟合，缩小气象修正表待拟合区域的范围，选取面积最大的区域作为一个分区。然后继续根据贪心策略遍历所有剩余连续矩形区域，将贪心搜索所得的下一个精度合格的最大区域作为一个分区，进一步缩小问题范围。

单次迭代中，为避免出现碎片区域，每个被搜索子区域的距离和高度的上限或下限均与区域M的边界重合，即从区域M的距离边界和高度边界开始搜索，在每一整块子区域Pi内遍历计算每个高度也连续的子区域Qi(Qi⊆P)的拟合残差，如果被搜索的区域Qi面积小于本次迭代中已搜索到的最大区域，则放弃计算，选择下一个可处理的子区域Qi+1或Pi+1，直到遍历整个区域M，精度满足要求的最大区域作为一个分区。

贪心搜索最大区域时，每次移动距离或高度游标，可同时计算两个游标到各自起始点共4个区域的拟合残差，移动游标，实现遍历搜索。

贪心搜索程序设计如下：

GetDataArea//区域M

for i←startD to stopD//子区域Pi

for j←startH to stopH//子区域Qi

//距离上的前进搜索高度上的前进搜索

GetDataAndFit(startD, i , startH , j)

//距离上的前进搜索高度上的后退搜索

GetDataAndFit( startD , i , j , stopH )

//距离上的后退搜索高度上的前进搜索

GetDataAndFit( i , stopD , startH , j )

//距离上的后退搜索高度上的后退搜索

GetDataAndFit( i , stopD , j , stopH )

end

end

return maxArea

1.5 算法流程

在上述关键处理环节的基础上，设计如图1所示的完整算法流程。

计算开始后进行参数初始化。首先加载气象修正表数据，将其存储于可以读取任意距离和高度范围内所有修正量数据的二维数组中；建立一个映射表，该表记录气象修正表定义域上的已经拟合过的数据，下次迭代前会通过该表查询未拟合的区域；选择当前要拟合的参数类型，例如气温变化引起的距离修正量；选择基函数形式和函数阶数；根据修正表数据精度设置贪心策略。然后开始迭代，自动搜索低于残差门限的最大区域进行分区。如果某次迭代的区域长度或宽度不足，则无法保证自变量多项式的线性无关性，可能导致逼近函数不连续、不平滑的错误结果，对于这种情况需增加区域长度或宽度，必要时可将已拟合的区域包含进来。如果某次迭代的数据定义域与上次相同，表明上次没有搜索到满足精度要求的区域，需退出程序，对该区域人工拟合后检查残差分布。

2 算法验证

在分区结束后，修正表所有的修正量都可通过逼近函数计算得出，可验算全部分区的拟合函数系数。由于有贪心策略控制，气象修正表中的所有值域数据的拟合残差都已低于残差门限，验算时遍历全部值域。通过值域数据对应的水平距离和高度查找其所在分区的拟合系数，用逼近函数结合该组系数计算所得的修正量计算出拟合残差。除残差外，还需检查逼近函数的连续性[13]。调用逼近函数，通过距离不变、高度步进和高度不变、距离步进的方法，计算修正表中未给出的离散点数据，检查步进过程中的修正量增量是否存在不合理波动。如果自变量超出定义域，则增量置为0。观察增量数据曲线，曲线应随距离或高度的增加平滑变化。虽然各区域均为闭区间，相邻分区的边界重合，但气象修正表本身存在误差，区域间结合处附近的函数值仍可能存在跳点，如出现米级以上的增量变化，必定是区域间结合处的增量不连续，则需滑动两区域的分界线，重新分区，直到新分区的拟合系数可实现增量连续。

图2、3为残差分布图，一般情况下，受试验现场气象测量精度、弹道测量精度和随机因素的影响，在气象修正表定义域最大距离或最大高度附近的数值随机误差较大，残差也较大，可适当提高门限，也可能存在拟合残差总是超过门限的点或极个别残差明显偏大的“贼点”数据，可将“贼点”附近的数据做平滑处理，用插值数据替代“贼点”数据[6]。

图4为5条在高度不变条件下的气温距离修正量随目标水平距离步进的增量曲线，曲线表示目标分别为高度H=500,600,700,800,900 m时，逼近函数计算出的目标在水平距离为sd时，水平距离每增加0.1 m时的距离修正量的增量，即Δsd=f(sd+0.1,H)-f(sd,H)，sd为修正表水平距离定义域内的任一数值。

图5为5条在距离不变条件下的气温距离修正量随目标高度步进的增量曲线，曲线表示目标水平距离分别为SD=0.5，1.0，1.5，2.0，2.5 km时，逼近函数计算出的目标在高度为h时，高度每增加0.1 m时的距离修正量的增量，即Δsd=f(SD,h+0.1)-f(SD,h)，h为修正表高度定义域内的任一数值。

图6为5条在高度不变条件下的气温高度修正量随水平距离步进的增量曲线，图7为5条在距离不变条件下的气温高度修正量随高度步进的部分增量曲线。图6和图7中的每条曲线代表随水平距离或高度按0.1 m步进的修正量增量走势。

由图4～7可知，每条曲线连续平滑，没有拐点、不合理波动或“毛刺”。在定义域内如若所有增量曲线均连续平滑，则表明拟合系数正确可用，分区方式合理，逼近函数可用于非标条件的修正。

3 结束语

基于贪心算法的自动分区拟合算法可以根据残差门限自动搜索满足精度要求的区域自动进行分区，通过对离散数据的拟合计算出逼近函数系数，从而实现对气象修正量的快速计算。验算结果表明，在修正表的定义域内，逼近函数可以准确逼近修正表已给出的离散修正量数据，对于修正表未给出的任意距离和高度上的修正量可正确计算。