激光雷达方程应用（使用MATLAB和Python语言得到回波强度、距离平方矫正信号、消光系数随距离的变化曲线）

---

前言

上学期刚上完《激光大气探测》这门儿课程，课程报告要求是对激光雷达得到的数据进行处理，最终得到回波强度、距离平方矫正信号和消光系数随距离的变化曲线。报告是用MATLAB写的，假期又用Python尝试了一下，写出来还有点儿成就感，在此分享。对于每一步在此不做具体介绍，相信大家看一下注释就可以明白。

---

程序中所使用的txt文件数据格式为：
第一列：距离
第二列：回波信号强度

一、MATLAB语言程序内容

%% 本程序用于绘制出激光雷达Pr(回波强度)、Pr_rr(距离平方矫正信号)、alpha_r(消光系数)三者分别随距离的变化曲线
%   为了减少因噪声产生的误差，根据自行判断的r_m数值，可在程序的数据读取部分以及(3)中自行设置r_m以及rs(range_start)的值
close all;
clear;
clc;
Data = load("F:\Study\20181222105152.txt");           %数据第一列为距离/Km 第二列为回波强度
Prmin = min(Data(:,2));                      %得到因噪声造成的回波信号的最大负值数据
r_m = 1.095;                                 %设置一个合适的r_m，需是0.0075的倍数（由于本程序第三部分采用求和计算积分的方...
                                             ...法，所以对于点数据，在设置时需要满足r_m在数据中，即0.0075的倍数）
i= round(r_m/0.0075);                        %得出规定的r_m所对应的序号
Data(1:i,2) = Data(1:i,2)+abs(Prmin);        %从初始到r_m的回波信号均加上负最大回波信号的绝对值

%% （1）得到回波信号Pr随距离的变化曲线
figure(1)
plot(Data(1:i,1),Data(1:i,2));
title("回波强度随距离的变化曲线");
xlabel("Distance/Km");
ylabel("Pr");
%% （2）得到距离矫正信号随距离的变化曲线
figure(2)
Pr_rr = Data(1:i,1).^2.*Data(1:i,2);        %距离平方矫正信号:回波信号乘以对应距离的平方
plot(Data(1:i,1),Pr_rr);
title("距离矫正信号随距离的变化曲线");
xlabel("Distance/Km");
ylabel("距离平方矫正信号Pr*r^2");

%% （3）计算消光系数alpha_r
Pr_r_m = Data(i,2);                         % r_m所对应的回波强度
c1 = r_m^2*Pr_r_m;
%% 公式一 alpha_r = exp((S-Sm)/k)/(alfa_m^(-1)+2*int(exp(S-Sm)/k),r,r_m)/k)
     %% S(r) = log(Pr*r^2)
     %% Sm = S(r_m）
     %% alpha_m = alpha(r_m)
     %% k = 1   
     %% alpha_m = 10e-5
%% 原公式一可化简为 alpha_r = r^2*Pr/(c1*(10e5+2*int((r^2*Pr/c1),r,r_m)))
 rs = 0.0075;                               %自行设置初始rs的值，此处从第一个开始  
 rss = round(rs/0.0075);
 r0 = Data(rss:1:i,1);                           
 r = r0(end:-1:rss);                        %为方便下面cumtrpz的使用，将r0中的数据倒置
 Pr0 = Data(rss:1:i,2);
 Pr = Pr0(end:-1:rss);                      %同样将Pr0中的数据倒置
 ji = r.*r.*Pr;
 jic = ji/c1;                               %(r^2*Pr/r_m^2*Pr)
 C = cumtrapz(jic);                         %得到从第i个向前积分的每一个积分值
 c3 = 2*C(end:-1:1);                        %得到正序的积分值
 
 
Alpha_r = @(r,Pr,c3) r.^2.*Pr./(c1.*(10e5+c3));       %将Alpha_r表示为r和Pr的函数

%% （4）得到消光系数随距离的变化曲线
figure(3)
alpha_r = Alpha_r(r,Pr,c3);          
plot(r,alpha_r);
title("消光系数随距离的变化曲线");
xlabel("Distance/Km");
ylabel("消光系数");

得到的三张结果图如下：

二、Python语言程序内容

## 本程序用于绘制出激光雷达Pr（回波强度）、Pr_rr（距离平方矫正信号）、alpha_r(消光系数) 三者分别随距离的变化曲线
import numpy as np
import matplotlib.pyplot as plt
import pylab as pl
r=[]
Pr=[]
with open('F:\\Study\\20181222105152.txt','r') as data:
    for line in data:
        ch = line.split()
        r.append(float(ch[0]))
        Pr.append(float(ch[1]))
Prmin = min(Pr)                                         #得到回波强度最小值
Pr = (np.array(Pr)-(Prmin-0.0001)).tolist()             #由于所有回波强度实际上都应该是正数，所以针对所有的回波强度值都减去最小值（此处的-0.0001，可以使得数据中不出现回波强度为零的情况）
c = (np.array(r)*np.array(Pr)).tolist()
Pr_rr = (np.array(r)*np.array(r)*np.array(Pr)).tolist() #得到所有距离平方矫正信号

r_m = 1.095                                             #设置一个参考距离，视从雷达到该距离内的信噪比满足数据要求
num =  r.index(r_m)                                     #得到该距离处所对应的位置
Pr_r_m = Pr[num]                                        #得到参考位置出的回波强度
c1 = r_m*r_m*Pr_r_m                                     #得到参考位置出的距离平方矫正信号

r = r[0:num]                                            #选取雷达到参考位置
Pr = Pr[0:num]
Pr_rr = Pr_rr[0:num]

## 公式一 alpha_r = exp((S-Sm)/k)/((alpha_m)^(-1)+2*int((exp(S-Sm)/k),r,r,m)/k)
 # S(r) = log(Pr*r^2)
 # Sm = S(r_m)
 # alpha_m = alpha(r_m)
 # k = 1
 # alpha_m = 10e-5
## 原公式一化简为 alpha_r = r^2*Pr/(c1*(10e5+2*int((r^2*Pr/c1),r,r_m)))

a = ((np.array(Pr_rr))/c1).tolist()
a = a[::-1]                                            #为了便于使用cunsum函数，cunsum可以求list任意一个数和其之前所有数值之和（化积分为求和）
c2 = np.cumsum(a)
c3 = c1*(10e5+2*c2)

alpha_r = Pr_rr/c3
 
figure = plt.figure(1)
plt.plot(r,Pr,'red')
plt.title('Pr-r')
plt.xlabel('r/Km')
plt.ylabel('Pr')

ax = pl.gca()                                         # 获取当前图像的坐标轴信息
ax.xaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位，1是指科学计数法时的位数
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位，1是指科学计数法时的位数


figure = plt.figure(2)
plt.plot(r,Pr_rr,'blue')
plt.title('Pr_rr-r')
plt.xlabel('r/Km')
plt.ylabel('Pr*r^2')

ax = pl.gca()  
ax.xaxis.get_major_formatter().set_powerlimits((0,1)) 
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) 


figure = plt.figure(3)
plt.plot(r,alpha_r,'green')
plt.title('alpha_r-r')
plt.xlabel('r/Km')
plt.ylabel('alpha_r')

ax = pl.gca()  
ax.xaxis.get_major_formatter().set_powerlimits((0,1)) 
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) 


figure.show()

得到的三张结果图如下：

总结

对比两种方法，MATLAB处理这种数据会更加方便，结果图也更加美观。Python绘图title不能包含中文字符这就很难受了。两个程序内容的重点就在于 cumtrapz和cumsum两个函数的使用，使得积分化和简单了不少。okk，如果有朋友发现上述程序的问题以及有什么更好的方法来处理，超级欢迎一起交流沟通！！！