位置: IT常识 - 正文

开源代码 | FMCW-MIMO雷达仿真MATLAB(开源代码网站github)

编辑:rootadmin
开源代码 | FMCW-MIMO雷达仿真MATLAB

推荐整理分享开源代码 | FMCW-MIMO雷达仿真MATLAB(开源代码网站github),希望有所帮助,仅作参考,欢迎阅读内容。

文章相关热门搜索词:开源代码在哪找,开源代码网站,开源代码公布意味着什么,开源代码与组件使用情况说明怎么写,开源代码商用是否违法,开源代码是什么意思,开源代码与组件使用情况说明怎么写,开源代码公布意味着什么,内容如对您有帮助,希望把文章链接给更多的朋友!

本文编辑:调皮哥的小助理

本程序来源:https://github.com/ekurtgl/FMCW-MIMO-Radar-Simulation,作者是阿拉巴马大学博士生艾库特格尔,研究方向主要是雷达信号处理人类活动识别以及雷达数据的机器学习应用,这份比较新的开源雷达仿真代码,值得大家学习。

下面主要分析代码的主要内容,方便大家解读。

程序目录如下:

图片

FMCW_simulation.m是创建点目标并估计其范围、速度和角度信息的主脚本,首先研究这个脚本程序。

一、雷达参数

雷达参数的设置,属于是老生常谈了,之前的文章已经谈到很多了,不再详细重复论述。不过在本程序中,需要注意PRI默认为等于Chirp持续时长,实际上还要考虑空闲时间。这里的带宽指的是有效带宽,而不是发射信号带宽。程序中一共设置了10帧、1T8R,这些参数都是可以修改的。

```c> %% Radar parametersc = physconst('LightSpeed'); %speed of lightBW = 150e6; %bandwidth 有效fc = 77e9; % carrier frequencynumADC = 256; % # of adc samplesnumChirps = 256; % # of chirps per framenumCPI = 10;T = 10e-6; % PRI,默认不存在空闲时间PRF = 1/T;Fs = numADC/T; % sampling frequencydt = 1/Fs; % sampling intervalslope = BW/T;lambda = c/fc;N = numChirps*numADC*numCPI; % total # of adc samplest = linspace(0,T*numChirps*numCPI,N); % time axis, one frame 等间隔时间/点数t_onePulse = 0:dt:dt*numADC-dt; %单chirp时间numTX = 1;numRX = 8; %等效后Vmax = lambda/(T*4); % Max Unamb velocity m/sDFmax = 1/2*PRF; % = Vmax/(c/fc/2); % Max Unamb Dopp FreqdR = c/(2*BW); % range resolRmax = Fs*c/(2*slope); % TI's MIMO Radar docRmax2 = c/2/PRF; % lecture 2.3dV = lambda/(2*numChirps*T); % velocity resol, lambda/(2*framePeriod)d_rx = lambda/2; % dist. between rxsd_tx = 4*d_rx; % dist. between txsN_Dopp = numChirps; % length of doppler FFTN_range = numADC; % length of range FFTN_azimuth = numTX*numRX;R = 0:dR:Rmax-dR; % range axisV = linspace(-Vmax, Vmax, numChirps); % Velocity axisang_ax = -90:90; % angle axis二、目标参数

这里的目标参数就采用极坐标系转直角坐标系,然后分别用X方向的速度和Y方向上的速度乘上时间,带入后续的计算。

> %% 目标参数r1_radial = 50; v1_radial = 10; % velocity 1tar1_angle = -10;r1_y = cosd(tar1_angle)*r1_radial;r1_x = sind(tar1_angle)*r1_radial;v1_y = cosd(tar1_angle)*v1_radial;v1_x = sind(tar1_angle)*v1_radial;r1 = [r1_x r1_y 0];r2_radial = 100;v2_radial = -15; % velocity 2tar2_angle = 10;r2_y = cosd(tar2_angle)*r2_radial;r2_x = sind(tar2_angle)*r2_radial;v2_y = cosd(tar2_angle)*v2_radial;v2_x = sind(tar2_angle)*v2_radial;r2 = [r2_x r2_y 0];%发射天线位置tx_loc = cell(1,numTX);for i = 1:numTX tx_loc{i} = [(i-1)*d_tx 0 0]; scatter3(tx_loc{i}(1),tx_loc{i}(2),tx_loc{i}(3),'b','filled') hold onend% 接收天线位置rx_loc = cell(1,numRX);for i = 1:numRX rx_loc{i} = [tx_loc{numTX}(1)+d_tx+(i-1)*d_rx 0 0]; scatter3(rx_loc{i}(1),rx_loc{i}(2),rx_loc{i}(3),'r','filled')endtar1_loc = zeros(length(t),3);tar1_loc(:,1) = r1(1) + v1_x*t;tar1_loc(:,2) = r1(2) + v1_y*t;tar2_loc = zeros(length(t),3);tar2_loc(:,1) = r2(1) + v2_x*t;tar2_loc(:,2) = r2(2) + v2_y*t;

同时,这里还绘制了发射天线和接收天线的位置图,蓝色为发射天线位置,红色为接收天线位置。这个位置信息后面会用于计算目标的时间延迟。

三、发射信号建模

利用收发天线的位置以及目标参数中雷达的位置信息,先求目标到雷达的2-范数(也就是空间中两点的直线距离),然后转化为目标的延迟时间τ,如此以来得到的信号模型精度更高!这样的方法,与往常采用的不同,可以说以前的方法有点粗暴,虽同为远场条件,但是没有考虑阵列的位置对信号的影响!

> %% TX siganldelays_tar1 = cell(numTX,numRX);delays_tar2 = cell(numTX,numRX);r1_at_t = cell(numTX,numRX);r2_at_t = cell(numTX,numRX);tar1_angles = cell(numTX,numRX);tar2_angles = cell(numTX,numRX);tar1_velocities = cell(numTX,numRX);tar2_velocities = cell(numTX,numRX);for i = 1:numTX for j = 1:numRX delays_tar1{i,j} = (vecnorm(tar1_loc-repmat(rx_loc{j},N,1),2,2)+vecnorm(tar1_loc-repmat(tx_loc{i},N,1),2,2))/c; delays_tar2{i,j} = (vecnorm(tar2_loc-repmat(rx_loc{j},N,1),2,2)+vecnorm(tar2_loc-repmat(tx_loc{i},N,1),2,2))/c; endend四、接收信号模型

这里,接收信号没有采用发射与接收混频的形式,而是相位直接做差,分别计算两个目标的中频信号相加,此法等效为混频。其中@(tx,fx)相当于传参函数,用于后面的等式计算。

中频信号除了目标的距离表征的频率之外,还有目标运动速度引起的多普勒频移,是两部分的叠加。

> %% Complex signalphase = @(tx,fx) 2*pi*(fx.*tx+slope/2*tx.2); % transmittedphase2 = @(tx,fx,r,v) 2*pi*(2*fx*r/c+tx.*(2*fx*v/c + 2*slope*r/c)); % downconverted% f_oneChirp = slope*t(1:sum(t<=T));% f_t = repmat(f_oneChirp,1,numChirps*numCPI)-(BW/2); % transmit freq% f_t = BW/2*sawtooth(t/T*2*pi); fr1 = 2*r1(2)*slope/c; fr2 = 2*r2(2)*slope/c;fd1 = 2*v1_radial*fc/c; % doppler freqfd2 = 2*v2_radial*fc/c;f_if1 = fr1 + fd1; % beat or IF freqf_if2 = fr2 + fd2;% mixed1 = cell(numTX,numRX);% mixed2 = cell(numTX,numRX);mixed = cell(numTX,numRX);for i = 1:numTX for j = 1:numRX disp(['Processing Channel: ' num2str(j) '/' num2str(numRX)]); for k = 1:numChirps*numCPI phase_t = phase(t_onePulse,fc); phase_1 = phase(t_onePulse-delays_tar1{i,j}(k*numADC),fc); % received phase_2 = phase(t_onePulse-delays_tar2{i,j}(k*numADC),fc); signal_t((k-1)*numADC+1:k*numADC) = exp(1j*phase_t); signal_1((k-1)*numADC+1:k*numADC) = exp(1j*(phase_t - phase_1)); signal_2((k-1)*numADC+1:k*numADC) = exp(1j*(phase_t - phase_2)); end mixed{i,j} = signal_1 + signal_2; endend

绘制出的局部结果如下,若要观察全部的信号,则需要修改X轴的范围限制。

五、2D-FFT

上面我一句雷达信号处理原理方面的内容都没有提及,主要是这些基础的内容,我默认各位都已经清楚了。如果不清楚,那就先去搞清楚再回头来看(调皮连续波公众号里有),不然可能会遇到理解困难。

首先生成一个三维的DataCube,然后每个CPI做一次2D-FFT就OK了,没啥难度,顶多就是数据组成或者格式,理解起来有些麻烦,不过这都不是核心难题。

> RDC = reshape(cat(3,mixed{:}),numADC,numChirps*numCPI,numRX*numTX); % radar data cubeRDMs = zeros(numADC,numChirps,numTX*numRX,numCPI);for i = 1:numCPI RD_frame = RDC(:,(i-1)*numChirps+1:i*numChirps,:); RDMs(:,:,:,i) = fftshift(fft2(RD_frame,N_range,N_Dopp),2);endfigureimagesc(V,R,20*log10(abs(RDMs(:,:,1,1))/max(max(abs(RDMs(:,:,1,1))))));colormap(jet(256))% set(gca,'YDir','normal')clim = get(gca,'clim');caxis([clim(1)/2 0])xlabel('Velocity (m/s)');ylabel('Range (m)');

处理结果如下:

化成三维的可以标注数据,方便查看解算的结果以及误差,可见距离和速度都基本吻合,存在些许误差。

六、CA-CFAR开源代码 | FMCW-MIMO雷达仿真MATLAB(开源代码网站github)

代码一读就懂,老生常谈的东西了,没啥可讲的。

> numGuard = 2; % # of guard cellsnumTrain = numGuard*2; % # of training cellsP_fa = 1e-5; % desired false alarm rate SNR_OFFSET = -5; % dBRDM_dB = 10*log10(abs(RDMs(:,:,1,1))/max(max(abs(RDMs(:,:,1,1)))));[RDM_mask, cfar_ranges, cfar_dopps, K] = ca_cfar(RDM_dB, numGuard, numTrain, P_fa, SNR_OFFSET);figureh=imagesc(V,R,RDM_mask);xlabel('Velocity (m/s)')ylabel('Range (m)')title('CA-CFAR')

处理结果:

七、角度估计(一)3D-FFT

常规操作,基本内容,没啥难度。

```c> rangeFFT = fft(RDC(:,1:numChirps,:),N_range);angleFFT = fftshift(fft(rangeFFT,length(ang_ax),3),3);range_az = squeeze(sum(angleFFT,2)); % range-azimuth mapfigurecolormap(jet)% imagesc(ang_ax,R,20*log10(abs(range_az)./max(abs(range_az(:))))); mesh(ang_ax,R,20*log10(abs(range_az)./max(abs(range_az(:))))); xlabel('Azimuth Angle')ylabel('Range (m)')title('FFT Range-Angle Map')set(gca,'clim', [-35, 0])doas = zeros(K,181); % direction of arrivalsfigurehold on; grid on;for i = 1:K doas(i,:) = fftshift(fft(rangeFFT(cfar_ranges(i),cfar_dopps(i),:),181)); plot(ang_ax,10*log10(abs(doas(i,:))))endxlabel('Azimuth Angle')ylabel('dB')

运行结果如下: 可见,角度并不是很准确,读者可以验证,目标角度越偏离发现,误差越大。这其实与雷达的角度分辨率有关,用下文的超分辨算法可以改善。

(二)MUSIC算法

有读者常问,说MUSIC在网上找到的代码都是仿真的,很少看到有对实际目标的数据进行处理的,这不就来了嘛。不过,说到底还是见的太少。

MUSIC算法原理默认都清楚,不清楚的自己查阅相关资料。

> d = 0.5;M = numCPI; % # of snapshotsfor k=1:length(ang_ax) a1(:,k)=exp(-1i*2*pi*(d*(0:numTX*numRX-1)'*sin(ang_ax(k).'*pi/180)));endfor i = 1:K Rxx = zeros(numTX*numRX,numTX*numRX); for m = 1:M A = squeeze(RDMs(cfar_ranges(i),cfar_dopps(i),:,m)); Rxx = Rxx + 1/M * (A*A'); end [Q,D] = eig(Rxx); % Q: eigenvectors (columns), D: eigenvalues [D, I] = sort(diag(D),'descend'); Q = Q(:,I); % Sort the eigenvectors to put signal eigenvectors first Qs = Q(:,1); % Get the signal eigenvectors Qn = Q(:,2:end); % Get the noise eigenvectors for k=1:length(ang_ax) music_spectrum(i,k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*(Qn*Qn')*a1(:,k)); endend

运行结果漂亮的很,不要管幅度差异,幅度不代表功率。

(三)点云生成

坐标转换即可,若想要点云数量更多,可以降低CFAR门限,放更多的目标点进来。

> [~, I] = max(music_spectrum(2,:));angle1 = ang_ax(I);[~, I] = max(music_spectrum(1,:));angle2 = ang_ax(I);coor1 = [cfar_ranges(2)*sind(angle1) cfar_ranges(2)*cosd(angle1) 0];coor2 = [cfar_ranges(1)*sind(angle2) cfar_ranges(1)*cosd(angle2) 0];figurehold on;title('3D Coordinates (Point Cloud) of the targets')scatter3(coor1(1),coor1(2),coor1(3),100,'m','filled','linewidth',9)scatter3(coor2(1),coor2(2),coor2(3),100,'b','filled','linewidth',9)xlabel('Range (m) X')ylabel('Range (m) Y')zlabel('Range (m) Z')

处理结果:

(四)MUSIC 距离-AOA谱

算法没变,绘图方式不同而已。

> rangeFFT = fft(RDC);for i = 1:N_range Rxx = zeros(numTX*numRX,numTX*numRX); for m = 1:M A = squeeze(sum(rangeFFT(i,(m-1)*numChirps+1:m*numChirps,:),2)); Rxx = Rxx + 1/M * (A*A'); end% Rxx = Rxx + sqrt(noise_pow/2)*(randn(size(Rxx))+1j*randn(size(Rxx))); [Q,D] = eig(Rxx); % Q: eigenvectors (columns), D: eigenvalues [D, I] = sort(diag(D),'descend'); Q = Q(:,I); % Sort the eigenvectors to put signal eigenvectors first Qs = Q(:,1); % Get the signal eigenvectors Qn = Q(:,2:end); % Get the noise eigenvectors for k=1:length(ang_ax) music_spectrum2(k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*(Qn*Qn')*a1(:,k)); end range_az_music(i,:) = music_spectrum2;endfigurecolormap(jet)imagesc(ang_ax,R,20*log10(abs(range_az_music)./max(abs(range_az_music(:))))); xlabel('Azimuth')ylabel('Range (m)')title('MUSIC Range-Angle Map')clim = get(gca,'clim');

(五)压缩感知

原理捎带复杂,但主要还是一个凸优化问题。

> figurehold on; grid on;title('Angle Estimation with Compressed Sensing')xlabel('Azimuth')ylabel('dB')for i = 1:K A = squeeze(RDMs(cfar_ranges(i),cfar_dopps(i),:,1)); cvx_begin variable s(numTheta) complex; %alphax(numTheta,1) phix(numTX*numRX,numTheta)... %cap_theta(numTX*numRX,numTheta) %B(numTX*numRX,numTheta)%psix(numTheta,numTheta) %A(numRX*numTX,1) % A is the initial measurement % cap_theta == phix * psix; % minimize(norm(alphax,1)) % pow_p(norm(A-cap_theta*alphax,2),2) <= 1; % norm(A-cap_theta*alphax,2) <= 1; % minimize(norm(A-cap_theta*alphax,1)) minimize(norm(s,1)) norm(A-B*s,2) <= 1; cvx_end cvx_status cvx_optval plot(ang_ax,10*log10(abs(s)))end

第一个脚本的内容就说完了,下面是第二个脚本。

FMCW_sim_v2.m是读取 Kinect v2 设备捕获的人体骨骼关节三维坐标,提取原始雷达数据立方体 (RDC) 并播放距离多普勒图、输出微多普勒频谱图的主脚本。

Kinect v2里人体骨骼结构如下:

距离多普勒图:

微多普勒频谱图:

> rBin = 1:256;nfft = 212;window = 256;noverlap = 200;shift = window - noverlap;sx = myspecgramnew(sum(RDC(rBin,:,:)),window,nfft,shift); % mti filter and IQ correctionsx2 = abs(flipud(fftshift(sx,1)));timeAxis = [1:numCPI]*frameDuration; % TimefreqAxis = linspace(-PRF/2,PRF/2,nfft); % Frequency Axisfig = figure('visible','on');colormap(jet(256));% set(gca,'units','normalized','outerposition',[0,0,1,1]);doppSignMTI = imagesc(timeAxis,[-PRF/2 PRF/2],20*log10(abs(sx2/max(sx2(:)))));% axis xy% set(gca,'FontSize',10)title('micro-Doppler Spectrogram');% title(fOut(end-22:end-4))xlabel('Time (sec)');ylabel('Frequency (Hz)');caxis([-45 0]) % 40set(gca, 'YDir','normal')set(gcf,'color','w');

对该开源仿真代码感兴趣的读者,可以访问链接:https://github.com/ekurtgl/FMCW-MIMO-Radar-Simulation。

修改后的代码放到云盘了,修改的不多,如有需要,后台回复“2022”获得下载链接。

【本期结束】

本文链接地址:https://www.jiuchutong.com/zhishi/299086.html 转载请保留说明!

上一篇:【Tensorflow深度学习】实现手写字体识别、预测实战(附源码和数据集 超详细)(tensorflow gan)

下一篇:安装nvm,并使用nvm安装nodejs及配置环境变量(nvme安装win10教程)

  • 华硕x43s驱动(华硕x43s配置)(华硕x42j显卡驱动)

    华硕x43s驱动(华硕x43s配置)(华硕x42j显卡驱动)

  • 虾头变黑是因为重金属超重吗(虾头变黑是因为重金属超标吗)(虾头变黑因为什么)

    虾头变黑是因为重金属超重吗(虾头变黑是因为重金属超标吗)(虾头变黑因为什么)

  • wifi开启ap隔离什么意思(无线网设置开启ap隔离什么意思)

    wifi开启ap隔离什么意思(无线网设置开启ap隔离什么意思)

  • 华为mate30pro怎么重启手机(华为mate30pro怎么刷机)

    华为mate30pro怎么重启手机(华为mate30pro怎么刷机)

  • 微信朋友圈发送失败轻触以再次发送(微信朋友圈发送失败是什么原因)

    微信朋友圈发送失败轻触以再次发送(微信朋友圈发送失败是什么原因)

  • 32位操作系统,基于x64的处理器是什么意思(32位操作系统基于x6处理器装系统)

    32位操作系统,基于x64的处理器是什么意思(32位操作系统基于x6处理器装系统)

  • 华硕路由器恢复出厂(华硕路由器恢复出厂设置)

    华硕路由器恢复出厂(华硕路由器恢复出厂设置)

  • 全网通充话费能退钱吗(全网通话费充值卡可以全部通用吗)

    全网通充话费能退钱吗(全网通话费充值卡可以全部通用吗)

  •   淘宝如何刷单(淘宝如何刷人气)

      淘宝如何刷单(淘宝如何刷人气)

  • 如何消除抖音号水印(如何消除抖音视频里的抖音号)

    如何消除抖音号水印(如何消除抖音视频里的抖音号)

  • 为什么投影仪和电脑不能同时显示视频(为什么投影仪和电脑有色差)

    为什么投影仪和电脑不能同时显示视频(为什么投影仪和电脑有色差)

  • 抖音怎么显示好友在线(抖音怎么显示好友在线时间)

    抖音怎么显示好友在线(抖音怎么显示好友在线时间)

  • 手机4G旁边突然出来个加号(手机4g旁边突然没信号)

    手机4G旁边突然出来个加号(手机4g旁边突然没信号)

  • qq闺蜜标识怎么弄(qq闺蜜标识怎么不显示)

    qq闺蜜标识怎么弄(qq闺蜜标识怎么不显示)

  • 数位屏不连电脑可以用吗(数位屏不连电脑会坏吗)

    数位屏不连电脑可以用吗(数位屏不连电脑会坏吗)

  • 小米传送门长按用不了(小米传送门长按识别无法启用)

    小米传送门长按用不了(小米传送门长按识别无法启用)

  • ipad pro可以接鼠标吗(ipad pro可以接鼠标键盘吗)

    ipad pro可以接鼠标吗(ipad pro可以接鼠标键盘吗)

  • 手机vsco怎么注册账号(vsco怎么用手机号注册)

    手机vsco怎么注册账号(vsco怎么用手机号注册)

  • 美图t9参数配置(美图t9手机参数配置)

    美图t9参数配置(美图t9手机参数配置)

  • wps如何删除表格中的线(wps如何删除表格横线)

    wps如何删除表格中的线(wps如何删除表格横线)

  • 小米9多少倍变焦(小米九几倍光学变焦)

    小米9多少倍变焦(小米九几倍光学变焦)

  • win10固态128g够用吗(win10系统128g固态硬盘)

    win10固态128g够用吗(win10系统128g固态硬盘)

  • 华为怎样设置台后运行(华为手机怎么调台)

    华为怎样设置台后运行(华为手机怎么调台)

  • DXEnum.exe是什么进程  DXEnum进程是什么病毒(电脑dwm.exe是什么)

    DXEnum.exe是什么进程 DXEnum进程是什么病毒(电脑dwm.exe是什么)

  • 免税收入对应的成本费用可以扣除吗
  • 增值税发票作废了税钱退还吗
  • 增值税明细账用三栏式还是多栏式
  • 租赁办公场所的请示
  • 增值税税控系统专用设备
  • 销售清单要装订吗
  • 咨询服务费预收率怎么算
  • 合并报表成本法调整为权益法
  • 固定资产投资会计处理
  • 公司租房水电费如何入帐
  • 发票能加盖公章吗
  • 企业所得税季度预缴怎么计算
  • 材料款付款流程
  • 融资租赁出租人和承租人的区别
  • 单独计价入账的土地是无形资产吗
  • 固定资产交换的账务处理
  • 机械设备租赁征求意见稿
  • 应付账款多付了
  • 季度所得税申报表营业收入怎么填
  • 简化丢失专用发票的处理流程
  • 供应商把价格算错了怎么办
  • 行政事业单位绿化工程完工应结转资产吗
  • 如何接收银行承兑
  • 公司员工结婚礼金规定
  • win10禁止使用网络
  • win11如何启用远程访问
  • 减少子公司合并的风险
  • Win10升级win11是激活状态吗?
  • 印花税的会计处理是什么
  • mac电脑遇到问题而重新启动
  • laravel批量更新数据
  • 预付款多长时间
  • php字符串变量
  • tim模块
  • 电脑椅什么牌子质量好
  • 购货折让会计分录
  • 增值税减免后按3%记提税款吗
  • ipad最长使用时间
  • 应付债券的应付利息怎么计算
  • 融资租赁的固定资产所有权归谁
  • 注册教育培训机构需要哪些手续
  • 基于Python实现五大常用分类算法(原理+代码)
  • chatGPT身份指令
  • 反射机制python
  • 库存盘点差异会计分录
  • 坏账的确认条件税务会查吗
  • 月末制造费用可以有余额吗
  • 实收资本确认时间什么入账
  • 小企业会计准则和企业会计准则的区别
  • 即征即退相关证明材料
  • 企业固定资产账户核算的固定资产其所有权均属于本企业
  • 公允价值变动是一级科目吗
  • 金税四期上线后对个人的影响
  • 工业企业新产品收入的调研报告
  • 企业广告费以后会涨吗
  • 积分换物品是真的吗
  • 固定资产清理主动
  • 免税设备是什么意思
  • 已认证进项税转出口退税怎么处理
  • 教育费附加免征文件
  • 会计可以用收据做账吗
  • 缴纳的车辆购置税需要计提吗
  • centos6.9关闭防火墙命令
  • linux在服务器领域的应用状况
  • xwizard.exe是什么
  • ssms.exe是什么进程
  • xp的命令提示符在哪
  • win1020h2版本是正式版吗
  • Linux系统下的Samba客户端配置
  • win10registry
  • 跑酷角色左右移动怎么弄
  • linux中切换目录命令符
  • 获取nodejs命令行信息
  • 你应该知道的2000个地理常识
  • Unity destructor Or OnDestory
  • python放歌
  • 国家税务局发票验证查询系统
  • 常规巡察和专项巡察相结合
  • 国家税务总局江西省
  • 土地增值税发票加计扣除5%年限
  • 免责声明:网站部分图片文字素材来源于网络,如有侵权,请及时告知,我们会第一时间删除,谢谢! 邮箱:opceo@qq.com

    鄂ICP备2023003026号

    网站地图: 企业信息 工商信息 财税知识 网络常识 编程技术

    友情链接: 武汉网站建设