本发明涉及电力系统，尤其涉及一种基于最小特征值的电力系统振荡稳定性分析方法及系统。
背景技术：
1、为实现“双碳”目标，新能源机组正逐步替代传统火电机组，成为主要的能源供应主体。新能源机组通常分为两类，一类以直驱风电机组和光伏机组为代表，通过变流器与电网相连，另一类以双馈风电机组为代表，通过变流器控制异步电机与电网相连。其中前一类机组可能与弱电网相互作用，从而引发电磁振荡现象；后一类机组可能与电网中的容性设备(例如，串补电容，并联电容)相连，产生电磁振荡现象。振荡的发生会严重危害电力系统安全稳定运行，因此需要提出一种方法能够准确定量分析振荡稳定性。
2、现有振荡稳定性分析方法常见有特征值分析方法和阻抗法。特征值分析方法首先建立系统的微分代数方程组，然后将其转化为状态空间方程的形式，通过计算矩阵特征值来得到系统各振荡模式。现有阻抗法大多基于(广义)nyquist方法来进行稳定性分析。该方法的首先建立设备侧和系统侧的阻抗模型，然后计算系统开环传递函数，最后通过(广义)nyquist判据来实现稳定性判断。对于特征值分析法，在实际工程中，由于系统电力电子设备数量多且拓扑复杂，导致整体系统阶数高，计算量极大。对于阻抗法，需要建立设备侧和系统侧的阻抗模型，该阻抗模型建立过程需要系统节点导纳矩阵的求逆操作。对于大规模电力系统而言，系统节点数目多，节点导纳矩阵的阶数会很高，导致矩阵求逆操作的计算量很大，系统振荡稳定性分析过慢。因此需要提出一种方法能够准确定量分析振荡稳定性，并且满足计算速度快的条件。
技术实现思路
1、本发明提供一种基于最小特征值的电力系统振荡稳定性分析方法及系统，能够准确定量分析大规模电力系统的振荡稳定性，同时相比于现有阻抗法和特征值分析法，本发明具有更快的分析速度。
2、有益效果：
3、第一方面，本发明实施例提供一种基于最小特征值的电力系统振荡稳定性分析方法，包括：
4、获取电力系统中各个电力设备的运行数据，所述运行数据包括各个电力设备的工作频段以及工作时所产生的小信号电压、电流相量；
5、根据所述运行数据构建各个电力设备各自对应的阻抗模型；
6、使用电力系统的拓展节点-支路关联矩阵对所述各个阻抗模型进行稀疏化与耦合，得到采用稀疏矩阵表示的频率耦合阻抗网络模型；
7、计算获得所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，其中k为预设值；
8、对所述各个稀疏矩阵各自对应的最小k个特征值进行特征值曲线分析，确定所述电力系统的振荡稳定性。
9、本发明实施例提供一种基于最小特征值的电力系统振荡稳定性分析方法，首先根据各个电力设备的运行数据构建各个对应的阻抗模型，然后根据电力系统各个节点的拓扑连接方式即拓展节点-支路关联矩阵将各个阻抗模型进行耦合并稀疏化，得到频率耦合阻抗网络模型，相比于传统的阻抗法，本发明所提供的方法在建模过程中对整体的频率耦合阻抗网络模型进行了稀疏化，有效降低了后续稳定性分析过程中的计算量，提高了系统的计算速度；进一步的，本发明还利用了最小特征值分析方法，计算所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，并对所述各个稀疏矩阵各自对应的最小k个特征值进行特征值曲线分析，进而确定所述电力系统的振荡稳定性，通过最小特征值代替传统的节点导纳矩阵来表征电力系统动态特性，使得计算效率大大提高。而现有技术不足以在系统节点数目极大，例如超过10000节点时，系统的振荡稳定性分析速度明显变慢，不能及时得到系统振荡稳定性结果。
10、进一步的，所述根据所述运行数据构建各个电力设备各自对应的阻抗模型，具体为：
11、
12、其中，ω表示频率，取值在ωmin～ωmax范围内，ωmin和ωmax分别为某个电力设备的工作频段中的最低频率和最高频率；δup和δip分别为小信号电压、电流相量，δupc和δipc分别为δup和δip对应的频率耦合小信号电压、电流相量，上标*表示取共轭；y11、y12、y21、y22为某个电力设备对应的阻抗模型2×2矩阵的四个元素。
13、进一步的，所述使用电力系统的拓展节点-支路关联矩阵对所述各个阻抗模型进行稀疏化与耦合，得到采用稀疏矩阵表示的频率耦合阻抗网络模型，具体为：
14、
15、其中，为频率ω下的对角矩阵，其对角元为各个电力设备对应的阻抗模型，采用稀疏矩阵表示，为拓展节点-支路关联矩阵，i2为二维单位矩阵，为kronecker乘积符号，a为节点-支路关联矩阵。
16、在一种可能实现的方式中，当所述频率为第一频率时，所述计算获得所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，包括：
17、初始化获得初始降维向量，并根据所述初始降维向量构建第一降维矩阵；
18、根据所述第一降维矩阵对所述第一频率对应的第一稀疏矩阵进行降维，得到降维稀疏矩阵；
19、计算获得所述降维稀疏矩阵的所有特征值并按照各个特征值的幅值对所述各个特征值进行排序，得到幅值最小的k个临时特征值；
20、根据预设更新公式更新所述初始降维向量，并根据更新后的初始降维向量更新所述k个临时特征值，不断迭代更新所述初始降维向量和所述k个临时特征值，直到各个所述k个临时特征值均满足第一预设条件，此时的k个临时特征值即为所述第一稀疏矩对应的最小k个特征值。
21、本发明实施例提供一种计算稀疏矩阵最小k个特征值的方法，通过迭代更新的方式不断更新初始降维向量和k个临时特征值，直到各个所述k个临时特征值均满足第一预设条件时，结束迭代，输出k个临时特征值。其中，在每一次迭代时，通过初始降维向量构建第一降维矩阵并对频率耦合阻抗网络模型中的第一稀疏矩阵进行降维，得到降维稀疏矩阵，对于大规模电力系统而言，由于系统节点数目多，节点导纳矩阵的阶数会很高，若使用传统的特征值分析方法会导致特征值的计算量很大，进而影响系统振荡稳定性分析速度，而本方法在每次计算特征值前，先对稀疏矩阵进行降维，减少后续求特征值的计算量，同时，为了提高振荡稳定性分析的准确性，通过多次迭代的方式不断更新初始降维向量和k个临时特征值，直到所计算出来的k个临时特征值满足分析条件才停止迭代，输出最小k个特征值用于后续的振荡稳定性分析，提高了振荡稳定性分析的准确性与速度。
22、在一种可能实现的方式中，所述根据所述第一降维矩阵对所述第一频率对应的第一稀疏矩阵进行降维，得到降维稀疏矩阵，包括：
23、使用gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系；
24、根据预设的迭代求解公式迭代求解所述第一降维矩阵和所述降维稀疏矩阵中的各个元素，得到所述降维稀疏矩阵，其中，所述第一降维矩阵中的第一列列向量由所述初始降维向量计算获得，所述降维稀疏矩阵中的各个元素为gram-schmidt系数。
25、本发明实施例提供一种对稀疏矩阵进行降维的方法，通过gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系，然后逐步求解第一降维矩阵和降维稀疏矩阵中的各个元素，得到降维稀疏矩阵，实现了对第一稀疏矩阵的降维，提高了振荡稳定性分析的速度。
26、进一步的，所述使用gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系，具体公式为：
27、
28、其中，a为所述第一稀疏矩阵，qm+1为所述第一降维矩阵，为所述降维稀疏矩阵，m+1为所述第一降维矩阵的列向量数量，m为所述降维稀疏矩阵的列向量数量；
29、所述根据预设的迭代求解公式迭代求解所述第一降维矩阵和所述降维稀疏矩阵中的各个元素，得到所述降维稀疏矩阵，具体公式为：
30、
31、
32、
33、
34、其中，x为所述初始降维向量，q1为所述第一降维矩阵的第一列列向量，qi+1为所述第一降维矩阵的第i+1列列向量，hji为所述gram-schmidt系数，在所述迭代求解过程中，令i的取值从1逐步增加到m，不断求解qi+1和hji，最后求解出hm中的所有元素，得到所述降维稀疏矩阵。
35、在一种可能实现的方式中，所述根据预设更新公式更新所述初始降维向量，包括：
36、使用qr分解法对所述降维稀疏矩阵进行分解，得到所述降维稀疏矩阵对应的正交矩阵，具体公式为：
37、hm-λii＝qiri
38、其中，hm为所述降维稀疏矩阵去掉最后一行所得到的m×m矩阵，λi为按照幅值升序排序后的第i个特征值，i可取k+1至m，i为单位向量，qi和ri分别为qr分解得到的正交矩阵和上三角矩阵；
39、根据所述正交矩阵和预设更新公式更新所述初始降维向量，具体公式为：
40、
41、
42、
43、其中，v为本次迭代更新后的初始降维向量，v为本次迭代更新前的初始降维向量；qk+1为所述第一降维矩阵的第k+1列，为每次qr分解得到的各个qi的乘积，为的第m行第k列的元素；hm为所述降维稀疏矩阵去掉最后一行的m×m维矩阵，为的第k行第k+1列的元素。
44、本发明实施例提供一种更新所述初始降维向量的方法，在求出降维稀疏矩阵的所有特征值并排序后，利用最小k个临时特征值以外的其他特征值，通过qr分解法对所述降维稀疏矩阵进行分解，得到对应的正交矩阵和上三角矩阵，然后根据所述正交矩阵和预设更新公式更新所述初始降维向量，用于下一轮迭代更新各个特征值，实现初始降维向量与特征值的迭代更新，搜索到满足预设条件的最小k个特征值，后续使用各个频率下的最小k个特征值进行系统振荡稳定性分析，提高了振荡稳定性分析的准确性。
45、进一步的，所述直到各个所述k个临时特征值均满足第一预设条件，具体为：
46、
47、其中，tol为设定的阈值，ui为第i个临时特征值λi对应的特征向量，i可取1至k，hm+1,m为所述降维稀疏矩阵中的最后一行最后一列的元素，qm+1为所述第一降维矩阵中的第m+1列列向量。
48、在一种可能实现的方式中，所述对所述各个稀疏矩阵各自对应的最小k个特征值进行特征值曲线分析，确定所述电力系统的振荡稳定性，包括：
49、结合各个频率各自对应的最小k个特征值，构建k条频率-特征值曲线；
50、分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线；
51、对各个对数导数曲线进行振荡模式识别，如果第一对数导数曲线的实部曲线存在极值点，并且所述第一对数导数曲线对应的虚部曲线在所述极值点处的斜率为负，则所述第一对数导数在所述极值点处存在主导振荡模式；
52、检索所有识别到的主导振荡模式，若存在实部大于0的主导振荡模式，则确定所述电力系统不稳定，否则，确定所述电力系统稳定。
53、本发明实施例提供一种特征值曲线分析方法，首先结合各个频率各自对应的最小k个特征值，构建k条频率-特征值曲线，由于之前计算获得了各个频率各自对应的k个特征值，因此每个频率下都有k个按幅值由小到大排序的特征值，因此可以构建k条频率-特征值曲线，在每条频率-特征值曲线中，每个频率对应一个特征值；分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线并根据各个对数导数曲线进行振荡模式识别，只需要判断对数导数曲线的实部曲线是否存在极值点以及虚部曲线在该极值点处的斜率是否为负，即可识别出所有的主导振荡模式；最后判断各个主导振荡模式的实部是否大于0，即可得到振荡稳定性分析的结果。本发明实施例通过构建特征值曲线，并对特征值曲线进行曲线分析的方式，可以快速且准确地评估电力系统的振荡稳定性，提高了振荡稳定性分析的速度。
54、进一步的，所述分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线，具体公式为：
55、
56、其中，dl(λi)为特征值第i条频率-特征值曲线对应的对数导数曲线，λi(ω)为第i条频率-特征值曲线，i取值为1到k，ω为当前频率，δω为频率步长，j为虚数符号；
57、所述第一对数导数在所述极值点处存在主导振荡模式，具体公式为：
58、z＝α+jω
59、
60、其中，z＝α+jω为所述主导振荡模式的表达式，α为所述主导振荡模式的实部，-α为振荡模式阻尼，ω为振荡模式频率，re(dl(λi))为所述对数导数曲线的实部。
61、第二方面，相应的，本发明提供一种基于最小特征值的电力系统振荡稳定性分析系统，包括获取模块、模型构建模块、耦合模块、特征值计算模块以及分析模块；
62、其中，所述获取模块用于获取电力系统中各个电力设备的运行数据，所述运行数据包括各个电力设备的工作频段以及工作时所产生的小信号电压、电流相量；
63、所述模型构建模块用于根据所述运行数据构建各个电力设备各自对应的阻抗模型；
64、所述耦合模块用于使用电力系统的拓展节点-支路关联矩阵对所述各个阻抗模型进行稀疏化与耦合，得到采用稀疏矩阵表示的频率耦合阻抗网络模型；
65、所述特征值计算模块用于计算获得所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，其中k为预设值；
66、所述分析模块用于对所述各个稀疏矩阵各自对应的最小k个特征值进行特征值曲线分析，确定所述电力系统的振荡稳定性。
67、进一步的，所述模型构建模块根据所述运行数据构建各个电力设备各自对应的阻抗模型，具体为：
68、
69、其中，ω表示频率，取值在ωmin～ωmax范围内，ωmin和ωmax分别为某个电力设备的工作频段中的最低频率和最高频率；δup和δip分别为小信号电压、电流相量，δupc和δipc分别为δup和δip对应的频率耦合小信号电压、电流相量，上标*表示取共轭；y11、y12、y21、y22为某个电力设备对应的阻抗模型2×2矩阵的四个元素。
70、进一步的，所述耦合模块使用电力系统的拓展节点-支路关联矩阵对所述各个阻抗模型进行稀疏化与耦合，得到采用稀疏矩阵表示的频率耦合阻抗网络模型，具体为：
71、
72、
73、其中，为频率ω下的对角矩阵，其对角元为各个电力设备对应的阻抗模型，采用稀疏矩阵表示，为拓展节点-支路关联矩阵，i2为二维单位矩阵，为kronecker乘积符号，a为节点-支路关联矩阵。
74、在一种可能实现的方式中，当所述频率为第一频率时，所述特征值计算模块计算获得所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，包括：
75、初始化获得初始降维向量，并根据所述初始降维向量构建第一降维矩阵；
76、根据所述第一降维矩阵对所述第一频率对应的第一稀疏矩阵进行降维，得到降维稀疏矩阵；
77、计算获得所述降维稀疏矩阵的所有特征值并按照各个特征值的幅值对所述各个特征值进行排序，得到幅值最小的k个临时特征值；
78、根据预设更新公式更新所述初始降维向量，并根据更新后的初始降维向量更新所述k个临时特征值，不断迭代更新所述初始降维向量和所述k个临时特征值，直到各个所述k个临时特征值均满足第一预设条件，此时的k个临时特征值即为所述第一稀疏矩对应的最小k个特征值。
79、在一种可能实现的方式中，所述根据所述第一降维矩阵对所述第一频率对应的第一稀疏矩阵进行降维，得到降维稀疏矩阵，包括：
80、使用gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系；
81、根据预设的迭代求解公式迭代求解所述第一降维矩阵和所述降维稀疏矩阵中的各个元素，得到所述降维稀疏矩阵，其中，所述第一降维矩阵中的第一列列向量由所述初始降维向量计算获得，所述降维稀疏矩阵中的各个元素为gram-schmidt系数。
82、进一步的，所述使用gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系，具体公式为：
83、
84、其中，a为所述第一稀疏矩阵，qm+1为所述第一降维矩阵，为所述降维稀疏矩阵，m+1为所述第一降维矩阵的列向量数量，m为所述降维稀疏矩阵的列向量数量；
85、所述根据预设的迭代求解公式迭代求解所述第一降维矩阵和所述降维稀疏矩阵中的各个元素，得到所述降维稀疏矩阵，具体公式为：
86、
87、
88、
89、
90、其中，x为所述初始降维向量，q1为所述第一降维矩阵的第一列列向量，qi+1为所述第一降维矩阵的第i+1列列向量，hji为所述gram-schmidt系数，在所述迭代求解过程中，令i的取值从1逐步增加到m，不断求解qi+1和hji，最后求解出hm中的所有元素，得到所述降维稀疏矩阵。
91、在一种可能实现的方式中，所述根据预设更新公式更新所述初始降维向量，包括：
92、使用qr分解法对所述降维稀疏矩阵进行分解，得到所述降维稀疏矩阵对应的正交矩阵，具体公式为：
93、hm-λii＝qiri
94、其中，hm为所述降维稀疏矩阵去掉最后一行所得到的m×m矩阵，λi为按照幅值升序排序后的第i个特征值，i可取k+1至m，i为单位向量，qi和ri分别为qr分解得到的正交矩阵和上三角矩阵；
95、根据所述正交矩阵和预设更新公式更新所述初始降维向量，具体公式为：
96、
97、
98、
99、其中，为本次迭代更新后的初始降维向量，v为本次迭代更新前的初始降维向量；qk+1为所述第一降维矩阵的第k+1列，为每次qr分解得到的各个qi的乘积，为的第m行第k列的元素；hm为所述降维稀疏矩阵去掉最后一行的m×m维矩阵，为的第k行第k+1列的元素。
100、进一步的，所述直到各个所述k个临时特征值均满足第一预设条件，具体为：
101、
102、其中，tol为设定的阈值，ui为第i个临时特征值λi对应的特征向量，i可取1至k，hm+1,m为所述降维稀疏矩阵中的最后一行最后一列的元素，qm+1为所述第一降维矩阵中的第m+1列列向量。
103、在一种可能实现的方式中，所述分析模块包括频率-特征值曲线构建单元、对数导数曲线构建单元、振荡模式识别单元以及系统稳定性判断单元：
104、其中，所述频率-特征值曲线构建单元用于结合各个频率各自对应的最小k个特征值，构建k条频率-特征值曲线；
105、所述对数导数曲线构建单元用于分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线；
106、所述振荡模式识别单元用于对各个对数导数曲线进行振荡模式识别，如果第一对数导数曲线的实部曲线存在极值点，并且所述第一对数导数曲线对应的虚部曲线在所述极值点处的斜率为负，则所述第一对数导数在所述极值点处存在主导振荡模式；
107、所述系统稳定性判断单元用于检索所有识别到的主导振荡模式，若存在实部大于0的主导振荡模式，则确定所述电力系统不稳定，否则，确定所述电力系统稳定。
108、进一步的，所述对数导数曲线构建单元分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线，具体公式为：
109、
110、其中，dl(λi)为特征值第i条频率-特征值曲线对应的对数导数曲线，λi(ω)为第i条频率-特征值曲线，i取值为1到k，ω为当前频率，δω为频率步长，j为虚数符号；
111、所述第一对数导数在所述极值点处存在主导振荡模式，具体公式为：
112、z＝α+jω
113、
114、其中，z＝α+jω为所述主导振荡模式的表达式，α为所述主导振荡模式的实部，-α为振荡模式阻尼，ω为振荡模式频率，re(dl(λi))为所述对数导数曲线的实部。
1.一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，包括：
2.如权利要求1所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述根据所述运行数据构建各个电力设备各自对应的阻抗模型，具体为：
3.如权利要求1所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述使用电力系统的拓展节点-支路关联矩阵对所述各个阻抗模型进行稀疏化与耦合，得到采用稀疏矩阵表示的频率耦合阻抗网络模型，具体为：
4.如权利要求1所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，当所述频率为第一频率时，所述计算获得所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，包括：
5.如权利要求4所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述根据所述第一降维矩阵对所述第一频率对应的第一稀疏矩阵进行降维，得到降维稀疏矩阵，包括：
6.如权利要求5所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述使用gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系，具体公式为：
7.如权利要求4所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述根据预设更新公式更新所述初始降维向量，包括：
8.如权利要求4所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述直到各个所述k个临时特征值均满足第一预设条件，具体为：
9.如权利要求1所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述对所述各个稀疏矩阵各自对应的最小k个特征值进行特征值曲线分析，确定所述电力系统的振荡稳定性，包括：
10.如权利要求9所述的一种基于最小特征值的电力系统振荡稳定性分析方法，其特征在于，所述分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线，具体公式为：
11.一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，包括获取模块、模型构建模块、耦合模块、特征值计算模块以及分析模块；
12.如权利要求11所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述模型构建模块根据所述运行数据构建各个电力设备各自对应的阻抗模型，具体为：
13.如权利要求11所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述耦合模块使用电力系统的拓展节点-支路关联矩阵对所述各个阻抗模型进行稀疏化与耦合，得到采用稀疏矩阵表示的频率耦合阻抗网络模型，具体为：
14.如权利要求11所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，当所述频率为第一频率时，所述特征值计算模块计算获得所述频率耦合阻抗网络模型中各个频率对应的稀疏矩阵各自的最小k个特征值，包括：
15.如权利要求14所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述根据所述第一降维矩阵对所述第一频率对应的第一稀疏矩阵进行降维，得到降维稀疏矩阵，包括：
16.如权利要求15所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述使用gram-schmidt法建立第一稀疏矩阵、第一降维矩阵和降维稀疏矩阵之间的等式关系，具体公式为：
17.如权利要求14所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述根据预设更新公式更新所述初始降维向量，包括：
18.如权利要求14所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述直到各个所述k个临时特征值均满足第一预设条件，具体为：
19.如权利要求11所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述分析模块包括频率-特征值曲线构建单元、对数导数曲线构建单元、振荡模式识别单元以及系统稳定性判断单元：
20.如权利要求19所述的一种基于最小特征值的电力系统振荡稳定性分析系统，其特征在于，所述对数导数曲线构建单元分别求解各个频率-特征值曲线的对数导数，得到k条对数导数曲线，具体公式为：
