From fdb402561da56b49ded15d106bfc6e1bea87e8e0 Mon Sep 17 00:00:00 2001 From: lvjunjie Date: Thu, 27 Aug 2026 11:08:07 +0800 Subject: [PATCH] =?UTF-8?q?=E6=A0=B9=E6=8D=AE=E7=9B=AE=E6=A0=87=E6=9B=B2?= =?UTF-8?q?=E7=BA=BF=E5=88=92=E5=88=86=E6=97=B6=E9=97=B4=E6=AE=B5=EF=BC=8C?= =?UTF-8?q?=E6=A0=B9=E6=8D=AE=E6=97=B6=E9=97=B4=E6=AE=B5=E7=9A=84=E8=AF=AF?= =?UTF-8?q?=E5=B7=AE=E8=B0=83=E6=95=B4=E5=8F=82=E6=95=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../nmCalculation/nmCalculationAutoFitLM.h | 11 +- .../nmCalculation/nmCalculationAutoFitLM.cpp | 643 +++++++++++------- 2 files changed, 392 insertions(+), 262 deletions(-) diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h index 71e22de..3cbbe88 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h @@ -13,14 +13,21 @@ class nmDataWellBase; -// 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据， -// 其余诊断量用于有限差分灵敏度分析和信赖域选参。 +// 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据；时间窗口误差 +// 用于细搜选参，原有上下、左右和形状量仅保留为过程诊断。 struct AutoFitObjectiveBreakdownLM { bool valid; double total; double pressureLoss; double derivativeLoss; QVector residualVector; + QVector timeWindowBegin; + QVector timeWindowEnd; + QVector timeWindowStart; + QVector timeWindowEndTime; + QVector timeWindowPressureLoss; + QVector timeWindowDerivativeLoss; + QVector timeWindowLoss; double verticalCommonBias; double verticalLoss; bool verticalReliable; diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp index e0a15e6..571097a 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp @@ -169,12 +169,14 @@ static double fromTrustRegionCoordinate(double coordinate, return lower + coordinate * (upper - lower); } -enum TrustRegionErrorComponent +enum TrustRegionTimeWindow { - TRUST_REGION_VERTICAL_COMPONENT = 0, - TRUST_REGION_HORIZONTAL_COMPONENT, - TRUST_REGION_SHAPE_COMPONENT, - TRUST_REGION_TOTAL_COMPONENT + TRUST_REGION_WELLBORE_STORAGE_WINDOW = 0, + TRUST_REGION_TRANSITION_WINDOW, + TRUST_REGION_IARF_WINDOW, + TRUST_REGION_LATE_WINDOW, + TRUST_REGION_TIME_WINDOW_COUNT, + TRUST_REGION_TOTAL_WINDOW }; // 一次真实求解的完整快照。除了参数和总误差，还保存内部坐标、诊断分量和 @@ -202,7 +204,13 @@ static bool trustRegionResidualsValid( { // 损失函数固定使用 80 个压力点和 80 个导数点。严格校验长度，避免 // Jacobian 沿用旧维度后访问另一候选的短残差向量。 - if(!breakdown.valid || breakdown.residualVector.size() != 160) { + if(!breakdown.valid || breakdown.residualVector.size() != 160 || + breakdown.timeWindowBegin.size() != + TRUST_REGION_TIME_WINDOW_COUNT || + breakdown.timeWindowEnd.size() != + TRUST_REGION_TIME_WINDOW_COUNT || + breakdown.timeWindowLoss.size() != + TRUST_REGION_TIME_WINDOW_COUNT) { return false; } @@ -211,6 +219,17 @@ static bool trustRegionResidualsValid( return false; } } + for(int window = 0; + window < TRUST_REGION_TIME_WINDOW_COUNT; + ++window) { + if(breakdown.timeWindowBegin[window] < 0 || + breakdown.timeWindowEnd[window] <= + breakdown.timeWindowBegin[window] || + breakdown.timeWindowEnd[window] > 80 || + !isFiniteNumber(breakdown.timeWindowLoss[window])) { + return false; + } + } return true; } @@ -224,68 +243,167 @@ static double trustRegionSquaredNorm(const QVector& values) return sum; } -// 计算同维向量内积；维度不一致表示局部模型无效，返回零让调用方放弃修正。 -static double trustRegionDotProduct(const QVector& left, - const QVector& right) +// trace 和运行日志使用稳定的英文标识，便于离线过程分析按窗口筛选。 +static QString trustRegionWindowName(int window) { - if(left.size() != right.size()) { - return 0.0; + if(window == TRUST_REGION_WELLBORE_STORAGE_WINDOW) { + return "w0_storage"; } - - double sum = 0.0; - for(int i = 0; i < left.size(); ++i) { - sum += left[i] * right[i]; + if(window == TRUST_REGION_TRANSITION_WINDOW) { + return "w1_transition"; } - return sum; + if(window == TRUST_REGION_IARF_WINDOW) { + return "w2_iarf"; + } + if(window == TRUST_REGION_LATE_WINDOW) { + return "w3_late"; + } + return "total"; } -// trace 和运行日志使用稳定的英文标识，便于现有离线脚本继续按字段筛选。 -static QString trustRegionComponentName(int component) +// 四个窗口在目标曲线的固定对数时间网格上识别。井储段依据早期单位斜率， +// IARF 依据压力导数平台；特征不清楚时按固定比例分段，保证每次评价维度一致。 +static void buildTrustRegionTimeWindows( + const QVector& time, + const QVector& logPressure, + const QVector& logDerivative, + QVector* windowBegin, + QVector* windowEnd) { - if(component == TRUST_REGION_VERTICAL_COMPONENT) { - return "vertical"; + if(!windowBegin || !windowEnd) { + return; } - if(component == TRUST_REGION_HORIZONTAL_COMPONENT) { - return "horizontal"; + + const int pointCount = qMin( + time.size(), qMin(logPressure.size(), logDerivative.size())); + windowBegin->clear(); + windowEnd->clear(); + if(pointCount < 16) { + return; } - if(component == TRUST_REGION_SHAPE_COMPONENT) { - return "shape"; + + QVector logTime(pointCount, 0.0); + QVector derivativeSlope( + pointCount, std::numeric_limits::quiet_NaN()); + for(int i = 0; i < pointCount; ++i) { + logTime[i] = qLn(time[i]); } - return "total"; -} -// 三类损失量纲一致，直接选择当前最大的可靠分量；都很小时退回总残差梯度。 -static int trustRegionDominantComponent( - const AutoFitObjectiveBreakdownLM& breakdown, - double diagnosisThreshold) -{ - int component = TRUST_REGION_TOTAL_COMPONENT; - double largestLoss = diagnosisThreshold; - - if(breakdown.verticalReliable && - isFiniteNumber(breakdown.verticalLoss) && - breakdown.verticalLoss > largestLoss) { - component = TRUST_REGION_VERTICAL_COMPONENT; - largestLoss = breakdown.verticalLoss; + // 五点局部回归比相邻两点差分更不容易把目标曲线的小幅波动误判成流态边界。 + for(int center = 0; center < pointCount; ++center) { + int begin = qMax(0, center - 2); + int end = qMin(pointCount, center + 3); + double meanX = 0.0; + double meanY = 0.0; + for(int i = begin; i < end; ++i) { + meanX += logTime[i]; + meanY += logDerivative[i]; + } + const int count = end - begin; + meanX /= count; + meanY /= count; + double covariance = 0.0; + double variance = 0.0; + for(int i = begin; i < end; ++i) { + double offsetX = logTime[i] - meanX; + covariance += offsetX * (logDerivative[i] - meanY); + variance += offsetX * offsetX; + } + if(variance > 1.0e-14) { + derivativeSlope[center] = covariance / variance; + } } - if(breakdown.horizontalReliable && - !breakdown.registrationAmbiguous && - isFiniteNumber(breakdown.horizontalLoss) && - breakdown.horizontalLoss > largestLoss) { - component = TRUST_REGION_HORIZONTAL_COMPONENT; - largestLoss = breakdown.horizontalLoss; + + const int minimumWindowPoints = qMax(5, pointCount / 16); + int storageEnd = qBound( + minimumWindowPoints, + static_cast(pointCount * 0.25 + 0.5), + pointCount - 3 * minimumWindowPoints); + bool storageIdentified = false; + int lastStoragePoint = -1; + int nonStorageCount = 0; + const int maximumStorageSearch = qMin( + pointCount - 3 * minimumWindowPoints, + static_cast(pointCount * 0.45)); + for(int i = 0; i < maximumStorageSearch; ++i) { + bool storageLike = isFiniteNumber(derivativeSlope[i]) && + derivativeSlope[i] >= 0.65 && + derivativeSlope[i] <= 1.35 && + qAbs(logPressure[i] - logDerivative[i]) <= 1.0; + if(storageLike) { + storageIdentified = true; + lastStoragePoint = i; + nonStorageCount = 0; + } else if(storageIdentified) { + ++nonStorageCount; + if(nonStorageCount >= 3) { + break; + } + } else if(i >= minimumWindowPoints) { + break; + } } - if(isFiniteNumber(breakdown.shapeLoss) && - breakdown.shapeLoss > largestLoss) { - component = TRUST_REGION_SHAPE_COMPONENT; + if(storageIdentified && lastStoragePoint + 1 >= minimumWindowPoints) { + storageEnd = qBound( + minimumWindowPoints, + lastStoragePoint + 1, + pointCount - 3 * minimumWindowPoints); } - return component; + int platformBegin = -1; + int platformEnd = -1; + int runBegin = -1; + int longestRun = 0; + const int platformSearchBegin = storageEnd + minimumWindowPoints; + const int platformSearchEnd = pointCount - minimumWindowPoints; + for(int i = platformSearchBegin; i <= platformSearchEnd; ++i) { + bool platformLike = i < platformSearchEnd && + isFiniteNumber(derivativeSlope[i]) && + qAbs(derivativeSlope[i]) <= 0.15; + if(platformLike && runBegin < 0) { + runBegin = i; + } + if((!platformLike || i == platformSearchEnd) && runBegin >= 0) { + int runEnd = i; + if(runEnd - runBegin > longestRun) { + longestRun = runEnd - runBegin; + platformBegin = runBegin; + platformEnd = runEnd; + } + runBegin = -1; + } + } + + if(longestRun < minimumWindowPoints) { + // 没有清晰平台时使用固定对数时间比例，不对目标曲线强行作流态解释。 + storageEnd = qBound( + minimumWindowPoints, + static_cast(pointCount * 0.25 + 0.5), + pointCount - 3 * minimumWindowPoints); + platformBegin = qBound( + storageEnd + minimumWindowPoints, + static_cast(pointCount * 0.45 + 0.5), + pointCount - 2 * minimumWindowPoints); + platformEnd = qBound( + platformBegin + minimumWindowPoints, + static_cast(pointCount * 0.75 + 0.5), + pointCount - minimumWindowPoints); + } + + // 边界两侧各重叠三个采样点，降低流态边界轻微识别偏差造成的选参跳变。 + const int overlapPoints = 3; + *windowBegin << 0 + << qMax(0, storageEnd - overlapPoints) + << qMax(0, platformBegin - overlapPoints) + << qMax(0, platformEnd - overlapPoints); + *windowEnd << qMin(pointCount, storageEnd + overlapPoints) + << qMin(pointCount, platformBegin + overlapPoints) + << qMin(pointCount, platformEnd + overlapPoints) + << pointCount; } -// 求解选中参数对应的阻尼正规方程。上下和左右诊断量保留方向；形状没有 -// 天然正负，因此使用 shapeLoss 对参数的局部导数。参数最多八维，使用带 -// 部分主元的高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。 +// 求解选中参数对应的阻尼正规方程。参数最多八维，使用带部分主元的 +// 高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。 static bool solveTrustRegionLinearSystem( QVector > matrix, QVector rightHandSide, @@ -419,33 +537,6 @@ static void updateTrustRegionJacobian( } } -// 对上下偏差、左右偏差和形状损失的梯度执行同样的割线秩一修正，使诊断 -// 选参模型与完整残差 Jacobian 保持在同一个已接受工作点。 -static void updateTrustRegionScalarGradient( - QVector* gradient, - double oldValue, - double newValue, - const QVector& coordinateStep) -{ - if(!gradient || gradient->size() != coordinateStep.size() || - !isFiniteNumber(oldValue) || !isFiniteNumber(newValue)) { - return; - } - - double denominator = trustRegionSquaredNorm(coordinateStep); - if(denominator <= 1.0e-12) { - return; - } - - double predictedChange = trustRegionDotProduct( - *gradient, coordinateStep); - double correction = - (newValue - oldValue - predictedChange) / denominator; - for(int i = 0; i < gradient->size(); ++i) { - (*gradient)[i] += correction * coordinateStep[i]; - } -} - static QStringList traceParameterNames() { QStringList names; @@ -849,10 +940,30 @@ void nmCalculationAutoFitLM::writeTraceHeader() << "solver_success" << "elapsed_ms" << "decision" - << "enabled_param_indices" - << "pressure_loss" - << "derivative_loss" - << "vertical_common_bias" + << "enabled_param_indices" + << "pressure_loss" + << "derivative_loss" + << "w0_begin_time" + << "w0_end_time" + << "w0_pressure_loss" + << "w0_derivative_loss" + << "w0_loss" + << "w1_begin_time" + << "w1_end_time" + << "w1_pressure_loss" + << "w1_derivative_loss" + << "w1_loss" + << "w2_begin_time" + << "w2_end_time" + << "w2_pressure_loss" + << "w2_derivative_loss" + << "w2_loss" + << "w3_begin_time" + << "w3_end_time" + << "w3_pressure_loss" + << "w3_derivative_loss" + << "w3_loss" + << "vertical_common_bias" << "vertical_loss" << "vertical_reliable" << "horizontal_physical_shift" @@ -901,7 +1012,7 @@ void nmCalculationAutoFitLM::writeTraceMetaFile() QTextStream out(&metaFile); out << "{\n"; - out << " \"schema_version\": 1,\n"; + out << " \"schema_version\": 2,\n"; out << " \"trace_type\": \"finite_difference_lm_trust_region\",\n"; out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n"; out << " \"created_at\": " @@ -969,8 +1080,29 @@ void nmCalculationAutoFitLM::writeTraceRow( if(objectiveBreakdown && objectiveBreakdown->valid) { cols << traceNumber(objectiveBreakdown->pressureLoss) - << traceNumber(objectiveBreakdown->derivativeLoss) - << traceNumber(objectiveBreakdown->verticalCommonBias) + << traceNumber(objectiveBreakdown->derivativeLoss); + for(int window = 0; + window < TRUST_REGION_TIME_WINDOW_COUNT; + ++window) { + cols << (window < objectiveBreakdown->timeWindowStart.size() + ? traceNumber(objectiveBreakdown->timeWindowStart[window]) + : QString()) + << (window < objectiveBreakdown->timeWindowEndTime.size() + ? traceNumber(objectiveBreakdown->timeWindowEndTime[window]) + : QString()) + << (window < objectiveBreakdown->timeWindowPressureLoss.size() + ? traceNumber( + objectiveBreakdown->timeWindowPressureLoss[window]) + : QString()) + << (window < objectiveBreakdown->timeWindowDerivativeLoss.size() + ? traceNumber( + objectiveBreakdown->timeWindowDerivativeLoss[window]) + : QString()) + << (window < objectiveBreakdown->timeWindowLoss.size() + ? traceNumber(objectiveBreakdown->timeWindowLoss[window]) + : QString()); + } + cols << traceNumber(objectiveBreakdown->verticalCommonBias) << traceNumber(objectiveBreakdown->verticalLoss) << QString::number(objectiveBreakdown->verticalReliable ? 1 : 0) << traceNumber(objectiveBreakdown->horizontalPhysicalShift) @@ -983,7 +1115,7 @@ void nmCalculationAutoFitLM::writeTraceRow( << QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0) << traceNumber(objectiveBreakdown->coverage); } else { - for(int i = 0; i < 14; ++i) { + for(int i = 0; i < 34; ++i) { cols << QString(); } } @@ -1633,7 +1765,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() const double minimumTrustRadius = 2.0e-3; const double maximumTrustRadius = 0.30; const double columnCorrelationLimit = 0.995; - const double diagnosisThreshold = 1.0e-5; // 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。 // 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。 const double effectiveRelativeImprovement = 2.0e-3; @@ -1654,12 +1785,9 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() bool stagnationConfirmationRequested = false; StopReasonLM stopReason = LM_MAX_ITERATIONS; - // jacobian 的行对应固定 160 维残差，列对应用户勾选的参数。 - // 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。 + // jacobian 的行对应固定 160 维残差，列对应用户勾选的参数。窗口选参直接 + // 使用对应残差行和 Jacobian，不再维护额外的上下、左右或形状梯度。 QVector > jacobian; - QVector verticalGradient(dimensions, 0.0); - QVector horizontalGradient(dimensions, 0.0); - QVector shapeGradient(dimensions, 0.0); QVector jacobianColumnValid(dimensions, false); // 参数向量的顺序始终与 m_enabledParamIndices 一致，不能按完整参数索引 @@ -1828,9 +1956,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() jacobian = QVector >( residualCount, QVector(dimensions, 0.0)); - verticalGradient.fill(0.0, dimensions); - horizontalGradient.fill(0.0, dimensions); - shapeGradient.fill(0.0, dimensions); jacobianColumnValid.fill(false, dimensions); TrustRegionEvaluation bestProbe; @@ -1916,27 +2041,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() base.breakdown.residualVector[row]) / delta; } - // 有符号诊断量只有在基点和试算点都可靠时才能计算方向梯度； - // shapeLoss 无方向可靠性标志，始终记录其局部变化率。 - if(base.breakdown.verticalReliable && - probe.breakdown.verticalReliable && - !base.breakdown.registrationAmbiguous && - !probe.breakdown.registrationAmbiguous) { - verticalGradient[column] = - (probe.breakdown.verticalCommonBias - - base.breakdown.verticalCommonBias) / delta; - } - if(base.breakdown.horizontalReliable && - probe.breakdown.horizontalReliable && - !base.breakdown.registrationAmbiguous && - !probe.breakdown.registrationAmbiguous) { - horizontalGradient[column] = - (probe.breakdown.horizontalPhysicalShift - - base.breakdown.horizontalPhysicalShift) / delta; - } - shapeGradient[column] = - (probe.breakdown.shapeLoss - - base.breakdown.shapeLoss) / delta; jacobianColumnValid[column] = true; columnBuilt = true; @@ -1972,27 +2076,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() base.breakdown.residualVector, bestProbe.breakdown.residualVector, acceptedStep); - if(base.breakdown.verticalReliable && - bestProbe.breakdown.verticalReliable) { - updateTrustRegionScalarGradient( - &verticalGradient, - base.breakdown.verticalCommonBias, - bestProbe.breakdown.verticalCommonBias, - acceptedStep); - } - if(base.breakdown.horizontalReliable && - bestProbe.breakdown.horizontalReliable) { - updateTrustRegionScalarGradient( - &horizontalGradient, - base.breakdown.horizontalPhysicalShift, - bestProbe.breakdown.horizontalPhysicalShift, - acceptedStep); - } - updateTrustRegionScalarGradient( - &shapeGradient, - base.breakdown.shapeLoss, - bestProbe.breakdown.shapeLoss, - acceptedStep); current = bestProbe; publishAcceptedPoint(current); @@ -2069,21 +2152,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() } } - // 先确定当前最突出的可靠诊断误差，用其梯度回答“哪些参数最能改善 - // 当前问题”；实际 LM 方向仍由完整残差梯度和 Jacobian 共同计算。 - int dominantComponent = trustRegionDominantComponent( - current.breakdown, diagnosisThreshold); - const QVector* componentGradient = nullptr; - if(dominantComponent == TRUST_REGION_VERTICAL_COMPONENT) { - componentGradient = &verticalGradient; - } else if(dominantComponent == TRUST_REGION_HORIZONTAL_COMPONENT) { - componentGradient = &horizontalGradient; - } else if(dominantComponent == TRUST_REGION_SHAPE_COMPONENT) { - componentGradient = &shapeGradient; - } - // 主目标采用 0.5*||r||^2，其对参数的梯度为 J^T*r。这里不再叠加 - // vertical/horizontal/shape，保证诊断分量不会改变真实接受目标。 + // 窗口损失，保证分段只改变选参，不改变真实接受目标和 LM 方向。 QVector totalGradient(dimensions, 0.0); for(int column = 0; column < dimensions; ++column) { if(!jacobianColumnValid[column]) { @@ -2096,57 +2166,11 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() } } - // 每轮最多联合调整三个灵敏参数。按当前诊断梯度绝对值由大到小选取， - // 并剔除 Jacobian 响应过度共线的列，降低弱可辨识参数互相补偿的风险。 - QVector selectedColumns; - QVector alreadyConsidered(dimensions, false); - for(int selection = 0; selection < qMin(3, dimensions); ++selection) { - int bestColumn = -1; - double bestScore = 0.0; - for(int column = 0; column < dimensions; ++column) { - if(alreadyConsidered[column] || - !jacobianColumnValid[column]) { - continue; - } - - double score = componentGradient - ? qAbs((*componentGradient)[column]) - : qAbs(totalGradient[column]); - if(!isFiniteNumber(score) || score <= bestScore) { - continue; - } - - bool excessivelyCorrelated = false; - for(int selectedIndex = 0; - selectedIndex < selectedColumns.size(); - ++selectedIndex) { - if(trustRegionJacobianColumnCorrelation( - jacobian, - column, - selectedColumns[selectedIndex]) > - columnCorrelationLimit) { - excessivelyCorrelated = true; - break; - } - } - if(!excessivelyCorrelated) { - bestColumn = column; - bestScore = score; - } - } - if(bestColumn < 0 || bestScore <= 1.0e-12) { - break; - } - selectedColumns.append(bestColumn); - alreadyConsidered[bestColumn] = true; - } - - // 诊断梯度接近零时，说明该分量在当前局部无法可靠选参，退回完整残差 - // 梯度，但接受标准仍然只有真实 total，诊断值不会重复计入目标函数。 - if(selectedColumns.isEmpty() && componentGradient) { - dominantComponent = TRUST_REGION_TOTAL_COMPONENT; - componentGradient = nullptr; - alreadyConsidered.fill(false, dimensions); + // 每轮最多联合调整三个灵敏参数，并继续剔除全曲线响应过度共线的列。 + auto selectColumnsByScore = [&](const QVector& scores) + -> QVector { + QVector selected; + QVector alreadyConsidered(dimensions, false); for(int selection = 0; selection < qMin(3, dimensions); ++selection) { @@ -2154,21 +2178,16 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() double bestScore = 0.0; for(int column = 0; column < dimensions; ++column) { if(alreadyConsidered[column] || - !jacobianColumnValid[column]) { - continue; - } - double score = qAbs(totalGradient[column]); - if(score <= bestScore) { + !jacobianColumnValid[column] || + column >= scores.size() || + !isFiniteNumber(scores[column]) || + scores[column] <= bestScore) { continue; } bool excessivelyCorrelated = false; - for(int selectedIndex = 0; - selectedIndex < selectedColumns.size(); - ++selectedIndex) { + for(int i = 0; i < selected.size(); ++i) { if(trustRegionJacobianColumnCorrelation( - jacobian, - column, - selectedColumns[selectedIndex]) > + jacobian, column, selected[i]) > columnCorrelationLimit) { excessivelyCorrelated = true; break; @@ -2176,15 +2195,104 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() } if(!excessivelyCorrelated) { bestColumn = column; - bestScore = score; + bestScore = scores[column]; } } if(bestColumn < 0 || bestScore <= 1.0e-12) { break; } - selectedColumns.append(bestColumn); + selected.append(bestColumn); alreadyConsidered[bestColumn] = true; } + return selected; + }; + + QVector windowOrder; + for(int window = 0; + window < TRUST_REGION_TIME_WINDOW_COUNT; + ++window) { + windowOrder.append(window); + } + std::sort(windowOrder.begin(), windowOrder.end(), + [&](int left, int right) -> bool { + double leftLoss = left < current.breakdown.timeWindowLoss.size() + ? current.breakdown.timeWindowLoss[left] : -1.0; + double rightLoss = right < current.breakdown.timeWindowLoss.size() + ? current.breakdown.timeWindowLoss[right] : -1.0; + return leftLoss > rightLoss; + }); + + QVector selectedColumns; + int dominantWindow = TRUST_REGION_TOTAL_WINDOW; + const int residualPointCount = jacobian.size() / 2; + bool windowsValid = jacobian.size() % 2 == 0 && + current.breakdown.timeWindowBegin.size() == + TRUST_REGION_TIME_WINDOW_COUNT && + current.breakdown.timeWindowEnd.size() == + TRUST_REGION_TIME_WINDOW_COUNT && + current.breakdown.timeWindowLoss.size() == + TRUST_REGION_TIME_WINDOW_COUNT; + + // score 是该参数在窗口内按一维最小二乘预计可消除的残差能量。 + // 它只用于回答“调整哪些参数”，正负方向仍由后面的完整 LM 方程决定。 + if(windowsValid) { + for(int order = 0; + order < windowOrder.size() && selectedColumns.isEmpty(); + ++order) { + int window = windowOrder[order]; + int begin = qBound( + 0, current.breakdown.timeWindowBegin[window], + residualPointCount); + int end = qBound( + begin, current.breakdown.timeWindowEnd[window], + residualPointCount); + int pointCount = end - begin; + if(pointCount <= 0 || + !isFiniteNumber( + current.breakdown.timeWindowLoss[window])) { + continue; + } + + QVector windowScores(dimensions, 0.0); + for(int column = 0; column < dimensions; ++column) { + if(!jacobianColumnValid[column]) { + continue; + } + double gradient = 0.0; + double curvature = 0.0; + for(int point = begin; point < end; ++point) { + int pressureRow = point; + int derivativeRow = residualPointCount + point; + gradient += + jacobian[pressureRow][column] * + current.breakdown.residualVector[pressureRow] + + jacobian[derivativeRow][column] * + current.breakdown.residualVector[derivativeRow]; + curvature += + jacobian[pressureRow][column] * + jacobian[pressureRow][column] + + jacobian[derivativeRow][column] * + jacobian[derivativeRow][column]; + } + windowScores[column] = + gradient * gradient / + (curvature + 1.0e-20) / pointCount; + } + selectedColumns = selectColumnsByScore(windowScores); + if(!selectedColumns.isEmpty()) { + dominantWindow = window; + } + } + } + + // 所有窗口都缺少可用灵敏方向时退回完整残差梯度，避免分段识别异常 + // 阻断细搜；候选仍然只按全曲线 total 是否下降来接受。 + if(selectedColumns.isEmpty()) { + QVector totalScores(dimensions, 0.0); + for(int column = 0; column < dimensions; ++column) { + totalScores[column] = qAbs(totalGradient[column]); + } + selectedColumns = selectColumnsByScore(totalScores); } // 当前局部没有可用方向时先缩小半径并重建灵敏度；只有已经在最小 @@ -2370,37 +2478,11 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() consecutiveSolverFailures = 0; // 有效候选即使最终被拒绝，也提供了一条真实割线，可用于修正下一轮 // 局部模型；是否成为新工作点仍只由下面的 total 严格比较决定。 - const AutoFitObjectiveBreakdownLM oldBreakdown = current.breakdown; updateTrustRegionJacobian( &jacobian, - oldBreakdown.residualVector, + current.breakdown.residualVector, candidate.breakdown.residualVector, coordinateStep); - if(oldBreakdown.verticalReliable && - candidate.breakdown.verticalReliable && - !oldBreakdown.registrationAmbiguous && - !candidate.breakdown.registrationAmbiguous) { - updateTrustRegionScalarGradient( - &verticalGradient, - oldBreakdown.verticalCommonBias, - candidate.breakdown.verticalCommonBias, - coordinateStep); - } - if(oldBreakdown.horizontalReliable && - candidate.breakdown.horizontalReliable && - !oldBreakdown.registrationAmbiguous && - !candidate.breakdown.registrationAmbiguous) { - updateTrustRegionScalarGradient( - &horizontalGradient, - oldBreakdown.horizontalPhysicalShift, - candidate.breakdown.horizontalPhysicalShift, - coordinateStep); - } - updateTrustRegionScalarGradient( - &shapeGradient, - oldBreakdown.shapeLoss, - candidate.breakdown.shapeLoss, - coordinateStep); // reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确； // 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。 @@ -2409,7 +2491,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() candidate.fitness * candidate.fitness); double reductionRatio = actualReduction / predictedReduction; bool accepted = candidate.fitness < current.fitness; - QString componentName = trustRegionComponentName(dominantComponent); + QString windowName = trustRegionWindowName(dominantWindow); if(accepted) { // 真实总误差下降后才正式替换 current，并同步发布参数、曲线和诊断。 @@ -2466,27 +2548,29 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting() candidate.parameters, candidate.fitness, true, - candidate.elapsedMs, - accepted - ? QString("accepted_%1").arg(componentName) - : QString("rejected_%1").arg(componentName), - &candidate.breakdown); - - QString componentDisplayName = componentName; - if(componentName == "vertical") { - componentDisplayName = tr("vertical deviation"); - } else if(componentName == "horizontal") { - componentDisplayName = tr("horizontal deviation"); - } else if(componentName == "shape") { - componentDisplayName = tr("shape deviation"); - } else if(componentName == "total") { - componentDisplayName = tr("total error"); + candidate.elapsedMs, + accepted + ? QString("accepted_%1").arg(windowName) + : QString("rejected_%1").arg(windowName), + &candidate.breakdown); + + QString windowDisplayName = windowName; + if(dominantWindow == TRUST_REGION_WELLBORE_STORAGE_WINDOW) { + windowDisplayName = tr("wellbore-storage interval"); + } else if(dominantWindow == TRUST_REGION_TRANSITION_WINDOW) { + windowDisplayName = tr("early transition interval"); + } else if(dominantWindow == TRUST_REGION_IARF_WINDOW) { + windowDisplayName = tr("IARF interval"); + } else if(dominantWindow == TRUST_REGION_LATE_WINDOW) { + windowDisplayName = tr("late-time interval"); + } else { + windowDisplayName = tr("total error"); } emit logMessageGenerated( tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5") .arg(iteration + 1) - .arg(componentDisplayName) + .arg(windowDisplayName) .arg(selectedColumns.size()) .arg(candidate.fitness, 0, 'e', 4) .arg(accepted ? tr("accepted") : tr("rejected"))); @@ -3653,6 +3737,45 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError( derivativeResidual[i]); } + // 窗口只由目标曲线确定，因此所有候选使用完全相同的时间分段。每个窗口 + // 同时统计压力、导数和二者等权联合 RMS，供细搜选择当前关注区间。 + buildTrustRegionTimeWindows( + commonX, + targetLogPressure, + targetLogDerivative, + &breakdown.timeWindowBegin, + &breakdown.timeWindowEnd); + if(breakdown.timeWindowBegin.size() != + TRUST_REGION_TIME_WINDOW_COUNT || + breakdown.timeWindowEnd.size() != + TRUST_REGION_TIME_WINDOW_COUNT) { + return invalidLoss; + } + for(int window = 0; + window < TRUST_REGION_TIME_WINDOW_COUNT; + ++window) { + int begin = breakdown.timeWindowBegin[window]; + int end = breakdown.timeWindowEnd[window]; + if(begin < 0 || end <= begin || end > numPoints) { + return invalidLoss; + } + double windowPressureLoss = rmse( + pressureResidual, begin, end); + double windowDerivativeLoss = rmse( + derivativeResidual, begin, end); + if(!isFiniteNumber(windowPressureLoss) || + !isFiniteNumber(windowDerivativeLoss)) { + return invalidLoss; + } + breakdown.timeWindowStart.append(commonX[begin]); + breakdown.timeWindowEndTime.append(commonX[end - 1]); + breakdown.timeWindowPressureLoss.append(windowPressureLoss); + breakdown.timeWindowDerivativeLoss.append(windowDerivativeLoss); + breakdown.timeWindowLoss.append(qSqrt( + 0.5 * windowPressureLoss * windowPressureLoss + + 0.5 * windowDerivativeLoss * windowDerivativeLoss)); + } + const double logGridStep = (targetLogMaxX - targetLogMinX) / (numPoints - 1);