基于光纤传感器的平面曲线重建算法建模

Modeling of Plane Curve Reconstruction Algorithm Based on Optical Fiber Sensor

摘要: 为推动光纤传感技术的多领域应用，解决其平面曲线重建建模的核心难题，本文构建了基于四阶龙格–库塔法的平面曲线重构模型，提出了一套完整的曲率计算与曲线拟合方案。针对离散测点连续曲率分布的获取问题，首先依据曲率与波长的关联式，借助传感器波长测量数据，计算出两个初始状态下FBG1-FBG6的曲率；采用三次样条插值法拟合曲率曲线，经SPSS软件验证，模型相关系数R²均近似为1，精准明确了曲率与横坐标x的量化关系，并实现对应位置曲率的可靠估算。为完成平面曲线的精准重构，基于上述曲率数据构建微分方程，代入曲率函数后采用四阶龙格–库塔法求解微分方程近似解，通过Python编程迭代得到重构平面曲线，并分析了重构曲线的核心特征。此外，为优化曲线重建的采样策略并探究误差来源，以平面曲线y = x³ +x (0 ≤ x ≤ 1)为研究对象，通过Python程序计算其弧长，设定采样点数量n = 10并确定各采样点的横坐标与弧长间隔；求解各采样点的一、二阶导数，结合曲率公式计算曲率值，再通过迭代法求解弧长、切线方向变化量等参数以补充坐标点，最终绘制重构曲线与原曲线的对比图像。分析表明，重构过程中的误差主要来源于采样设计、离散数据处理、曲率近似计算及数值求解等环节。

Abstract: To promote the multi-field applications of fiber optic sensing technology and address the core challenge of planar curve reconstruction modeling, this paper constructs a planar curve reconstruction model based on the fourth-order Runge-Kutta method and proposes a complete set of solutions for curvature calculation and curve fitting. To obtain the continuous curvature distribution from discrete measurement points, first, based on the correlation between curvature and wavelength, the curvatures of FBG1-FBG6 in two initial states are calculated using the sensor’s wavelength measurement data. Then, the cubic spline interpolation method is employed to fit the curvature curve. Verified by SPSS software, the model’s correlation coefficient R² is approximately 1, accurately defining the quantitative relationship between curvature and the horizontal coordinate x, and enabling reliable estimation of curvature at corresponding positions. To achieve precise reconstruction of the planar curve, a differential equation is constructed based on the above curvature data. After substituting the curvature function, the fourth-order Runge-Kutta method is used to solve the approximate solution of the differential equation. Through Python programming iteration, the reconstructed planar curve is obtained, and its core characteristics are analyzed. Additionally, to optimize the sampling strategy for curve reconstruction and explore error sources, the plane curve y = x³ + x (0 ≤ x ≤ 1) is taken as the research object. Using a Python program, its arc length is calculated, the number of sampling points n = 10 is set, and the horizontal coordinates and arc length intervals of each sampling point are determined. The first and second derivatives at each sampling point are solved, and curvature values are calculated in combination with the curvature formula. Then, iterative methods are used to solve parameters such as arc length and tangent direction changes to supplement coordinate points, ultimately drawing a comparative image of the reconstructed curve and the original curve. Analysis shows that the errors in the reconstruction process mainly originate from sampling design, discrete data processing, approximate curvature calculation, and numerical solution.

1. 引言

当前光纤传感技术在结构形变监测中的应用需求与重建精度不足的矛盾日益突出[1]，如何基于离散传感数据实现平面曲线精准重建成为工程检测的核心挑战。传统方法多依赖经验或低阶拟合模型，难以应对连续曲率的平滑性需求[2]。尤其在高精度场景下，传感数据的离散性对曲线拟合影响显著，亟需精细化的数值方法[3]。三次样条插值与四阶龙格–库塔法作为高效计算工具，在数据平滑、微分方程求解领域展现出优势[4]，为曲线重构提供了新思路。然而现有研究多聚焦单一方法应用，对“曲率估算–曲线重构”全流程的耦合优化仍须深化。

在工业检测、医疗仪器等领域对形变精度要求持续提升的背景下[5]，优化光纤传感曲线重建算法、提升重构可靠性，是推动该技术深度应用的重要组成部分。因此，本文试图从光纤传感器平面曲线重建的角度出发，在三次样条插值法与四阶龙格–库塔法的基础上，分析光纤传感技术在曲线重建中面临的实际问题，针对如何实现平面曲线精准重建探讨解决思路与方法。详细数据见第十六届“华中杯”学生数学建模挑战赛C题[6]。

2. 模型建立

2.1. 目标横坐标x处曲率插值估算模型

2.1.1. 求不同状态下各传感器处的曲率

本文所用波长测量数据来源于第十六届“华中杯”大学生数学建模挑战赛C题[6]，包含FBG1-FBG6六个传感点在两组初始状态下受力前后的波长信息，具体为初始状态1、测试1、初始状态2、测试2四组数据。

波长$\lambda$ 与曲线曲率$k$ 之间的关系近似为

$k=\frac{c\left(\lambda -{\lambda }_0\right)}{{\lambda }_0},$ (1)

其中，${\lambda }_0$ 表示初始波长，$\lambda$ 表示受力后的波长，$c$ 是给定的常数4200。

利用此公式并结合提供的波长数据可以将各个传感点测得的波长数据转化为曲率数据。对于状态$j\left(j=1,2\right)$ ，分别计算各个传感点$i\left(i=1,2,\cdots ,6\right)$ 的曲率$k_{ij}$ ：

$k_{ij}=\frac{4200\left({\lambda }_{ij}-{\lambda }_j\right)}{{\lambda }_j},$ (2)

其中，${\lambda }_{ij}$ 是第$i$ 个传感点在状态$j$ 下的受力后波长(表1中“测试$j$ ”列数据)，${\lambda }_j$ 表示在状态$j$ 下的初始波长(${\lambda }_1=1529\text{ nm},{\lambda }_2=1540\text{ nm}$ )。

计算出每个传感点的曲率，结果如表1传感点(FBG1-FBG6)的曲率所示：

Table 1. Curvatures of discrete sensing points FBG1-FBG6

表1. 离散传感点FBG1-FBG6的曲率

	FBG1	FBG2	FBG3	FBG4	FBG5	FBG6
初始状态1	2.21948986	2.21674297	2.23322433	2.23047744	2.23597122	2.22223676
初始状态2	2.98636364	2.97818182	2.97272727	2.98090909	2.98363636	2.97545455

2.1.2. 建立目标横坐标x处曲率插值估算模型

假设点FBG1的坐标为$\left(0,0\right)$ ，且过该点的切线与水平方向的夹角为45˚，假设有$n+1$ 个已知数据点$\left(x_i,k_i\right)$ ，其中$i\left(i=0,1,\cdots ,n\right)$ ，三次样条插值将数据拟合三次多项式，三次多项式可以表示为：

$S_i\left(x\right)=a_i+b_i\left(x-x_i\right)+c_i{\left(x-x_i\right)}^2+d_i{\left(x-x_i\right)}^3,$ (3)

其中，$x\in \left[x_i,x_{i+1}\right]$ ，$a_i=k_i$ ，需满足连续，一、二阶导数连续条件，确保曲线平滑性[7]。

2.2. 基于四阶龙格–库塔法的平面曲线重构模型

2.2.1. 推导微分方程

根据曲率公式考虑曲线整体呈上升趋势，即得一阶常微分方程组：

$\frac{\text{d}p}{\text{d}x}=k\left(x\right)\cdot \sqrt{{\left(1+p^2\right)}^3},$ (4)

代入前面得到的曲率函数得到两组完整的微分方程组。

2.2.2. 利用四阶龙格–库塔法重构平面曲线

通过四次斜率估计实现高精度数值积分，针对上述一阶微分方程组，定义两个待求函数

$f_1\left(x,p\right)=k\left(x\right)\cdot \sqrt{{\left(1+p^2\right)}^3}$ 和$f_2\left(x,p\right)=p$ 斜率估计值公式如下[8]：

对于$\frac{\text{d}p}{\text{d}x}$ 的斜率估计：

$\{\begin{array}{l}{k}_{1p}=h\cdot f_1\left(x_n,p_n\right)\\ k_{2p}=h\cdot f_1\left(x_n+\frac{h}{2},p_n+\frac{k_{1p}}{2}\right)\\ k_{3p}=h\cdot f_1\left(x_n+\frac{h}{2},p_n+\frac{k_{2p}}{2}\right)\\ k_{4p}=h\cdot f_1\left(x_n+h,p_n+k_{3p}\right)\end{array}.$ (5)

对于$\frac{\text{d}y}{\text{d}x}$ 的斜率估计：

$\{\begin{array}{l}{k}_{1y}=h\cdot f_1\left(x_n,p_n\right)\\ k_{2y}=h\cdot f_2\left(x_n+\frac{h}{2},p_n+\frac{k_{1p}}{2}\right)\\ k_{3y}=h\cdot f_2\left(x_n+\frac{h}{2},p_n+\frac{k_{2p}}{2}\right)\\ k_{4y}=h\cdot f_2\left(x_n+h,p_n+k_{3p}\right)\end{array}.$ (6)

迭代更新：

$\{\begin{array}{l}{p}_{n+1}=p_n+\frac{1}{6}\left(k_{1p}+2k_{2p}+2k_{3p}+k_{4p}\right)\\ y_{n+1}=y_n+\frac{1}{6}\left(k_{1p}+2k_{2p}+2k_{3p}+k_{4p}\right)\\ x_{n+1}=x_n+h\end{array},$ (7)

其中，$h$ 为步长，$x_n$ 、$p_n$ 分别为第$n$ 步的横坐标和切线斜率，通过分步计算斜率，有效捕捉函数在微小区间内的变化趋势，降低截断误差。

2.3. 基于等间距弧长采样与曲率反推的平面曲线重构模型

2.3.1. 求所有采样点的曲率值

要求出在平面曲线$y=x^3+x\left(0\le x\le 1\right)$ 上的等间距弧长采样点的曲率值[9]，首先要计算出此平面曲线方程的弧长；然后定义采样点的个数进而可求出等间距弧长的间隔，找出各个采样点位置，即知晓各个采样点的$x$ 值；最后根据曲率计算公式$K=\frac{\left|y^{″}\right|}{\sqrt{{\left(1+{y^{\prime }}^2\right)}^3}}$ ，通过求出各个采样点的一阶导数和二阶导数值，代入公式即可求出曲率值。

Part 1：计算弧长

由平面曲线方程$y=x^3+x\left(0\le x\le 1\right)$ ，根据弧长公式$s=\sqrt{1+{y^{\prime }}^2}$ 可求出平面曲线的弧长。

Part 2：求出各个采样点横坐标值

首先定义采样点的个数为10，则等间距弧长为10，然后利用编程公式求出每段等间距弧长，最后可得出所有采样点横坐标为0.145，0.281，0.404，0.513，0.608，0.692，0.767，0.834，0.895，0.951。

Part 3：计算各个采样点的曲率值

根据计算曲率公式$K=\frac{\left|y^{″}\right|}{\sqrt{{\left(1+{y^{\prime }}^2\right)}^3}}$ 得，求曲率值要计算出各个采样点的一阶导数值和二阶导数值，再

代入公式即可求出。

2.3.2. 重构曲线

首先，根据平面曲线方程$y=x^3+x\left(0\le x\le 1\right)$ ，利用求出的各个采样点的横坐标值，求出各个采样点的$y$ 值。然后累计弧长，可先通过计算相邻点之间的欧几里得距离作为弧长增量的近似值，再进行弧长的累计。其次计算切线方向，利用求出的曲率值，用弧长乘以曲率即可计算出切线方向的变化量，实现切线方向的更新。最后通过计算出的坐标进行曲线的重构。

Part 1：计算各个采样点的$y$ 值

已知平面曲线方程$y=x^3+x\left(0\le x\le 1\right)$ ，将求出的各个采样点的横坐标值代入平面曲线中，可求出各个采样点的$y$ 值。

Part 2：累计弧长

依题意得起始点的坐标为$\left(0,0\right)$ ，首先计算相邻点之间的欧几里得距离作为弧长增量的近似值，然后对弧长进行累计。

$\text{d}x=x_i-x_{i-1},$ (8)

$\text{d}y=\text{d}{y}_i-\text{d}{y}_{i-1},$ (9)

$\text{d}s=\sqrt{\text{d}{x}^2-\text{d}{y}^2}.$ (10)

Part 3：计算切线方向

利用求出的曲率值，然后用弧长乘以曲率计算出切线方向的变化量$\Delta \theta =\text{d}s\cdot k_{i-1}$ ，由此实现切线方向的更新。

Part 4：曲线重构

先用切线方向和弧长增量来计算点的位置变化，然后将计算出来的坐标进行对曲线的重构。

$\text{d}x=\cos \theta \cdot \text{d}s,$ (11)

$\text{d}y=\sin \theta \cdot \text{d}s.$ (12)

3. 模型求解与分析

3.1. 目标横坐标x处曲率插值估算模型

利用$$ \mathrm{P}\mathrm{y}\mathrm{t}\mathrm{h}\mathrm{o}\mathrm{n}$$程序计算出横坐标$x$ 轴相应位置处的曲率，结果如表2所示：

Table 2. Curvatures from cubic spline interpolation

表2. 三次样条插值曲率

横坐标$x$ (米)	0.3	0.4	0.5	0.6	0.7
测试1曲率$k$	2.21198169	2.21257402	2.21425324	2.21674297	2.21976681
测试2曲率$k$	2.98317045	2.98158474	2.97989127	2.97818182	2.97654812

利用$$ \mathrm{S}\mathrm{P}\mathrm{S}\mathrm{S}$$软件进行拟合，并用$R^2$ 值检验三次样条插值的拟合效果，见图1、图2：

Figure 1. Verification of cubic spline interpolation for test 1

图1. 测试1三次样条插值检验

Figure 2. Verification of cubic spline interpolation for test 2

图2. 测试2三次样条插值检验

由图1、图2可知：在测试1中，对三次样条插值模型$y=2.22-0.05x+0.11x^2-0.05x^3$ 进行检验得到的相关系数$R^2$ 近似为1；在测试2中，对三次样条插值模型$y=2.99-4.9\times {10}^{-3}x-0.02x^2+0.02x^3$ 进行检验得到的相关系数$R^2$ 近似为1；三次样条插值在保证平滑性的前提下，拟合精度高。

3.2. 基于四阶龙格–库塔法的平面曲线重构模型

通过四阶龙格–库塔法对测试1和测试2平面曲线重构，图中红色线表示三次样条插值曲线，蓝色线表示四阶龙格–库塔法重构曲线，结果见图3、图4。

测试1：基于曲率函数$k\left(x\right)=2.22-0.05x+0.11x^2-0.05x^3$ 。

Figure 3. Reconstruction of plane curve based on fourth-order Runge-Kutta method for test 1

图3. 测试1基于四阶龙格–库塔法的平面曲线重构

测试2：基于曲率函数$k\left(x\right)=2.99-0.0049-0.02x^2+0.02x^3$ 。

Figure 4. Reconstruction of plane curve based on fourth-order Runge-Kutta method for test 2

图4. 测试2基于四阶龙格–库塔法的平面曲线重构

综上两组曲线共性为单调上升、平滑连续，个性差异源于曲率大小：测试2曲率更大，垂直位移更大、弯曲更复杂，符合光纤形变物理规律。数学层面，四阶龙格–库塔法全局截断误差为$O\left(\text{d}{x}^4\right)$ ，步长0.001下误差可控在${10}^{-12}$ 量级，保障拟合精度；物理层面，重构曲线形态与光纤受力形变逻辑一致，曲率越大形变越显著。对比可知，初始状态2的初始波长大于状态1，导致测试2曲率更大，验证了模型对初始波长变化的准确响应。

3.3. 基于等间距弧长采样与曲率反推的平面曲线重构模型

用Python软件求出各采样点曲率值为0.280，0.419，0.420，0.358，0.287，0.227，0.181，0.147，0.120，0.100。最后，利用SPSS软件对采样点的横坐标及其曲率绘制散点图(见图5)。

Figure 5. Curvature scatter plot

图5. 曲率散点图

图5为曲率散点图，呈现了等间距弧长采样点的横坐标及其对应的曲率值分布。从图中可以看到，曲率值随采样点横坐标的变化呈现出一定的波动特征。在横坐标较小的区域，曲率值有上升趋势，而在横坐标较大的区域，曲率值又有下降趋势，反映出曲线在该段的弯曲程度变化情况。这种曲率的波动体现了原曲线$y=x^3+x\left(0\le x\le 1\right)$ 弯曲特性的非均匀性，即曲线不同位置的弯曲程度存在差异，等弧长采样能够捕捉到这种曲率的变化规律。

重构平面曲线后，将重构曲线图像与原平面曲线图像进行对比，如图6。

Figure 6. Comparison diagram of the original plane curve and the reconstructed plane curve

图6. 原平面曲线与重构平面曲线的对比图

图6原平面曲线与重构平面曲线的对比图，展示了原平面曲线$y=x^3+x\left(0\le x\le 1\right)$ 与通过等弧长采样、曲率计算和坐标重构得到的重构曲线的对比。可以观察到，重构曲线(橙色线)与原曲线(蓝色线)整体趋势较为一致，都呈现出从原点开始单调上升的形态，这表明基于等弧长采样和曲率的曲线重构方法能够大致还原原曲线的走势。然而，重构曲线与原曲线之间也存在一定的偏差，采样点(红色点)分布在两条曲线之间，反映出重构过程中存在误差。

3.4. 重构曲线与原曲线之间的误差分析

(1) 采样误差：等弧长采样本身就是一个近似过程，实际采样点可能并不完全满足等弧长条件。

(2) 离散数据误差：选取的是离散点集，而曲线本身是连续的。这种离散化本身就引入了误差，因为原始曲线在离散点之间的具体形状并未被完全捕获。

(3) 曲率近似误差：在实际计算中，我们假设在相邻点之间的弧长上曲率是恒定的，这种假设引入了曲率近似误差。

(4) 数值计算误差：由于计算机内部表示浮点数的方式和三角函数的计算精度限制，这些操作可能引入微小的误差。

4. 结束语

本文主要研究基于光纤传感器的平面曲线重建算法建模问题，构建了以三次样条插值与四阶龙格–库塔法为核心的技术模型，为光纤传感数据的曲线重构提供了切实可行的方案。通过整合离散波长数据，依托数值方法实现连续曲率估算与曲线重建，有效提升了平面曲线重构结果与实际形态的贴合度，强化了光纤传感技术在形变监测场景的应用可靠性。该模型还可拓展至多芯光纤的三维形态重构、复杂结构的形变轨迹追踪等领域，为航空航天、精密制造、土木工程、医疗仪器设计等多个场景提供精准的形状重建支持，为相关技术的实际应用与发展奠定基础。

基金项目

2025年度湖南省大学生创新训练计划一般项目(S202510551066)。

NOTES

^*通讯作者。