前言
上学期刚上完《激光大气探测》这门儿课程，课程报告要求是对激光雷达得到的数据进行处理，最终得到回波强度、距离平方矫正信号和消光系数随距离的变化曲线。报告是用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，如果有朋友发现上述程序的问题以及有什么更好的方法来处理，超级欢迎一起交流沟通！！！
扫一扫
3306
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
