计算机工程
Computer Engineering
文章编号:1000#428(2018)07-0279-06
文献标志码:A
2018年7月
July 2018中图分类号:TP391
基于边缘卡尔曼滤波的GM-PHD多目标被动跟踪算法
曲
长
文
冯
奇
1,毛
宇
2
!
周强1
(1.海军航空工程学院电子信息工程系,山东烟台2001; 2.中国人民91431,海南文昌571300)
摘要:针对杂波干扰条件下,非线性、个数时变的多目标被动跟踪问题,提出一种基于边缘卡尔曼滤波的高斯混
合概率假设密度(PHD)滤波算法。采用边缘化变换计算目标状态的概率分布特性,获得目标状态及其协方差矩阵 估计的闭式解,解决目标模型非线性问题。利用量测信息生成新生目标强度,使滤波器具备对观测空间任意位置 随机出现新目标的跟踪能力。实验结果表明,与扩展卡尔曼算法、无迹卡尔曼算法和容积卡尔曼FK。 算法相比,该算法在生成目标轨迹、目标个数估计和跟踪精度等方面有更好的性能。关键词\"多目标跟踪;随机有限集;边缘卡尔曼滤波;概率假设密度;量测驱动
中文引用格式:曲长文,冯奇,毛宇,等.基于边缘卡尔曼滤波的GM-PHD多目标被动跟踪算法[J].计算机工 程,2018,44(7):279-284.英文引用格式:QU Changwen,FENG Qi,MAO Yu,et al. GM-PHD multi-objective passive tracking algorithm based on marg^in Kalman filtering* J] . Computer Engineering,2018 &44 (7) ,279-284.
GM-PHD Multi-objective Passive Tracking Algorithm
Based on Margin Kalman Filtering
QU Changwen1,FENG Qi1,MAO Yu2,ZHOU Qiang1
(1. Department of Electronic and Information Engineering,Naval Aeronautical and Astronautical University,Yantai,Shandong 2001,China;2.91431 Troops of PLA,Wenchang,Hainan 571300,China)
[Abstract] For theproblem of passive multi-objective tracking with thetime varying number in clutter, aGaussianmixture Probability Hypothesis Density ( GM-PHD) based on the marginalized Kalman filter for nonlinear Gaussiansystem is proposed. The marginalized transformation is applied to calculate the prediction and update distribution of target states for the nonlinear multi-objective models,and then the close-form solution of the target state and its covariance can be got. The measurements are used to generate thedensity of new born targets that appear randomly anywhere in the observation space. Experimental results show that compared with the Extended Kalman PHD ( EK-PHD),UnscentedKalman PHD ( UK-PHD) and Centralized Kalman PHD ( CK-PHD),the margin Kalman GM-PHD ( MK-PHD) is better in aspects of generating target trajectory,estimation of target number and tracking precision.
[Keywords] multi-objective tracking ; Random Finite Sets ( RFS); margin Kalman filtering; Probability Hypothesis Density( PHD); measurement driving DOI:10. 19678/j. issn. 1000-3428.0047524
〇概述
在实际的多目标跟踪场景中,目标个数未知且
时变、量测信息不确定(杂波干扰、虚警)等因素,为 多目标跟踪带来巨大困难。传统的多目标跟踪方 法[1-]是先关联后跟踪,但是当目标个数及杂波干扰 较多情况下,数据关联算法将会出现组合爆炸问题, 计算量剧增。为此,文献[4]基于随机有限集 (Random Finite Sets,RFS)理论,提出了概率假设密
度(Probability Hypothesis Density,PHD)滤波算法。该滤波算法将复杂的多目标状态空间运算投影到单
目标状态空间,然后通过传递全局后验概率密度的 一阶矩得到各目标状态估计。
与先关联后跟踪的传统多目标跟踪算法相比, PHD滤波算法较好地避开了数据关联运算,在最大 程度保持了信息完整性的同时降低了计算量。但是 PHD算法在递推过程中存在集合积分运算,计算上 难以实现。粒子 PHD(Particle PHD,P-PHD)[5-]和
作者简介:曲长文(1963—),男,教授、博士生导师,主研方向为无源定位与跟踪技术、合成孔径雷达成像、目标检测与识别;冯奇,博士
研究生;毛宇,学士;周强,助理研究员、博士。收稿日期:2017-06-08 修回日期:2017-07-14 E-mail:3793391+qq. com
280计算机工程2018年7月15日
高斯混合 PHD (Gaussian Mixture PHD, GM-PHD)是当前PHD算法的2种主要实现方法。顾名思义P-
PHD具有粒子滤波的优点,可用于非线性非高斯情 况下的多目标跟踪,但其不足也很明显,计算复杂度 高且目标状态提取精度严重依赖于聚类算法的稳定 性。文献[7]在线性高斯模型假设下提出了 GM-
PHD滤波器,在运算效率和状态提取方面优于P- PHD,且能给出后验概率密度的解析表达式。
在非线性情况下,GM-PHD滤波器无法得到闭 式解。针对此问题,文献[7-8]提出相应的扩展卡尔 曼 PHD (Extended Kalman PHD,EK-PHD)滤波器、 无迹卡尔曼 PHD ( Unscented Kalman PHD,UK- PHD)滤波器。而后,容积卡尔曼PHD ( Centralized Kalman PHD,CK-PHD)滤波器[(]也相继提出,此类 滤波器属于数值积分近似高维积分问题的范畴。
本文结合边缘卡尔曼滤波(Margin Kalman Filtering,MKF)[1051]提出一种 MK-PHD 滤波算法。 边缘化变换通过合适的先验分布矩阵可得到积分的 闭式解,对协方差矩阵的估计精度更高,且能保证协 方差矩阵的半正定性。此外,利用量测生成新生目 标强度[12 ],使新生目标强度覆盖所有量测,并采用 文献[13]方法解决权重归一化失衡问题。1
概率假设密度滤波
在RFS理论下,多目标目标状态与量测都可以 建模为如下随机有限集合K =丨)11,),2,…,\" 汽1),[=|^,1,^1,0,^,|\"以[)。其中,)1表 示1时刻第3个目标的状态向量,Ay表示1时刻 第J'个观测向量,K和[表示1时刻目标状态集合 和量测集合,和,为目标数和量测数,K和Z分 别表示目标状态空间和量测空间,S(K)和S(Z)表 示所有有限子集构成的集合。
PHD滤波器利用PHD代替多目标全局后验概 率密度在Bayes递推式中传递。假设1 - 1时刻 V_&( •)表示目标一阶矩,则PHD滤波器递推过程 可表示如下:
1) 预测
V 7 ~1()) =* + & (31 i-1() 7) +
.8())A7-1()7))V-1(.喏
(1)
2) 更新V ()) =( 1
()) ) V | t -1 ()) +#
!,:());(z|))
|-1())
(2)
其中,.8t())为目标生存概率,)表示目标状态向
量,=|-1( • ^)为单目标的状态转移密度,3||-1 (• | •)表示衍生目标的概率密度,*|为I时刻新生
目标的强度函数,.7,|())为I时刻目标检测概率,
9(a)为I时刻杂波强度,z表示观测向量,;t(z7) 为单目标似然函数。为方便讨论,本文假设没有衍 生目标。基于量测驱动的GM-PHD实现方法在文 献[13-15]已具体描述,本文不再赘述。&
边缘卡尔曼滤波
MKF滤波算法的基础是边缘化变换(Marginalized Transformation,MT),MT也是一种计算随机变量经
非线性变换后统计量的方法。假设随机向量I- *(/,>),其中,/和>为已知的均值和协方差。记 经过非线性变换;:^1 — ^\"后的随机向量为(尽 管;为确定性函数,MT中将其建模为具有先验分布 &(;)的随机过程。除了先验分布外,可利用的信息 为sigma点集;==[I0,1,…,I24]和经过非线性变换后 构成的取值向量z = [;(10),;(1),…,;(124)],则非 线性函数;的后验分布可记为.(;7,;=)。首先关于 ;的均值和协方差矩阵表达式可以表示为:
((;)=j *(,
;/,!);(1)d1(3)
!(;)=
=&N(t ■./,,! )[ ;(1 ) --■(;)] x[;(1) --\"(;)] Td1
(!)
对;进行边缘化积分可得J的均值:(=6[ ((;) |z,^l =&\"( ;).(g\\z,x) d; (5)
协方差矩阵为:
!'&(<) =6[!(;) |z,;==
&&(;)P(;\\z,x)d; ())
为了避开数值近似,需要找到合适的先验分布 &(;)使得式(5)和式(6)具有解析表达式。假设; 可以由Hermite多项式的线性组合近似表示。因为 Hermite多项式的计算是针对服从标准正态分布的 随机变量,需要随机解耦处理:引人标准随机向量
1 ~*(0,$4)并定义向量少=;(1 =;( I + 槡&1),
则其与J = ;(1) I1 ~ *( 1,!J服从相同的分布。
将;表示为P阶Hermite多项式线性组合的 形式:
;(1 ;!)=!?<( 1)
(7)
其中,<(1 ) = [E0,E1 ( -1 ),…,EP ( -1 ),E1 ( !2 ),…,Ep( !),…,E( ?„),…,Ep( !„)]?为 Hermite 多项 式组成的基函数,4表示状态向量的维数,P为多项 式的阶数。P维列向量!3和标量#构成了加权向 量!,其中,!3( 3 = 1,2,…,4,j' = 1,2,…,m)描述了! 到'的变换:
#0
…#… #
…
!4,1……!
第44卷第7期曲长文,冯奇,毛宇,等:基于边缘卡尔曼滤波的GM-PHD多目标被动跟踪算法281
因此描述了《维向量I通过基函数<(I转 换为量测'的变换,即:
'=(!)?<(〇
$8)
因此,合适的先验分布!3〜*(〇,,!‘),是求 解随机变量^的均值和方差的关键。值得说明的 是,对于!只需要知道其先验分布即可,其具体值并 不重要,其影响可通过边缘化积分消除。3 MK-PHD滤波器实现
本节结合MKF算法将GM-PHD滤波器的应用 扩展至非线性目标模型。目标状态转移方程和观测 方程为:
X1 = \"t ( X1-1,Vk-l )
(9)
(t= 初始化:令%3% = %,,!% }@%由观测量集 合及观测方程求得,其中, 表示新生高斯项均值 向量,!*%表示新生高斯项协方差矩阵。 选取先验协方差矩阵1,1为.阶对角阵,2 ( .(5,有最少2个非零元素。生成sigma点集;== [)%,)&,…,)2\"],其中,4为目标状态维数。当2 ( .(5时,选用UT变换sigma点集;当2(.(3时,选 用容积规则8gma点集。令!3 = 1,构造! = diag丨 (%/,,!/‘,…,!\"丨,其中,(#%/,取任意值,并计算 MKF算法中的不变参量,公式如下: 3 = [1,〇,…,〇]T/(„. + i) (11) C=diag {0,1!,…,.!,1!,…,.!}(„. + &) (12) #?$)) =[<()〇),...,<()2\")](„. + 1)x(2„ + 1) $13) P# , ( = ( / - P##T [ HP##T] 1 #)P #,D = [ 0,L] (14) 其中,L {\"…,\"….,\"=f1,0,…,0]1x.。在初始化结束之后,1时刻已知5 , Pt-1丨丨=11,观测量集合[及1时刻新生目标权重 ,其中,Qt-1表示存活高斯项均值向量,Pt-1表示 存活高斯项协方差矩阵。步骤1存活目标预测1) 选定sigma点集后,传递sigma点集可得: (= (Qt-1 + 槡?^0 ),…’<\"(qH + 槡―?)] (15) 2) 计算状态向量预测值Ql|)i-1 :f/# Z=P#H *hp0H]-1 (T (16)Qtl)l-1 =/3 (17) 3)估计比例因子:a:i =(2«+1) +2Z ( HP#H?) &|=1 &2& (18) (=[\"t(qI-1 + 槡P-1 )0) _QtJi)t-1,…, \"(Q-1 + yp|)2K -1] (19) 其中,(表示矩阵A的第j行。计算状态预测值的 协方差矩阵: Pf\\k-1 =/\\,〇/#\\z +diag|,,,,…,《„| x trace( p#, ZC) +* (20) 4)计算权值:1 -1 =.8t %1|1 (21)步骤&新生目标预测 由观测量集合[及观测方程求解得新生目标 状态向量预测值i i -1和协方差矩阵P*| i i -1,权重 | t-1 = %*,t。 步骤3存活目标更新 1) 考虑所有高斯项均为被检测到的情况,权重,状态向量均值和协方差矩阵: =( 1 -.7,1 ) @1 -1 (22)q |)=q|I)i-1 (23)Pi =P||t -1 (24) 2) 用Qt|t-1代替Q|-1,Pt|t-1代替P|-1,重复步 骤1中MKF处理方法,可得1时刻观测向量预测值 及其协方差矩阵: =M“C//#| z + diag{,1,,2 ’ …,,,} X !a#(P#|zC) +*t (25)3) 状态预测值与观测量预测值的互协方阵为: P& =槡!@1-10/#| [ (26)4) 卡尔曼增益矩阵为: =P& ( S|>) -1 (27) 5) 考虑能够被检测到的情况,利用观测量对所有高斯项更新均值,协方差矩阵和权重: ■100) =Q||k -1 +5|> ( Hk Z -Hk@ )( 28 %rP k( zZk - 1++)=Pk@k_1 G|)(5k@ ) ?( 29 %% ( Zk -1++) =.7,t %i)k-1 *( Hi;H@,sk@) (30%%k Zk _1++) =%kZk-1++)/2⑴,|:1,2,…,1 (31%2(z =9(z +#1%r*z+〇4■ %*,t 11 -1 (32% 步骤4 新生目标更新 按照一个目标对应一组观测量的假设R,t = *,对 每一个新生目标高斯项更新权重,均值和协方差矩阵: 2( ( =9t ( ( \"1 -1+\" +W*^,t |t' -1 (33) 更新 差矩 集合 282 计算机工程 2018年7月15日 当预测与更新结束之后,计算1 + 1时刻新生目标权重为%*,i + i 步骤M裁剪与合并 裁剪与合并方法与标准GM-PHD滤波算法相 同,存活目标进行裁剪与合并处理后再整合新生目 标得{%1@,q1@,!1@ };1i。步骤.多目标状态提取 令目标状态提取门限为0. 5,提取目标编号集合 6 . 〇.1,2,…,>0.5},则目标个数估计* =*疆(%1\" %,/\" 6,目标状态估计 Q3 . ,,\" 6,3.1,2,…,*1。 需要说明的是步骤3和步骤4中2((的表达式 应用了文献* 13 ]中的改进处理方法,避免了权重求解 中的归一化失衡问题。此外,针对新生目标权重的选 取,初始化阶段,令%@ . %〇,;• . 1,2,…,R*,〇。经预测、 更新处理后,得1 - 1时刻新生目标权重为%*@_1,则1时刻新生目标权重取% i . 1_1。 4 仿真结果与分析 本文在被动跟踪场景下,采用ek-phd、uk- PHD、CK-PHD及MK-PHD 4种方法对杂波环境下个数时变的多目标进行跟踪。定位采用双站 TDOA/DOA联合定位,与文献* 16 ]中收发分离的 多基地雷达相比,实现完全意义上被动定位与跟踪。 假设观测空间为二维平面区域*0 ',2 000 m] / *0 m,1 000 m],目标存活概率& =0.(8,目标检测 概率.7 =0.(9。观测站1位于m .*0,0]?,观测站 2位于m . *2 000 ',0 m]T,定位观测量为两站协 同估计的TDOA与观测站1估计得到的DOA。观 测时间为60 8采样周期O.1 8仿真中MKF算法 积分精度..5,先验分布矩阵为1 . diag 5 1,1 x 10-1,1 x10 2,1 x10-3,1 x10 46。 假设1时刻目标状态矢量为).* &,& 1,'1, '1]?,目标的状态方程为, ' .:)1-1 +0-1 (37%其中,01-1为过程噪声,且是同分布高斯白噪声, 即01-1〜*(0,P%,:为目标的状态转移矩阵。假设 目标模型为匀速(Constant Velocity,CV %模型,即: : . CV 1 [ :cvl , 其中 : _「1 O _「〇/4 O/2-:cv = [0 1],Pc叫 O/2 O2. (%为过程噪声标准差,取(% =0. 1 m/s2。量测 模型可描述为: (=< (&) +-1 (39 % <())=「2arcta1_rn('1/$11 )J] %3 .槡(&1 -&%2 + ('1 -'%2 其中,%3表示1时刻目标到观测站/的距离,-为观测 噪声,且是同分布高斯白噪声,即-〜*(0,*,* . diag5(2,(2#6,其中,(表示距离差估计标准差,(#表示 方位角估计标准差,仿真中取(% . 1 ',(# =0. 1町/ 180 rad。 假设杂波模型为泊松随机过程,其强度为:9(A) (a) (40%其中,为单位体积内杂波个数,5为观测空间的体 积,n( •)为均匀分布的概率密度函数,表示杂波出现 位置在观测区域内服从均匀分布。仿真中每平方米 取.5 x10—6,观测区域面积为2 x106 m2,也就是 说杂波平均数为10个。仿真中新生目标权重初始值 %。=0.01,裁剪门限O.1 x10—4,合并门限U.10,最 大高斯项Rmax.100。图1为目标真实运动轨迹示意 图。目标1在第1s新生,在第30 s消失;目标2在第 10 s新生,在第40 s消失;目标3在第20 s新生,在第 60 s消失;目标4在第5 s新生,在第50 s消失。 图1观测区域及目标的真实运动轨迹 图 2 为分别采用 EK-PHD、UK-PHD、CK-PHD 及MK-PHD的目标轨迹跟踪结果。由图2可得, EK-PHD算法对多目标轨迹跟踪效果最差,几乎不 能辨别目标轨迹,这是由于定位观测量TDOA与 DOA非线性较强导致的。UK-PHD与CK-PHD对 多目标跟踪效果相近,且明显由于EK-PHD。这是 因为UK-PHD和CK-PHD在计算随机变量通过非 线性变换后统计量问题上能精确到二阶矩,而EK- PHD只能精确到一阶矩。此外,由文献*16]关于 UK-PHD及CK-PHD性能对比分析可知,当目标状 态维数4(3时,UK-PHD与CK-PHD的稳定因子都 等于1,两者性能相当。在本文中只对目标的状态进 第44卷第7期曲长文,冯奇,毛宇,等:基于边缘卡尔曼滤波的GM-PHD多目标被动跟踪算法283 行估计,即n = 2,因此UK-PHD与CK-PHD跟踪性能相近,在下面的目标数估计及OSPA距离对比中 也能得到相同的结论。显然,4种方法中MK-PHD 跟踪效果最好,目标轨迹最清晰、最准确。 图3从目标个数估计方面对4种算法进行对 比。由图3可得,EK-PHD算法对目标个数漏估情 况严重。UK-PHD与CK-PHD对目标个数估计情况 相近在目标数最大时大约比实际目标个数少〇. 3。 MK-PHD对目标个数估计最准确,在目标个数最大 时大约比实际目标数少0.1。值得注意的是,所有算 法均在目标数目增长时刻表现出目标个数估计的滞 后性。这是因为目标的个数估计是通过处理存活目 标的权重得到的,所以当前的新生目标必须要等到 下一时刻才可以被估计到。 +真—实MK-PHD 目标数 CK-PHD........UK-PHD/ >« EK-PHD 图3目标个数估计对比 图4采用最优子模式指派$ Optimal Subpattern Assignment,OSPA)距离来度量4种算法的跟踪误差。 对2个有限集K .{&,&,…,&}和L . {'&,'$,…,',},〇SPA距离定义如下: Qospa $.,#,K,L). (/I( min #7# $4 &,'&$ +#( m-n)) )\\ . 丄 ,4(, (41) QOSPA(.,#,K,L) = Q OSPA ( P,#, L,K),m 其中,#是一个水平参数,1 (.< 5是距离敏感性参 数。仿真中取.=2,# = 1〇〇。 284 计算机工程2018年7月15日 由图4可得与前面仿真类似的结论,EK-PHD 算法的跟踪精度最差,UK-PHD、CK-PHD及MK- PHD算法跟踪效果明显优于EK-PHD,且MK-PHD 跟踪效果优于UK-PHD与CK-PHD算法。在第5 8、 第10 8和第20 S目标个数增加时刻4种算法的 OSPA距离均出现尖峰,这是由于个数估计的滞后性 导致的,与图3相对应。 最后,给出4种算法时间复杂度的对比结果。 仿真硬件为CPU主频2.5 GHz,内存3 GB的台式 机。仿真软件为MATLAB。通过500次Mon* Carlo仿真实验,计算得到各算法完成单次完整跟踪 平均耗时如表1所示。可以看出,-K-PHD算法耗 时最少,CK-PHD算法次之,UK-PHD算法耗时略大 于CK-PHD,这是因为每次递推计算CK-PHD所需 sigma点比UK-PHD少1个。MK-PHD耗时最多, 这是由于仿真中采用UK-PHD所用sigma点,以及解耦处理与估计比例因子,等运算导致的。 表1跟踪算法运行时间对比 跟踪算法 单次跟踪平均耗时/s 相对计算强度 EK-PHD2.93 1.00 UK-PHD3.431. 17CK-PHD3.361.14MK-PHD 3.60 1.23 5结束语 本文结合MKF滤波算法,将GM-PHD算法应 用扩展至非线性目标模型,给出MK-PHD算法的实 现步骤。通过数据仿真在目标个数估计、OSPA距离 及运算耗时等方面与EK-PHD、UK-PHD和CK-PHD 等算法进行对比分析。仿真结果表明,MK-PHD算 法生成目标轨迹更清晰,目标个数估计更准确,并且 跟踪精度更高。但在相同的软、硬件运行环境下, MK-PHD算法耗时最多,下一步将针对该问题展开 研究,从而提高MK-PHD算法的运算速度。 参考文献 [1 ] PASHA S A,TUAN HD,VO B N. Nonlinear Bayesian filtering using the unscented linear fractional transformation model [ J ]. IEEE Transactions on Signal Processing,2010,58 (2) :477-4.[2 ] BAR-SHALOM Y,FORTMAN T E. Tracking and data association [ M ] . San Diego, USA: Academic Press,1988.[3]何友,修建娟,关欣.雷达数据处理及应用[M]. 3版.北京:电子工业出版社,2013. [4 ] MAHLER R. Multi-target Bayes filtering via first-order multi-target moments [ J ] . IEEE Transactions on Aerospace and Electronic Systems,2003,39 (4 ): 1152- 1178.[5] VOB-N,SINGH S,DOUCET A. Sequential monte Carlo methods for multi-target filtering with random finite sets [ J]. IEEE Transactions on Aerospace and Electronic Systems,2005,41 (4) :1224-1245.[6 ] ZAJIC T,MAHLER R. A particle-systems implementation of the PHD multitarget tracking filter[C]//Proceedings of SPIE4 03. Washington D. C.,USA: IEEE Press,2003: 291-299.[7 ] VO B-N,MA W-K. The Gaussian mixture probability hypothesis density filter [ J ] $IEEE Transactions on Signal Processing,2006,54(11):4091-4104.[8 ] WANG Xiaoying,WANG Jiacun. Simulation analysis of EKF and UKF implementations in PHD filter [ C ] // Proceedings of IEEE International Conference on Networking,Sensing,and Control. Washington D. C., USA:IEEE Press,2016 :28-30.[9 ]王品,谢维信,刘宗香,等.一种非线性GM-PHD滤 波新方法[J].电子学报,2012,40(8):1597-1602.[10] SANDBLOMF,SVENSSON L . Moment estimation using a mar^ginalized transform [ J ]. IEEE Transactions on Signal Processing,2012,60(12):6138-6150.[11] 徐征,曲长文,王昌海.多站无源跟踪边缘卡尔曼滤 波算法[J].信号处理,2013,29(8):949-955.[12] 徐从安,刘瑜,熊伟,等.新生目标强度未知的双 门限粒子PHD滤波器[J].航空学报,2015,36( 12): 3957-3969$[13] 欧阳成,华云,高尚伟.改进的自适应新生目标强度 PHD滤波[J],系统工程与电子技术,2013,35 ( 12): 2452-2458$[14] RISTIC B,CLARK D,VO B-N,et al. Adaptive target birth intensity forPHD and CPHD filters [ J ]. IEEE Transactions on Aerospace and Electronic Systems, 2012,48 (2) :1656-1668.[15] CANG Yan,CHEN Di,SUN Weijin. Adaptive target birth intensity for Gaussian mixture probability hypothesis density filter [ C ] //Proceedings of IEEE International Conference on Control Science and Systems Engineering. Washington D. C.,USA : IEEE Press, 2014:29-30.[16] 王海环,王俊.基于容积卡尔曼的粒子PHD多目标 跟踪算法[J].系统工程与电子技术,2015,37 (9): 1960-1966. 编辑刘冰
因篇幅问题不能全部显示,请点此查看更多更全内容
Copyright © 2019- 7swz.com 版权所有 赣ICP备2024042798号-8
违法及侵权请联系:TEL:199 18 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务