SAR合成孔径雷达图像点目标仿真报告(附matlab代码)(共23页).doc

上传人:飞****2 文档编号:13826927 上传时间:2022-05-01 格式:DOC 页数:23 大小:775KB
返回 下载 相关 举报
SAR合成孔径雷达图像点目标仿真报告(附matlab代码)(共23页).doc_第1页
第1页 / 共23页
SAR合成孔径雷达图像点目标仿真报告(附matlab代码)(共23页).doc_第2页
第2页 / 共23页
点击查看更多>>
资源描述

《SAR合成孔径雷达图像点目标仿真报告(附matlab代码)(共23页).doc》由会员分享,可在线阅读,更多相关《SAR合成孔径雷达图像点目标仿真报告(附matlab代码)(共23页).doc(23页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、精选优质文档-倾情为你奉上SAR图像点目标仿真报告徐一凡1 SAR原理简介合成孔径雷达(Synthetic Aperture Radar .简称SAR)是一种高分辨率成像雷达技术。它利用脉冲压缩技术获得高的距离向分辨率.利用合成孔径原理获得高的方位向分辨率.从而获得大面积高分辨率雷达图像。SAR回波信号经距离向脉冲压缩后.雷达的距离分辨率由雷达发射信号带宽决定:.式中表示雷达的距离分辨率.表示雷达发射信号带宽.表示光速。同样.SAR回波信号经方位向合成孔径后.雷达的方位分辨率由雷达方位向的多谱勒带宽决定:.式中表示雷达的方位分辨率.表示雷达方位向多谱勒带宽.表示方位向SAR平台速度。在小斜视角

2、的情况下.方位分辨率近似表示为.其中为方位向合成孔径的长度。2 SAR的几何关系雷达位置和波束在地面覆盖区域的简单几何模型如图1所示。此次仿真考虑的是正侧视的条带式仿真.也就是说倾斜角为零.SAR波束中心和SAR平台运动方向垂直的情况。图1 雷达数据获取的几何关系建立坐标系XYZ如图2所示.其中XOY平面为地平面;SAR平台距地平面高H.以速度V沿X轴正向匀速飞行;P点为SAR平台的位置矢量.设其坐标为(x,y,z); T点为目标的位置矢量.设其坐标为;由几何关系.目标与SAR平台的斜距为: (1)由图可知:;令.其中为平台速度.s为慢时间变量(slow time).假设.其中表示SAR平台的

3、x坐标为的时刻;再令.表示目标与SAR的垂直斜距.重写(1)式为: (2) 就表示任意时刻时.目标与雷达的斜距。一般情况下.于是通过傅里叶技术展开.可将(2)式可近似写为: (3)可见.斜距是的函数.不同的目标.也不一样.但当目标距SAR较远时.在观测带内.可近似认为不变.即。 图2:空间几何关系 (a)正视图 (b)侧视图图2(a)中.表示合成孔径长度.它和合成孔径时间的关系是。(b)中.为雷达天线半功率点波束角.为波束轴线与Z轴的夹角.即波束视角.为近距点距离.为远距点距离.W为测绘带宽度.它们的关系为: (4)3 SAR的回波信号模型SAR在运动中以一定的周期()发射和接收信号.具体过程

4、如图3所示。发射机以的时间发射啁啾脉冲.然后切换天线开关接收回波信号。图3 雷达发射脉冲串的时序当雷达不处于发射状态时.它接收3反射回波。发射和接收回波的时间序列如图4所示。在机载情况下.每个回波可以在脉冲发射间隔内直接接收到。但是在星载情况下.由于距离过大.某个脉冲的回波要经过610个脉冲间隔才能接收到。这里仿真为了方便.默认为机载情况。图4 脉冲雷达的发射与接收周期假设为chirp信号持续时间.下标表示距离向;PRF为重复频率.PRT为重复周期.等于。接收序列中.表示发射第个脉冲时.目标回波相对于发射序列的延时。雷达的发射序列数学表达式为式(5): (5)式中.表示矩形信号.为距离向的ch

5、irp信号调频率.为载频。雷达回波信号由发射信号波形.天线方向图.斜距.目标RCS.环境等因素共同决定.若不考虑环境因素.则单点目标雷达回波信号可写成式(6)所示: (6)其中.表示点目标的雷达散射截面.表示点目标天线方向图双向幅度加权.表示载机发射第n个脉冲时.电磁波再次回到载机时的延时.带入式(6)中得: (7)式(7)就是单点目标回波信号模型.其中.是chirp分量.它决定距离向分辨率;为多普勒分量.它决定方位向分辨率。对于任意一个脉冲.回波信号可表示为式(8)所示:(8)我们知道.由于随慢时间s的变化而变化.所以计算机记录到的回波数据存储形式如图5所示:图5 目标照射时间内.单个点目标

6、回波能量在信号处理器的二维存储器中的轨迹4 距离徙动及校正根据图2可知.在倾斜角为零或很小的时候.目标与雷达的瞬时距离为.根据几何关系可知.根据泰勒级数展开可得: (9)由式(9)可知.不同慢时间对应着不同的.并且是一个双曲线形式或者近似为一个二次形式。如图5所示.同一目标的回波存储在计算机里不在同一直线上.存在距离徙动。从而定义距离徙动量: (10) 为了进行方位向的压缩.方位向的回波数据必须在同一条直线上.也就是说必须校正距离徙动。由式(10)可知.不同的最近距离r对应着不同的.因此在时域处理距离徙动会非常麻烦。因此.对方位向进行傅里叶变换.对距离向不进行变换.得到新的域。由于方位向的频率

7、即为多普勒频率.所以这个新的域也称为距离多普勒域。将斜距R写成多普勒fa的函数.即。众所周知.对最近距离为r的点目标P.回波多普勒是倾斜角的函数.即.斜距.于是 (11)所以距离多普勒域中的我距离徙动为=.可发现它不随慢时间变换.同一最短距离r对应着相同大小的距离徙动。因此在距离多普勒域对一个距离徙动校正就是对一组具有相同最短距离的点目标的距离徙动校正.这样可以节省运算量。为了对距离徙动进行校正.需要得到距离徙动单元.即距离徙动体现在存储单元中的移动数值.距离徙动单元可以表示为.这个值通常为一个分数.由于存储单元都是离散的.所以不同通过在存储单元简单的移动得到准确的值。为了得到准确的徙动校正值

8、.通常需要进行插值运算。本仿真采用了两种插值方法最近邻点插值和sinc插值.下面分别进行介绍。最近邻点插值法的优点是简单而快速.缺点是不够精确。.其中N为整数部分.n为小数部分.整数部分徙动可以直接通过平移消除.对于小数部分则通过四舍五入的方法变为0或者1.这样就可以得到较为精确的插值。Sinc插值原理如下:在基带信号下.卷积核是sinc函数(12)插值信号为(13)即为所有输入样本的加权平均。可通过频域来理解.如图6所示.采样信号的频谱等于以采样率重复的信号频谱。为了重建信号.只需要一个周期频谱(如基带周期).因此需要理想矩形低通滤波器在频域中提取基带频谱(如图6)所示。已知该理想滤波器在时

9、域中是sinc函数。由于频域相乘相当于时域卷积.故插值可以通过与sinc核的卷积来实现。图6 理想低通滤波器怎样对采样信号进行插值5 点目标成像matlab仿真5.1距离多普勒算法距离多普勒算法(RDA)是在1976年至1978年为民用星载SAR提出的.它兼顾了成熟、简单、高效和精确等因素.至今仍是使用最广泛的成像算法。它通过距离和方位上的频域操作.到达了高效的模块化处理要求.同时又具有了一维操作的简便性。图7示意了RDA的处理流程。这里主要讨论小倾斜角及短孔径下的基本RDA处理框图。1.当数据处在方位时域时.可通过快速卷积进行距离压缩。也就是说.距离FFT后随即进行距离向匹配滤波.再利用距离

10、IFFT完成距离压缩。回波信号为:(14)距离向压缩后的信号为: (15) (16)2.通过方位FFT将数据变换至距离多普勒域.多普勒中心频率估计以及大部分后续操作都在该域进行。方位向傅里叶变换后信号为: (17)3.在距离多普勒域进行随距离时间及方位频率变化的RCMC.该域中同一距离上的一组目标轨迹相互重合。RCMC将距离徙动曲线拉直到与方位频率轴平行的方向。这里可以采用最近邻点插值法或者sinc插值法.具体插值方法见前面。假设RCMC插值是精确的.信号变为: (18)4.通过每一距离门上的频域匹配滤波实现方位压缩。为进行方位压缩.将RCMC后的乘以频域匹配滤波器。(19) (20)5.最后

11、通过方位IFFT将数据变换回时域.得到压缩后的复图像。复原后的图像为:(21)图8 距离多普勒算法流程图5.2 Chirp Scaling算法距离多普勒算法具有诸多优点.但是距离多普勒算法有两点不足:首先.当用较长的核函数提高距离徙动校正(RCMC)精度时.运算量较大;其次.二次距离压缩(SRC)对方位频率的依赖性问题较难解决.从而限制了其对某些大斜视角和长孔径SAR的处理精度。Chirp Scaling算法避免了RCMC中的插值操作.通过对Chirp信号进行频率调制.实现了对该信号的尺度变换或平移。图8显示了Chirp Scaling算法处理流程。这里主要讨论小倾斜角及短孔径下的基本CSA处

12、理框图。主要步骤包括四次FFT和三次相位相乘。1.通过方位向FFT将数据变换到距离多普勒域。2.通过相位相乘实现Chirp Scaling操作.使所有目标的距离徙动轨迹一致化。这是第一步相位相乘。用以改变线调频率尺度的Chirp Scaling二次相位函数为: (22)3.通过距离向FFT将数据变到二维频域。4.通过与参考函数进行相位相乘.同时完成距离压缩、SRC和一致RCMC。这是第二步相位相乘。用于距离压缩.距离徙动校正的相位函数写为: (23)5.通过距离向IFFT将数据变回到距离多普勒域。6.通过与随距离变化的匹配滤波器进行相位相乘.实现方位压缩。此外.由于步骤2中的Chirp Sca

13、ling操作.相位相乘中还需要附加一项相位校正。这是第三步相位相乘。补偿由Chirp Scaling引起的剩余相位函数是:(24)7.最后通过方位向IFFT将数据变回到二维时域.即SAR图像域。图8Chirp Scaling算法流程图简而言之.R-D算法是将徙动曲线逐一校正.CS算法是以某一徙动曲线为参考.在Doppler域内消除不同距离门的徙动曲线的差异.令这些曲线成为一组相互“平行”的曲线.然后在二维频率域内统一的去掉距离徙动。通俗一点就是.RD算法是将弯曲的信号一根根掰直.而CS算法是先把所有信号都掰得一样弯.然后再统一掰直。6 仿真结果6.1 使用最近邻点插值的距离多普勒算法仿真结果本

14、文首先对5个点目标的回波信号进行了仿真.5个点目标构成了矩形的4个顶点和中心.其坐标分别如下.格式为(方位向,距离向,后向反射系数): 0 9750 1 100 9750 1 50 10000 1 0 10250 1 100 10250 1图9的上图是距离向压缩后的图像.从图中可以看到5条回波信号(其中有几条部分重合.但仍能看出来)目标回波信号存在明显的距离徙动.需要进行校正。图9的下图是通过最近邻点插值法校正后的图像.可以看出图像基本被校正为直线。图9距离向压缩后最近邻点插值的结果图10为进行方位向压缩后形成的图像.可以明显看出5个点目标.并且5个点目标构成了矩形的四个顶点及其中心。图10

15、通过最近邻点插值生成的点目标图像6.2 使用最近邻点插值的距离多普勒算法仿真结果图11上图为通过距离压缩后的图像.图11的下图为通过sinc插值法校正后的图像。图11 距离向压缩后sinc插值的结果图12为进行方位向压缩后形成的图像.可以明显看出5个点目标.并且5个点目标构成了矩形的四个顶点及其中心。图12 通过sinc插值生成的点目标图像6.3 Chirp Scaling算法仿真结果同样.在Chirp Scaling中.对5个点目标的回波信号进行了仿真.5个点目标构成了矩形的4个顶点和中心.其坐标分别如下.格式为(方位向,距离向,后向反射系数): 1200 0 1 1250 -50 1 12

16、50 50 1 1150 -50 1 1150 50 1图13是仿真的雷达回波信号图图13 仿真出来的SAR回波信号图14是经过第一次相位校正之后.通过距离向压缩后的距离时域-方位时域信号图(Chirp Scaling算法的七个步骤中并不包含该信号.该信号是将步骤2之后的信号通过方位向傅里叶逆变换.再进行距离向压缩得到的.只为了验证原理)。按照理论.该图中所有点的距离徙动都应该一样。从图中大致可看出.五个点的距离徙动是差不多的。图14 Chirp Scaling第2步之后、经过距离向压缩得到的图图15为步骤5之后.信号距离压缩.距离徙动校正之后的距离多普勒域中的信号图。图15 距离徙动校正之后

17、的图图16为步骤6之后.消除相位偏移的图。图16 消除相位偏移的图图17为通过Chirp Scaling算法生成的点目标图像图17通过Chirp Scaling算法生成的点目标图像6.4 几种算法比较本文讨论了距离多普勒算法和Chirp Scaling算法.其中距离多普勒算法考虑了最近邻点插值和sinc插值两种插值方法。距离多普勒算法兼顾了成熟、简单、高效和精确等因素.至今仍被广泛使用.但是距离多普勒算法有两点不足:首先.当用较长的核函数提高距离徙动校正(RCMC)精度时.运算量较大;其次.二次距离压缩(SRC)对方位频率的依赖性问题较难解决.从而限制了其对某些大斜视角和长孔径SAR的处理精度

18、。最近邻点插值的优点是速度快.该插值的运行时间为2.秒.缺点是不够精确;sinc插值的优点是精确.该方法的运行时间为29.秒.缺点是速度慢;Chirp Scaling算法避免了插值运算.提高了速度.运行时间为0.秒.但是其算法较为复杂。%=%文件名:NearSAR.m%作者:徐一凡%功能:合成孔径雷达距离多普勒算法点目标成像 %=clear;clc;close all;%=%常数定义C=3e8; %光速%雷达参数Fc=1e9; %载频1GHzlambda=C/Fc; %波长%目标区域参数Xmin=0; %目标区域方位向范围Xmin,XmaxXmax=50; Yc=10000; %成像区域中线Y

19、0=500; %目标区域距离向范围Yc-Y0,Yc+Y0 %成像宽度为2*Y0%轨道参数V=100; %SAR的运动速度100 m/sH=5000; %高度 5000 mR0=sqrt(Yc2+H2); %最短距离%天线参数D=4; %方位向天线长度Lsar=lambda*R0/D; %SAR合成孔径长度.合成孔径雷达成像算法与实现P.100Tsar=Lsar/V; %SAR照射时间%慢时间域参数Ka=-2*V2/lambda/R0; %多普勒频域调频率P.93Ba=abs(Ka*Tsar); %多普勒频率调制带宽PRF=Ba; %脉冲重复频率.PRF其实为多普勒频率的采样率.又为复频率.所以

20、等于Ba.P.93PRT=1/PRF; %脉冲重复时间ds=PRT; %慢时域的时间步长Nslow=ceil(Xmax-Xmin+Lsar)/V/ds); %慢时域的采样数.ceil为取整函数.结合P.76的图理解Nslow=2nextpow2(Nslow); %nextpow2为最靠近2的幂次函数.这里为fft变换做准备sn=linspace(Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢时间域的时间矩阵PRT=(Xmax-Xmin+Lsar)/V/Nslow; %由于Nslow改变了.所以相应的一些参数也需要更新.周期减小了PRF=1/PRT;ds=PRT

21、;%快时间域参数设置Tr=5e-6; %脉冲持续时间5usBr=30e6; %chirp频率调制带宽为30MHzKr=Br/Tr; %chirp调频率Fsr=2*Br; %快时域采样频率.为3倍的带宽dt=1/Fsr; %快时域采样间隔Rmin=sqrt(Yc-Y0)2+H2);Rmax=sqrt(Yc+Y0)2+H2+(Lsar/2)2);Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%快时域的采样数量Nfast=2nextpow2(Nfast); %更新为2的幂次.方便进行fft变换tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);

22、 %快时域的离散时间矩阵dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; %更新间隔Fsr=1/dt;%分辨率参数设置DY=C/2/Br; %距离向分辨率DX=D/2; %方位向分辨率%点目标参数设置Ntarget=5; %点目标的数量%点目标格式x,y,反射系数sigmaPtarget=Xmin,Yc-50*DY,1 %点目标位置.这里设置了5个点目标.构成一个矩形以及矩形的中心 Xmin+50*DX,Yc-50*DY,1 Xmin+25*DX,Yc,1 Xmin,Yc+50*DY,1 Xmin+50*DX,Yc+50*DY,1; disp(Parameters:) %参数显

23、示disp(Sampling Rate in fast-time domain);disp(Fsr/Br)disp(Sampling Number in fast-time domain);disp(Nfast)disp(Sampling Rate in slow-time domain);disp(PRF/Ba)disp(Sampling Number in slow-time domain);disp(Nslow)disp(Range Resolution);disp(DY)disp(Cross-range Resolution);disp(DX) disp(SAR integration

24、 length);disp(Lsar) disp(Position of targets);disp(Ptarget)%=%生成回波信号K=Ntarget; %目标数目N=Nslow; %慢时域的采样数M=Nfast; %快时域的采样数T=Ptarget; %目标矩阵Srnm=zeros(N,M); %生成零矩阵存储回波信号for k=1:1:K %总共K个目标 sigma=T(k,3); %得到目标的反射系数 Dslow=sn*V-T(k,1); %方位向距离.投影到方位向的距离 R=sqrt(Dslow.2+T(k,2)2+H2); %实际距离矩阵 tau=2*R/C; %回波相对于发射波

25、的延时 Dfast=ones(N,1)*tm-tau*ones(1,M); %(t-tau).其实就是时间矩阵.ones(N,1)和ones(1,M)都是为了将其扩展为矩阵 phase=pi*Kr*Dfast.2-(4*pi/lambda)*(R*ones(1,M);%相位.公式参见P.96 Srnm=Srnm+sigma*exp(j*phase).*(0Dfast&DfastTr).*(abs(Dslow)Lsar/2)*ones(1,M);%由于是多个目标反射的回波.所以此处进行叠加end%=%距离-多普勒算法开始%距离向压缩tic;tr=tm-2*Rmin/C;Refr=exp(j*pi

26、*Kr*tr.2).*(0tr&trM %判断是否超出边界 Sa_RD(n,m)=Sa_RD(n,M/2); else if delta_RMC=0.5 %五入 Sa_RD(n,m)=Sa_RD(n,m+round(RMC)+1); else %四舍 Sa_RD(n,m)=Sa_RD(n,m+round(RMC); end end end waitbar(n/N)endclose(h)%=Sr_rmc=iftx(Sa_RD); %距离徙动校正后还原到时域Ga = abs(Sr_rmc);%方位向压缩ta=sn-Xmin/V;Refa=exp(j*pi*Ka*ta.2).*(abs(ta)Tsa

27、r/2);Sa=iftx(ftx(Sr_rmc).*(conj(ftx(Refa).*ones(1,M);Gar=abs(Sa);toc;%=%绘图colormap(gray);figure(1)subplot(211);row=tm*C/2-2008;col=sn*V-26;imagesc(row,col,255-Gr); %距离向压缩.未校正距离徙动的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel(距离向),ylabel(方位向),title(距离向压缩.未校正距离徙动的图像),subplot(212);imagesc(row,col,

28、255-Ga); %距离向压缩.校正距离徙动后的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel(距离向),ylabel(方位向),title(距离向压缩.校正距离徙动后的图像),figure(2)colormap(gray);imagesc(row,col,255-Gar); %方位向压缩后的图像axis(Yc-Y0,Yc+Y0,Xmin-Lsar/2,Xmax+Lsar/2);xlabel(距离向),ylabel(方位向),title(方位向压缩后的图像),%=%文件名:SincSAR.m%作者:徐一凡%功能:合成孔径雷达距离多普勒算法点

29、目标成像 %=clear;clc;close all;%=%常数定义C=3e8; %光速%雷达参数Fc=1e9; %载频1GHzlambda=C/Fc; %波长%目标区域参数Xmin=0; %目标区域方位向范围Xmin,XmaxXmax=50; Yc=10000; %成像区域中线Y0=500; %目标区域距离向范围Yc-Y0,Yc+Y0 %成像宽度为2*Y0%轨道参数V=100; %SAR的运动速度100 m/sH=5000; %高度 5000 mR0=sqrt(Yc2+H2); %最短距离%天线参数D=4; %方位向天线长度Lsar=lambda*R0/D; %SAR合成孔径长度.合成孔径雷

30、达成像算法与实现P.100Tsar=Lsar/V; %SAR照射时间%慢时间域参数Ka=-2*V2/lambda/R0; %多普勒频域调频率P.93Ba=abs(Ka*Tsar); %多普勒频率调制带宽PRF=Ba; %脉冲重复频率.PRF其实为多普勒频率的采样率.又为复频率.所以等于Ba.P.93PRT=1/PRF; %脉冲重复时间ds=PRT; %慢时域的时间步长Nslow=ceil(Xmax-Xmin+Lsar)/V/ds); %慢时域的采样数.ceil为取整函数.结合P.76的图理解Nslow=2nextpow2(Nslow); %nextpow2为最靠近2的幂次函数.这里为fft变换

31、做准备sn=linspace(Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow);%慢时间域的时间矩阵PRT=(Xmax-Xmin+Lsar)/V/Nslow; %由于Nslow改变了.所以相应的一些参数也需要更新.周期减小了PRF=1/PRT;ds=PRT;%快时间域参数设置Tr=5e-6; %脉冲持续时间5usBr=30e6; %chirp频率调制带宽为30MHzKr=Br/Tr; %chirp调频率Fsr=2*Br; %快时域采样频率.为3倍的带宽dt=1/Fsr; %快时域采样间隔Rmin=sqrt(Yc-Y0)2+H2);Rmax=sqrt(Yc+Y0)2+H

32、2+(Lsar/2)2);Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%快时域的采样数量Nfast=2nextpow2(Nfast); %更新为2的幂次.方便进行fft变换tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast); %快时域的离散时间矩阵dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; %更新间隔Fsr=1/dt;%分辨率参数设置DY=C/2/Br; %距离向分辨率DX=D/2; %方位向分辨率%点目标参数设置Ntarget=5; %点目标的数量%点目标格式x,y,反射系数sigmaPtarget=Xmin,Y

33、c-50*DY,1 %点目标位置.这里设置了5个点目标.构成一个矩形以及矩形的中心 Xmin+50*DX,Yc-50*DY,1 Xmin+25*DX,Yc,1 Xmin,Yc+50*DY,1 Xmin+50*DX,Yc+50*DY,1; disp(Parameters:) %参数显示disp(Sampling Rate in fast-time domain);disp(Fsr/Br)disp(Sampling Number in fast-time domain);disp(Nfast)disp(Sampling Rate in slow-time domain);disp(PRF/Ba)disp(Sampling Number in slow-time domain);disp(Nslow)disp(Range Resolution);disp(DY)disp(Cross-range Resolution);disp(DX) disp(SAR integration length);disp(Lsar) disp(Position of targets);disp(Ptarget)%=%生成回波信号K=Ntar

展开阅读全文
相关资源
相关搜索

当前位置:首页 > 教育专区 > 教案示例

本站为文档C TO C交易模式,本站只提供存储空间、用户上传的文档直接被用户下载,本站只是中间服务平台,本站所有文档下载所得的收益归上传人(含作者)所有。本站仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私,请立即通知淘文阁网,我们立即给予删除!客服QQ:136780468 微信:18945177775 电话:18904686070

工信部备案号:黑ICP备15003705号© 2020-2023 www.taowenge.com 淘文阁