CN109507495B - 一种变窗长的准同期并网参数测量方法 - Google Patents

Abstract

本发明提供了一种变窗长的准同期并网参数测量方法，通过汉宁卷积窗抑制频谱泄露；通过变窗长相位差法测量并网点两侧电压信号基波频率，根据所测基波频率调整卷积窗长度减小整周期采样的时间误差；采用离线生成不同窗长卷积窗、离线计算卷积窗与正余弦因子乘积并建立加权系数表的方法，加快算法的实时性；推导了方便快捷的汉宁卷积窗幅值校正公式，易于实现，计算量小。

Description

一种变窗长的准同期并网参数测量方法

技术领域

本发明涉及并网参数测量领域，尤其涉及一种变窗长的准同期并网参数测量方法。

背景技术

电网运行过程中常需要把发电机或微电网与大电网进行并列，这种将小系统通过断路器等开关设备并入大系统的操作称为同期操作。所谓同期即开关设备两侧电压大小相等、频率相等、相位相同，同期装置的作用是用来判断路器两侧是否达到同期条件，从而决定能否执行合闸并网的专用装置，其关键在于对并网双方电压信号频率、幅值以及相位的快速准确测量。为便于描述，开关设备的两侧分别称为s侧和g侧，s侧电压记为u_s，g侧电压记为u_g，。目前，同期参数(并网双方电压信号频率、幅值以及相位)测量方法主要分为硬件法和软件法。

硬件法是利用硬件电路实现对并网两侧电压信号的频率差、幅值差及相位差的测量，其运算简单，但需要专门的频率、相位检测电路，且存在抗干扰能力弱、易受高次谐波影响等缺点。

软件法主要是对并网两侧电压进行交流采样并通过数值方法计算电压参数，该方法不需专门的频率和相位测量电路，其测量精度及同期性能主要取决于频率、相位及幅值的微机算法。常用的同期并网参数测量算法是傅立叶变换算法，但在非同步采样时受频谱泄露的影响会出现较大的误差，尤其当微电网电压谐波含量较大时影响更为突出。对信号进行时域加窗处理，可以抑制频谱泄露。常用窗函数主要有矩形窗、三角窗、hanning窗等组合余弦窗，以及矩形自卷积窗、Hanning自卷积窗、Rife-Vincent自卷积窗等一系列自卷积窗。但并网装置测量的2个电压信号(u_s与u_g)频率不等，同时实现2个信号的整周期采样十分困难，加窗虽在一定程度上抑制频谱泄露，但当偏离整周期采样较多时，频谱泄露仍然比较大。中国发明专利CN201710647804.4“基于自适应序列全相位DFT的同期并网参数测量方法”通过计算两段采样序列的相关系数动态改变窗长以接近整周期采样，虽然其在一定程度上降低了频谱泄漏，提高了采样精度和参数的测量精度，但是其整个方法的计算量大，且相关系数的阀值选择依赖主观经验，测量精度并不稳定；而且其卷积窗采用在线生成的方式，进一步的加大了算法的计算量。

发明内容

本发明的目的在于提供一种变窗长的准同期并网参数测量方法，能够解决现有的算法偏离整周期采样较多时，频谱泄露仍然比较大的问题，而且整个计算量大大降低。

为了实现上述目的，本发明采用以下技术方案：

一种变窗长的准同期并网参数测量方法，包括以下步骤：

步骤A：根据频率的波动范围，确定L个周期的采样点数范围，对范围内每个偶数N，离线生成对应N点的卷积窗，再离线生成正余弦因子

然后将对应N点的卷积窗和正余弦因子相乘直接得到不同点数卷积窗对应的一维加权系数表；N为对应点数，j为虚数的单位符号，n为从0至N-1的变量，这里k＝L；

步骤B：启动同期检测装置，以采样频率f_s分别同时对s侧电压u_s与g侧电压u_g进行采样Lf_s/f_e+1个点，f_e表示电网额定频率50Hz；然后对u_s的Lf_s/f_e+1个采样点，去掉最后一个采样点剩下Lf_s/f_e个点进行一次单谱线DFT，结果记为X_s00，再去掉第一个采样点剩下Lf_s/f_e个点进行一次单谱线DFT，结果记为X_s01，显然前后两次DFT分析的数据段长度相同且相差一个采样间隔，利用常规相位差测频法初步计算u_s的频率f_us0，根据f_us0初始化变量N_s0＝N_s1＝even(Lf_s/f_us0)，其中even(x)为取最接近x的偶数，L为采样周期数，f_s为采样频率；同样对u_g的Lf_s/f_e+1采样点，进行同样的处理初步计算出u_g的频率f_ug0，根据f_ug0初始化变量N_g0＝N_g1＝even(Lf_s/f_ug0)，其中even(x)为取最接近x的偶数；

步骤C：继续以采样频率f_s分别同时对s侧电压u_s与g侧电压u_g进行采样，并将采集到的数据分别放入对应的缓冲区sb和gb中，利用变量crt，记录当前最新采样点存储的位置，变量num记录采样的数量，当采样的数量为整数Y时，然后进入下一步，其中缓冲区sb和gb的存储长度记为X，则Y<X；所述的Y＝even(Lf_s/f_0min)，其中even(x)为取最接近x的偶数，L为采样周期数，f_s为采样频率，f_0min为并网断路器两侧电压的最小可能频率；

步骤D：当采样数据个数达到设定值Y个后，从缓冲区sb、gb中分别截取一段数据u_s0和u_g0进行单线谱DFT分析，得到分析结果；其中数据u_s0的长度为N_s0，u_g0的长度为N_g0；对u_s0数据段进行单线谱DFT分析结果记为X_s0L，对u_g0数据段进行单线谱DFT分析结果记为X_g0L；

步骤E：再分别采样P个点并存入对应的缓冲区sb和gb中，则最新采样点在缓冲区中的位置自然就变成了i+P，再从sb和gb中分别截取一段数据u_s1和u_g1进行单谱线DFT分析，得到分析结果，其中，数据u_s1的长度为N_s1，数据u_g1的长度为N_g1；对u_s1数据段进行单线谱DFT分析结果记为X_s1L，对u_g1数据段进行单线谱DFT分析结果记为X_g1L；

步骤F：分别利用X_s0L、X_s1L和X_g0L、X_g1L根据变窗长的相位差法分别计算电压u_s和电压u_g的频率f_us和f_ug、幅值A_s和A_g、相位

和

并判断电压u_s、u_g是否满足并网条件，满足则发出合闸指令，结束测量；不满足，则进入下一步；

步骤G：令X_s0L＝X_s1L、N_s0＝N_s1,X_g0L＝X_g1L、N_g0＝N_g1，为下次测频做准备，同时利用所测频率f_us和f_ug、计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1,即更新每周期采样点数N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数，用于让下一次DFT分析截取的数据段接近整数个周期；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；

同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2；

步骤H：由上步骤得到的需要截取的数据段的长度N_s1和N_g1，重复步骤E-G，直至满足发出并网指令，结束测量。

所述的步骤C中，当数据存满缓冲区sb或缓冲区gb，新采样的数据会覆盖最老的采样数据；其中，当某时刻采样值存入sb(Ln-1)，缓冲区会被填满之后，新的采样值uy被存入sb(0)，覆盖了sb(0)中老的采样数据，这样缓冲区内的数据总是最新的采样数据；针对采样数据的检索访问方式为sb(j_d)＝sb(mod(j_d,Ln))，mod(x，y)表示求x/y的余数；缓冲区gb也采用同样的存储和访问方式。

所述的从缓冲区sb中截取长度为N的一段数据的方法具体为：以当前最新采样点为末尾，向前数N个点；假设某时刻最新采样点存储的位置为i，此位置的数据为sb(i),则单线谱DFT分析应该截取的数据段为从sb(i)开始向前数N个数据，即应该是由sb(i-N+1)至sb(i)存储单元中的采样数据构成，因为这两个点之间的数据个数为i-(i-N+1)+1＝N个，N即为数据的长度；所述的从缓冲区gb中截取长度为N的一段数据的方法与从缓冲区sb中截取长度为N的一段数据的方法相同，这里N表示N_s0、N_s1或N_g0、N_g1中的任何一个参数。

所述的单谱线DFT的具体过程为：分别对缓冲区sb和gb中最近采样的N个点数据与N个点对应的一维加权系数表对应相乘并累加，得到第L条谱线的DFT结果。

所述步骤A中N点卷积窗对应的一维加权系数表具体包括如下步骤：

步骤S1：定义N1＝N/2，步骤S1.1：离线生成N1点Hanning窗；

步骤S1.2：对N1点Hanning窗进行自卷积得到2N1-1点卷积窗，对2N1-1点卷积窗末端补1个零得到2N1点卷积窗，即N点卷积窗；

步骤S2：由N点卷积窗离线生成N点正余弦因子；

步骤S3：将N点卷积窗与N点正余弦因子相乘即可得到N点卷积窗对应的N点一维加权系数表；

步骤S4：建立数组addr，存储步骤S3得到的N点卷积窗对应的N点一维加权系数表的首地址；以频率上限对应的采样点数N_min为基准，N_min点一维加权系数表的首地址存在addr第一个元素位置，其他采样点数的一维加权系数表首地址按点数递增顺序存放，则N点一维加权系数表的首地址为addr((N-N_min)/2)。

所述利用变窗长相位差法计算电压频率f_us的方法为：

记数据段u_s1为最近一次进行DFT分析的数据段，N_s1表示u_s1数据段的长度，u_s0为上一次进行DFT分析的数据段，N_s0表示u_s0数据段的长度；u_s1数据段是u_s0数据段又采样P个点后截取的，因此s侧两次DFT分析的采样序列起始点对应的间隔为N_s0-N_s1+P，因此u_s1和u_s0的采样起始点的相位差值满足如下关系：

φ_s1-φ_s0＝2πf_us×(N_s0-N_s1+P)T_s； (19)

其中，φ_s1为u_s1的采样起始点的相位，φ_s0为u_s0的采样起始点的相位，T_s为采样周期，f_s为采样频率，T_s＝1/f_s；f_us为s侧的电压信号频率；根据傅立叶变换，采样起始点的相位θ又和谱线X(k)有如下关系：

其中，arg为求复数的相角，k₀＝Nf₀/f_s，N为DFT分析数据段的长度，f₀为信号频率；W(2πx/N)为汉宁卷积窗的频谱；

为减小计算量，本发明采用单谱线DFT，即只计算第L条谱线的值X_s0L、X_s1L，X_g0L、X_g1L，结合汉宁卷积窗的频谱公式有：

因此：

因为相角是以2π为周期的，因此根据式22求出的相角差和式19等式右侧计算出的相角差可能有2π整数倍的差距，因此可得：

其中，q为一整数，进一步的：

2πf_us×(N_s0-N_s1+P)T_s-πf_usT_s(N_s0-N_s1)＝arg(X_s1L)-arg(X_s0L)±2πq；

进一步的：

πf_us(N_s0-N_s1+2P)T_s＝arg(X_s1L)-arg(X_s0L)±2πq；

令f_us＝50+df，df∈[-5Hz,5Hz]，反映了频率偏离50Hz的程度，则上式进一步推导如下：

则：

进而s侧的电压信号频率f_us＝50+df：

同理：

g侧的电压频率f_ug的方法为：

f_ug＝50+df。

所述利用变窗长相位差法计算s侧的电压相位

和g侧的电压相位

通过如下公式计算：

计算信号幅值A_s、A_g的公式分别为：

其中，

Δf＝f_s/N_s1,

其中，

Δf＝f_s/N_g1,

所述步骤F中利用所测频率f_us和f_ug计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1，计算公式为N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2。

所述的

的计算方法为：

令x在(-0.5,0.5)之间以0.0001的速度递增变化，对应的计算F(x)，然后利用matlab多项式逼近函数polyfit(x,F(x),n)，n为多项式最高幂次，得到F(x)的近似多项式：

F(x)＝8.00000000+1.28987025x²+0.12451951x⁴+0.00994809x⁶。

本发明的有益效果：

本发明所述的一种变窗长的准同期并网参数测量方法，通过汉宁卷积窗抑制频谱泄露；通过变窗长相位差法测量并网点两侧电压信号基波频率，根据所测基波频率调整卷积窗长度减小整周期采样的时间误差；采用离线生成不同窗长卷积窗、离线计算卷积窗与正余弦因子乘积并建立加权系数表的方法，加快算法的实时性；推导了方便快捷的汉宁卷积窗幅值校正公式，易于实现，计算量小。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的方法流程图；

图2为本发明所述卷积窗加权系数表生成原理图；

图3为本发明所述变窗长相位差测频起始点间隔示意图。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，一种变窗长的准同期并网参数测量方法，包括以下步骤：

步骤A：根据频率的波动范围，确定L个周期的采样点数范围，对范围内每个偶数N，离线生成对应N点的卷积窗，再离线生成正余弦因子

然后将对应N点的卷积窗和正余弦因子相乘直接得到不同点数卷积窗对应的一维加权系数表；N为对应点数，j为虚数的单位符号，n为从0到N-1的变量,这里K＝L。

所述步骤A中N点卷积窗对应的一维加权系数表具体包括如下步骤：

步骤S1：定义N1＝N/2，步骤S1.1：离线生成N1点Hanning窗；

步骤S1.2：对N1点Hanning窗进行自卷积得到2N1-1点卷积窗，对2N1-1点卷积窗末端补1个零得到2N1点卷积窗，即N点卷积窗；

步骤S2：由N点卷积窗离线自动生成N点正余弦因子；

步骤S3：将N点卷积窗与N点正余弦因子相乘即可得到N点卷积窗对应的N点一维加权系数表；

步骤S4：建立数组addr，存储步骤S3得到的N点卷积窗对应的N点一维加权系数表的首地址；以频率上限对应的采样点数N_min为基准，N_min点一维加权系数表的首地址存在addr第一个元素位置，其他采样点数的一维加权系数表首地址按点数递增顺序存放，则N点一维加权系数表的首地址为addr((N-N_min)/2)。

步骤B：启动同期检测装置，以采样频率f_s分别同时对s侧电压u_s与g侧电压u_g进行采样Lf_s/f_e+1个点，f_e表示电网额定频率50Hz；然后对u_s的Lf_s/f_e+1个采样点，去掉最后一个采样点剩下Lf_s/f_e个点进行一次单谱线DFT，结果记为X_s00，再去掉第一个采样点剩下Lf_s/f_e个点进行一次单谱线DFT，结果记为X_s01，显然前后两次DFT分析的数据段长度相同且相差一个采样间隔，可利用常规相位差测频法初步计算u_s的频率f_us0；根据f_us0初始化变量N_s0＝N_s1＝even(Lf_s/f_us0)，其中even(x)为取最接近x的偶数，L为采样周期数，f_s为采样频率；同样对u_g的Lf_s/f_e+1采样点，进行同样的处理初步计算出u_g的频率f_ug0，根据f_ug0初始化变量N_g0＝N_g1＝even(Lf_s/f_ug0)，其中even(x)为取最接近x的偶数。常规相位差测频法是现有技术，不再赘述。

步骤C：继续以采样频率f_s分别同时对电压u_s与u_s进行采样，并把采集到的数据分别放入对应的缓冲区sb和gb中，利用变量crt，记录当前最新采样点存储的位置，变量num记录采样的数量，当采样的数量为整数Y时，然后进入下一步，其中缓冲区sb和gb的存储长度记为X，则Y<X；所述的Y＝even(Lf_s/f_0min)，其中even(x)为取最接近x的偶数，其中L为采样周期数，f_s为采样频率，f_0min为并网断路器两侧电压的最小可能频率；

开始时，缓冲区中的采样数据为0个，而DFT分析需要一定量的数据才能进行，因此步骤C执行时需要边采样边等待，等到采样数据个数达到一定数值后才开始进行DFT分析，这个数值被初始化为最接近Lf_s/f_0min的偶数，L为采样周期数，f_s为采样频率，f_0min为并网断路器两侧电压的最小可能频率。以采样频率f_s对电压u_s、u_g进行同时采样，每采样一个点，就存入对应的缓冲区sb或gb中，并用一个变量crt，记录当前最新采样点存储的位置；以存储u_s的采样值为例，sb(0)sb(1)…sb(Ln-1)代表采样数据缓冲区sb中的一个个存储单元，Ln是缓冲区长度，某时刻采样值u0存入sb(0)单元，则crt＝0,下一时刻采样值u1存入sb(1)单元，则crt＝1，依次类推，i时刻采样值ui存入sb(i),则crt＝i，crt始终存储当前时刻对应的最新采样值在缓冲区中的位置；

所述的步骤C中，当数据存满缓冲区sb或缓冲区gb，新采样的数据会覆盖最老的采样数据；即当某时刻采样值存入sb(Ln-1)，缓冲区会被填满，之后，新的采样值u_y被存入sb(0)，覆盖了sb(0)中老的采样数据，这样缓冲区内的数据总是最新的采样数据；针对采样数据的检索访问方式为sb(j_d)＝sb(mod(j_d,Ln))mod(x，y)表示求x/y的余数。缓冲区gb具有和sb类似的存储和访问方式，不再赘述。

步骤D：当采样数据个数达到设定值Y个后，从对应的缓冲区sb和gb中分别截取一段数据u_s0和u_g0进行单线谱DFT分析，得到分析结果；其中数据u_s0的长度为N_s0，u_g0的长度为N_g0；对u_s0数据段进行单线谱DFT分析结果记为X_s0L，对u_g0数据段进行单线谱DFT分析结果记为X_g0L；所述的单线谱DFT的具体过程为：分别对缓冲区sb最近采样的N_s0个点数据与N_s0个点对应的一维加权系数表对应相乘并累加，得到u_s0第L条谱线的DFT结果,同理对缓冲区gb最近采样的N_g0个点数据与N_g0个点对应的一维加权系数表对应相乘并累加，得到u_g0第L条谱线的DFT结果；

所述的从缓冲区sb中截取长度为N_s0的一段数据的方法具体为：以当前最新采样点为末尾，向前数N_s0个点；假设某时刻最新采样点存储的位置为i，此位置的数据为sb(i),则单线谱DFT分析应该截取的数据段为从sb(i)开始向前数N_s0个数据，即应该是由sb(i-N_s0+1)至sb(i)存储单元中的采样数据构成，因为这两个点之间的数据个数为i-(i-N_s0+1)+1＝N_s0个，N_s0即为数据的长度。同样的方式可以从缓冲区gb中截取长度为N_g0的一段数据。

步骤E：再分别采样P个点并存入对应的缓冲区sb和gb中，则最新采样点在缓冲区中的位置自然就变成了i+P，再从sb和gb中分别截取一段数据u_s1和u_g1进行单谱线DFT分析，其中u_s1数据段的长度为N_s1，u_g1数据段的长度为N_g1；对u_s1数据段的单谱线DFT分析结果记为X_s1L，对u_g1数据段的单谱线DFT分析结果记为X_g1L；

截取数据的方法：

以最新采样点为末尾，针对缓冲区sb，向前数N_s1个点，则本段数据在缓冲区中的位置为sb(i+P-N_s1+1)至sb(i+P)，它们之间的数据个数为(i+P)-(i+P-N_s1+1)+1＝N_s1，对u_s1数据段进行DFT分析结果记为X_s1L；

针对缓冲区gb，向前数N_g1个点，则本段数据在缓冲区中的位置为gb(i+P-N_g1+1)至gb(i+P)，它们之间的数据个数为(i+P)-(i+P-N_g1+1)+1＝N_g1，对u_g1数据段进行DFT分析结果记为X_g1L；

步骤F：分别利用X_s0L、X_s1L和X_g0L、X_g1L根据变窗长的相位差法分别计算电压u_s和电压u_g的频率f_us和f_ug、幅值A_s和A_g、末尾点相位

和

并判断电压u_s、u_g是否满足并网条件，满足则发出合闸指令，结束测量；不满足，则进入下一步；步骤G：令X_s0L＝X_s1L、N_s0＝N_s1,X_g0L＝X_g1L、N_g0＝N_g1，为下次测频做准备，同时利用所测频率f_us和f_ug、计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1,即更新每周期采样点数N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数，便于让下一次DFT分析截取的数据段接近整数个周期；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2。

步骤H：重复步骤E-G，直至满足发出并网指令，结束测量。

所述利用变窗长相位差法计算电压频率f_us的方法为：

记u_s1为最近一次进行DFT分析的数据段，u_s0为上一次进行DFT分析的数据段，u_s1数据段是u_s0数据段又采样P个点后截取的，因此s侧两次DFT分析的采样序列起始点对应的间隔为N_s0-N_s1+P，因此u_s0和u_s1的起始点的相位差值满足如下关系：

φ_s1-φ_s0＝2πf_us×(N_s0-N_s1+P)T_s； (19)

其中，φ_s1为u_s1的起始点的相位，φ_s0为u_s0的起始点的相位，T_s为采样周期，f_s为采样频率，T_s＝1/f_s；f_us为s侧的电压信号频率；

根据傅立叶变换，采样起始点的相位θ又和谱线X(k)有如下关系：

其中，arg为求复数的相角，k₀＝Nf₀/f_s，N为DFT分析数据段的长度，f₀为信号频率。

W(2πx/N)为汉宁卷积窗的频谱，

为减小计算量，本发明采用单谱线DFT，即只计算第L条谱线的值X_s0L、X_s1L，X_g0L、X_g1L，结合汉宁卷积窗的频谱公式有：

因此：

因为相角是以2π为周期的，因此根据式22求出的相角差和式19等式右侧计算出的相角差可能有2π整数倍的差距，因此可得：

其中，q为一整数，进一步的：

2πf_us×(N_s0-N_s1+P)T_s-πf_usT_s(N_s0-N_s1)＝arg(X_s1L)-arg(X_s0L)±2πq；

进一步的：

πf_us(N_s0-N_s1+2P)T_s＝arg(X_s1L)-arg(X_s0L)±2πq；

令f_us＝50+df，df∈[-5Hz,5Hz]，反映了频率偏离50Hz的程度，则上式进一步推导如下：

则：

进而s侧的电压信号频率f_us＝50+df：

上式中，分母一般不等于零，但如果窗长在变动过程中偶然出现分母为零的情况则会无法测频，因此在步骤G确定窗长时应判断N_s0-N_s1+2P是否等于0，如果等于0则令N_s1＝N_s1-2，以避免在变窗长相位差测频时出现分母为0的情况。

同理：

g侧的电压频率f_ug的方法为：

f_ug＝50+df。

上式中，分母一般不等于零，但如果窗长在变动过程中偶然出现分母为零的情况则会无法测频，因此在步骤G确定窗长时应判断N_g0-N_g1+2P是否等于0，如果等于0则令N_g1＝N_g1-2，以避免在变窗长相位差测频时出现分母为0的情况。

所述利用变窗长相位差法计算s侧的电压末尾采样点相位

和g侧的电压末尾采样点相位

通过如下公式计算：

计算信号幅值A_s、A_g的公式分别为：

其中，

Δf＝f_s/N_s1,

其中，

Δf＝f_s/N_g1,

所述的

的计算方法为：

令x在(-0.5,0.5)之间以0.0001的速度递增变化，对应的计算F(x)，然后利用matlab多项式逼近函数polyfit(x,F(x),n)，n为多项式最高幂次，得到F(x)的近似多项式：

F(x)＝8.00000000+1.28987025x²+0.12451951x⁴+0.00994809x⁶。

所述步骤G中利用所测频率f_us和f_ug、计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1,计算公式为

N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数；；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2。

为了方便本领域技术人员更清楚的了解本发明的内容，下面将进一步介绍单一频率信号x(t)加卷积窗后的频谱分析方法：

首先，定义频率为f₀、幅值为A、初相位为θ的单一频率信号x(t)，经过采样率为f_s的模数变换后得到离散信号x(n)：

其次，定义加窗函数的时域形式为w(n)，其连续频谱为W(2πf)，则加窗后x(n)的连续傅立叶变换为：

若忽略x(n)负频点-f₀处频峰的旁瓣影响，则x(n)在正频点f₀附近的连续频谱函数X(f)可以表达为：

对式(1)进行离散抽样，即可得到X(f)的离散傅立叶变换的表达式X(kΔf)为：

其中，Δf为离散频率间隔，Δf＝f_s/N，N为数据截断长度，k为谱线坐标，k＝0,1,2,…,N-1；进一步：

其中，k₀为峰值频率对应的谱线坐标，k₀＝f₀/Δf＝Nf₀/f_s，X(kΔf)简记为X(k)；由于x(n)的峰值频率f₀很难正好位于离散谱线的频点上，因此，k₀一般不是整数；假设距离峰值频率最近的谱线为第k_i条，则有：

进一步的：信号幅值A：

信号初相位θ：

其中，| |表示求复数的模，arg表示求复数的相位；

为了较好地抑制非同步采样下的频谱泄露，采用二阶汉宁卷积窗对信号进行时域加窗处理，N点汉宁卷积窗由两个N/2点汉宁窗卷积得到。N/2点的汉宁窗的计算式为：

因为M点汉宁窗频谱W_hann(ω)可表示为：

其中：

因此，N/2点汉宁窗的频谱可表示为：

其中：

根据卷积定理，函数在时域卷积等效于频域相乘，因此汉宁卷积窗的频率响应为原汉宁窗频谱的平方；则因此汉宁卷积窗的频谱W(ω)可表示为：

为方便进一步推导，令ω＝2πx/N对W(ω)进行变量代换：

为简化频谱表达式，分别对式中

进行化简：

将式12、13带入式11，可得：

因为一般N很大，则x<<N时，上式可简化为：

结合式4、5，信号的幅值A可按下式计算：

其中，

本发明所述的一种变窗长的准同期并网参数测量方法，通过汉宁卷积窗抑制频谱泄露；通过变窗长相位差法测量并网点两侧电压信号基波频率，根据所测基波频率调整卷积窗长度减小整周期采样的时间误差；采用离线生成不同窗长卷积窗、离线计算卷积窗与正余弦因子乘积并建立加权系数表的方法，加快算法的实时性；推导了方便快捷的汉宁卷积窗幅值校正公式，易于实现，计算量小。

如图2所示：优选方案为：所述通过数组addr检索到N点一维加权系数表的首地址的方法为：

步骤S1：定义N1＝N/2，

步骤S1.1：离线生成N1点Hanning窗；

步骤S1.2：对N1点Hanning窗进行自卷积得到2N1-1点卷积窗，对2N1-1点卷积窗末端补1个零得到2N1点卷积窗，即N点卷积窗；

其中，步骤S1.1、步骤S1.2的方法在背景技术中已做详细描述，这里不再赘述；

步骤S2：由N点卷积窗离线自动生成N点正余弦因子；

步骤S3：将N点卷积窗与N点正余弦因子相乘即可得到N点卷积窗对应的N点一维加权系数表；

步骤S4：建立数组addr，存储步骤S3得到的N点卷积窗对应的N点一维加权系数表的首地址；以频率上限对应的采样点数N_min为基准，N_min点一维加权系数表的首地址存在addr第一个元素位置，其他采样点数的一维加权系数表首地址按点数递增顺序存放，则N点一维加权系数表的首地址为addr((N-N_min)/2)。

需要注意的是：加窗插值算法虽然在一定程度上抑制了频谱泄露，但并网装置测量的两个电压信号频率不等，同时实现两个信号的整周期采样十分困难，当偏离整周期采样较多时，频谱泄露仍然比较大；为了尽量降低频谱泄露，本发明采用变窗长方法，根据信号的频率调整采样数据的时间窗长度使其尽量接近整数个周期，以降低频谱泄露；但随着采样时间窗长度的变化，对应的卷积窗长度也必须随之变化，由于卷积窗是由余弦窗自卷积得到的，如果采用在线生成卷积窗，然后再按下式进行DFT：

其中，k＝0～N-1，w(n)为汉宁卷积窗；

上述过程非常耗时，难以满足准同期实时性要求，同时由于采样时间窗长度可变，不能满足采样点数恒为2的整次幂，因此快速傅立叶变换也无法使用。

所以本发明采用离线生成汉宁卷积窗的方法，但随着卷积窗长度的变换，DFT变换所需的正余弦因子也会变化，在线生成正余弦因子需要大量三角函数运算，同样费时；因此本发明在离线生成卷积窗后，再离线生成正余弦因子

然后将对应点数的卷积窗和正余弦因子相乘直接得到不同点数卷积窗对应的加权系数表

式中k为DFT的谱线坐标，k＝0,1,2…,N-1；

同时，考虑到准同期并网只需考虑基波信号，因此只需计算距离基波频率最近的谱线即可，可根据基波频率、采样频率和采样点数确定距离基波频率最近的谱线坐标，即k₀＝Nf₀/f_s，由于本发明采用变窗长方法，使采样窗长尽可能接近L个基波周期(即

因此k₀＝Nf₀/f_s≈L，即第L根谱线是距离基波频率最近的谱线，只需计算第L根谱线的DFT结果即可,相应地只需计算第L根谱线的DFT对应的系数表即可，即xsb_N_k中只需计算k＝L时对应的系数表；实施例L＝4，这样一来就只需建立每个N值下第4根谱线对应的一维加权系数表，所需存储空间小。

上述方法在加卷积窗的同时也完成了DFT变换所需的乘法运算，进行DFT时直接将采样数据和对应的加权系数表相乘再累加即可，大大节省了计算量。

以下以具体例子进行说明：因为同期装置投入运行时，并网点两侧的电压频率波动一般处于45-55Hz范围内，45Hz为频率下限，55Hz为频率上限，以采样频率f_s＝4000Hz为例，则4个周期的采样点数可能在：

4×f_s/f＝4×4000/(45～55)＝290.9～355.5之间；

对应的生成290,292,294,296,…,356点汉宁自卷积窗，以及290,292,294,296,…,356点正余弦因子，其中290为频率上限(55Hz)对应的采样点数，356为频率下限(45Hz)对应的采样点数。

其中N(N＝290,292,294,296,…,356)点汉宁卷积窗由N/2点汉宁窗自卷积后末尾补1个零得到；为简化计算及节省存储空间，不采用奇数点的汉宁卷积窗。

将不同点数的汉宁卷积窗与对应点数的正余弦因子相乘得到不同点数的汉宁卷积窗对应的一维加权系数表。

变窗长DFT能够根据测得的信号频率调整采样数据的时间窗长度，以使其尽量接近整数个周期；进一步的，可利用相位差法测量电网频率，方便实用且精度高，但由于窗长是可变的，因此固定窗长的相位差测频不能适用，本发明提出了一种变窗长的相位差测频法，与固定窗长的相位差测频方法不同，变窗长的相位差测频法窗长不再固定不变如图3所示：图中u_s1为本次s侧的采样序列，u_s0为s测上一次的采样序列，u_s1是u_s0又采样P个点后截取的；其中，P值越大，意味着进行下一次测频需要间隔的时间越长，实时性就越差，但P值太小，测频间隔的时间很小，处理器要足够快，否则难以在测频间隔时间内完成测频及其它功能所需的计算量，综合考虑，本实施例中P＝8；当窗长固定不变时，N_s1＝N_s0＝定值，此时前后两段数据的起始数据总是相差8个间隔，当窗长不固定时，N_s、N_s0是变化的且二者并不一定相等，前后两段数据的起始数据相差的间隔也不再固定不变，而是随着N_s1、N_s0的变化而变化，对照图3可以推算出两次采样序列对应的间隔为N_s0-N_s1+P；因此前后两段数据的起始点相位差值满足如下关系：

φ_s1-φ_s0＝2πf_us×(N_s0-N_s1+P)T_s； (19)

其中，φ_s1为u_s1的起始点的相位，φ_s0为u_s0的起始点的相位，T_s为采样周期，f_s为采样频率，T_s＝1/f_s；f_us为s侧的电压信号频率；

根据傅立叶变换，采样起始点的相位θ又和谱线X(k)有如下关系：

其中，arg为求复数的相角，k₀＝Nf₀/f_s，N为DFT分析数据段的长度，f₀为信号频率。

W(2πx/N)为汉宁卷积窗的频谱。

为减小计算量，本发明采用单谱线DFT，即只计算第L条谱线的值X_s0L、X_s1L，X_g0L、X_g1L，结合汉宁卷积窗的频谱公式有：

因此：

因为相角是以2π为周期的，因此根据式22求出的相角差和式19等式右侧计算出的相角差可能有2π整数倍的差距，因此可得：

其中，q为一整数，进一步的：

2πf_us×(N_s0-N_s1+P)T_s-πf_usT_s(N_s0-N_s1)＝arg(X_s1L)-arg(X_s0L)±2πq；

进一步的：

πf_us(N_s0-N_s1+2P)T_s＝arg(X_s1L)-arg(X_s0L)±2πq； (30)

令f_us＝50+df，df∈[-5Hz,5Hz]，反映了频率偏离50Hz的程度，则式30进一步推导如下：

考虑条件df∈[-5Hz,5Hz]，N_s0,N_s1∈[290,356]，P＝8，T_s＝1/4000，可得方程(32)左侧取值范围在[-π，π]，因此通过合适的q值调整方程(32)右侧的值也在[-π，π]范围内，则：

f_us＝50+df；

因此，可得：

步骤A：根据频率的波动范围，确定4个周期的采样点数范围，对范围内每个偶数N，离线生成对应N点的卷积窗，再离线生成正余弦因子

然后将对应N点的卷积窗和正余弦因子相乘直接得到不同点数卷积窗对应的一维加权系数表；N为对应点数，j为虚数的单位符号，n为从0至N-1的变量，这里k＝4；

步骤B：启动同期检测装置，以采样频率f_s分别同时对s侧电压u_s与g侧电压u_g进行采样320+1点；然后对u_s的321采样点，去掉最后一个采样点剩下320点进行一次单谱线DFT，结果记为X_s00，再去掉第一个采样点剩下320点进行一次单谱线DFT，结果记为X_s01，利用常规相位差测频法初步计算u_s的频率f_us0；根据f_us0初始化变量N_s0＝N_s1＝even(4f_s/f_us0)，其中even(x)为取最接近x的偶数，4为采样周期数，f_s为采样频率；同样对u_g的321个采样点，进行同样的处理初步计算出u_g的频率f_ug0，根据f_ug0初始化变量N_g0＝N_g1＝even(4f_s/f_ug0)，其中even(x)为取最接近x的偶数。

步骤C：继续利用采集装置获取采样数据开辟两个400点存储区sb、gb分别存储s侧电压采样数据、g电压采样数据，将采样数据按循环存储的方法存入sb、gb；

当采样数据达到356点时进入下一步；

步骤D：通过数组addr检索到N_s0点一维加权系数表的首地址，对sb最近采样的N_s0点数据u_s0与N_s0点一维加权系数表对应相乘并累加得到第4条谱线的DFT结果记为X_s0L，同理对gb最近采样的N_g0点数据u_g0与N_g0点一维加权系数表对应相乘并累加得到第4条谱线的DFT结果X_g0L；

步骤E：再分别采样8个点并存入从对应的缓冲区sb和gb中，通过数组addr检索到N_s1点一维加权系数表的首地址，对sb存储区最新N_s1点采样数据与N_s1点一维加权系数表对应相乘并累加得到第4条谱线的DFT结果,记为X_s1L；

通过数组addr检索到N_g1点一维加权系数表的首地址，对gb存储区最新N_g1点采样数据与N_g1点一维加权系数表对应相乘并累加得到第4条谱线的DFT结果,记为X_g1L；

步骤F：利用变窗长相位差法计算s侧电压信号频率f_us：

f_us＝50+df；

同理，利用变窗长相位差法计算g侧电压信号频率f_ug：

f_ug＝50+df；

进而：利用变窗长相位差法计算s侧的电压末尾采样点的相位

和g侧的电压末尾采样点

相位：

最后，利用变窗长相位差法计算s侧的电压幅值A_s和g侧的电压幅值A_g：

其中，

其中，

最终，计算s侧电压信号和g侧电压信号的频率差、幅值差和相位差等参数，判断是否满足并网条件，若满足则发出并网指令，不满足则进入下一步。

步骤G：令X_s0L＝X_s1L、N_s0＝N_s1,X_g0L＝X_g1L、N_g0＝N_g1，为下次测频做准备，同时利用所测频率f_us和f_ug计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1,即更新每周期采样点数N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数，便于让下一次DFT分析截取的数据段接近整数个周期；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2。

重复步骤E-G，直至满足并网条件发出并网指令，结束测量。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (10)

1.一种变窗长的准同期并网参数测量方法，其特征在于，包括以下步骤：

步骤A：根据频率的波动范围，确定L个周期的采样点数范围，对范围内每个偶数N，离线生成对应N点的卷积窗，再离线生成正余弦因子

然后将对应N点的卷积窗和正余弦因子相乘直接得到不同点数卷积窗对应的一维加权系数表；N为对应点数，j为虚数的单位符号，n为从0至N-1的变量，这里k＝L；

步骤B：启动同期检测装置，以采样频率f_s分别同时对s侧电压u_s与g侧电压u_g进行采样Lf_s/f_e+1个点，f_e表示电网额定频率50Hz；然后对u_s的Lf_s/f_e+1个采样点，去掉最后一个采样点剩下Lf_s/f_e个点进行一次单谱线DFT，结果记为X_s00，再去掉第一个采样点剩下Lf_s/f_e个点进行一次单谱线DFT，结果记为X_s01，显然前后两次DFT分析的数据段长度相同且相差一个采样间隔，利用常规相位差测频法初步计算u_s的频率f_us0，根据f_us0初始化变量N_s0＝N_s1＝even(Lf_s/f_us0)，其中even(x)为取最接近x的偶数，L为采样周期数，f_s为采样频率；同样对u_g的Lf_s/f_e+1采样点，进行同样的处理初步计算出u_g的频率f_ug0，根据f_ug0初始化变量N_g0＝N_g1＝even(Lf_s/f_ug0)，其中even(x)为取最接近x的偶数；

步骤C：继续以采样频率f_s分别同时对s侧电压u_s与g侧电压u_g进行采样，并将采集到的数据分别放入对应的缓冲区sb和gb中，利用变量crt，记录当前最新采样点存储的位置，变量num记录采样的数量，当采样的数量为整数Y时，然后进入下一步，其中缓冲区sb和gb的存储长度记为X，则Y<X；所述的Y＝even(Lf_s/f_0min)，其中even(x)为取最接近x的偶数，L为采样周期数，f_s为采样频率，f_0min为并网断路器两侧电压的最小可能频率；

步骤D：当采样数据个数达到设定值Y个后，从缓冲区sb、gb中分别截取一段数据u_s0和u_g0进行单线谱DFT分析，得到分析结果；其中数据u_s0的长度为N_s0，u_g0的长度为N_g0；对u_s0数据段进行单线谱DFT分析结果记为X_s0L，对u_g0数据段进行单线谱DFT分析结果记为X_g0L；

步骤E：再分别采样P个点并存入对应的缓冲区sb和gb中，某时刻最新采样点存储的位置为i，则最新采样点在缓冲区中的位置自然就变成了i+P，再从sb和gb中分别截取一段数据u_s1和u_g1进行单谱线DFT分析，得到分析结果，其中，数据u_s1的长度为N_s1，数据u_g1的长度为N_g1；对u_s1数据段进行单线谱DFT分析结果记为X_s1L，对u_g1数据段进行单线谱DFT分析结果记为X_g1L；

步骤F：分别利用X_s0L、X_s1L和X_g0L、X_g1L根据变窗长的相位差法分别计算电压u_s和电压u_g的频率f_us和f_ug、幅值A_s和A_g、相位

和

并判断电压u_s、u_g是否满足并网条件，满足则发出合闸指令，结束测量；不满足，则进入下一步；

步骤G：令X_s0L＝X_s1L、N_s0＝N_s1,X_g0L＝X_g1L、N_g0＝N_g1，为下次测频做准备，同时利用所测频率f_us和f_ug、计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1,即更新每周期采样点数N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数，用于让下一次DFT分析截取的数据段接近整数个周期；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；

同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2；

步骤H：由上步骤得到的需要截取的数据段的长度N_s1和N_g1，重复步骤E-G，直至满足发出并网指令，结束测量。

2.根据权利要求1所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述的步骤C中，当数据存满缓冲区sb或缓冲区gb，新采样的数据会覆盖最老的采样数据；其中，当某时刻采样值存入sb(Ln-1)，缓冲区会被填满之后，新的采样值uy被存入sb(0)，覆盖了sb(0)中老的采样数据，这样缓冲区内的数据总是最新的采样数据；针对采样数据的检索访问方式为sb(j_d)＝sb(mod(j_d,Ln))，mod(x，y)表示求x/y的余数，Ln是缓冲区长度；缓冲区gb也采用同样的存储和访问方式。

3.根据权利要求2所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述的从缓冲区sb中截取长度为N的一段数据的方法具体为：以当前最新采样点为末尾，向前数N个点；假设某时刻最新采样点存储的位置为i，此位置的数据为sb(i),则单线谱DFT分析应该截取的数据段为从sb(i)开始向前数N个数据，即应该是由sb(i-N+1)至sb(i)存储单元中的采样数据构成，因为这两个点之间的数据个数为i-(i-N+1)+1＝N个，N即为数据的长度；所述的从缓冲区gb中截取长度为N的一段数据的方法与从缓冲区sb中截取长度为N的一段数据的方法相同，这里N表示N_s0、N_s1或N_g0、N_g1中的任何一个参数。

4.根据权利要求3所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述的单谱线DFT的具体过程为：分别对缓冲区sb和gb中最近采样的N个点数据与N个点对应的一维加权系数表对应相乘并累加，得到第L条谱线的DFT 结果。

5.根据权利要求1所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述步骤A中N点卷积窗对应的一维加权系数表具体包括如下步骤：

步骤S1：定义N1＝N/2，步骤S1.1：离线生成N1点Hanning窗；

步骤S1.2：对N1点Hanning窗进行自卷积得到2N1-1点卷积窗，对2N1-1点卷积窗末端补1个零得到2N1点卷积窗，即N点卷积窗；

步骤S2：由N点卷积窗离线生成N点正余弦因子；

步骤S3：将N点卷积窗与N点正余弦因子相乘即可得到N点卷积窗对应的N点一维加权系数表；

步骤S4：建立数组addr，存储步骤S3得到的N点卷积窗对应的N点一维加权系数表的首地址；以频率上限对应的采样点数N_min为基准，N_min点一维加权系数表的首地址存在addr第一个元素位置，其他采样点数的一维加权系数表首地址按点数递增顺序存放，则N点一维加权系数表的首地址为addr((N-N_min)/2)。

6.根据权利要求1所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述利用变窗长相位差法计算电压频率f_us的方法为：

记数据段u_s1为最近一次进行DFT分析的数据段，N_s1表示u_s1数据段的长度，u_s0为上一次进行DFT分析的数据段，N_s0表示u_s0数据段的长度；u_s1数据段是u_s0数据段又采样P个点后截取的，因此s侧两次DFT分析的采样序列起始点对应的间隔为N_s0-N_s1+P，因此u_s1和u_s0的采样起始点的相位差值满足如下关系：

φ_s1-φ_s0＝2πf_us×(N_s0-N_s1+P)T_s；(19)

其中，φ_s1为u_s1的采样起始点的相位，φ_s0为u_s0的采样起始点的相位，T_s为采样周期，f_s为采样频率，T_s＝1/f_s；f_us为s侧的电压信号频率；根据傅立叶变换，采样起始点的相位θ又和谱线X(k)有如下关系：

其中，arg为求复数的相角，k₀＝Nf₀/f_s，N为DFT分析数据段的长度，f₀为信号频率；W(2πx/N)为汉宁卷积窗的频谱；

为减小计算量，采用单谱线DFT，即只计算第L条谱线的值X_s0L、X_s1L，X_g0L、X_g1L，结合汉宁卷积窗的频谱公式有：

因此：

因为相角是以2π为周期的，因此根据式22求出的相角差和式19等式右侧计算出的相角差可能有2π整数倍的差距，因此可得：

其中，q为一整数，进一步的：

2πf_us×(N_s0-N_s1+P)T_s-πf_usT_s(N_s0-N_s1)＝arg(X_s1L)-arg(X_s0L)±2πq；

进一步的：

πf_us(N_s0-N_s1+2P)T_s＝arg(X_s1L)-arg(X_s0L)±2πq；

令f_us＝50+df，df∈[-5Hz,5Hz]，反映了频率偏离50Hz的程度，则上式进一步推导如下：

则：

进而s侧的电压信号频率f_us＝50+df：

同理：

g侧的电压频率f_ug的方法为：

f_ug＝50+df。

7.根据权利要求6所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述利用变窗长相位差法计算s侧的电压相位

和g侧的电压相位

通过如下公式计算：

8.根据权利要求7所述的一种变窗长的准同期并网参数测量方法，其特征在于：

计算信号幅值A_s、A_g的公式分别为：

其中，

Δf＝f_s/N_s1,

其中，

Δf＝f_s/N_g1,

9.根据权利要求8所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述步骤F中利用所测频率f_us和f_ug计算下一次DFT分析所需要截取的数据段的长度N_s1和N_g1，计算公式为N_s1＝even(Lf_s/f_us)、N_g1＝even(Lf_s/f_ug)，其中even(x)为取最接近x的偶数；如果N_s0-N_s1+2P＝0，则令N_s1＝N_s1-2，以避免在变窗长相位差测频过程出现分母为零的情况；同理如果N_g0-N_g1+2P＝0，则令N_g1＝N_g1-2。

10.根据权利要求9所述的一种变窗长的准同期并网参数测量方法，其特征在于：所述的

的计算方法为：

令x在(-0.5,0.5)之间以0.0001的速度递增变化，对应的计算F(x)，然后利用matlab多项式逼近函数polyfit(x,F(x),n)，n为多项式最高幂次，得到F(x)的近似多项式：

F(x)＝8.00000000+1.28987025x²+0.12451951x⁴+0.00994809x⁶。

Images