一种采用听音定位技术实时监测目标物体运动轨迹的方法

文档序号:1002365 发布日期:2020-10-23 浏览:5次 >En<

阅读说明:本技术 一种采用听音定位技术实时监测目标物体运动轨迹的方法 (Method for monitoring motion trail of target object in real time by adopting listening positioning technology ) 是由 陈文兵 车文刚 于 2020-07-14 设计创作,主要内容包括:一种采用听音定位技术实时监测目标物体运动轨迹的方法,涉及军事、无人驾驶、智能机器人等技术领域。包括如下步骤:(1)用空间几何数学模型、传播过程中音强与声音传播距离的平方成反比、以及能量守恒定理等三种不同的方法建立微分方程来求解声源点(p)到各支麦克风(M&lt;Sub&gt;i&lt;/Sub&gt;)的距离d&lt;Sub&gt;i&lt;/Sub&gt;;(2)设计麦克风数量少且运算量小的克风分布树Tr;(3)采用最小二乘法从误差解集合中确定最优声源点位置p&lt;Sub&gt;i&lt;/Sub&gt;;(4)由惠更斯原理采用回声定位方法确定各支麦克风的相对位置;(5)随机采样已解声源点的位置来拟合目标物体的运动曲线。(A method for monitoring a motion trail of a target object in real time by adopting a listening positioning technology relates to the technical fields of military affairs, unmanned driving, intelligent robots and the like. The method comprises the following steps: (1) three different methods, namely a space geometric mathematical model, the inverse proportion of the sound intensity to the square of the sound propagation distance in the propagation process, the energy conservation theorem and the like are used for establishing a differential equation to solve the problem from a sound source point (p) to each microphone (M) i ) Distance d of i (ii) a (2) Designing a small number of microphones and a small amount of calculation for a microphone distribution tree Tr; (3) determining the optimal sound source point position p from the error solution set by adopting a least square method i (ii) a (4) Determining the relative position of each microphone by adopting an echo positioning method according to a Huygens principle; (5) the positions of the de-noised source points are randomly sampled to fit the motion curve of the target object.)

一种采用听音定位技术实时监测目标物体运动轨迹的方法

技术领域

本发明涉及目标物体运动轨迹监测技术领域,特别涉及一种采用听音定位技术实时监测目标物体运动轨迹的方法。

背景技术

听音定位模型在军事、商业上都有很重要的用途。我们希望把听音定位模型用在人工智能方向,由于图像在拍摄过程受到多种不定因素的影响,并且对图像的处理和分析技术难度高,数据量较大,即使近来比较火热的深度学习方法也受到数据多样性和计算量大等特点的限制。因此,听音定位模型可以简化一些在三维空间上难于分析的问题,如自动驾驶技术、机器人工作等方面。同时此项技术和医学技术相结合可以改善残疾人的生活方式。

和已经存在的回声定位方法相比,听音定位技术不用靠探测物体主动发出声音再接受回声来计算目标物体的位置,其原理和动物利用耳朵分辨声源方位一样,靠两个和两个以上不同位置的麦克风接受到的声音来判断声源的方位。本专利主要从几何数学、运动学公式、能量守恒定理等不同理论方向建立微分方程式分析问题,利用最小二乘法改进了DTOA算法的误差优化方式。最后,通过麦克风树随机时间点采样的声音来拟合出目标物体的运动轨迹,节约资源,减少运算量。

发明内容

本发明的目的是提供一种采用听音定位技术实时监测目标物体运动轨迹的方法,可以实现任意密闭空间中物体运动状态的实时监测。

为达到上述目的,本发明采用如下技术方案:

一种采用听音定位技术实时监测目标物体运动轨迹的方法,包括以下步骤:(1)采用三种不同的方法建立微分方程来求解声源点p到各支麦克风Mi的距离di,1)、空间几何数学模型,2)、传播过程中音强与声音传播距离的平方成反比,3)能量守恒定理;(2)设计麦克风数量少且运算量小的克风分布树Tr;(3)采用最小二乘法从误差解集合中确定最优声源点位置pi;(4)由惠更斯原理采用回声定位方法确定各支麦克风的相对位置;(5)随机采样已解声源点的位置来拟合目标物体的运动曲线以便做预测减少计算量。

进一步具体包括如下步骤:

(1)求解声源点p到各支麦克风Mi的距离di

在空间直角坐标系中已知若干点的坐标Mi(xi,yi,zi)求未知点p(x,y,0),利用已知点M到点p距离di建立一个方程组来解出两个未知量x和y,未知点p就是两个方程曲线的交点;

1)用空间几何数学模型求解di,1

设n只麦克风M1,M2,…,Mn分别接收到同一声源信号的时刻为ti,则Mi和M1接收到同一声源发出信号的时间差:

ti,1=|ti-t1| (9)

声音在空气中的传播速度为定值c,声源p分别到麦克风M1和麦克风Mi的距离差为d1、di

di,1=c.ti,1 (10)

di,1=di-d1 (11)

由公式(10)和(11)可以得出:

Figure BDA0002583926240000031

和公式(1)联立可以得出:

Figure BDA0002583926240000032

公式(13)中用上式减下式有:

Figure BDA0002583926240000033

其中I=1,2,…,n;令xi,1=xi-x1,yi,1=yi-y1,目标物位于地平面上,因此z=0,化简(14)有:

Figure BDA0002583926240000035

线性方程组(15)中x,y,d1为未知量,结合公式(1),再建立两个关系式可以求解x,y;当麦克风数量为3时,可以得到两个TDOA测量值,先假定d1已知,则声源的位置(x,y,0)可以由线性方程组的矩阵相乘形式表述:

两边同时乘以公式(16)中的第一个矩阵的逆可以求出

Figure BDA0002583926240000037

的代数式,将之带入公式(1)

Figure BDA0002583926240000038

中,可以得到一个关于d1的二次方程,将其正根带入(16)就可以求出声源p的位置。

2)由传播过程中音强与声音传播距离的平方成反比求解di

声音在不变介质中传播时,麦克风处的音强V与到音源的距离平方d2成反比,有:

V=V0/d2 (17)

因此,测出声源在两只麦克风处的音强V1,Vi,联立方程组可以求解出d1,di代会公式(1),可以求解出声源点坐标p(x,y,0)。

Figure BDA0002583926240000041

3)由能量守恒定理求解di

声音在传播过程中遵循能量守恒定理,θ为声音在空气中传播造成能量减少的因素集,Δd表示声音的直线传播位移,建立微分方程:

Δw=-θ.Δd (19)

对微分方程(19)两边同时积分有:

∫dw=∫-θdx (20)

对(20)求解我们可以得到声音在传播过程中随位移变化的曲线:

w=w0-θx (21)

其中w0是声源的初始能量,x是当前方向上的声信号的直线位移;由于声音的发出时间、和w0未知,只能测到同一声源传到麦克风处的能量wi,现假设测量到三只麦克风声音的能量w1,wi

Figure BDA0002583926240000042

计算机中矩阵运算方式:

Figure BDA0002583926240000051

其中公式(22)的第三式可以由公式(9)推出,解方程组求出x1,xi(di)带入公式(1)可以求出声源点坐标p(x,y,0)。

(2)设计麦克风数量少且运算量小的克风分布树Tr

如果不考虑误差分析和优化,在已知物体在水平地面情况下,同时综合考虑对麦克风相对位置的求解复杂程度比,只需要M1,M2,M3三只麦克就能求解此问题,当三只麦克风相互垂直时最节约成本且可以达到使用要求的。

(3)采用最小二乘法从误差解集合中确定最优声源点位置pi

在实际情况中,声音在传播的过程中会受到噪音和混响音的影响,声音传播到两只不同的麦克风处存在着时延误差,故 把模型放到平面上来分析,其中理想情况下所有曲线会交于一点p,但是所有曲线大概率情况下会有一个误差交空间,其中所有的点构成一个解集合H;

实际情况下,当麦克风的数量越多(>3)误差空间会越小,直到出现最优解,但是考虑成本效率问题,适当的增加麦克风构建解空间,然后把误差处理转化到在解集合H中求解最优的点问题,通过先验知识把距离三条曲线距离和最小的点作为最优解点位置,令W为解线性方程组参数的矩阵,X为因变量矩阵,h为自变量,则理想解的线性方程组可简化为:

WX+b=h (2)

点p对应的实际误差解方程组可以表示成:

e=h-(WX+b) (3)

由于外因素和测量误差的影响,公式(3)的解会形成一个解区域;因此,一种新的降低误差的方法,从解集点空间H中寻找一个到三条曲线的距离和最小的点(最优点),研究发现,解集点空间必定在三条曲线同时向外凸出的区域,并且最优点一定会在由三个交点直接相连组成的三角形里面;

通过求解方程组可以解出三角形的三个顶点坐标,转换可以得到三角形三条边的直线方程,解集合中点p(xi,yi)的遍历用以下方式,当x从X1到X3以每次K cm的间隔跳动取值xi,同时y在f1(xi)和f2(xi)或f1(xi)和f3(xi)之间以K cm的间隔取值,如果间隔取得越小误差越低,但是计算量也会随之增大,用这种方式遍历所有的误差点;

由以上模型,p(xi,yi),(i=1,2,…,n),用εi,k表示第i个点到第k条曲线的距离,可以由点到直线的距离公式求解。

(4)由惠更斯原理采用回声定位方法确定各支麦克风的相对位置

在求解麦克风相对于大厅的坐标时,利用麦克风处发出的子声源来建立回声定位模型;在封闭空间中,子声源接收到垂直反射波和斜角反射波两种反射波;基于先验知识和实体模型推导,从麦克风M出扩散出的声波反射轨迹可以有以下结论;

A.斜角反射波不能返回到M处,即使能返回也会先在空间中做多次往返运动;

B.直线碰撞的两个反射波和垂直碰撞的两个反射波的“混响音”波普不同;

C.麦克风会最先接收到6个垂直于墙面的垂直反射波,取垂直反射波路径可以确定麦克风相对于大厅的位置;

当麦克风Mk已经接收到6组1/2k频率回音(k为麦克风编号),往返时间为ΔT1,ΔT1,…,ΔT6,通过“混响音”的平频谱分析得出ΔT1和ΔT3、ΔT2和ΔT4、ΔT4和ΔT6为三组直线碰撞音的往返时间差,且ΔT1、ΔT2和ΔT4为每组碰撞中较先返回到M的回声的往返时间差,因此,由回声的垂直往返线路构成的空间角就是距麦克风树最接近的墙角;

求解麦克风M1,M2,…,Mn的位置,先要利用一个麦克风来回声定位确定一个墙角为坐标原点,其他麦克风Mi的位置都以第一个麦克风确定的空间坐标系为基础去定位;

由以上公式可以确定M1(x1,y1,z1)的坐标和坐标系;所有麦克风Mi都是在同一麦克风支架上,只要所有麦克风不在同一条线上,最后不会影响求解声源位置p,因此,为了提高计算效率可以改进麦克风树,其中M1的垂直上空必须放置一只麦克风M2,其它任意麦克风可以必须和麦克风M1(或M2)在同一水平面上;M2的坐标可以确定为:Di,j为两麦克风间的距离;

Mi的垂直坐标为z1,因此求麦克风Mi(xi,yi,z2)坐标,只需要求解xi,yi两个未知数;

由方程组(5)可以解出xi,yi从而确定Mi(xi,yi,z2)的坐标;并且z1和z2就是麦克风相对于水平地面的高度;前面已经求出声源p的坐标和麦克风M坐标的关系,把麦克风Mi(xi,yi,z2)回代可以确定目标物体在大厅中的位置。

(5)随机采样已解声源点的位置来拟合目标物体的运动曲线以便做预测减少计算量;

取随机时间点求解一遍目标物体的位置,由于目标物体在水平面上运动,通过上面的步骤求出m个样本点pi(xi,yi),(i=1,2,…,m),求解出物体的运动轨迹曲线函数f(x),再通过物体的平均速度v转化为一个关于时间的f(t)函数,由麦克风首次接收到声音的初始时间和初始点位置就能用时间t预测物体接下来的行驶位置;

由泰勒公式可知,任何n阶可导函数都能在x=0处展开:

Figure BDA0002583926240000082

公式(6)可以化简为:

f(x)=c0+c1x+c2x2+c3x3+…+cnxn (7)

由公式(7)求解出C(c0,c1,c2,…,cn)就能求解出f(x),当采样点数m>n时,求解方程组就能得到f(x);确定C(c0,c1,c2,…,cn)采用最小二乘准则,使m个点pi(xi,yi)与曲线f(x)的距离δi的平方和最小;

记:

Figure BDA0002583926240000084

物体在运动过程中存在着惯性,一般都会是直线行驶(这时候也可以使用线性回归模拟)f(x)为一次函数,假定人在曲线行驶,这时候其运动轨迹最复杂为一元三次函数,因此,把m个点带入(7)构成的方程组中,一共有m-n(n≤3)组C的解,把这m-n组参数分别带入f(x),再代入公式(8)找到一组最优的C(c0,c1,c2,c3);

一元二次方程为:f(x)=c0+c1x+c2x2

一元三次方程为:f(x)=c0+c1x+c2x2+c3x3

本发明与现有技术相比,具有以下优点:使用三种不同的方法来确定目标物体相对于麦克风的位置关系,合理设计的麦克风树降低了对麦克风坐标求解的复杂程度。用加权最小二乘法对由“混响音”和噪音造成的时延误差解集合做了优化,在一定程度上降低了误差。最后再用曲线拟合预测人的运动轨迹减少整个模型的计算量。

附图说明

图1是本发明的流程图;

图2是麦克风树的最优设计方案;

图3是理想解情况的示意图;

图4是误差解空间集合H。

具体实施方式

以下结合附图对本发明作进一步详细说明。

实施例一,一种采用听音定位技术实时监测目标物体运动轨迹的方法,参照图1-4,包括以下步骤:(1)求解声源点p到各支麦克风Mi的距离di;(2)设计麦克风数量少且运算量小的克风分布树Tr;(3)采用最小二乘法从误差解集合中确定最优声源点位置pi;(4)由惠更斯原理采用回声定位方法确定各支麦克风的相对位置;(5)随机采样已解声源点的位置来拟合目标物体的运动曲线以便做预测减少计算量。

进一步具体包括如下步骤:

(1)求解声源点p到各支麦克风Mi的距离di

在空间直角坐标系中已知若干点的坐标Mi(xi,yi,zi)求未知点p(x,y,0),利用已知点M到点p距离di建立一个方程组来解出两个未知量x和y,未知点p就是两个方程曲线的交点;

Figure BDA0002583926240000101

用空间几何数学模型求解di,1

设n只麦克风M1,M2,…,Mn分别接收到同一声源信号的时刻为ti,则Mi和M1接收到同一声源发出信号的时间差:

ti,1=|ti-t1| (9)

声音在空气中的传播速度为定值c,声源p分别到麦克风M1和麦克风Mi的距离差为d1、di

di,1=c.ti,1 (10)

di,1=di-d1 (11)

由公式(10)和(11)可以得出:

Figure BDA0002583926240000102

和公式(1)联立可以得出:

Figure BDA0002583926240000111

公式(13)中用上式减下式有:

Figure BDA0002583926240000112

其中

Figure BDA0002583926240000113

I=1,2,…,n;令xi,1=xi-x1,yi,1=yi-y1,目标物位于地平面上,因此z=0,化简(14)有:

Figure BDA0002583926240000114

线性方程组(15)中x,y,d1为未知量,结合公式(1),再建立两个关系式可以求解x,y;当麦克风数量为3时,可以得到两个TDOA测量值,先假定d1已知,则声源的位置(x,y,0)可以由线性方程组的矩阵相乘形式表述:

Figure BDA0002583926240000115

两边同时乘以公式(16)中的第一个矩阵的逆可以求出的代数式,将之带入公式(1)

Figure BDA0002583926240000117

中,可以得到一个关于d1的二次方程,将其正根带入(16)就可以求出声源p的位置。

(2)设计麦克风数量少且运算量小的克风分布树Tr

如果不考虑误差分析和优化,在已知物体在水平地面情况下,同时综合考虑对麦克风相对位置的求解复杂程度比,只需要M1,M2,M3三只麦克就能求解此问题,当三只麦克风相互垂直时最节约成本且可以达到使用要求的;

(3)采用最小二乘法从误差解集合中确定最优声源点位置pi

在实际情况中,声音在传播的过程中会受到噪音和混响音的影响,声音传播到两只不同的麦克风处存在着时延误差,故 把模型放到平面上来分析,其中理想情况下所有曲线会交于一点p,但是所有曲线大概率情况下会有一个误差交空间,其中所有的点构成一个解集合H;

实际情况下,当麦克风的数量越多(>3)误差空间会越小,直到出现最优解,但是考虑成本效率问题,适当的增加麦克风构建解空间,然后把误差处理转化到在解集合H中求解最优的点问题,通过先验知识把距离三条曲线距离和最小的点作为最优解点位置,令W为解线性方程组参数的矩阵,X为因变量矩阵,h为自变量,则理想解的线性方程组可简化为:

WX+b=h (2)

点p对应的实际误差解方程组可以表示成:

e=h-(WX+b) (3)

由于外因素和测量误差的影响,公式(3)的解会形成一个解区域;

因此,一种新的降低误差的方法,从解集点空间H中寻找一个到三条曲线的距离和最小的点(最优点),研究发现,解集点空间必定在三条曲线同时向外凸出的区域,并且最优点一定会在由三个交点直接相连组成的三角形里面;

通过求解方程组可以解出三角形的三个顶点坐标,转换可以得到三角形三条边的直线方程,解集合中点p(xi,yi)的遍历用以下方式,当x从X1到X3以每次K cm的间隔跳动取值xi,同时y在f1(xi)和f2(xi)或f1(xi)和f3(xi)之间以K cm的间隔取值,如果间隔取得越小误差越低,但是计算量也会随之增大,用这种方式遍历所有的误差点;

由以上模型,p(xi,yi),(i=1,2,…,n),用εi,k表示第i个点到第k条曲线的距离,可以由点到直线的距离公式求解;

下面是求解最优点p(xi,yi)的伪代码:

(4)由惠更斯原理采用回声定位方法确定各支麦克风的相对位置

在求解麦克风相对于大厅的坐标时,利用麦克风处发出的子声源来建立回声定位模型;在封闭空间中,子声源接收到垂直反射波和斜角反射波两种反射波;基于先验知识和实体模型推导,从麦克风M出扩散出的声波反射轨迹可以有以下结论;

A.斜角反射波不能返回到M处,即使能返回也会先在空间中做多次往返运动;

B.直线碰撞的两个反射波和垂直碰撞的两个反射波的“混响音”波普不同;

C.麦克风会最先接收到6个垂直于墙面的垂直反射波,取垂直反射波路径可以确定麦克风相对于大厅的位置;

当麦克风Mk已经接收到6组1/2k频率回音(k为麦克风编号),往返时间为ΔT1,ΔT1,…,ΔT6,通过“混响音”的平频谱分析得出ΔT1和ΔT3、ΔT2和ΔT4、ΔT4和ΔT6为三组直线碰撞音的往返时间差,且ΔT1、ΔT2和ΔT4为每组碰撞中较先返回到M的回声的往返时间差,因此,由回声的垂直往返线路构成的空间角就是距麦克风树最接近的墙角;

求解麦克风M1,M2,…,Mn的位置,先要利用一个麦克风来回声定位确定一个墙角为坐标原点,其他麦克风Mi的位置都以第一个麦克风确定的空间坐标系为基础去定位;

Figure BDA0002583926240000141

由以上公式可以确定M1(x1,y1,z1)的坐标和坐标系;所有麦克风Mi都是在同一麦克风支架上,只要所有麦克风不在同一条线上,最后不会影响求解声源位置p,因此,为了提高计算效率可以改进麦克风树,其中M1的垂直上空必须放置一只麦克风M2,其它任意麦克风可以必须和麦克风M1(或M2)在同一水平面上;M2的坐标可以确定为:Di,j为两麦克风间的距离;

Mi的垂直坐标为z1,因此求麦克风Mi(xi,yi,z2)坐标,只需要求解xi,yi两个未知数;

Figure BDA0002583926240000151

由方程组(5)可以解出xi,yi从而确定Mi(xi,yi,z2)的坐标;并且z1和z2就是麦克风相对于水平地面的高度;前面已经求出声源p的坐标和麦克风M坐标的关系,把麦克风Mi(xi,yi,z2)回代可以确定目标物体在大厅中的位置;

(5)随机采样已解声源点的位置来拟合目标物体的运动曲线以便做预测减少计算量;

取随机时间点求解一遍目标物体的位置,由于目标物体在水平面上运动,通过上面的步骤求出m个样本点pi(xi,yi),(i=1,2,…,m),求解出物体的运动轨迹曲线函数f(x),再通过物体的平均速度v转化为一个关于时间的f(t)函数,由麦克风首次接收到声音的初始时间和初始点位置就能用时间t预测物体接下来的行驶位置;

由泰勒公式可知,任何n阶可导函数都能在x=0处展开:

Figure BDA0002583926240000152

公式(6)可以化简为:

f(x)=c0+c1x+c2x2+c3x3+…+cnxn (7)

由公式(7)求解出C(c0,c1,c2,…,cn)就能求解出f(x),当采样点数m>n时,求解方程组就能得到f(x);确定C(c0,c1,c2,…,cn)采用最小二乘准则,使m个点pi(xi,yi)与曲线f(x)的距离δi的平方和最小;

记:

Figure BDA0002583926240000153

物体在运动过程中存在着惯性,一般都会是直线行驶(这时候也可以使用线性回归模拟)f(x)为一次函数,假定人在曲线行驶,这时候其运动轨迹最复杂为一元三次函数,因此,把m个点带入(7)构成的方程组中,一共有m-n(n≤3)组C的解,把这m-n组参数分别带入f(x),再代入公式(8)找到一组最优的C(c0,c1,c2,c3);

一元二次方程为:f(x)=c0+c1x+c2x2

一元三次方程为:f(x)=c0+c1x+c2x2+c3x3

实施例二,实施例二与实施例一相同的技术方案不再赘述,其不同的技术方案在于,在步骤(1)求解声源点p到各支麦克风Mi的距离di时,由传播过程中音强与声音传播距离的平方成反比求解di

声音在不变介质中传播时,麦克风处的音强V与到音源的距离平方d2成反比,有:

V=V0/d2 (17)

因此,测出声源在两只麦克风处的音强V1,Vi,联立方程组可以求解出d1,di代会公式(1),可以求解出声源点坐标p(x,y,0)。

Figure BDA0002583926240000161

实施例三,实施例三与实施例一相同的技术方案不再赘述,其不同的技术方案在于,在步骤(1)求解声源点p到各支麦克风Mi的距离di时,由能量守恒定理求解di

声音在传播过程中遵循能量守恒定理,θ为声音在空气中传播造成能量减少的因素集,Δd表示声音的直线传播位移,建立微分方程:

Δw=-θ.Δd (19)

对微分方程(19)两边同时积分有:

∫dw=∫-θdx (20)

对(20)求解我们可以得到声音在传播过程中随位移变化的曲线:

w=w0-θx (21)

其中w0是声源的初始能量,x是当前方向上的声信号的直线位移;由于声音的发出时间、和w0未知,只能测到同一声源传到麦克风处的能量wi,现假设测量到三只麦克风声音的能量w1,wi

计算机中矩阵运算方式:

其中公式(22)的第三式可以由公式(9)推出,解方程组求出x1,xi(di)带入公式(1)可以求出声源点坐标p(x,y,0)。

代入数据进行测试:

Pi是声源的位置坐标,M1、M2和M3分别是三个麦克风的位置,d2,1=d2-d1表示距离差,di表示每支麦克风到声源的距离,Δt两两麦克风之间接收到同一声源信号的时间差,c表示音速。

表1 测试数据

Figure BDA0002583926240000173

Figure BDA0002583926240000181

因为此问题中物体是在水平地面上的,所以我们带入最后一组测试用例得到的结果为:

result:[(1.00000000000000,1.00000000000000,1.00000000000000)],括号中前两个数据和声源(1,1,0)的x和y对应,括号中第三个数据时M1到声源的距离。

最后在模型的测试中,如果不考虑外观因素的情况下,利用此模型计算出的点就是物体的实际位置。考虑外观因素时,如果不断增加麦克风的数量可以降低误差使计算点无限接近物体的实际位置,误差分析中利用最小二乘法从误差解集合中求解的最优点位置能以0.94的概率接近物体的实际位置。并且拟合出的光滑曲线能以较高的仿真度接近人的实际运动轨迹。

本具体实施例仅仅是对本发明的解释,其并不是对本发明的限制,本领域技术人员在阅读完本说明书后可以根据需要对本实施例做出没有创造性贡献的修改,但只要在本发明的权利要求范围内都受到专利法的保护。

19页详细技术资料下载
上一篇:一种医用注射器针头装配设备
下一篇:基于机器视觉和全息方法的声场测试分析方法及系统

网友询问留言

已有0条留言

还没有人留言评论。精彩留言会获得点赞!

精彩留言,会给你点赞!