企业官网建设流程全解析

很多人第一次接触随机潮流，基本都是从IEEE34节点这个算例入手的。我自己也是这样：一开始看文献里“半不变量”“Gram-Charlier级数”这些词，觉得门槛很高，等真正把基于半不变量的概率潮流计算在IEEE34节点系统上完整跑通之后，才发现核心逻辑并不复杂，难的是把细节做对。这篇内容就把我怎么处理数据、怎么写Matlab代码、怎么校验结果、以及踩过的坑一次说清楚。它适合正在做配电网规划、新能源并网分析或者电压质量评估的同学，尤其是想绕开蒙特卡洛、追求计算速度的工程场景。

随机潮流的应用场景现在越来越多。分布式光伏和电动车大规模接入后，负荷曲线和出力曲线都不是确定的一条线，调度和规划要回答的问题也从“某个工况下电压合不合格”变成了“一年里电压越限的概率有多大”。这种问题用确定性潮流没法回答，只能靠随机潮流来求节点电压和支路潮流的概率分布。半不变量法是目前工程上兼顾速度和精度的主流选择，它不需要成千上万次采样，一次确定性潮流加若干次矩阵运算就能出结果。下面我按“为什么这么做-系统怎么建模-代码怎么实现-出了问题怎么排查”的顺序来写。

1. 为什么要用半不变量做概率潮流

1.1 确定性潮流的局限与随机潮流要解决什么问题

常规潮流计算给定一组负荷和发电机出力，求解得到一组节点电压和支路功率，这是确定性的。但实际系统里，负荷预测误差、光伏波动、风电随机性、电动汽车充电行为都会让“真实运行点”偏离“计算工况点”。以前配电网裕度大，按最大最小场景估算一下就够用；现在分布式电源渗透率上来，很多馈线末端电压波动很大，只算一个工况点根本看不出风险。

举个例子：某个节点电压在基准工况下是0.98 pu，看着没问题。但考虑负荷波动后，这个节点电压可能有30%的概率低于0.95 pu。确定性潮流给不出这30%，但随机潮流可以。随机潮流的输入是各个节点注入功率的分布（均值、方差、偏度等），输出是节点电压、支路潮流的概率密度函数PDF和累积分布函数CDF，进而得到越限概率和分位数。并网导则里经常要求提供P95、P99这类指标，这些指标必须依赖概率手段才能算出来。

随机潮流的方法大体分三类：蒙特卡洛模拟法、解析法、近似法。蒙特卡洛精度高、实现简单，但计算量太大；一次算例做一万次潮流抽样，单次潮流哪怕只要几十毫秒，累计起来也到分钟级，如果做数百个场景的规划比选，时间根本扛不住。解析法里面最典型的就是半不变量法，它把输入随机量的概率特征通过线性化潮流方程映射到输出量，计算速度非常快，非常适合配电网这种需要批量评估的场景。

1.2 半不变量法的数学骨架：矩、半不变量与Gram-Charlier展开

半不变量，也叫累积量（cumulant），本质上是随机变量的另一种“概率指纹”。我们常用均值、方差、偏度、峰度来描述分布，这些指标都可以由半不变量组合出来。半不变量最大的优势是它有很好的代数性质：

• 若多个独立随机变量相加，则和的各阶半不变量等于各个变量对应阶半不变量之和；
• 若随机变量做线性变换 (Y=aX+b)，则一阶半不变量满足 (\kappa_{Y,1}=a\kappa_{X,1}+b)，二阶及以上满足 (\kappa_{Y,k}=a^k\kappa_{X,k})。

这两条性质对概率潮流来说非常关键。潮流方程在基准运行点附近做泰勒展开并忽略二阶以上项，可以写成：

[ \Delta X = J^{-1} \Delta W ]

其中 (\Delta W) 是节点注入功率扰动，(\Delta X) 是节点电压幅值和相角扰动，(J) 是最后一次牛顿-拉夫逊迭代的雅可比矩阵。也就是说，输出状态量可以写成输入注入量的线性组合。借助半不变量的可加性和齐次性，不需要做卷积或数值积分，直接就能把输出状态量的各阶半不变量算出来。

但半不变量本身不是概率分布，要还原成我们熟悉的PDF和CDF，常用的工具是Gram-Charlier级数。它的思路是以标准正态分布为基础，用Hermite多项式逐项修正偏度、峰度等非正态特征。标准化变量 (z=(x-\mu)/\sigma) 后，密度函数写成：

[ f(x)=\phi(z)\left[1+\frac{\lambda_3}{6}H_3(z)+\frac{\lambda_4}{24}H_4(z)+\frac{\lambda_5}{120}H_5(z)+\cdots\right] ]

其中 (\lambda_k) 是标准化后的半不变量，(H_k) 是Hermite多项式。阶数越高，对偏度和尾部特征的刻画越精细。这套组合拳就是半不变量概率潮流的核心：先求半不变量，再用级数展开还原分布。

2. 方法选型：半不变量法与蒙特卡洛的取舍

2.1 蒙特卡洛为什么“慢”，半不变量为什么“快”

蒙特卡洛的原理很朴素：从输入变量的概率分布里抽样，每抽一组样本就做一次确定性潮流，最后统计所有潮流结果的分布特征。这个方法的好处是几乎不受系统非线性限制，也能处理任意分布和相关性。代价就是计算量线性增长。以IEEE34节点系统为例，一次牛顿-拉夫逊潮流在普通PC上大约需要20到50毫秒，1万次抽样就是几分钟到十几分钟。如果系统规模扩大到几百上千个节点，抽样次数通常还得增加，耗时更快。

半不变量法完全不同，它的计算量主要由一次基准潮流和几个矩阵乘法决定。基准潮流只做一次，雅可比矩阵只求逆一次，之后无论输入分布多复杂，都不会再做潮流迭代。从实测来看，IEEE34节点算例在Matlab里跑完整个半不变量计算加Gram-Charlier级数还原，通常不到0.1秒，比蒙特卡洛快两到三个数量级。

2.2 我的选型建议与适用边界

选方法不能只看快慢，还要看场景。下面这个表是我在实际项目里做选型时的对比思路：

我的习惯是：先用半不变量法做快速筛查，找出电压越限概率较高的节点，再对这些高风险区域用蒙特卡洛做精细化校验。这样既保证速度，又不牺牲对关键节点的精度。如果系统里调压器或者OLTC动作对结果影响很大，半不变量法就要非常小心，因为分接头切换是离散事件，线性化模型天然刻画不了。

3. IEEE34节点系统建模与随机输入设置

3.1 IEEE34节点测试馈线基本情况

IEEE34节点测试馈线是IEEE PES配电网测试馈线系列里非常有代表性的一个算例。它原型在美国亚利桑那州，是一条实际运行过的配电馈线，主要电压等级是24.9 kV，末端有一段4.16 kV的降压支路。系统包含两条主要的调压器支路、电容器组、种类丰富的负荷，以及单相、两相、三相混合的线路结构。相比于IEEE 33节点这类纯三相平衡算例，IEEE34节点更接近真实配电网的复杂情况。

不过有一点要提醒：原版IEEE34节点数据是三相不平衡模型，直接拿来做单相或三相对称的随机潮流，需要对数据做一些等值处理。在我这次实现里，为了把重点放在半不变量方法本身，我是把三相线路按单相等值（取单位长度正序阻抗）来处理的。如果你想做完整的三相不平衡随机潮流，思路完全一样，只是把雅可比矩阵从 (2n \times 2n) 扩成 (3n\times3n) 的复合矩阵，代码量会增加不少，但核心还是半不变量那套流程。

3.2 负荷与分布式电源的不确定性建模

随机潮流的效果很大程度上取决于输入分布设得合不合理。常规负荷我习惯用正态分布，均值取基准潮流里的负荷值，标准差取均值的5%到20%。这个范围基本覆盖了负荷预测误差的典型水平。正态分布有个偷懒的好消息：它的三阶及以上半不变量都是0，计算时只需要考虑一阶和二阶。

光伏出力和风电出力就不能简单用正态了。光伏出力在晴空条件下接近Beta分布，在云层遮挡时会出现大量零值和低出力段；风电出力则更接近Weibull分布。处理这类分布有两种方式：已知分布参数时，直接推导半不变量解析式；不知道分布参数但有历史数据时，就从样本数据估算中心矩再递推半不变量。后面这种更通用，我的代码里就是用样本估算的方式实现的。

分布式电源并网的位置也需要注意。如果多个光伏电站处在同一个区域，天气条件相近，它们的出力具有很强的正相关性。这时候如果按独立变量处理，会低估系统电压波动的范围。处理相关性的一个可行方案是先对相关正态变量做Cholesky分解生成相关样本，再从样本估算半不变量；或者用Copula函数建模尾部相关性。不过这会大幅增加实现复杂度，建议在基础版本跑通之后再加。

3.3 MatLab数据组织与输入文件设计

Matlab实现里最容易被忽视的是数据结构设计。IEEE34节点系统有34个节点和33条支路，如果全部裸写在脚本里，后面调试会非常痛苦。我的做法是用结构体集中管理：

% 节点数据 bus_data = struct( ... 'id', (1:34)', ... 'type', type_array, ... % 1-PQ节点 2-PV节点 3-平衡节点 'P_load', P_load, ... 'Q_load', Q_load, ... 'V_base', V_base); % 支路数据 branch_data = struct( ... 'from', from_idx, ... 'to', to_idx, ... 'R', R_pu, ... 'X', X_pu, ... 'B', B_pu, ... 'tap', tap_ratio);

这种结构的好处是：后续做基准潮流、算灵敏度矩阵、做随机注入设置，所有函数都只与结构体交互，不会出现索引对不上的问题。原始数据文件可以从IEEE PES测试馈线页面下载，然后在Matlab里写一个转换脚本，把标准的bus.txt和branch.txt转成上述结构体，顺便统一单位成标幺值。单位不统一是后面出错的一大来源，建议在数据加载阶段就全部转成pu。

4. 半不变量概率潮流计算流程与Matlab实现

4.1 完整计算流程

我在代码里把整个流程分成了十个步骤，每一步都单独封装函数，方便排查：

1. 加载IEEE34节点数据，转标幺值，形成bus_data和branch_data。
2. 取负荷期望值作为确定性潮流输入，运行牛顿-拉夫逊潮流。
3. 保存基准状态量 (X_0)（电压幅值和相角）和最后一次雅可比矩阵 (J)。
4. 对各节点注入功率设置随机分布，计算各阶半不变量。
5. 对雅可比矩阵求逆，提取与PQ节点对应的灵敏度子矩阵 (S)。
6. 合成节点电压幅值和相角的各阶半不变量。
7. 计算支路潮流的灵敏度矩阵，合成支路潮流的各阶半不变量。
8. 对每个输出量做标准化，用Gram-Charlier级数还原PDF和CDF。
9. 计算越限概率和指定分位数。
10. 与蒙特卡洛结果对比，输出误差指标。

第2步和第3步是关键中的关键。基准潮流的质量直接决定线性化误差。我在实际测试中发现，如果基准潮流本身不收敛或者雅可比矩阵病态，后面算出来的概率分布会非常离谱，而且很难排查，因为问题不在随机算法，而在最基础的潮流求解。

4.2 基准潮流与雅可比矩阵获取

对配电网这种辐射状结构，直接用牛顿-拉夫逊容易因为初值差而不收敛。我习惯先用前推回代法跑一遍，得到一组更接近真实解的初值，再切到牛顿-拉夫逊求精确解和雅可比矩阵。这样做还有个额外好处：前推回代对配电网收敛性更稳，牛顿-拉夫逊则能自然地给出雅可比矩阵，两者互补。

% 牛顿-拉夫逊潮流核心函数（示意） function [V, theta, J] = nr_powerflow(bus_data, branch_data, pq_idx) % 初始化 V = bus_data.V_base; theta = zeros(size(V)); Ybus = build_ybus(bus_data, branch_data); % 迭代求解 for iter = 1:20 [P_calc, Q_calc] = calc_injection(V, theta, Ybus); dP = bus_data.P_spec - P_calc; dQ = bus_data.Q_spec - Q_calc; dW = [dP(pq_idx); dQ(pq_idx)]; J = build_jacobian(V, theta, Ybus, pq_idx); dX = J \ dW; % 更新电压和相角 ... if norm(dW, inf) < 1e-8 break; end end end

这里有个容易踩的坑：雅可比矩阵的列顺序、行顺序必须和后面灵敏度矩阵的索引完全对应。否则半不变量合成时，矩阵乘法的维度对不上，结果全是NaN。我在代码里统一用pq_idx来标记PQ节点集合，所有涉及注入功率的向量和矩阵都从这些索引取出，确保顺序一致。

4.3 输入随机变量半不变量计算

前面说了，正态分布的高阶半不变量都是0，所以如果只考虑负荷波动，输入侧的半不变量计算非常简单。但光伏和风电往往带偏度，这时候就要从样本或者已知PDF中求半不变量。下面这个函数是我从历史样本估算中心矩再递推半不变量的通用实现：

function kappa = cumulant_from_samples(x, order) % 从样本估算各阶半不变量 % 输入x为列向量，order为最高阶数 mu = mean(x); mu2 = mean((x - mu).^2); mu3 = mean((x - mu).^3); mu4 = mean((x - mu).^4); kappa = zeros(1, order); kappa(1) = mu; kappa(2) = mu2; kappa(3) = mu3; kappa(4) = mu4 - 3 * mu2^2; % 更高阶可以根据需要继续递推 end

计算样本半不变量时要注意数值稳定性。如果样本量太小，三阶和四阶中心矩的估计误差会被放大，尤其是四阶半不变量计算里有减去 (3\mu_2^2) 的操作，很容易因为两个大数相减出现明显误差。所以样本量建议至少几千条，并且最好对原始数据进行去趋势处理。

4.4 输出状态量半不变量合成

有了灵敏度矩阵 (S) 和输入半不变量 (\kappa_W)，输出状态量的半不变量可以用下面的代码合成。这里默认各节点注入功率相互独立；如果考虑了相关性，需要在矩阵乘法里额外加入协方差项。

function kappa_X = combine_cumulants(X0, S, kappa_W, order) % X0：基准状态量列向量 % S：灵敏度矩阵，行对应状态量，列对应输入注入量 % kappa_W：输入半不变量矩阵，行对应阶数，列对应输入变量 % kappa_X：输出半不变量矩阵，行对应阶数，列对应状态量 n_state = size(S, 1); kappa_X = zeros(order, n_state); % 一阶半不变量 = 基准值 + 灵敏度矩阵乘输入均值 mu_W = kappa_W(1, :)'; kappa_X(1, :) = X0' + S * mu_W; % 二阶及以上半不变量：逐阶做矩阵乘法 for k = 2:order S_k = S .^ k; % 每个元素取k次幂 kappa_X(k, :) = (S_k * kappa_W(k, :).')'; end end

这段代码是整个算例的核心。理解它的关键在于：半不变量法的二阶以上部分不涉及基准值，只与灵敏度系数的 (k) 次方有关。举个例子，如果某个节点对输入功率的灵敏度系数是0.3，输入功率二阶半不变量是0.01，那么输出电压的二阶半不变量贡献就是 (0.3^2 \times 0.01 = 0.0009)。这就是“线性变换下 (a^k) 倍”的直接体现。

实际运行中，我只关心电压幅值的概率分布，不关心相角，所以在合成之后把相角对应的行丢掉就行。但注意灵敏度矩阵中相角部分对后续支路潮流计算是有用的，不能一开始就删。

4.5 Gram-Charlier级数还原概率分布

得到输出半不变量后需要还原成PDF和CDF。理论公式已经放在前面了，这里给出一个可直接用的Matlab函数骨架：

function [pdf_y, cdf_y, y] = gram_charlier_pdf(kappa_X, mu, sigma, n_points) % kappa_X：状态量各阶半不变量（行向量） % mu、sigma：该状态量的均值和标准差 % n_points：横轴采样点数 lambda3 = kappa_X(3) / sigma^3; lambda4 = kappa_X(4) / sigma^4; y = linspace(mu - 4*sigma, mu + 4*sigma, n_points); z = (y - mu) / sigma; pdf_y = normpdf(z); pdf_y = pdf_y .* (1 + lambda3/6 * (z.^3 - 3*z) ... + lambda4/24 * (z.^4 - 6*z.^2 + 3)); % 防止数值上出现负概率 pdf_y = max(pdf_y, 0); cdf_y = cumtrapz(y, pdf_y) / max(cumtrapz(y, pdf_y)); end

这里有几个细节要提醒。Gram-Charlier级数在概率密度的尾部容易出现振荡，甚至出现负概率；我在代码里做了max(pdf_y,0)的截断，但截断会影响CDF归一化，所以后面又用cumtrapz重新归一化了一次。如果你只是想求分位数，我更推荐用Cornish-Fisher展开，它对尾部的数值稳定性比Gram-Charlier好，尤其当偏度比较大、需要算P95/P99时，Cornish-Fisher的误差更可控。

另外，标准化用的均值应该是输出量的一阶半不变量，也就是合成后的kappa_X(1)，标准差是二阶半不变量的平方根。这个对应关系看着简单，但我在初版代码里犯过错：直接用蒙特卡洛样本的均值去替换合成后的均值，导致偏度对齐不上。正确做法是全部从半不变量体系内部取，不要混用两套统计口径。

5. 结果分析与精度校验

5.1 电压概率分布与越限概率怎么看

跑完Gram-Charlier级数还原，你会得到每个节点的电压PDF和CDF。在IEEE34节点算例里，末端节点（比如840以后）的电压分布通常更宽，对负荷波动更敏感。我习惯把每个节点的CDF画在同一张图上，然后叠加0.95 pu和1.05 pu这两条越限线，这样系统里哪些节点电压越限风险高，一眼就能看出来。

越限概率的计算公式很简单：

[ P_{\text{low}}=F(0.95), \quad P_{\text{high}}=1-F(1.05) ]

其中 (F(\cdot)) 是该节点电压的累积分布函数。需要提醒的是，配电网电压允许范围在不同标准下不一样，有些场景要求0.93~1.07 pu，有些对分布式电源接入后的电压要求更严格到0.95~1.05 pu。这个范围不要写死在代码里，最好作为参数传入，方便切换不同考核标准。

5.2 与蒙特卡洛结果的对比

半不变量法本质是线性化方法，所以它和蒙特卡洛的差异会随着系统非线性增强而变大。以我的算例为例，我在同样输入分布下跑了10000次蒙特卡洛，再和半不变量法结果对比，主要看三件事：均值、标准差、CDF曲线最大偏差。

实测下来，正常负荷波动范围内（标准差为均值10%），IEEE34节点末端电压均值的偏差通常在0.0005 pu以内，标准差偏差在2%左右，CDF曲线的最大绝对偏差在0.02到0.05之间。如果负荷波动比例继续加大，或者系统已经接近电压失稳边界，偏差会明显增大。这时候不要急着增加Gram-Charlier级数阶数，而是要考虑系统是否已经偏离了基准运行点的线性化区间，可能需要分段线性化。

我建议把蒙特卡洛结果作为“标尺”，保留在代码里作为回归测试项。以后改数据、改代码，只要跑一下对比脚本，看关键节点的均值误差是不是还在可接受范围内，就知道有没有改坏东西。这个方法帮我在后期迭代中省了很多排查时间。

5.3 支路潮流的概率特征

支路潮流的处理方式和节点电压类似，只是灵敏度矩阵的来源不同。支路有功和无功对节点注入的偏导，需要通过支路两端电压和相角推导出来。相比节点电压，支路潮流的非线性更强，尤其当支路重载时，Gram-Charlier级数的截断误差会更明显。

在IEEE34节点算例里，越靠近根节点（电源点）的支路潮流波动通常越大，因为上游支路承载了所有下游负荷和光伏出力的叠加波动。这个特征可以在支路潮流概率分布的均值变化曲线里清楚看到。如果你要做线路容量评估，推荐重点关注根节点附近支路的P95值，它代表了重载场景下的真实压力。

6. 常见问题与排查技巧

6.1 潮流不收敛或雅可比奇异

我在调试半不变量法时，遇到最多的不是随机计算问题，而是基准潮流本身不收敛。IEEE34节点系统里包含调压器和变压器，分接头设置不当或者PV节点处理不对，都可能让牛顿-拉夫逊迭代发散。

排查思路：先把所有负荷和光伏出力的随机波动清零，只用期望值做确定性潮流，看能不能收敛。如果确定性潮流都算不出来，问题就不在随机部分，而在系统建模。配电网结构弱环或者辐射状，推荐先用前推回代法求一组初值，再交给牛顿法，能解决大部分收敛问题。如果是雅可比矩阵奇异，检查是否有节点注入功率没有对应的方程，比如PV节点的无功方程被错误引入。

6.2 半不变量出现NaN或负方差

出现NaN通常有两个原因：一是样本量太少，高阶中心矩不稳定，导致四阶半不变量出现负值；二是灵敏度矩阵里有Inf，说明雅可比矩阵求逆出了问题。负方差则基本可以锁定二阶半不变量算错了。

排查时先检查输入半不变量矩阵，打印出每一阶的值，看看符号和量级是否合理。比如正态分布的三阶以上半不变量理论上应该是0，如果你的结果里出现很大的非零值，说明样本处理或者分布参数设置有误。另外，从样本估算半不变量时，我建议用中心矩公式而不是原始矩递推公式，后者在数值上更容易出现大数相减的精度丢失。

6.3 Gram-Charlier级数尾部振荡

Gram-Charlier级数在阶数有限的情况下，尾部会像多项式逼近那样出现振荡，甚至给出负概率密度。这个在偏度较大的光伏出力场景里特别常见。如果只是计算越限概率，可以用CDF的补集来做；如果是求分位数，换成Cornish-Fisher展开会更稳。

还有一种做法是增加级数阶数到6阶或8阶，但并不是阶数越高越好，过高阶数可能引入新的数值振荡。我的建议是：先看CDF曲线，如果尾部单调性破坏，就降低阶数或者换展开方式。不要迷信“高阶等于高精度”，工程计算里稳定比理论精度更重要。

6.4 离散控制设备和相关性怎么处理

调压器分接头和电容器组投切是离散控制，半不变量法本质上很难直接处理，因为在线性化模型里它们要么被固定，要么被忽略。一个可行的折中方案是：把分接头运行范围分为几个典型档位，每个档位下单独做一次半不变量潮流，最后按档位概率加权合成总分布。这个方法会增加计算量，但比完全忽略要可靠得多。

相关性处理则是另一个大坑。忽略节点注入之间的相关性会低估系统电压波动的空间同步性，尤其是同一馈线上多个光伏电站同时出力波动时。严格做法是在输入侧用Cholesky分解生成相关正态样本，再估算相关半不变量；或者用协方差矩阵扩展灵敏度合成公式。这些内容已经超出本文基础范围，但如果你做的是高渗透率光伏接入场景，建议尽早把相关性模型考虑进去。

下面这个表是我把排查经验整理成的速查表，实际调试时可以直接对照：

7. 实操心得与后续扩展

7.1 我踩过的几个坑

第一个坑是最开始做线性化时，我把负荷波动直接叠加到了初始电压上，而不是叠加到节点注入功率、再通过灵敏度矩阵映射到电压。结果一阶半不变量和蒙特卡洛对不上，排查了很久才发现是把基准状态量和注入扰动的因果关系搞反了。正确的逻辑是：基准潮流先解出 (X_0)，然后注入功率扰动 (\Delta W) 引起状态量扰动 (\Delta X=S\Delta W)，一阶半不变量是 (X_0+S\mu_W)，不是直接把 (\mu_W) 加到 (X_0) 上。

第二个坑是Matlab自带的moment函数在计算中心矩时，默认采用样本方差的无偏修正，这对高阶矩不一定适用。为了避免这种隐性差异，我干脆自己写中心矩计算函数，固定用mean((x-mu).^k)的方式，虽然少了统计上的无偏性，但在半不变量递推里保持了数值口径一致。做算法验证时，宁可“不标准”也不能“不一致”。

第三个坑是调压器的影响。IEEE34节点系统里的调压器会把末端电压强行维持在一个范围内，导致电压概率分布在接近调压阈值时出现明显的截断。如果不管调压器，直接让末端电压随负荷波动自由下降，算出来的越限概率会比实际情况严重很多。我后来是先把调压器固定在典型档位，再评估剩余波动，这样结果才合理一些。

7.2 可以扩展的方向

半不变量法这套代码框架其实可以做很多扩展。比如把输入从静态分布换成时间序列预测分布，就能做动态概率潮流；把输出CDF嵌入到配电网网架规划里，就能做考虑不确定性的方案比选；加上储能充放电策略后，还能用来评估储能对电压越限概率的改善效果。我自己的下一步打算是把相关性模型完善起来，让同一馈线下的多个光伏电站出力不再按独立变量处理。

最后分享一个小技巧：代码里把蒙特卡洛校验保留成“回归测试脚本”，每次改动算法或数据后都自动跑一遍，看关键节点的均值误差是否还在阈值内。半不变量法本身很快，配合少量蒙特卡洛校验，既能享受速度优势，又不至于对算法的线性化误差失去感知。这个方法我强烈建议你也试试。