Compare commits

...

1 Commits

@ -13,14 +13,21 @@
class nmDataWellBase; class nmDataWellBase;
// 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据 // 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据;时间窗口误差
// 其余诊断量用于有限差分灵敏度分析和信赖域选参 // 用于细搜选参,原有上下、左右和形状量仅保留为过程诊断
struct AutoFitObjectiveBreakdownLM { struct AutoFitObjectiveBreakdownLM {
bool valid; bool valid;
double total; double total;
double pressureLoss; double pressureLoss;
double derivativeLoss; double derivativeLoss;
QVector<double> residualVector; QVector<double> residualVector;
QVector<int> timeWindowBegin;
QVector<int> timeWindowEnd;
QVector<double> timeWindowStart;
QVector<double> timeWindowEndTime;
QVector<double> timeWindowPressureLoss;
QVector<double> timeWindowDerivativeLoss;
QVector<double> timeWindowLoss;
double verticalCommonBias; double verticalCommonBias;
double verticalLoss; double verticalLoss;
bool verticalReliable; bool verticalReliable;

@ -169,12 +169,14 @@ static double fromTrustRegionCoordinate(double coordinate,
return lower + coordinate * (upper - lower); return lower + coordinate * (upper - lower);
} }
enum TrustRegionErrorComponent enum TrustRegionTimeWindow
{ {
TRUST_REGION_VERTICAL_COMPONENT = 0, TRUST_REGION_WELLBORE_STORAGE_WINDOW = 0,
TRUST_REGION_HORIZONTAL_COMPONENT, TRUST_REGION_TRANSITION_WINDOW,
TRUST_REGION_SHAPE_COMPONENT, TRUST_REGION_IARF_WINDOW,
TRUST_REGION_TOTAL_COMPONENT TRUST_REGION_LATE_WINDOW,
TRUST_REGION_TIME_WINDOW_COUNT,
TRUST_REGION_TOTAL_WINDOW
}; };
// 一次真实求解的完整快照。除了参数和总误差,还保存内部坐标、诊断分量和 // 一次真实求解的完整快照。除了参数和总误差,还保存内部坐标、诊断分量和
@ -202,7 +204,13 @@ static bool trustRegionResidualsValid(
{ {
// 损失函数固定使用 80 个压力点和 80 个导数点。严格校验长度,避免 // 损失函数固定使用 80 个压力点和 80 个导数点。严格校验长度,避免
// Jacobian 沿用旧维度后访问另一候选的短残差向量。 // 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; return false;
} }
@ -211,6 +219,17 @@ static bool trustRegionResidualsValid(
return false; 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; return true;
} }
@ -224,68 +243,167 @@ static double trustRegionSquaredNorm(const QVector<double>& values)
return sum; return sum;
} }
// 计算同维向量内积;维度不一致表示局部模型无效,返回零让调用方放弃修正。 // trace 和运行日志使用稳定的英文标识,便于离线过程分析按窗口筛选。
static double trustRegionDotProduct(const QVector<double>& left, static QString trustRegionWindowName(int window)
const QVector<double>& right)
{ {
if(left.size() != right.size()) { if(window == TRUST_REGION_WELLBORE_STORAGE_WINDOW) {
return 0.0; return "w0_storage";
} }
if(window == TRUST_REGION_TRANSITION_WINDOW) {
double sum = 0.0; return "w1_transition";
for(int i = 0; i < left.size(); ++i) {
sum += left[i] * right[i];
} }
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<double>& time,
const QVector<double>& logPressure,
const QVector<double>& logDerivative,
QVector<int>* windowBegin,
QVector<int>* windowEnd)
{ {
if(component == TRUST_REGION_VERTICAL_COMPONENT) { if(!windowBegin || !windowEnd) {
return "vertical"; 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<double> logTime(pointCount, 0.0);
QVector<double> derivativeSlope(
pointCount, std::numeric_limits<double>::quiet_NaN());
for(int i = 0; i < pointCount; ++i) {
logTime[i] = qLn(time[i]);
} }
return "total";
}
// 三类损失量纲一致,直接选择当前最大的可靠分量;都很小时退回总残差梯度。 // 五点局部回归比相邻两点差分更不容易把目标曲线的小幅波动误判成流态边界。
static int trustRegionDominantComponent( for(int center = 0; center < pointCount; ++center) {
const AutoFitObjectiveBreakdownLM& breakdown, int begin = qMax(0, center - 2);
double diagnosisThreshold) int end = qMin(pointCount, center + 3);
{ double meanX = 0.0;
int component = TRUST_REGION_TOTAL_COMPONENT; double meanY = 0.0;
double largestLoss = diagnosisThreshold; for(int i = begin; i < end; ++i) {
meanX += logTime[i];
if(breakdown.verticalReliable && meanY += logDerivative[i];
isFiniteNumber(breakdown.verticalLoss) && }
breakdown.verticalLoss > largestLoss) { const int count = end - begin;
component = TRUST_REGION_VERTICAL_COMPONENT; meanX /= count;
largestLoss = breakdown.verticalLoss; 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 && const int minimumWindowPoints = qMax(5, pointCount / 16);
isFiniteNumber(breakdown.horizontalLoss) && int storageEnd = qBound(
breakdown.horizontalLoss > largestLoss) { minimumWindowPoints,
component = TRUST_REGION_HORIZONTAL_COMPONENT; static_cast<int>(pointCount * 0.25 + 0.5),
largestLoss = breakdown.horizontalLoss; pointCount - 3 * minimumWindowPoints);
bool storageIdentified = false;
int lastStoragePoint = -1;
int nonStorageCount = 0;
const int maximumStorageSearch = qMin(
pointCount - 3 * minimumWindowPoints,
static_cast<int>(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) && if(storageIdentified && lastStoragePoint + 1 >= minimumWindowPoints) {
breakdown.shapeLoss > largestLoss) { storageEnd = qBound(
component = TRUST_REGION_SHAPE_COMPONENT; 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<int>(pointCount * 0.25 + 0.5),
pointCount - 3 * minimumWindowPoints);
platformBegin = qBound(
storageEnd + minimumWindowPoints,
static_cast<int>(pointCount * 0.45 + 0.5),
pointCount - 2 * minimumWindowPoints);
platformEnd = qBound(
platformBegin + minimumWindowPoints,
static_cast<int>(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( static bool solveTrustRegionLinearSystem(
QVector<QVector<double> > matrix, QVector<QVector<double> > matrix,
QVector<double> rightHandSide, QVector<double> rightHandSide,
@ -419,33 +537,6 @@ static void updateTrustRegionJacobian(
} }
} }
// 对上下偏差、左右偏差和形状损失的梯度执行同样的割线秩一修正,使诊断
// 选参模型与完整残差 Jacobian 保持在同一个已接受工作点。
static void updateTrustRegionScalarGradient(
QVector<double>* gradient,
double oldValue,
double newValue,
const QVector<double>& 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() static QStringList traceParameterNames()
{ {
QStringList names; QStringList names;
@ -849,10 +940,30 @@ void nmCalculationAutoFitLM::writeTraceHeader()
<< "solver_success" << "solver_success"
<< "elapsed_ms" << "elapsed_ms"
<< "decision" << "decision"
<< "enabled_param_indices" << "enabled_param_indices"
<< "pressure_loss" << "pressure_loss"
<< "derivative_loss" << "derivative_loss"
<< "vertical_common_bias" << "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_loss"
<< "vertical_reliable" << "vertical_reliable"
<< "horizontal_physical_shift" << "horizontal_physical_shift"
@ -901,7 +1012,7 @@ void nmCalculationAutoFitLM::writeTraceMetaFile()
QTextStream out(&metaFile); QTextStream out(&metaFile);
out << "{\n"; out << "{\n";
out << " \"schema_version\": 1,\n"; out << " \"schema_version\": 2,\n";
out << " \"trace_type\": \"finite_difference_lm_trust_region\",\n"; out << " \"trace_type\": \"finite_difference_lm_trust_region\",\n";
out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n"; out << " \"run_id\": " << jsonEscape(m_traceRunId) << ",\n";
out << " \"created_at\": " out << " \"created_at\": "
@ -969,8 +1080,29 @@ void nmCalculationAutoFitLM::writeTraceRow(
if(objectiveBreakdown && objectiveBreakdown->valid) { if(objectiveBreakdown && objectiveBreakdown->valid) {
cols << traceNumber(objectiveBreakdown->pressureLoss) cols << traceNumber(objectiveBreakdown->pressureLoss)
<< traceNumber(objectiveBreakdown->derivativeLoss) << traceNumber(objectiveBreakdown->derivativeLoss);
<< traceNumber(objectiveBreakdown->verticalCommonBias) 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) << traceNumber(objectiveBreakdown->verticalLoss)
<< QString::number(objectiveBreakdown->verticalReliable ? 1 : 0) << QString::number(objectiveBreakdown->verticalReliable ? 1 : 0)
<< traceNumber(objectiveBreakdown->horizontalPhysicalShift) << traceNumber(objectiveBreakdown->horizontalPhysicalShift)
@ -983,7 +1115,7 @@ void nmCalculationAutoFitLM::writeTraceRow(
<< QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0) << QString::number(objectiveBreakdown->registrationAmbiguous ? 1 : 0)
<< traceNumber(objectiveBreakdown->coverage); << traceNumber(objectiveBreakdown->coverage);
} else { } else {
for(int i = 0; i < 14; ++i) { for(int i = 0; i < 34; ++i) {
cols << QString(); cols << QString();
} }
} }
@ -1633,7 +1765,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
const double minimumTrustRadius = 2.0e-3; const double minimumTrustRadius = 2.0e-3;
const double maximumTrustRadius = 0.30; const double maximumTrustRadius = 0.30;
const double columnCorrelationLimit = 0.995; const double columnCorrelationLimit = 0.995;
const double diagnosisThreshold = 1.0e-5;
// 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。 // 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。
// 更小的下降仍保留为最佳解,但不能反复清除停滞状态、延长拟合时间。 // 更小的下降仍保留为最佳解,但不能反复清除停滞状态、延长拟合时间。
const double effectiveRelativeImprovement = 2.0e-3; const double effectiveRelativeImprovement = 2.0e-3;
@ -1654,12 +1785,9 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
bool stagnationConfirmationRequested = false; bool stagnationConfirmationRequested = false;
StopReasonLM stopReason = LM_MAX_ITERATIONS; StopReasonLM stopReason = LM_MAX_ITERATIONS;
// jacobian 的行对应固定 160 维残差,列对应用户勾选的参数。 // jacobian 的行对应固定 160 维残差,列对应用户勾选的参数。窗口选参直接
// 三个 gradient 单独描述诊断分量对参数的局部变化,只用于本轮选参 // 使用对应残差行和 Jacobian不再维护额外的上下、左右或形状梯度
QVector<QVector<double> > jacobian; QVector<QVector<double> > jacobian;
QVector<double> verticalGradient(dimensions, 0.0);
QVector<double> horizontalGradient(dimensions, 0.0);
QVector<double> shapeGradient(dimensions, 0.0);
QVector<bool> jacobianColumnValid(dimensions, false); QVector<bool> jacobianColumnValid(dimensions, false);
// 参数向量的顺序始终与 m_enabledParamIndices 一致,不能按完整参数索引 // 参数向量的顺序始终与 m_enabledParamIndices 一致,不能按完整参数索引
@ -1828,9 +1956,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
jacobian = QVector<QVector<double> >( jacobian = QVector<QVector<double> >(
residualCount, QVector<double>(dimensions, 0.0)); residualCount, QVector<double>(dimensions, 0.0));
verticalGradient.fill(0.0, dimensions);
horizontalGradient.fill(0.0, dimensions);
shapeGradient.fill(0.0, dimensions);
jacobianColumnValid.fill(false, dimensions); jacobianColumnValid.fill(false, dimensions);
TrustRegionEvaluation bestProbe; TrustRegionEvaluation bestProbe;
@ -1916,27 +2041,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
base.breakdown.residualVector[row]) / delta; 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; jacobianColumnValid[column] = true;
columnBuilt = true; columnBuilt = true;
@ -1972,27 +2076,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
base.breakdown.residualVector, base.breakdown.residualVector,
bestProbe.breakdown.residualVector, bestProbe.breakdown.residualVector,
acceptedStep); 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; current = bestProbe;
publishAcceptedPoint(current); publishAcceptedPoint(current);
@ -2069,21 +2152,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
} }
} }
// 先确定当前最突出的可靠诊断误差,用其梯度回答“哪些参数最能改善
// 当前问题”;实际 LM 方向仍由完整残差梯度和 Jacobian 共同计算。
int dominantComponent = trustRegionDominantComponent(
current.breakdown, diagnosisThreshold);
const QVector<double>* 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。这里不再叠加 // 主目标采用 0.5*||r||^2其对参数的梯度为 J^T*r。这里不再叠加
// vertical/horizontal/shape保证诊断分量不会改变真实接受目标 // 窗口损失,保证分段只改变选参,不改变真实接受目标和 LM 方向。
QVector<double> totalGradient(dimensions, 0.0); QVector<double> totalGradient(dimensions, 0.0);
for(int column = 0; column < dimensions; ++column) { for(int column = 0; column < dimensions; ++column) {
if(!jacobianColumnValid[column]) { if(!jacobianColumnValid[column]) {
@ -2096,57 +2166,11 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
} }
} }
// 每轮最多联合调整三个灵敏参数。按当前诊断梯度绝对值由大到小选取, // 每轮最多联合调整三个灵敏参数,并继续剔除全曲线响应过度共线的列。
// 并剔除 Jacobian 响应过度共线的列,降低弱可辨识参数互相补偿的风险。 auto selectColumnsByScore = [&](const QVector<double>& scores)
QVector<int> selectedColumns; -> QVector<int> {
QVector<bool> alreadyConsidered(dimensions, false); QVector<int> selected;
for(int selection = 0; selection < qMin(3, dimensions); ++selection) { QVector<bool> alreadyConsidered(dimensions, false);
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);
for(int selection = 0; for(int selection = 0;
selection < qMin(3, dimensions); selection < qMin(3, dimensions);
++selection) { ++selection) {
@ -2154,21 +2178,16 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
double bestScore = 0.0; double bestScore = 0.0;
for(int column = 0; column < dimensions; ++column) { for(int column = 0; column < dimensions; ++column) {
if(alreadyConsidered[column] || if(alreadyConsidered[column] ||
!jacobianColumnValid[column]) { !jacobianColumnValid[column] ||
continue; column >= scores.size() ||
} !isFiniteNumber(scores[column]) ||
double score = qAbs(totalGradient[column]); scores[column] <= bestScore) {
if(score <= bestScore) {
continue; continue;
} }
bool excessivelyCorrelated = false; bool excessivelyCorrelated = false;
for(int selectedIndex = 0; for(int i = 0; i < selected.size(); ++i) {
selectedIndex < selectedColumns.size();
++selectedIndex) {
if(trustRegionJacobianColumnCorrelation( if(trustRegionJacobianColumnCorrelation(
jacobian, jacobian, column, selected[i]) >
column,
selectedColumns[selectedIndex]) >
columnCorrelationLimit) { columnCorrelationLimit) {
excessivelyCorrelated = true; excessivelyCorrelated = true;
break; break;
@ -2176,15 +2195,104 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
} }
if(!excessivelyCorrelated) { if(!excessivelyCorrelated) {
bestColumn = column; bestColumn = column;
bestScore = score; bestScore = scores[column];
} }
} }
if(bestColumn < 0 || bestScore <= 1.0e-12) { if(bestColumn < 0 || bestScore <= 1.0e-12) {
break; break;
} }
selectedColumns.append(bestColumn); selected.append(bestColumn);
alreadyConsidered[bestColumn] = true; alreadyConsidered[bestColumn] = true;
} }
return selected;
};
QVector<int> 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<int> 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<double> 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<double> 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; consecutiveSolverFailures = 0;
// 有效候选即使最终被拒绝,也提供了一条真实割线,可用于修正下一轮 // 有效候选即使最终被拒绝,也提供了一条真实割线,可用于修正下一轮
// 局部模型;是否成为新工作点仍只由下面的 total 严格比较决定。 // 局部模型;是否成为新工作点仍只由下面的 total 严格比较决定。
const AutoFitObjectiveBreakdownLM oldBreakdown = current.breakdown;
updateTrustRegionJacobian( updateTrustRegionJacobian(
&jacobian, &jacobian,
oldBreakdown.residualVector, current.breakdown.residualVector,
candidate.breakdown.residualVector, candidate.breakdown.residualVector,
coordinateStep); 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 表示预测准确; // reductionRatio 衡量局部线性模型的可信度:接近 1 表示预测准确;
// 值较小表示虽然可能下降,但模型低估了非线性,需要收紧下一步。 // 值较小表示虽然可能下降,但模型低估了非线性,需要收紧下一步。
@ -2409,7 +2491,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
candidate.fitness * candidate.fitness); candidate.fitness * candidate.fitness);
double reductionRatio = actualReduction / predictedReduction; double reductionRatio = actualReduction / predictedReduction;
bool accepted = candidate.fitness < current.fitness; bool accepted = candidate.fitness < current.fitness;
QString componentName = trustRegionComponentName(dominantComponent); QString windowName = trustRegionWindowName(dominantWindow);
if(accepted) { if(accepted) {
// 真实总误差下降后才正式替换 current并同步发布参数、曲线和诊断。 // 真实总误差下降后才正式替换 current并同步发布参数、曲线和诊断。
@ -2466,27 +2548,29 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
candidate.parameters, candidate.parameters,
candidate.fitness, candidate.fitness,
true, true,
candidate.elapsedMs, candidate.elapsedMs,
accepted accepted
? QString("accepted_%1").arg(componentName) ? QString("accepted_%1").arg(windowName)
: QString("rejected_%1").arg(componentName), : QString("rejected_%1").arg(windowName),
&candidate.breakdown); &candidate.breakdown);
QString componentDisplayName = componentName; QString windowDisplayName = windowName;
if(componentName == "vertical") { if(dominantWindow == TRUST_REGION_WELLBORE_STORAGE_WINDOW) {
componentDisplayName = tr("vertical deviation"); windowDisplayName = tr("wellbore-storage interval");
} else if(componentName == "horizontal") { } else if(dominantWindow == TRUST_REGION_TRANSITION_WINDOW) {
componentDisplayName = tr("horizontal deviation"); windowDisplayName = tr("early transition interval");
} else if(componentName == "shape") { } else if(dominantWindow == TRUST_REGION_IARF_WINDOW) {
componentDisplayName = tr("shape deviation"); windowDisplayName = tr("IARF interval");
} else if(componentName == "total") { } else if(dominantWindow == TRUST_REGION_LATE_WINDOW) {
componentDisplayName = tr("total error"); windowDisplayName = tr("late-time interval");
} else {
windowDisplayName = tr("total error");
} }
emit logMessageGenerated( emit logMessageGenerated(
tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5") tr("Iteration %1: focus=%2, parameters=%3, error=%4, result=%5")
.arg(iteration + 1) .arg(iteration + 1)
.arg(componentDisplayName) .arg(windowDisplayName)
.arg(selectedColumns.size()) .arg(selectedColumns.size())
.arg(candidate.fitness, 0, 'e', 4) .arg(candidate.fitness, 0, 'e', 4)
.arg(accepted ? tr("accepted") : tr("rejected"))); .arg(accepted ? tr("accepted") : tr("rejected")));
@ -3653,6 +3737,45 @@ double nmCalculationAutoFitLM::calculateLogLogCurveError(
derivativeResidual[i]); 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 = const double logGridStep =
(targetLogMaxX - targetLogMinX) / (targetLogMaxX - targetLogMinX) /
(numPoints - 1); (numPoints - 1);

Loading…
Cancel
Save