光学陀螺仪标度因数的初步仿真（使用matlab）-CSDN博客

本文链接：https://blog.csdn.net/2403_89560913/article/details/147402506

0.写在前面

作者今年正在参加大学生科研训练（SRTP），我们的项目是关于光学陀螺仪的。最近做了一个光学陀螺仪标度因数的初步仿真，写这篇文章的目的就是为了再次整理思路，也是为了分享给我的队友们。由于作者现在还是大一，并且目前并不是相关专业的，所以仿真内容难免会有不准确甚至错误的地方，希望能够得到读到这篇文章的好心人的指正。

1.基本思路

光学陀螺仪的原理基于萨格纳克效应，就是在同一个环形光路中，两束传播方向相反的光会因为光路的旋转而产生相位差，监测此相位差就能够计算出旋转角速度。作者在此就不写具体的推导了，~~因为我不会，~~只给出计算的公式： $\Delta \varphi = \frac{4\pi A\Omega }{\lambda c}$ ，其中Δφ是两束光的相位差，A为光路面积，Ω是角速度，λ是光的波长（一般采用波长为1310nm和1550nm的光），c是光速。

在实际情况下，我们能检测到的相位差Δφ会受到噪声的影响，导致我们实际检测出来的相位差是Δφ(noise) = Δφ(ideal) + noise.（这里用英国话不是在放洋屁，是为了对接后面我代码中的写法，勿喷）

现在，我们需要把检测出来的相位差转化成信号输出，这里假设输出的信号是电压，那么对于这个输出信号，可以通过这个公式计算： $I=I_{0}cos\Delta \varphi$ ，这里 $I$ 代表输出信号， $I_{0}$ 代表光强。但是实际情况下，输出信号也会受到噪声干扰，因此结合上面的内容，我们最终得到的输出信号为： $I(output) = I_{0}cos(\Delta \varphi (ideal) + noise1) + noise2$ .

那么，所谓的标度因数(scale factor)，其实就是输出信号和角速度的比值，这里标度因数用K来表示，就可以写成 $I(output)=K\Omega$ .为什么要计算这个标度因数？因为如果我们没有掌握准确的标度因数，就算输出的信号再准确也很难求出精确的角速度。

现在，我们就可以理一下基本思路了：使用matlab模拟噪声，进而模拟出带有噪声的输出信号，对此信号进行滤波操作，再对滤波后的信号进行线性拟合，从而求出标度因数K。这个模拟可以做成百上千次，然后求出标度因数K的概率分布（也就是所谓的蒙特卡洛仿真），进而得到K的更具体的信息。

2.代码实现

首先是基本参数的设定：

%% 参数设置
lambda = 1550e-9;    %波长λ=1550nm
A = 0.01;    %光路面积A=0.01m^2
c = 3e8;    %光速c = 3e8m/s
Omega = linspace(-10,10,100);    %用一个向量来存储-10到10之间的不同角速度

其中linspace 是一个用于在指定的起点和终点之间生成均匀分布的数值序列的函数，返回一个向量，这里就储存了从-10到10rad/s之间100个均匀分布的样本点。

然后我们先对相位差Δφ进行噪声处理：

%% 相位差的噪声处理
delta_phi_ideal = ((4 * pi * A) / (lambda * c)) .* Omega;    %公式计算出的理想相位差
noise_random_walk = cumsum(0.01 * randn(size(Omega)));    %随机游走噪声，噪声强度0.01(可修改)
noise_white = 0.01 * randn(size(Omega));    %白噪声，噪声强度0.01(可修改)
delta_phi_noise = delta_phi_ideal + noise_white + noise_random_walk;    %最终得到的噪声影响后的相位差

这里作者对理想相位差施加了两种不同的噪声，一种是白噪声，通常来源于温度影响，一种是随机游走噪声，来源于电子器件的影响。可以先不管两种噪声的生成原理，只需要知道这个相位差已经被噪声干扰过了就行了！

接着我们计算噪声干扰后的输出信号：

%% 输出信号&噪声处理
I0 = 1;    %光强，设置为1
sigma_I = 0.05;    %信号输出噪声标准差(噪声强度?)
I_output = I0 * cos(delta_phi_noise) + sigma_I * randn(size(Omega));

这里的输出信号带有的噪声采用了白噪声的模型，也就是代码中的sigma_I * randn(size(Omega))，sigma_I可以看作是噪声强度（应该吧）。

现在，我们对得到的输出信号进行滤波处理：

%% 巴特沃斯滤波器
fc = 1;    %截止频率1Hz
fs = 10;    %采样频率10Hz
[b,a] = butter(2,fc / (fs / 2));    %二阶巴特沃斯低通滤波器
% delta_phi_filtered = filtfilt(b,a,delta_phi_noise);    %零相位滤波
I_filtered = filtfilt(b,a,I_output);

这里我们采用的是巴特沃斯滤波器，背后的原理作者自己也不是很明白了，总之滤波后的信号就是变得丝滑一些了。（也可以对噪声干扰后的相位差进行滤波处理，被注释掉的那一行就是）

最后，我们对滤波后的输出信号进行线性拟合：

%% 信号输出曲线拟合
p = polyfit(Omega,I_filtered,1);    %采用一阶多项式拟合
I_fitted = polyval(p,Omega);    %得到拟合后的信号输出
scale_factor = p(1);    %一阶多项式的系数->标度因数

其中，polyfit函数的用法如下：polyfit(x,y,n)，就是将y和x之间的函数关系用n阶多项式拟合，返回的是这个多项式的系数组成的向量。这里我们采用一阶多项式拟合，也就是线性拟合，polyval函数就是把自变量带入到拟合多项式中算出拟合结果，这里自变量是角速度Ω。显然，所求多项式的一次项系数就是我们要找的标度因数。

到这里，我们就可以对某一个仿真情况进行作图分析了：

%% 作图
figure;
plot(Omega,I_output,"r--","LineWidth",1.5,"DisplayName","实际信号输出");
hold on;
plot(Omega,I_filtered,"m-.","LineWidth",2,"DisplayName","滤波器处理后的信号输出");
hold on;
plot(Omega,I_fitted,"b","LineWidth",1.5,"DisplayName","拟合信号输出");
xlabel("角速度Ω/(rad/s)");
ylabel("信号强度");
title("不同角速度下的信号输出(某一次的测试结果)");
legend("show");
grid on;

得到的某一次仿真的结果如下：

通过拟合信号输出，我们可以计算这一次仿真下的标度因数为多少。

但是这样的数据我们需要计算成百上千次，然后才能得到它的概率分布，所以，我们需要采用一个循环来进行这一操作：

%% 用蒙特卡洛仿真求标度因数的概率分布
n_simu = 1000;    %进行1000次计算
scale_factors = zeros(n_simu,1);    %将每一次的计算结果储存在一个1000行，1列的矩阵中
for i = 1:n_simu
    %% 以下末尾带M的变量都是为了多次计算专门设置的，与上面的变量没有关系
    %% 相位差的噪声处理
    noise_random_walk_M = cumsum(0.01 * randn(size(Omega)));    %随机游走噪声，噪声强度0.01(可修改)
    noise_white_M = 0.01 * randn(size(Omega));    %白噪声，噪声强度0.01(可修改)
    delta_phi_noise_M = delta_phi_ideal + noise_white_M + noise_random_walk_M;

    %% 输出信号&噪声处理
    % I0 = 1;    %光强，设置为1
    % sigma_I = 0.05;    %信号输出噪声标准差(噪声强度?)
    I_output_M = I0 * cos(delta_phi_noise_M) + sigma_I * randn(size(Omega));

    %% 巴特沃斯滤波器
    % fc = 1;    %截止频率1Hz
    % fs = 10;    %采样频率10Hz
    % [b,a] = butter(2,fc / (fs / 2));    %二阶巴特沃斯低通滤波器
    % % delta_phi_filtered = filtfilt(b,a,delta_phi_noise);    %零相位滤波
    I_filtered_M = filtfilt(b,a,I_output_M);

    %% 信号输出曲线拟合
    p_M = polyfit(Omega,I_filtered_M,1);    %采用一阶多项式拟合
    % I_fitted = polyval(p,Omega);    %得到拟合后的信号输出
    scale_factor_M = p_M(1);    %一阶多项式的系数->标度因数
    scale_factors(i) = scale_factor_M;    %把单次计算的结果储存在scale_factors矩阵中

end

我在每个每个变量的末尾添加了M，代表这些都是参与到循环当中的变量，与上面的变量是没有关系的。这一块代码的原理和之前的一模一样，可以看出来我是把上面的代码CV下来的，然后做了注释掉了一些重复的部分，改了一下变量名。

然后再作图和计算，求出标度因数的概率分布，均值和标准差：

figure;
histogram(scale_factors,50);
xlabel("标度因数");
ylabel("概率分布");
title("标度因数的概率分布");
mean_scale_factor = mean(scale_factors);    %计算标度因数的均值
std_scale_factor = std(scale_factors);    %计算标度因数分布的标准差
fprintf("标度因数的均值为: %.6f\n",mean_scale_factor);
fprintf("标度因数分布的标准差为: %.6f\n",std_scale_factor);

这里histogram的作用是生成50个区间段的直方图。

某一次仿真中得到的结果如下：