YXS
/
nmWTAI-Platform
8
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

修改损失函数，增加上下，左右，形状，早中晚误差诊断

Browse Source
feature/MultiWellAutoFit-20260805
lvjunjie 2 weeks ago
parent 5539150b5b
commit 6abd1913c7
2 changed files with 460 additions and 118 deletions
Show Stats Download Patch File Download Diff File

64
Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h
Unescape Escape View File

@ -9,6 +9,7 @@
#include <QMutex>
#include <QDateTime>
#include <QFile>
#include <limits>
#include "nmCalculation_global.h"
@ -19,6 +20,67 @@ class nmDataWellBase;
class QTimer;
class QProcess;
// 双对数曲线误差分解，供误差诊断和后续参数调整读取。
//
// 所有 pressure/derivative/shape 数值均是在 log(value) 空间计算的无量纲误差。
// verticalBias* 保留正负号：正值表示模拟曲线整体高于目标，负值表示整体低于目标。
// horizontalPhysicalShift 是 log(time) 方向的等效平移量，正值表示模拟曲线相对目标偏右。
// total 仍是 PSO 当前使用的 fitness（用于粒子比较和收敛判断）；其余字段只描述误差
// 来源，不会在本次改动中直接修改粒子参数，避免损失诊断和参数更新策略互相耦合。
struct AutoFitObjectiveBreakdown {
// valid/total 是本次评价是否有效及其最终 fitness（用于排序和收敛判断）。
bool valid;
double total;
// pressureLoss/derivativeLoss 是压力和导数两条曲线的整体误差。
double pressureLoss;
double derivativeLoss;
// vertical* 描述整体上下偏移；保留 bias 的符号以判断偏高或偏低。
double verticalBiasPressure;
double verticalBiasDerivative;
double verticalLoss;
// horizontal* 描述等效的对数时间偏移；physicalShift 为正表示模拟曲线相对目标向右
// （时间延迟），为负表示向左。
double horizontalShift;
double horizontalPhysicalShift;
double horizontalLoss;
// shapeLoss 是去除整体上下和左右偏移后剩余的曲线形状差异。
double shapeLoss;
// 早、中、晚分段误差用于定位误差主要出现在哪个时间阶段。
double pressureEarlyLoss;
double pressureMiddleLoss;
double pressureLateLoss;
double derivativeEarlyLoss;
double derivativeMiddleLoss;
double derivativeLateLoss;
// coverage 是 50 点目标网格上的 min(有效点比例、连续 log-time 跨度比例)。
// coveragePenalty 是归一化覆盖缺口的平方惩罚，并以 0.1 权重加入 total。
double coverage;
double coveragePenalty;
// 无效评价使用 1e10 作为统一的“差解”标记；其他字段用 NaN 表示尚未得到诊断值。
AutoFitObjectiveBreakdown()
: valid(false)
, total(1.0e10)
, pressureLoss(std::numeric_limits<double>::quiet_NaN())
, derivativeLoss(std::numeric_limits<double>::quiet_NaN())
, verticalBiasPressure(std::numeric_limits<double>::quiet_NaN())
, verticalBiasDerivative(std::numeric_limits<double>::quiet_NaN())
, verticalLoss(std::numeric_limits<double>::quiet_NaN())
, horizontalShift(std::numeric_limits<double>::quiet_NaN())
, horizontalPhysicalShift(std::numeric_limits<double>::quiet_NaN())
, horizontalLoss(std::numeric_limits<double>::quiet_NaN())
, shapeLoss(std::numeric_limits<double>::quiet_NaN())
, pressureEarlyLoss(std::numeric_limits<double>::quiet_NaN())
, pressureMiddleLoss(std::numeric_limits<double>::quiet_NaN())
, pressureLateLoss(std::numeric_limits<double>::quiet_NaN())
, derivativeEarlyLoss(std::numeric_limits<double>::quiet_NaN())
, derivativeMiddleLoss(std::numeric_limits<double>::quiet_NaN())
, derivativeLateLoss(std::numeric_limits<double>::quiet_NaN())
, coverage(std::numeric_limits<double>::quiet_NaN())
, coveragePenalty(std::numeric_limits<double>::quiet_NaN())
{}
};
// PSO粒子结构
// 这里的 position / velocity / bestPosition 只保存“用户勾选参与拟合的参数”，
// 不是完整的 11 个储层/井筒参数。完整参数向量会在写 trace 或调用代理模型时
@ -96,6 +158,7 @@ public:
void stopFitting();
QVector<double> getBestSolution() const;
double getBestFitness() const;
AutoFitObjectiveBreakdown getLastObjectiveBreakdown() const;
QString getLastError() const;
bool isRunning() const;
int getCurrentIteration() const;
@ -288,6 +351,7 @@ private:
double m_previousBestFitness; // 上一轮全局最优误差，用于自适应参数更新。
QVector<QVector<double> > m_lastEvaluatedLogLogData; // 最近一次真实求解得到的 result log-log 曲线。
QVector<QVector<double> > m_globalBestLogLogData; // 当前全局最优对应的 result log-log 曲线。
mutable AutoFitObjectiveBreakdown m_lastObjectiveBreakdown; // 最近一次损失评价的误差分解。
QVector<QVector<double> > m_userInitialLogLogData; // 用户初始解对应的 result log-log 曲线，用于精英保护。
// ===== 优化配置 =====

486
Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp
Unescape Escape View File

@ -791,6 +791,12 @@ double nmCalculationAutoFitPSO::getBestFitness() const
return m_globalBestFitness;
}
AutoFitObjectiveBreakdown nmCalculationAutoFitPSO::getLastObjectiveBreakdown() const
{
// 返回最近一次损失评价的误差分解，供界面或后续优化逻辑读取。
return m_lastObjectiveBreakdown;
}
QString nmCalculationAutoFitPSO::getLastError() const
{
// 上一次失败的人类可读错误信息，主要给 UI 层弹窗或日志使用。
@ -809,6 +815,7 @@ void nmCalculationAutoFitPSO::resetOptimizer()
m_previousBestFitness = 1e10;
m_lastEvaluatedLogLogData.clear();
m_globalBestLogLogData.clear();
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown();
m_userInitialLogLogData.clear();
m_currentIteration = 0;
m_totalEvaluations = 0;
@ -4016,6 +4023,7 @@ double nmCalculationAutoFitPSO::evaluateFitness(const QVector<double>& parameter
static int callCount = 0;
callCount++;
m_lastEvaluatedLogLogData.clear();
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown();
try {
DEBUG_OUT(QString("%1: Call #%2 - Starting evaluation with %3 parameters")
@ -4989,170 +4997,440 @@ double nmCalculationAutoFitPSO::calculateLogLogCurveError(
const QVector<QVector<double> >& target,
const QVector<QVector<double> >& result) const
{
// 双对数曲线误差计算。
//
// target 通常来自目标井历史曲线，result 来自当前粒子参数下的模拟曲线。
// 两条曲线的时间点往往不完全一致，所以这里先取两者时间范围的重叠区间，
// 再在公共时间网格上插值对齐，最后分别计算压力曲线和压力导数曲线误差。
//
// 返回值越小表示拟合越好；返回 1e10 表示曲线无效或无法比较。
// 验证数据
// 这里只负责“曲线比较和误差诊断”，不根据诊断结果直接修改任何拟合参数。
// 调用方可以读取 m_lastObjectiveBreakdown 做诊断或展示；本函数本身不修改参数。
// 在统一的对数时间网格上计算压力和导数残差，并拆分为上下、左右、形状误差。
const double invalidLoss = 1.0e10;
const double valueFloor = 1.0e-12; // 导数接近零时的对数下限，避免 log(0)。
const double huberDelta = qLn(1.2); // 约对应 20% 的相对偏差拐点。
const double minimumCoverage = 0.95; // 点数比例和连续跨度比例都至少接近 95%。
const int numPoints = 50; // 固定网格使不同候选的损失具有可比性。
// 目标函数的主排序项为 0.5*pressureLoss + 0.5*derivativeLoss；
// coveragePenalty 只在接近覆盖边界时提供连续惩罚，上下、左右和形状分量
// 会写入 m_lastObjectiveBreakdown，供后续按误差类型选择参数。
m_lastObjectiveBreakdown = AutoFitObjectiveBreakdown();
if(!validateLogLogData(target) || !validateLogLogData(result)) {
return 1e10;
return invalidLoss;
}
try {
// 数据对齐：找到目标曲线与模拟曲线 time 轴的重叠区域。
// 不在重叠区域内的点不参与误差，避免外推导致误差失真。
double targetMinX = target[0][0];
double targetMaxX = target[0][0];
// 清洗曲线并拆成压力、导数两条曲线。导数可以为负，所以统一使用绝对值
// 进入双对数空间；时间和压力必须为正，否则无法进行对数插值。这里的清洗
// 只丢弃无法比较的采样点，不改变原始曲线或求解器输出。
auto prepareCurve = [valueFloor](const QVector<QVector<double> >& data,
QVector<QPointF>* pressure,
QVector<QPointF>* derivative) -> bool {
if(!pressure || !derivative || data.size() < 3 ||
data[0].size() != data[1].size() ||
data[0].size() != data[2].size()) {
return false;
}
for(int i = 1; i < target[0].size(); ++i) {
if(target[0][i] < targetMinX) targetMinX = target[0][i];
for(int i = 0; i < data[0].size(); ++i) {
if(!isFiniteNumber(data[0][i]) || !isFiniteNumber(data[1][i]) ||
!isFiniteNumber(data[2][i]) || data[0][i] <= 0.0 ||
data[1][i] <= 0.0) {
continue;
}
if(target[0][i] > targetMaxX) targetMaxX = target[0][i];
pressure->append(QPointF(data[0][i], data[1][i]));
derivative->append(QPointF(data[0][i],
qMax(qAbs(data[2][i]), valueFloor)));
}
double resultMinX = result[0][0];
double resultMaxX = result[0][0];
// 插值要求时间严格递增。重复时间点保留排序后的最后一个值，
// 避免重复横坐标导致对数插值分母为零。压力和导数分别去重，
// 这样即使某条曲线存在重复时间点，也不会污染另一条曲线的插值。
auto sortAndUnique = [](QVector<QPointF>* curve) {
std::stable_sort(curve->begin(), curve->end(),
[](const QPointF& left, const QPointF& right) {
return left.x() < right.x();
});
for(int i = 1; i < result[0].size(); ++i) {
if(result[0][i] < resultMinX) resultMinX = result[0][i];
QVector<QPointF> unique;
unique.reserve(curve->size());
if(result[0][i] > resultMaxX) resultMaxX = result[0][i];
for(int i = 0; i < curve->size(); ++i) {
if(unique.isEmpty() || curve->at(i).x() > unique.last().x()) {
unique.append(curve->at(i));
} else {
unique[unique.size() - 1] = curve->at(i);
}
}
double overlapMinX = qMax(targetMinX, resultMinX);
double overlapMaxX = qMin(targetMaxX, resultMaxX);
*curve = unique;
};
if(overlapMinX >= overlapMaxX) {
DEBUG_OUT("No overlap between target and result LogLog curves");
return 1e10;
sortAndUnique(pressure);
sortAndUnique(derivative);
return pressure->size() >= 3 && derivative->size() >= 3;
};
QVector<QPointF> targetPressure;
QVector<QPointF> targetDerivative;
QVector<QPointF> resultPressure;
QVector<QPointF> resultDerivative;
if(!prepareCurve(target, &targetPressure, &targetDerivative) ||
!prepareCurve(result, &resultPressure, &resultDerivative)) {
return invalidLoss;
}
// 生成公共 X 网格进行插值。使用对数均匀网格，是为了给早期时间段
// 更多分辨率；试井双对数曲线的早期形态通常对参数识别很敏感。
QVector<double> commonX;
int numPoints = 50;
// 在 log(time)-log(value) 空间做线性插值，而不是在线性坐标直接插值。
// 这样可以保持双对数曲线的时间尺度和数量级特征；返回值是
// log(abs(value))，后续残差因此可以直接解释为相对幅值误差。
auto interpolateLogValue = [valueFloor](const QVector<QPointF>& curve,
double x,
double* value) -> bool {
if(!value || curve.size() < 2 || x < curve.first().x() ||
x > curve.last().x() || x <= 0.0) {
return false;
}
int right = 1;
while(right < curve.size() && curve[right].x() < x) {
++right;
}
right = qMin(right, curve.size() - 1);
int left = qMax(0, right - 1);
double leftLogX = qLn(curve[left].x());
double rightLogX = qLn(curve[right].x());
double denominator = rightLogX - leftLogX;
double leftLogY = qLn(qMax(qAbs(curve[left].y()), valueFloor));
double rightLogY = qLn(qMax(qAbs(curve[right].y()), valueFloor));
if(qAbs(denominator) <= 1.0e-12) {
*value = leftLogY;
} else {
double ratio = (qLn(x) - leftLogX) / denominator;
*value = leftLogY + ratio * (rightLogY - leftLogY);
}
return isFiniteNumber(*value);
};
if(overlapMinX > 0 && overlapMaxX > 0) {
// 对数空间均匀分布
double logMin = qLn(overlapMinX);
double logMax = qLn(overlapMaxX);
const double targetMinX = targetPressure.first().x();
const double targetMaxX = targetPressure.last().x();
const double resultMinX = resultPressure.first().x();
const double resultMaxX = resultPressure.last().x();
if(targetMinX <= 0.0 || targetMaxX <= targetMinX ||
resultMinX <= 0.0 || resultMaxX <= resultMinX) {
return invalidLoss;
}
// 目标曲线的完整时间范围作为统一比较区间。若候选曲线覆盖不足，
// 后面的 coverage 检查会拒绝它，防止候选通过缩短时间范围来降低误差。
QVector<double> commonX(numPoints);
QVector<double> commonLogX(numPoints);
QVector<double> targetLogPressure(numPoints);
QVector<double> targetLogDerivative(numPoints);
const double targetLogMinX = qLn(targetMinX);
const double targetLogMaxX = qLn(targetMaxX);
for(int i = 0; i < numPoints; ++i) {
double logX = logMin + i * (logMax - logMin) / (numPoints - 1);
double x = qExp(logX);
double logX = targetLogMinX +
static_cast<double>(i) *
(targetLogMaxX - targetLogMinX) / (numPoints - 1);
commonLogX[i] = logX;
commonX[i] = qExp(logX);
if(!interpolateLogValue(targetPressure, commonX[i],
&targetLogPressure[i]) ||
!interpolateLogValue(targetDerivative, commonX[i],
&targetLogDerivative[i])) {
return invalidLoss;
}
}
// 残差采用“模拟减目标”，因此正值表示模拟曲线在对数幅值上高于目标，
// 负值表示模拟曲线偏低。NaN 表示该网格点不在模拟曲线支持范围内，
// 后续统计会自动跳过，但覆盖率检查仍会限制候选不能靠缺失数据降低损失。
QVector<double> pressureResidual(numPoints,
std::numeric_limits<double>::quiet_NaN());
QVector<double> derivativeResidual(numPoints,
std::numeric_limits<double>::quiet_NaN());
QVector<double> pressureSlope(numPoints, 0.0);
QVector<double> derivativeSlope(numPoints, 0.0);
int firstSupported = -1;
int lastSupported = -1;
int supportedCount = 0;
// 数值保护
if(!isFiniteNumber(x) || x <= 0) {
for(int i = 0; i < numPoints; ++i) {
if(commonX[i] < resultMinX || commonX[i] > resultMaxX) {
continue;
}
commonX.append(x);
double resultLogPressure = 0.0;
double resultLogDerivative = 0.0;
if(!interpolateLogValue(resultPressure, commonX[i],
&resultLogPressure) ||
!interpolateLogValue(resultDerivative, commonX[i],
&resultLogDerivative)) {
continue;
}
DEBUG_OUT("Using log-uniform grid for better early-time coverage");
pressureResidual[i] = resultLogPressure - targetLogPressure[i];
derivativeResidual[i] = resultLogDerivative - targetLogDerivative[i];
++supportedCount;
if(firstSupported < 0) {
firstSupported = i;
}
if(commonX.isEmpty()) {
DEBUG_OUT("Failed to generate common X grid");
return 1e10;
lastSupported = i;
}
// 同时使用点覆盖率和连续时间跨度覆盖率，避免只覆盖少数离散点也被判定为完整。
AutoFitObjectiveBreakdown breakdown;
breakdown.coverage = supportedCount > 0
? static_cast<double>(supportedCount) / numPoints
: 0.0;
if(firstSupported >= 0 && lastSupported >= firstSupported) {
double span = qMax(1.0e-12, targetLogMaxX - targetLogMinX);
double coveredSpan = commonLogX[lastSupported] -
commonLogX[firstSupported];
breakdown.coverage = qMin(breakdown.coverage,
qMax(0.0, coveredSpan / span));
}
// 插值目标曲线。target[1] 是压力，target[2] 是压力导数。
QVector<QPointF> targetCurve1, targetCurve2;
// 覆盖率越接近 1，惩罚越小；覆盖不足 minimumCoverage 时直接返回无效损失。
double coverageGap = qMax(0.0, 1.0 - breakdown.coverage);
breakdown.coveragePenalty =
qPow(coverageGap / (1.0 - minimumCoverage), 2.0);
for(int i = 0; i < target[0].size(); ++i) {
// 检查数据有效性
if(isFiniteNumber(target[0][i]) && isFiniteNumber(target[1][i]) &&
isFiniteNumber(target[2][i])) {
targetCurve1.append(QPointF(target[0][i], target[1][i]));
targetCurve2.append(QPointF(target[0][i], target[2][i]));
if(breakdown.coverage < minimumCoverage) {
breakdown.total = invalidLoss;
m_lastObjectiveBreakdown = breakdown;
return invalidLoss;
}
// 目标曲线斜率用于把“残差随时间的系统性变化”解释为左右平移。
// 斜率用对数坐标计算，与前面的插值空间保持一致；平坦区斜率接近零，
// 不会凭空制造水平偏移量。
for(int i = 0; i < numPoints; ++i) {
int left = i == 0 ? 0 : i - 1;
int right = i == numPoints - 1 ? numPoints - 1 : i + 1;
double denominator = commonLogX[right] - commonLogX[left];
if(qAbs(denominator) > 1.0e-12) {
pressureSlope[i] =
(targetLogPressure[right] - targetLogPressure[left]) /
denominator;
derivativeSlope[i] =
(targetLogDerivative[right] - targetLogDerivative[left]) /
denominator;
}
}
// Huber RMS 在小残差区域保持平方损失，在异常点区域转为线性增长，
// 避免少量求解器异常点完全主导候选排序。这里没有除以目标值，
// 因为残差已经是 log(value) 差值，本身就是相对误差的表达；返回值是
// Huber rho 均值的平方根，保持与 RMS 类似的尺度。
auto huberRms = [huberDelta](const QVector<double>& values,
int begin,
int end) -> double {
double sum = 0.0;
int count = 0;
for(int i = qMax(0, begin);
i < qMin(end, static_cast<int>(values.size())); ++i) {
if(!isFiniteNumber(values[i])) {
continue;
}
if(targetCurve1.isEmpty() || targetCurve2.isEmpty()) {
DEBUG_OUT("Target curves are empty after filtering");
return 1e10;
double absoluteValue = qAbs(values[i]);
double rho = absoluteValue <= huberDelta
? values[i] * values[i]
: 2.0 * huberDelta * absoluteValue -
huberDelta * huberDelta;
sum += rho;
++count;
}
QVector<QPointF> alignedTarget1 = interpolateData(targetCurve1, commonX);
QVector<QPointF> alignedTarget2 = interpolateData(targetCurve2, commonX);
return count > 0
? qSqrt(sum / count)
: std::numeric_limits<double>::quiet_NaN();
};
// 插值结果曲线。result 与 target 使用同一 commonX，保证逐点可比。
QVector<QPointF> resultCurve1, resultCurve2;
// 用 Huber 加权迭代估计残差中心，作为整体上下偏移。相比普通平均值，
// 它对局部尖峰更稳健，同时保留偏高/偏低的方向信息。迭代只用于诊断，
// 不会把残差“校正”后再写回求解器结果。
auto huberCenter = [huberDelta](const QVector<double>& values) -> double {
double center = 0.0;
int count = 0;
for(int i = 0; i < result[0].size(); ++i) {
// 检查数据有效性
if(isFiniteNumber(result[0][i]) && isFiniteNumber(result[1][i]) &&
isFiniteNumber(result[2][i])) {
resultCurve1.append(QPointF(result[0][i], result[1][i]));
resultCurve2.append(QPointF(result[0][i], result[2][i]));
for(int i = 0; i < values.size(); ++i) {
if(isFiniteNumber(values[i])) {
center += values[i];
++count;
}
}
if(resultCurve1.isEmpty() || resultCurve2.isEmpty()) {
DEBUG_OUT("Result curves are empty after filtering");
return 1e10;
if(count == 0) {
return std::numeric_limits<double>::quiet_NaN();
}
QVector<QPointF> alignedResult1 = interpolateData(resultCurve1, commonX);
QVector<QPointF> alignedResult2 = interpolateData(resultCurve2, commonX);
center /= count;
// 检查插值结果
if(alignedTarget1.isEmpty() || alignedTarget2.isEmpty() ||
alignedResult1.isEmpty() || alignedResult2.isEmpty()) {
DEBUG_OUT("LogLog interpolation failed");
return 1e10;
// 固定最多 8 次迭代，控制每个候选的计算开销并保持结果稳定。
for(int iteration = 0; iteration < 8; ++iteration) {
double weightedSum = 0.0;
double weightTotal = 0.0;
for(int i = 0; i < values.size(); ++i) {
if(!isFiniteNumber(values[i])) {
continue;
}
if(alignedTarget1.size() != alignedResult1.size() ||
alignedTarget2.size() != alignedResult2.size()) {
DEBUG_OUT("LogLog interpolation size mismatch");
return 1e10;
double distance = qAbs(values[i] - center);
// 距离接近零时直接取权重 1，避免除零并保持中心点不被放大。
double weight = distance <= huberDelta || distance < 1.0e-12
? 1.0
: huberDelta / distance;
weightedSum += weight * values[i];
weightTotal += weight;
}
// 计算两条曲线的误差。当前压力和导数各占 50%。
// 如果后续要让导数形态更重要，可以从这里调整权重。
double error1 = calculateCurveError(alignedTarget1, alignedResult1);
double error2 = calculateCurveError(alignedTarget2, alignedResult2);
if(weightTotal <= 1.0e-12) {
break;
}
// 检查个别误差是否有效
if(!isFiniteNumber(error1) || error1 > 1e9) {
DEBUG_OUT(QString("Curve1 error is invalid: %1").arg(error1));
error1 = 1e10;
double nextCenter = weightedSum / weightTotal;
if(qAbs(nextCenter - center) <= 1.0e-12) {
center = nextCenter;
break;
}
if(!isFiniteNumber(error2) || error2 > 1e9) {
DEBUG_OUT(QString("Curve2 error is invalid: %1").arg(error2));
error2 = 1e10;
center = nextCenter;
}
// 组合误差 - 添加保护
double combinedError;
return center;
};
if(error1 > 1e9 && error2 > 1e9) {
combinedError = 1e10;
} else if(error1 > 1e9) {
combinedError = error2;
} else if(error2 > 1e9) {
combinedError = error1;
} else {
combinedError = 0.5 * error1 + 0.5 * error2;
// 整体压力/导数误差用于排序；verticalLoss 主要用于解释整体上下偏移。
breakdown.pressureLoss = huberRms(pressureResidual, 0, numPoints);
breakdown.derivativeLoss = huberRms(derivativeResidual, 0, numPoints);
breakdown.verticalBiasPressure = huberCenter(pressureResidual);
breakdown.verticalBiasDerivative = huberCenter(derivativeResidual);
breakdown.verticalLoss =
0.5 * (qAbs(breakdown.verticalBiasPressure) +
qAbs(breakdown.verticalBiasDerivative));
// 去除上下中心后，把残差投影到目标曲线斜率上估计左右偏移。
// residual ~= -physicalShift * targetSlope，因此 physicalShift 取回归系数的相反数。
// 这是局部一阶近似，用于判断方向和大小，不等同于再次优化时间轴。
double horizontalNumerator = 0.0;
double horizontalDenominator = 0.0;
for(int i = 0; i < numPoints; ++i) {
if(!isFiniteNumber(pressureResidual[i]) ||
!isFiniteNumber(derivativeResidual[i])) {
continue;
}
DEBUG_OUT(QString("LogLog errors: Curve1=%1, Curve2=%2, Combined=%3")
.arg(error1, 0, 'e', 4).arg(error2, 0, 'e', 4).arg(combinedError, 0, 'e', 4));
double pressureCentered =
pressureResidual[i] - breakdown.verticalBiasPressure;
double derivativeCentered =
derivativeResidual[i] - breakdown.verticalBiasDerivative;
horizontalNumerator += pressureSlope[i] * pressureCentered +
derivativeSlope[i] * derivativeCentered;
horizontalDenominator += pressureSlope[i] * pressureSlope[i] +
derivativeSlope[i] * derivativeSlope[i];
}
return qMin(1e9, combinedError);
breakdown.horizontalShift = horizontalDenominator > 1.0e-12
? horizontalNumerator /
horizontalDenominator
: 0.0;
breakdown.horizontalPhysicalShift = -breakdown.horizontalShift;
breakdown.horizontalLoss = qAbs(breakdown.horizontalShift);
// 从原始残差中扣除“整体上下 + 等效左右”两部分，剩余项才作为形状误差。
// 因此 shapeLoss 较大而 vertical/horizontal 较小时，说明主要是曲率、拐点
// 或导数变化趋势不一致，而不是简单的整体平移。
QVector<double> shapePressure(numPoints,
std::numeric_limits<double>::quiet_NaN());
QVector<double> shapeDerivative(numPoints,
std::numeric_limits<double>::quiet_NaN());
for(int i = 0; i < numPoints; ++i) {
if(isFiniteNumber(pressureResidual[i])) {
shapePressure[i] = pressureResidual[i] -
breakdown.verticalBiasPressure -
breakdown.horizontalShift * pressureSlope[i];
}
if(isFiniteNumber(derivativeResidual[i])) {
shapeDerivative[i] = derivativeResidual[i] -
breakdown.verticalBiasDerivative -
breakdown.horizontalShift *
derivativeSlope[i];
}
}
breakdown.shapeLoss =
0.5 * (huberRms(shapePressure, 0, numPoints) +
huberRms(shapeDerivative, 0, numPoints));
// 将对数时间网格分成早、中、晚三段，用于定位误差集中出现的阶段。
// 网格本身按 log(time) 均匀分布，所以三段对应的是时间数量级，而非原始
// 线性时间长度，适合双对数试井曲线的早期/中期/晚期判读。
const int segment1 = numPoints / 3;
const int segment2 = (2 * numPoints) / 3;
breakdown.pressureEarlyLoss = huberRms(pressureResidual, 0, segment1);
breakdown.pressureMiddleLoss =
huberRms(pressureResidual, segment1, segment2);
breakdown.pressureLateLoss =
huberRms(pressureResidual, segment2, numPoints);
breakdown.derivativeEarlyLoss =
huberRms(derivativeResidual, 0, segment1);
breakdown.derivativeMiddleLoss =
huberRms(derivativeResidual, segment1, segment2);
breakdown.derivativeLateLoss =
huberRms(derivativeResidual, segment2, numPoints);
// 函数入口已将 m_lastObjectiveBreakdown 重置为无效状态，因此这里直接返回
// invalidLoss 时不会把上一候选的误差分解误报给调用方。
// 导数、压力或形状无法形成有效统计时，整个候选都视为无效，避免 NaN
// 进入粒子排序。
if(!isFiniteNumber(breakdown.pressureLoss) ||
!isFiniteNumber(breakdown.derivativeLoss) ||
!isFiniteNumber(breakdown.shapeLoss) ||
!isFiniteNumber(breakdown.verticalLoss)) {
return invalidLoss;
}
// 当前总损失作为 fitness 用于粒子比较、收敛/停止判断；上下、左右、形状和
// 分段分量先作为诊断输出，不在本次改动中直接参与参数更新。
breakdown.total = 0.5 * breakdown.pressureLoss +
0.5 * breakdown.derivativeLoss +
0.1 * breakdown.coveragePenalty;
breakdown.valid = isFiniteNumber(breakdown.total) &&
breakdown.total >= 0.0;
m_lastObjectiveBreakdown = breakdown;
DEBUG_OUT(QString("LogLog objective: pressure=%1, derivative=%2, vertical=%3, horizontal=%4, shape=%5, coverage=%6, total=%7")
.arg(breakdown.pressureLoss, 0, 'e', 4)
.arg(breakdown.derivativeLoss, 0, 'e', 4)
.arg(breakdown.verticalLoss, 0, 'e', 4)
.arg(breakdown.horizontalLoss, 0, 'e', 4)
.arg(breakdown.shapeLoss, 0, 'e', 4)
.arg(breakdown.coverage, 0, 'f', 4)
.arg(breakdown.total, 0, 'e', 4));
return breakdown.valid ? qMin(1.0e9, breakdown.total) : invalidLoss;
} catch(const std::exception& e) {
DEBUG_OUT(QString("Exception in LogLog error calculation: %1").arg(e.what()));
return 1e10;
return invalidLoss;
} catch(...) {
DEBUG_OUT("Unknown exception in LogLog error calculation");
return 1e10;
return invalidLoss;
}
}

Write Preview
Loading…
Cancel
Save