diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h index 3418b04..e65d1ee 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitPSO.h @@ -9,6 +9,7 @@ #include #include #include +#include #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::quiet_NaN()) + , derivativeLoss(std::numeric_limits::quiet_NaN()) + , verticalBiasPressure(std::numeric_limits::quiet_NaN()) + , verticalBiasDerivative(std::numeric_limits::quiet_NaN()) + , verticalLoss(std::numeric_limits::quiet_NaN()) + , horizontalShift(std::numeric_limits::quiet_NaN()) + , horizontalPhysicalShift(std::numeric_limits::quiet_NaN()) + , horizontalLoss(std::numeric_limits::quiet_NaN()) + , shapeLoss(std::numeric_limits::quiet_NaN()) + , pressureEarlyLoss(std::numeric_limits::quiet_NaN()) + , pressureMiddleLoss(std::numeric_limits::quiet_NaN()) + , pressureLateLoss(std::numeric_limits::quiet_NaN()) + , derivativeEarlyLoss(std::numeric_limits::quiet_NaN()) + , derivativeMiddleLoss(std::numeric_limits::quiet_NaN()) + , derivativeLateLoss(std::numeric_limits::quiet_NaN()) + , coverage(std::numeric_limits::quiet_NaN()) + , coveragePenalty(std::numeric_limits::quiet_NaN()) + {} +}; + // PSO粒子结构 // 这里的 position / velocity / bestPosition 只保存“用户勾选参与拟合的参数”， // 不是完整的 11 个储层/井筒参数。完整参数向量会在写 trace 或调用代理模型时 @@ -96,6 +158,7 @@ public: void stopFitting(); QVector 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 > m_lastEvaluatedLogLogData; // 最近一次真实求解得到的 result log-log 曲线。 QVector > m_globalBestLogLogData; // 当前全局最优对应的 result log-log 曲线。 + mutable AutoFitObjectiveBreakdown m_lastObjectiveBreakdown; // 最近一次损失评价的误差分解。 QVector > m_userInitialLogLogData; // 用户初始解对应的 result log-log 曲线，用于精英保护。 // ===== 优化配置 ===== diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp index 936d9fc..ccc05ea 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitPSO.cpp @@ -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& 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 >& target, const QVector >& 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 >& data, + QVector* pressure, + QVector* 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* curve) { + std::stable_sort(curve->begin(), curve->end(), + [](const QPointF& left, const QPointF& right) { + return left.x() < right.x(); + }); + + QVector unique; + unique.reserve(curve->size()); + + 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); + } + } - for(int i = 1; i < result[0].size(); ++i) { - if(result[0][i] < resultMinX) resultMinX = result[0][i]; + *curve = unique; + }; - if(result[0][i] > resultMaxX) resultMaxX = result[0][i]; - } + sortAndUnique(pressure); + sortAndUnique(derivative); + return pressure->size() >= 3 && derivative->size() >= 3; + }; - double overlapMinX = qMax(targetMinX, resultMinX); - double overlapMaxX = qMin(targetMaxX, resultMaxX); + QVector targetPressure; + QVector targetDerivative; + QVector resultPressure; + QVector resultDerivative; - if(overlapMinX >= overlapMaxX) { - DEBUG_OUT("No overlap between target and result LogLog curves"); - return 1e10; + if(!prepareCurve(target, &targetPressure, &targetDerivative) || + !prepareCurve(result, &resultPressure, &resultDerivative)) { + return invalidLoss; } - // 生成公共 X 网格进行插值。使用对数均匀网格，是为了给早期时间段 - // 更多分辨率；试井双对数曲线的早期形态通常对参数识别很敏感。 - QVector commonX; - int numPoints = 50; + // 在 log(time)-log(value) 空间做线性插值，而不是在线性坐标直接插值。 + // 这样可以保持双对数曲线的时间尺度和数量级特征；返回值是 + // log(abs(value))，后续残差因此可以直接解释为相对幅值误差。 + auto interpolateLogValue = [valueFloor](const QVector& curve, + double x, + double* value) -> bool { + if(!value || curve.size() < 2 || x < curve.first().x() || + x > curve.last().x() || x <= 0.0) { + return false; + } - if(overlapMinX > 0 && overlapMaxX > 0) { - // 对数空间均匀分布 - double logMin = qLn(overlapMinX); - double logMax = qLn(overlapMaxX); + int right = 1; - for(int i = 0; i < numPoints; ++i) { - double logX = logMin + i * (logMax - logMin) / (numPoints - 1); - double x = qExp(logX); + while(right < curve.size() && curve[right].x() < x) { + ++right; + } - // 数值保护 - if(!isFiniteNumber(x) || x <= 0) { - continue; - } + 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)); - commonX.append(x); + if(qAbs(denominator) <= 1.0e-12) { + *value = leftLogY; + } else { + double ratio = (qLn(x) - leftLogX) / denominator; + *value = leftLogY + ratio * (rightLogY - leftLogY); } - DEBUG_OUT("Using log-uniform grid for better early-time coverage"); + return isFiniteNumber(*value); + }; + + 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 commonX(numPoints); + QVector commonLogX(numPoints); + QVector targetLogPressure(numPoints); + QVector targetLogDerivative(numPoints); + const double targetLogMinX = qLn(targetMinX); + const double targetLogMaxX = qLn(targetMaxX); + + for(int i = 0; i < numPoints; ++i) { + double logX = targetLogMinX + + static_cast(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; + } } - if(commonX.isEmpty()) { - DEBUG_OUT("Failed to generate common X grid"); - return 1e10; - } + // 残差采用“模拟减目标”，因此正值表示模拟曲线在对数幅值上高于目标， + // 负值表示模拟曲线偏低。NaN 表示该网格点不在模拟曲线支持范围内， + // 后续统计会自动跳过，但覆盖率检查仍会限制候选不能靠缺失数据降低损失。 + QVector pressureResidual(numPoints, + std::numeric_limits::quiet_NaN()); + QVector derivativeResidual(numPoints, + std::numeric_limits::quiet_NaN()); + QVector pressureSlope(numPoints, 0.0); + QVector derivativeSlope(numPoints, 0.0); + int firstSupported = -1; + int lastSupported = -1; + int supportedCount = 0; + + for(int i = 0; i < numPoints; ++i) { + if(commonX[i] < resultMinX || commonX[i] > resultMaxX) { + continue; + } - // 插值目标曲线。target[1] 是压力，target[2] 是压力导数。 - QVector targetCurve1, targetCurve2; + double resultLogPressure = 0.0; + double resultLogDerivative = 0.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(!interpolateLogValue(resultPressure, commonX[i], + &resultLogPressure) || + !interpolateLogValue(resultDerivative, commonX[i], + &resultLogDerivative)) { + continue; } - } - if(targetCurve1.isEmpty() || targetCurve2.isEmpty()) { - DEBUG_OUT("Target curves are empty after filtering"); - return 1e10; + pressureResidual[i] = resultLogPressure - targetLogPressure[i]; + derivativeResidual[i] = resultLogDerivative - targetLogDerivative[i]; + ++supportedCount; + + if(firstSupported < 0) { + firstSupported = i; + } + + lastSupported = i; + } + + // 同时使用点覆盖率和连续时间跨度覆盖率，避免只覆盖少数离散点也被判定为完整。 + AutoFitObjectiveBreakdown breakdown; + breakdown.coverage = supportedCount > 0 + ? static_cast(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)); + } + + // 覆盖率越接近 1，惩罚越小；覆盖不足 minimumCoverage 时直接返回无效损失。 + double coverageGap = qMax(0.0, 1.0 - breakdown.coverage); + breakdown.coveragePenalty = + qPow(coverageGap / (1.0 - minimumCoverage), 2.0); + + 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; + } } - QVector alignedTarget1 = interpolateData(targetCurve1, commonX); - QVector alignedTarget2 = interpolateData(targetCurve2, commonX); + // Huber RMS 在小残差区域保持平方损失，在异常点区域转为线性增长， + // 避免少量求解器异常点完全主导候选排序。这里没有除以目标值， + // 因为残差已经是 log(value) 差值，本身就是相对误差的表达；返回值是 + // Huber rho 均值的平方根，保持与 RMS 类似的尺度。 + auto huberRms = [huberDelta](const QVector& values, + int begin, + int end) -> double { + double sum = 0.0; + int count = 0; - // 插值结果曲线。result 与 target 使用同一 commonX，保证逐点可比。 - QVector resultCurve1, resultCurve2; + for(int i = qMax(0, begin); + i < qMin(end, static_cast(values.size())); ++i) { + if(!isFiniteNumber(values[i])) { + continue; + } - 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])); + double absoluteValue = qAbs(values[i]); + double rho = absoluteValue <= huberDelta + ? values[i] * values[i] + : 2.0 * huberDelta * absoluteValue - + huberDelta * huberDelta; + sum += rho; + ++count; } - } - if(resultCurve1.isEmpty() || resultCurve2.isEmpty()) { - DEBUG_OUT("Result curves are empty after filtering"); - return 1e10; - } + return count > 0 + ? qSqrt(sum / count) + : std::numeric_limits::quiet_NaN(); + }; + + // 用 Huber 加权迭代估计残差中心，作为整体上下偏移。相比普通平均值， + // 它对局部尖峰更稳健，同时保留偏高/偏低的方向信息。迭代只用于诊断， + // 不会把残差“校正”后再写回求解器结果。 + auto huberCenter = [huberDelta](const QVector& values) -> double { + double center = 0.0; + int count = 0; + + for(int i = 0; i < values.size(); ++i) { + if(isFiniteNumber(values[i])) { + center += values[i]; + ++count; + } + } - QVector alignedResult1 = interpolateData(resultCurve1, commonX); - QVector alignedResult2 = interpolateData(resultCurve2, commonX); + if(count == 0) { + return std::numeric_limits::quiet_NaN(); + } - // 检查插值结果 - if(alignedTarget1.isEmpty() || alignedTarget2.isEmpty() || - alignedResult1.isEmpty() || alignedResult2.isEmpty()) { - DEBUG_OUT("LogLog interpolation failed"); - return 1e10; - } + center /= count; - if(alignedTarget1.size() != alignedResult1.size() || - alignedTarget2.size() != alignedResult2.size()) { - DEBUG_OUT("LogLog interpolation size mismatch"); - return 1e10; - } + // 固定最多 8 次迭代，控制每个候选的计算开销并保持结果稳定。 + for(int iteration = 0; iteration < 8; ++iteration) { + double weightedSum = 0.0; + double weightTotal = 0.0; - // 计算两条曲线的误差。当前压力和导数各占 50%。 - // 如果后续要让导数形态更重要，可以从这里调整权重。 - double error1 = calculateCurveError(alignedTarget1, alignedResult1); - double error2 = calculateCurveError(alignedTarget2, alignedResult2); + for(int i = 0; i < values.size(); ++i) { + if(!isFiniteNumber(values[i])) { + continue; + } - // 检查个别误差是否有效 - if(!isFiniteNumber(error1) || error1 > 1e9) { - DEBUG_OUT(QString("Curve1 error is invalid: %1").arg(error1)); - error1 = 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; + } - if(!isFiniteNumber(error2) || error2 > 1e9) { - DEBUG_OUT(QString("Curve2 error is invalid: %1").arg(error2)); - error2 = 1e10; - } + if(weightTotal <= 1.0e-12) { + break; + } - // 组合误差 - 添加保护 - double combinedError; + double nextCenter = weightedSum / weightTotal; + if(qAbs(nextCenter - center) <= 1.0e-12) { + center = nextCenter; + break; + } - 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; - } + center = nextCenter; + } - 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)); + return center; + }; + + // 整体压力/导数误差用于排序；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; + } + + 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]; + } + + breakdown.horizontalShift = horizontalDenominator > 1.0e-12 + ? horizontalNumerator / + horizontalDenominator + : 0.0; + breakdown.horizontalPhysicalShift = -breakdown.horizontalShift; + breakdown.horizontalLoss = qAbs(breakdown.horizontalShift); + + // 从原始残差中扣除“整体上下 + 等效左右”两部分，剩余项才作为形状误差。 + // 因此 shapeLoss 较大而 vertical/horizontal 较小时，说明主要是曲率、拐点 + // 或导数变化趋势不一致，而不是简单的整体平移。 + QVector shapePressure(numPoints, + std::numeric_limits::quiet_NaN()); + QVector shapeDerivative(numPoints, + std::numeric_limits::quiet_NaN()); + + for(int i = 0; i < numPoints; ++i) { + if(isFiniteNumber(pressureResidual[i])) { + shapePressure[i] = pressureResidual[i] - + breakdown.verticalBiasPressure - + breakdown.horizontalShift * pressureSlope[i]; + } - return qMin(1e9, combinedError); + 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; } }