2023年电赛F题——基于声传播的智能定位系统（基础部分（包括定位算法））

前言

       2023电赛也是刚刚结束，写这篇blog也是记录一下自己第一次电赛。由于时间和能力有限，奋战四天三夜暂时也是只做完基础部分，不过好在精度还可以，基础部分实现的效果也还不错。时间在方案选择上浪费了比较多，也试了很多错，和平时去做以往的电赛题有很大区别。

1 题目要求

 2 方案选择

2.1 总体方案

      系统总体框图如下图所示。

      各部分选择：

      电源：12V锂电池供电接DC-DC直流降压模块输出5V信号给单片机供电。

      发射装置：扬声器（喇叭）

      接收装置：压电陶瓷片

 

2.2 具体选型

      MCU：我们选择的是STM32F103ZET6的核心板，接口比C8T6要多，而且自带的ADC采样的频率足够我们采集声音信号，获取波峰时间。

     显示/交互模块：直接用的手头上有的LCD屏幕，正点原子的TFT LCD 4.3英寸带触控功能，本题选OLED即可，只是我们没有OLED。看个人习惯。

 

     下面是重头戏，也是我们一开始两天都在纠结的，用什么做发射，用什么做接收。

     发射装置：我们试过蜂鸣器/麦克风（带咪头）/超声波换能器/压电陶瓷片，通过亚克力板发射接收信号比较难以接收到，拿压电陶瓷片来说，单片机产生的步进信号需要经过双电源12V供电的大功率功放放大才能传入亚克力板中在接受端接收到信号。最终我们采用了大功率的广播喇叭的核心扬声器部分当作发射装置。

      接收装置：使用压电陶瓷片作为接收信号的装置，优点是耦合度较为一致，方便信号接收时进行采集与测量，缺点是压电陶瓷片作为接收装置时，对接收的信号有较大的能量衰减， 故接收信号的电压幅值较小，需要另行放大再进行采样。需要外加足够功率的功率放大器才能使用，而且最好不要直接通过硬性胶水连接。

 

      连接问题：选用耦合剂（我们采用的是超声波耦合剂，正常医院里做B超用的那种），凡士林效果也不错而且有一定粘性。

2.3 硬件搭建

 3 采集信号分析

3.1 单路信号

        通过压电陶瓷片接收信号再经功放放大，接入示波器进行分析，主要分析敲击产生的声音振动波形信号。我们实际采集到的单路信号的波形如下图所示。

         通过单片机对该信号做AD采样与分析即可获得该信号。该信号的处理有如下两种方法。

       （1）外部中断接收，处理记录当时的时间。

       （2）在定时器中断中处理。

        我们采用的是定时器中断记录时间进行处理。

3.2 多路信号

      同3.1，接入两路压电陶瓷片测试波形，由于敲击点距离两路陶瓷片距离不同，信号会有到达时间差，通过这个时间差，我们可以对应该距离的差值从而计算出具体敲击源的位置。

 

 4 软件算法及分析

4.1 定位模型分析——基于二维TDOA定位算法

       根据题目要求，设置四个区域Z的其中三个区域为信号接收端，各放置一个压电陶瓷片接收来自人为敲击的同一信号源，建立如下图所示模型。

        在上图模型中，设定 为设定的三个接收点， 为敲击信号注入点，由数学关系可以得到：

 

 

 联立上述方程组并化简得到下式（矩阵形式表示）：

       根据亚克力板本身固有的特性，点产生的敲击信号在亚克力板中传播的速度是稳定的，所以敲击信号传播到各点耗费的时间是不同的，故有信号到达亚克力板上任意两点的时间差与传播速度的乘积等于这两点到敲击信号点的距离差，故有下式成立：

                                                          

       其中为声音在亚克力板中的传播速度。

4.2 TDOA算法MATLAB仿真

      关于TDOA的算法网上有很多，找一份自己仿真一下就差不多能得出结果。

clc;
%% 初始化数值


x1 = 42;  y1 = 3;           % 传感器1

x2 = 42;  y2 = 42;          % 传感器2
x3 = 3;   y3 = 42;          % 传感器3
 

s1 = 0;    % 时延距离 s1 = r21
s2 = 0;    % 时延距离 s2 = r31


%% 求解未知点
k1 = x1^2 + y1^2;   % 中间值
k2 = x2^2 + y2^2;
k3 = x3^2 + y3^2;

p1_molecule = (y2 - y1)*s2^2 - (y3 - y1)*s1^2 + (y3 - y1)*(k2 - k1) - (y2 - y1)*(k3 - k1);
p1_denominator = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
p1 = p1_molecule / (p1_denominator*2);
q1 = ((y2 - y1)*s2 - (y3 - y1)*s1)/((x2 - x1)*(y3 - y1)-(x3 - x1)*(y2 - y1));

p2_molecule = (x2 - x1)*s2^2 - (x3 - x1)*s1^2 + (x3 - x1)*(k2 - k1) - (x2 - x1)*(k3 - k1);
p2_denominator = (x3 - x1)*(y2 - y1) - (x2 - x1)*(y3 - y1);
p2 = p2_molecule / (p2_denominator * 2);
q2 = ((x2 - x1)*s2 - (x3 - x1)*s1)/((x3 - x1)*(y2 - y1)-(x2 - x1)*(y3 - y1));
%% 求取方程
a = q1^2 + q2^2 -1;
b = -2*((x1 - p1)*q1 + (y1 - p2)*q2);
c = (x1 - p1)^2 + (y1 - p2)^2;

r1 = (-b + sqrt(b^2 - 4*a*c))/(2*a);
r2 = (-b - sqrt(b^2 - 4*a*c))/(2*a);

if r1 > 0
    x = p1 + q1*r1;
    y = p2 + q2*r1;
else
    x = p1 + q1*r2;
    y = p2 + q2*r2;
end

disp('未知点坐标: ')
disp(x), disp(y)

      根据上述分析，只要测得三传感器接收点的接收信号时间差，即可通过传播速度解得最终位置,将仿真代码移植为C语言程序，即可通过单片机接收时间差计算坐标获得结果。

     注意：传播速度的设定会影响计算，所以建议通过输入固定的最终位置坐标和测得的时间差计算得到该值。

 void calculate_loc(void) {
    double c=c1;     //  cm/us      亚克力板声速为2600-2700m/s
    double t3;
    double t4;

    double x1, y1, x2, y2, x3, y3;
    double s1;
    double s2,k1,k2,k3;
    double p1_molecule, p1_denominator, p1, q1;
    double p2_molecule, p2_denominator, p2, q2;
    double a, b;
    double discriminant;
    double r1, r2;
	  t3 = f_sec2-f_sec1;
    t4 = f_sec3-f_sec1;
    x1 = -19, y1 = -19;
    x2 = 19, y2 = -19;
    x3 = 19, y3 = 19;
    s1 = c*t3;
    s2 = c*t4;
    k1 = x1 * x1 + y1 * y1;                                                                                                                                                                                
    k2 = x2 * x2 + y2 * y2;
    k3 = x3 * x3 + y3 * y3;

    p1_molecule = (y2 - y1) * s2 * s2 - (y3 - y1) * s1 * s1 + (y3 - y1) * (k2 - k1) - (y2 - y1) * (k3 - k1);
    p1_denominator = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
    p1 = p1_molecule / (p1_denominator * 2);
    q1 = ((y2 - y1) * s2 - (y3 - y1) * s1) / ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1));

    p2_molecule = (x2 - x1) * s2 * s2 - (x3 - x1) * s1 * s1 + (x3 - x1) * (k2 - k1) - (x2 - x1) * (k3 - k1);
    p2_denominator = (x3 - x1) * (y2 - y1) - (x2 - x1) * (y3 - y1);
    p2 = p2_molecule / (p2_denominator * 2);
    q2 = ((x2 - x1) * s2 - (x3 - x1) * s1) / ((x3 - x1) * (y2 - y1) - (x2 - x1) * (y3 - y1));

    a = q1 * q1 + q2 * q2 - 1;
    b = -2 * ((x1 - p1) * q1 + (y1 - p2) * q2);
    c = (x1 - p1) * (x1 - p1) + (y1 - p2) * (y1 - p2);


    discriminant = b * b - 4 * a * c;
    r1 = (-b + sqrt(discriminant)) / (2 * a);
    r2 = (-b - sqrt(discriminant)) / (2 * a);


    if (r1 > 0) {
        x_1 = p1 + q1 * r1;
        y_1 = p2 + q2 * r1;
    }
    else {
        x_1 = p1 + q1 * r2;
        y_1 = p2 + q2 * r2;
    }

}

4.3 单片机扫频信号产生

4.3.1 定时器初始化

       

void TIM3_Int_Init(u16 arr, u16 psc)
{
//	
	TIM_TimeBaseInitTypeDef TIM_TimeBaseInitStrue; //定义一个定时中断的结构体
	NVIC_InitTypeDef NVIC_InitStrue; //定义一个中断优先级初始化的结构体
	
  RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM3, ENABLE); //使能通用定时器2时钟
	
	TIM_TimeBaseInitStrue.TIM_Period=arr; //计数模式为向上计数时，定时器从0开始计数，计数超过到arr时触发定时中断服务函数
	TIM_TimeBaseInitStrue.TIM_Prescaler=psc; //预分频系数，决定每一个计数的时长
	TIM_TimeBaseInitStrue.TIM_CounterMode=TIM_CounterMode_Up; //计数模式：向上计数
	TIM_TimeBaseInitStrue.TIM_ClockDivision=TIM_CKD_DIV1; //一般不使用，默认TIM_CKD_DIV1
	TIM_TimeBaseInit(TIM3, &TIM_TimeBaseInitStrue); //根据TIM_TimeBaseInitStrue的参数初始化定时器TIM2
	
	TIM_ITConfig(TIM3, TIM_IT_Update, ENABLE); //使能TIM2中断，中断模式为更新中断：TIM_IT_Update
	
	NVIC_InitStrue.NVIC_IRQChannel=TIM3_IRQn; //属于TIM2中断
	NVIC_InitStrue.NVIC_IRQChannelCmd=ENABLE; //中断使能
	NVIC_InitStrue.NVIC_IRQChannelPreemptionPriority=1; //抢占优先级为1级，值越小优先级越高，0级优先级最高
	NVIC_InitStrue.NVIC_IRQChannelSubPriority=1; //响应优先级为1级，值越小优先级越高，0级优先级最高
	NVIC_Init(&NVIC_InitStrue); //根据NVIC_InitStrue的参数初始化VIC寄存器，设置TIM2中断
	
	TIM_Cmd(TIM3, ENABLE); //使能定时器TIM2
}

        这部分内容简单，定时器初始化之后配置好分频即可产生对应频率的信号。

4.4 ADC采样

4.4.1 ADC初始化

void  Adc_Init(void)
{ 	
	ADC_InitTypeDef ADC_InitStructure; 
	GPIO_InitTypeDef GPIO_InitStructure;

	RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA |RCC_APB2Periph_ADC1 |RCC_APB2Periph_GPIOB	, ENABLE );	  //使能ADC1通道时钟
 

	RCC_ADCCLKConfig(RCC_PCLK2_Div6);   //设置ADC分频因子6 72M/6=12,ADC最大时间不能超过14M

	//PA1 作为模拟通道输入引脚                         
	GPIO_InitStructure.GPIO_Pin = GPIO_Pin_1|GPIO_Pin_3|GPIO_Pin_2;
	GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;		//模拟输入引脚
	GPIO_Init(GPIOA, &GPIO_InitStructure);	

//		//PA1 作为模拟通道输入引脚                         
//	GPIO_InitStructure.GPIO_Pin = ;
//	GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;		//模拟输入引脚
//	GPIO_Init(GPIOB, &GPIO_InitStructure);	
	
	ADC_DeInit(ADC1);  //复位ADC1 

	ADC_InitStructure.ADC_Mode = ADC_Mode_Independent;	//ADC工作模式:ADC1和ADC2工作在独立模式
	ADC_InitStructure.ADC_ScanConvMode = DISABLE;	//模数转换工作在单通道模式
	ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;	//模数转换工作在单次转换模式
	ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_None;	//转换由软件而不是外部触发启动
	ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;	//ADC数据右对齐
	ADC_InitStructure.ADC_NbrOfChannel = 1;	//顺序进行规则转换的ADC通道的数目
	ADC_Init(ADC1, &ADC_InitStructure);	//根据ADC_InitStruct中指定的参数初始化外设ADCx的寄存器   

  
	ADC_Cmd(ADC1, ENABLE);	//使能指定的ADC1
	
	ADC_ResetCalibration(ADC1);	//使能复位校准  
	 
	while(ADC_GetResetCalibrationStatus(ADC1));	//等待复位校准结束
	
	ADC_StartCalibration(ADC1);	 //开启AD校准
 
	while(ADC_GetCalibrationStatus(ADC1));	 //等待校准结束
 
//	ADC_SoftwareStartConvCmd(ADC1, ENABLE);		//使能指定的ADC1的软件转换启动功能

}	

4.4.2 ADC采样读值

u16 Get_Adc(u8 ch)   
{
  	//设置指定ADC的规则组通道，一个序列，采样时间
	ADC_RegularChannelConfig(ADC1, ch, 1, ADC_SampleTime_1Cycles5 );	//ADC1,ADC通道,采样时间为239.5周期	  			    
  
	ADC_SoftwareStartConvCmd(ADC1, ENABLE);		//使能指定的ADC1的软件转换启动功能	
	 
	while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC ));//等待转换结束

	return ADC_GetConversionValue(ADC1);	//返回最近一次ADC1规则组的转换结果
}

u16 Get_Adc_Average(u8 ch,u8 times)
{
	u32 temp_val=0;
	u8 t;
	for(t=0;t<times;t++)
	{
		temp_val+=Get_Adc(ch);
	}
	return temp_val/times;
} 

5 注意事项

    （1）调试时不要接电脑串口，要接电脑串口必须要和单片机共地，否则会影响ADC采样测量。

    （2）设置ADC采集判断信号到来的阈值电压信号是关键，要详细调整。

6  总结

       最后的作品不算完美，但这也是我们四天日夜奋战的结果，前期的各种探索经历过各种失败，最后完成了第一次比赛的大作品。如果本篇内容有不足之处，请大家多多包涵，提出建议，互相交流学习！