如何用OpenCV实现指数衰减曲线拟合并提取衰减常数？

一、OpenCV的非线性最小二乘拟合实现 #

OpenCV支持通过自定义代价函数实现非线性最小二乘拟合，可替代GSL的相关功能。针对指数衰减函数 y = A·exp(-x/τ) + B（τ为目标衰减常数），实现步骤如下：

1. 定义残差计算逻辑：计算每个数据点与拟合函数的偏差
2. 初始化参数猜测值：基于原始数据给出A、τ、B的初始估计（比如A取首点y值，B取末尾点均值，τ设为经验值）
3. 构建优化问题并求解：利用cv::MinProblemSolver执行最小化残差平方和的迭代优化

示例代码：

#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <vector>
#include <cmath>

// 计算残差：每个数据点的实际值与拟合值的差
void computeResiduals(const cv::Mat& params, const cv::Mat& x, const cv::Mat& y, cv::Mat& residuals) {
    const double A = params.at<double>(0);
    const double tau = params.at<double>(1);
    const double B = params.at<double>(2);
    residuals.create(x.rows, 1, CV_64F);
    for (int i = 0; i < x.rows; ++i) {
        const double xi = x.at<double>(i);
        const double yi = y.at<double>(i);
        residuals.at<double>(i) = yi - (A * exp(-xi / tau) + B);
    }
}

int main() {
    // 替换为你的600个数据点，转为CV_64F格式的Mat
    cv::Mat x_data = ...; 
    cv::Mat y_data = ...;

    // 初始化参数：[A, tau, B]
    cv::Mat params = (cv::Mat_<double>(3,1) << y_data.at<double>(0), 10.0, y_data.at<double>(y_data.rows-1));

    // 创建最小化求解器
    cv::Ptr<cv::MinProblemSolver> solver = cv::MinProblemSolver::create();
    // 定义代价函数：残差平方和
    auto costFunc = [&](const cv::Mat& p) {
        cv::Mat res;
        computeResiduals(p, x_data, y_data, res);
        return cv::norm(res, cv::NORM_L2SQR);
    };
    solver->setFunction(costFunc);
    // 设置终止条件：最大迭代1000次，精度1e-6
    solver->setTermCriteria(cv::TermCriteria(cv::TermCriteria::EPS + cv::TermCriteria::MAX_ITER, 1000, 1e-6));

    // 执行优化
    solver->minimize(params);

    // 提取衰减常数τ
    const double tau = params.at<double>(1);
    return 0;
}

二、fitline()的曲线拟合用法 #

cv::fitline()仅支持直线拟合，无法直接处理非线性的指数衰减曲线，但可通过线性变换将指数模型转为线性模型后使用：

1. 对指数衰减函数做变形：假设基线B已知（可通过末尾N个点均值估计），则 y - B = A·exp(-x/τ)，两边取自然对数得 ln(y-B) = ln(A) - x/τ
2. 令 Y = ln(y-B)、C = ln(A)、D = -1/τ，则模型转为线性模型 Y = C + D·x
3. 用fitline()拟合变换后的线性点集，再反推τ

示例代码：

#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <vector>
#include <cmath>

int main() {
    std::vector<cv::Point2d> raw_data; // 存储原始(x,y)数据点
    // 填充raw_data...

    // 估计基线B：取最后50个点的均值
    const int baseline_sample = 50;
    double B = 0.0;
    for (int i = raw_data.size() - baseline_sample; i < raw_data.size(); ++i) {
        B += raw_data[i].y;
    }
    B /= baseline_sample;

    // 构建线性拟合的点集，过滤无意义的负数值
    std::vector<cv::Point2d> linear_data;
    for (const auto& p : raw_data) {
        const double y_prime = p.y - B;
        if (y_prime > 1e-6) {
            linear_data.emplace_back(p.x, log(y_prime));
        }
    }

    // 拟合直线
    cv::Vec4f line_params;
    cv::fitline(linear_data, line_params, cv::DIST_L2, 0, 0.01, 0.01);
    // 计算斜率D = 方向向量y / 方向向量x
    const double D = line_params[1] / line_params[0];
    // 反推衰减常数τ
    const double tau = -1.0 / D;

    return 0;
}

注意：该方法依赖基线B的估计精度，对噪声敏感，适合噪声较小的数据集。

三、从SVR拟合结果中提取衰减常数 #

SVR是黑箱回归模型，不会直接输出指数衰减的参数。需先通过SVR得到拟合曲线，再对拟合曲线执行参数估计：

1. 用训练好的SVR模型对原始x范围的点进行预测，得到拟合后的y值序列
2. 对拟合后的(x,y)序列，使用第一部分的非线性最小二乘拟合，或第二部分的线性变换拟合，估计出τ

示例流程：

#include <opencv2/ml.hpp>
#include <vector>
#include <cmath>

int main() {
    // 加载训练好的SVR模型
    cv::Ptr<cv::ml::SVR> svr_model = cv::ml::SVR::load("trained_svr.xml");

    // 生成覆盖原始x范围的预测点
    const double x_min = ...; // 原始数据的最小x值
    const double x_max = ...; // 原始数据的最大x值
    std::vector<float> x_pred;
    for (double x = x_min; x <= x_max; x += 0.1) {
        x_pred.push_back(static_cast<float>(x));
    }
    cv::Mat x_mat(x_pred, true);
    x_mat = x_mat.reshape(1, x_pred.size());

    // 预测拟合y值
    cv::Mat y_pred_mat;
    svr_model->predict(x_mat, y_pred_mat);

    // 转换为(x,y)点集
    std::vector<cv::Point2d> fit_points;
    for (int i = 0; i < x_pred.size(); ++i) {
        fit_points.emplace_back(x_pred[i], y_pred_mat.at<float>(i));
    }

    // 对fit_points执行非线性最小二乘拟合（参考第一部分代码），提取τ
    // ...

    return 0;
}

注意：SVR拟合的曲线可能并非严格指数衰减，因此提取的τ是近似值，精度取决于SVR对指数模型的拟合程度。

内容的提问来源于stack exchange，提问作者Rolly Ng