网站公告设计,深圳市绿色建筑信息平台,本溪食品 中企动力提供网站建设,廊坊seo按天计费前一篇里#xff0c;获取了LORAwan的物理层波形#xff0c;并通过Octave查看了它的瞬时频率。LoRa是私有协议#xff0c;网上已经有了很不错的开源的实现#xff0c;如#xff1a;
S2_LoRa通信实验
LoRaPhy
以及GNU-Radio的Lora模块、LimeSDR的Lora实现。当我试图修改上…前一篇里获取了LORAwan的物理层波形并通过Octave查看了它的瞬时频率。LoRa是私有协议网上已经有了很不错的开源的实现如
S2_LoRa通信实验
LoRaPhy
以及GNU-Radio的Lora模块、LimeSDR的Lora实现。当我试图修改上述代码实现上一篇文章里获取的波形的还原时却无法得到正确的结果也找不到原因。为了知其所以然还是准备从头动手实现。由于是业余学习不打算解决所有LORA协议的接收处理问题只就事论事只做默认出厂配置SF-11 500K一类。
接收机是国产B205-mini或者B-210开发环境是开源的计算工具 GNU-Octave 。本文代码参考 a4_lora/octave 文章目录 0 .实验学习结论摘要1. 波形样本来源2. 粗略显示LoRa物理层突发的结构2.1 chirp的数学描述2.2 chirp的归一化离散形式 3. 在频差条件下检测头部3.1 钟差频差对头部检测几乎没有影响3.2 检测结果 4. 基于相关峰起伏周期的钟差估计算法4.1 低成本接收机的困境4.2 相关峰起伏周期的成因4.3 计算结果 5 基于 up-down chirp 边缘峰值交汇的时频差精确估计算法5.1 笨办法观察规律5.2 基于边缘最大值直线交点快速交汇的时频差估计 6 卡住位置补偿解调7 拉远距离加白噪声测试8 代码获取和运行 0 .实验学习结论摘要
虽然是站在LoRa-PHY、GNU-Radio巨人的肩膀上但是在近1个月的不断踩坑、尝试中还是非常的不顺利。最终得到正确结果时感觉就像攀过了一座陡峭的小山。
我们发现对于Lora这种使用Chirp的起始频率携带信息的低成本硬件其固有的时钟差、频率差对结果影响极大。尤其是使用另一个不靠谱的硬件来接收比如我的山寨 B205mini两者的钟差是不靠谱不靠谱离谱。能在Low SNR下解析出正确的数据需要对波形的深刻理解和认识每一个步骤都值得思考。对我学习的清华大学项目 LoRa-PHY里用于同步时钟误差的方法是假设主要的误差都来自LoRa设备用频差去折算钟差时差。在本实验里我们的接收设备也很不靠谱使得这个举措变得不稳定。我们通过观察相关峰值的最大值的周期变化发现了隐含在头部14组up-chirp里的低频周期提出了基于相关峰起伏周期的钟差估计算法并使用它对钟差进行补偿效果很好。对收发时钟、频率稳定度都很低的情形常见的二维时频相关算法耗时太大。提出基于 up-down chirp 边缘峰值交汇的时频差精确估计算法只在窗口的边缘检测峰值并两两连线计算交汇坐标。
本人已经退休不需要发Paper。谁如果论文里想引用这些idea只要标明参考出处即可。
1. 波形样本来源
样本来源就是上一节使用山寨B205mini/B210采集的IQ路数据。由于靠的很近几乎没有噪声。采样率是1MHz双字节 int16。LoRa的设备的带宽是500KHzSR11。
波形样本已经跟随代码一起放在Git服务器上。
2. 粗略显示LoRa物理层突发的结构
对于无噪声的波形可以采用相位差分的方式直接观察瞬时频率。读取一个突发数据 并利用复数求取角度差分就是相位差也就是瞬时频率。 上图表明我们购买的Lora模块默认开启了14个首部up-chirp2.25个down-chirp
%%%%%%%%%%%%%%%读取一个2倍速率采集的无噪声文件。
%% USRP B210 430MHz, 16bit IQ BW500KHz
%fid fopen(d:/lora/test.wav,rb);
fid fopen(data/iq_sf11_500KHz_spr1MHz.wav,rb);
raw_iq_1d fread(fid,int16);
fclose(fid);
%%组合复波形
raw_iq raw_iq_1d(1:2:end)1j*raw_iq_1d(2:2:end);
%幅度归一化
max_amp sum(abs(raw_iq))/length(raw_iq);
raw_iq raw_iq/max_amp;figure(1);
hold on
draw_freqs zeros(1,length(raw_iq));
draw_freqs(2:end) angle(raw_iq(2:end)) - angle(raw_iq(1:end-1));
draw_freqs(draw_freqs-pi ) draw_freqs(draw_freqs-pi ) 2*pi ;
draw_freqs(draw_freqspi ) draw_freqs(draw_freqspi ) - 2*pi ;
plot(draw_freqs,r);2.1 chirp的数学描述
根据LoRa的定义一个chirp 跨域的带宽是 B Hz持续的时间是 T 秒他们满足: B T 2 S F B T 2^{SF} BT2SF
在T秒内基带IQ向量旋转的频率从 -B/2 变化到 B/2。
由于频率 f ( t ) 2 π ( f 0 B T t ) f(t)2\pi(f_0 \frac{B}{T}t ) f(t)2π(f0TBt)
而向量的相位是频率的积分计算 p ( t ) ∫ f ( t ) d t 2 π ( f 0 B 2 T t ) t , t ⊂ [ 0 , T ) p(t)\int f(t)dt2\pi(f_0 \frac{B}{2T}t )t, t\subset [0,T) p(t)∫f(t)dt2π(f02TBt)t,t⊂[0,T)
此时样点的复平面向量可以表示为 s ( t ) e 2 π j ( f 0 B 2 T t ) t s(t)e^{2\pi j(f_0 \frac{B}{2T}t )t} s(t)e2πj(f02TBt)t
这个t的取值是 0 ~ T秒。
2.2 chirp的归一化离散形式
使用上述模拟波形的样式来计算chirp, 不容易出错。但如果换算为离散形式则更方便写程序。假设我们采用高于LoRa带宽B的K倍速率采样即 Sr BK样点为 n则参考本文开头的第一个式子里SF、B、T的关系一个chirp的样点数为 S Y M N ( S F , K ) 2 S F B S r 2 S F ⋅ K SYMN(SF,K)\frac {2^{SF}}{B} Sr2^{SF}\cdot K SYMN(SF,K)B2SFSr2SF⋅K 比如SF11, K2时N就是4096点。根据离散、连续的关系 t n S r n B K t\frac{n}{Sr}\frac{n}{BK} tSrnBKn p ( t ) 2 π ( f 0 B 2 T n B K ) n B K 2 π ( f 0 n 2 T ⋅ K ) n B K p(t)2\pi (f_0 \frac{B}{2T}\frac{n}{BK})\frac{n}{BK}2\pi (f_0 \frac{n}{2T\cdot K})\frac{n}{BK} p(t)2π(f02TBBKn)BKn2π(f02T⋅Kn)BKn
由于 B T 2 S F BT2^{SF} BT2SF, 有 p ( t ) 2 π ( f 0 B n 2 ⋅ 2 S F ⋅ K ) n K p(t)2\pi (\frac{f_0}{B} \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K} p(t)2π(Bf02⋅2SF⋅Kn)Kn
剩下一个 f0/B也可以兑换为比例其实f0取-B/2时他就是 -1/2。
经过上面的推导我们把 chirp 表示为和具体采样率、频率无关的值只和SF,K有关。 p u p − c h i r p ( t ) 2 π ( − 1 2 n 2 ⋅ 2 S F ⋅ K ) n K p_{up-chirp}(t)2\pi (-\frac{1}{2} \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K} pup−chirp(t)2π(−212⋅2SF⋅Kn)Kn S u p − c h i r p ( n ) e 2 π j ( − 1 2 n 2 ⋅ 2 S F ⋅ K ) n K S_{up-chirp}(n)e^{2\pi j (-\frac{1}{2} \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K}} Sup−chirp(n)e2πj(−212⋅2SF⋅Kn)Kn
我们可以使用Octave生成这样的chirp:
chirp_n 0:1:SYM-1;
chirp_up exp(2j*pi*(-1/2chirp_n/2./(2^SF)./K).*chirp_n./K);
plot(real(chirp_up));3. 在频差条件下检测头部
使用与上面的 chirp 对应的相关滤波器即可从噪声中检测头部的存在。也可以利用down-chirp 相乘后fft查看有无单音来检测那样速度更快。
检测这个步骤使用共轭能够和原始波形直接合成为很强的峰。 S u p − c h i r p ‾ ( n ) e − 2 π j ( − 1 2 n 2 ⋅ 2 S F ⋅ K ) n K \overline{S_{up-chirp}}(n)e^{-2\pi j (-\frac{1}{2} \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K}} Sup−chirp(n)e−2πj(−212⋅2SF⋅Kn)Kn
这个波形的瞬时频率为 f u p − c h i r p ‾ ( n ) − 2 π ( − 1 2 n 2 S F ⋅ K ) 1 K \overline{f_{up-chirp}}(n)-{2\pi (-\frac{1}{2} \frac{n}{2^{SF}\cdot K})\frac{1}{K}} fup−chirp(n)−2π(−212SF⋅Kn)K1
我们在第一阶段主要的工作就是用上面这个相关波形去和原始波形做相关两两相乘而后求和检测峰值。
3.1 钟差频差对头部检测几乎没有影响
即使在有频偏的情况下这种检测依然是有效的。只要频率偏移不是很离谱通过样点的移动总能用时间延迟的损失来抵消频偏
设存在一个时间偏移 n ′ δ n n\deltan n′δn, 以及一个固有的频率偏移 Δ \Delta Δ f u p − c h i r p ‾ ( n ′ ) f u p − c h i r p ‾ ( δ n ) − 2 π j ( Δ − 1 2 n δ 2 ⋅ 2 S F ⋅ K ) 1 K \overline{f_{up-chirp}}(n)\overline{f_{up-chirp}}(\deltan)-{2\pi j (\Delta-\frac{1}{2} \frac{n\delta}{2\cdot2^{SF}\cdot K})\frac{1}{K}} fup−chirp(n′)fup−chirp(δn)−2πj(Δ−212⋅2SF⋅Knδ)K1
当这个 δ \delta δ满足下式使得频偏 Δ \Delta Δ其恰好能被时间推移抵消 Δ δ 2 ⋅ 2 S F ⋅ K 0 \Delta\frac{\delta}{2\cdot2^{SF}\cdot K}0 Δ2⋅2SF⋅Kδ0
则结果变为 f u p − c h i r p ‾ ( n ′ ) − 2 π ( − 1 2 n ′ − δ 2 ⋅ 2 S F ⋅ K ) 1 K \overline{f_{up-chirp}}(n)-{2\pi (-\frac{1}{2} \frac{n-\delta}{2\cdot2^{SF}\cdot K})\frac{1}{K}} fup−chirp(n′)−2π(−212⋅2SF⋅Kn′−δ)K1
相关峰与无频偏时相比只是提前或者推迟了 δ \delta δ个样点幅度稍微有些降低而已可用的重叠频段不再是B Hz了。 3.2 检测结果
我们的头部检测结果如下 绿色的是瞬时频率蓝色的是相关峰黑色的是峰值标记位置。请格外注意这些蓝色峰值的起伏它很重要
4. 基于相关峰起伏周期的钟差估计算法
在第3节我们发现即使不去补偿时频差也能顺利检测到头部。但是LoRa靠频率的绝对起点f0来表示符号, 解扩、解调不做钟差、频差就得不到正确的结果。
4.1 低成本接收机的困境
学习资源 LoRa-PHY里认为钟差基本源自LoRa设备本身。这或许因为大学实力雄厚采集设备的时钟稳定而标准。此时可以沿用LoRa-PHY的算法使用频差/射频频率来兑换钟差。
然而我们的接收机成本很低本身的频差、钟差就很大和LoRa这种60多元的板子旗鼓相当。此时上述算法就不行了。LoRa-PHY抵消后存在一个微弱的残差使得在发送64字节数据时符号依旧会漂移3-4次。
4.2 相关峰起伏周期的成因
注意看3.2节的头部检测相关峰每个峰值的幅度并不一致 讲起来按照Lora的标准定义采用K倍速率采样不会有起伏。即使具有频偏相关位置的相位也是定值每次相关的幅度衰减是恒定的。为啥周期起伏还是个COS函数的样子那是因为时钟存在误差每次对齐的相位会发生漂移。
这个小细节蕴含着非常重要的一个参数正是LoRa设备的时钟和我们的接收机时钟之间的误差。这个误差使得原本应该 N 点1个符号的chirp接收到后可能变成 N0.12点如果不去补偿这个误差以N为周期去卡符号就会慢慢错位。错位后每次相关时相位初始值发生缓慢平移造成峰值周期性的起伏。
提取这个误差非常简单对这十几个点做FFT如对 14个点做1048576点的FFT并求取峰值峰值与窗口的比例偏移就是时钟偏移。我们把上图的顶部放大
可以看到这就是个余弦。用上面的数据做FFT求取补偿因子 cv_pay : vpk sqrt(var(prepeaks (:,2)));prepeaks (:,2) (prepeaks (:,2) - mpk)/vpk * 100;ck_fft prepeaks (:,2);cm_fft abs(fft(ck_fft,1048576));[ci_fft,cv_fft] max(cm_fft(1:1048576/2));cv_pay cv_fft / 1048576/SYM;cv_pay 的意思就是每个样点要补偿多少点是个很小的数。太大了就不对了。
4.3 计算结果
本例子里是 3.9700e-051个符号产生 3.9700e-05 * N 0.1626点偏移在K2时每6个chirp就会冒一个点每12个chip就会发生1次FFT峰值的移动造成错误。用LoRaPHY的方法计算出的是0.127二者还是存在相当的差异。
5 基于 up-down chirp 边缘峰值交汇的时频差精确估计算法
由于存在第三部分的原因无论是单独采用up chirp还是 down chirp都不能准确测定频差。但是Lora存在一个up-down chirp的符号反转构成了所谓的双线性调频(dual-chirp)就用这个波形进行匹配滤波。
5.1 笨办法观察规律
最笨的办法就是两层 for 循环把原始波形进行时间搬移、频率搬移而后求取各个搬移参数下的相关峰。
下图显示的是搜索的窗口。根据特定的频差首先把载波乘以 e 2 π j Δ n e^{2\pi j \Delta n} e2πjΔn,而后再平移 δ \delta δ样点计算共轭相乘的和。 在窗口内遍历所有的 Δ , δ \Delta,\delta Δ,δ, D ( m , n ) e 2 π j ( − 1 2 Δ m δ 2 ⋅ 2 S F ⋅ K ) m δ K ⋅ e − 2 π j ( − 1 2 n 2 ⋅ 2 S F ⋅ K ) n K D(m,n)e^{2\pi j (-\frac{1}{2}\Delta \frac{m\delta}{2\cdot2^{SF}\cdot K})\frac{m\delta}{K}}\cdot e^{-2\pi j (-\frac{1}{2} \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K}} D(m,n)e2πj(−21Δ2⋅2SF⋅Kmδ)Kmδ⋅e−2πj(−212⋅2SF⋅Kn)Kn
while(pos end_updownpos)x_ti pos - start_updownpos 1;for (x_fi 1:searchf)sumsig raw_iq(pos chirp_n2).*exp(1j.*2*pi*(x_fi-searchf/2)/2./(2^SF)./K./K.*chirp_n2) .* chirp_updown(1chirp_n2);sumv abs(sum(sumsig));updown_tf(x_ti,x_fi) sumv;endforpos pos 1;end %while结果如下
这个计算极其耗时实用性不高。不过这个图却很好的说明了up-down chirp做时频二维相关的特点。
中间的尖峰时频都对准了up,down同时相关上了产生双倍的扩频增益。两条相交的山脊由于时频可以互补发生了第三章阐述的平移部分相关上1条up或者down另一条没有增益。山脊是线性的、相交于尖峰。
5.2 基于边缘最大值直线交点快速交汇的时频差估计
上面的这个算法计算实在太慢了。注意到其实upchipdownchip各自的最大值在上述搜索窗口的边缘能够准确捕获只要在边缘4个边上进行计算而后通过连线交点估计位置即可。 if (pos end_updownpos)x_ti pos - start_updownpos 1;for (x_fi 1:searchf)if (posstart_updownpos pos1end_updownpos)if (x_fi1 x_fisearchf)continue;endifendifsumsig raw_iq(pos chirp_n2).*exp(1j.*2*pi*(x_fi-searchf/2)/2./(2^SF)./K./K.*chirp_n2) .* chirp_updown(1chirp_n2);sumv abs(sum(sumsig));updown_tf(x_ti,x_fi) sumv;endforpos pos 1;elsex_f zeros(1,4);x_t zeros(1,4);max_ct 0;raw_tf updown_tf;%找到4个峰值while (max_ct 4)[max_xtv,max_xti] max(updown_tf);[max_xfv,max_xfi] max(max(updown_tf));curr_f max_xfi;curr_t max_xti(curr_f);updown_tf(curr_t,curr_f) 0;check_pk 1;fake_peak 0;while (check_pk max_ct fake_peak0)distance sqrt((x_f(check_pk) - curr_f)^2 (x_t(check_pk) - curr_t)^2);if (distance 16 * K)fake_peak 1;endifcheck_pk check_pk 1;endwhileif (fake_peak1)continue;endifmax_ct max_ct 1;x_f(max_ct) curr_f;x_t(max_ct) curr_t;endwhile%求直线交点xrge size(updown_tf);[x_it1,y_it1] intersection_point(x_t(1),x_f(1),x_t(2),x_f(2),x_t(3),x_f(3),x_t(4),x_f(4));[x_it2,y_it2] intersection_point(x_t(1),x_f(1),x_t(3),x_f(3),x_t(2),x_f(2),x_t(4),x_f(4));[x_it3,y_it3] intersection_point(x_t(1),x_f(1),x_t(4),x_f(4),x_t(2),x_f(2),x_t(3),x_f(3));if (x_it1 0 x_it1 xrge(1) y_it1 0 y_it1 xrge(2))x_tc floor(x_it10.5);x_fc floor(y_it10.5);x_tf x_it1;x_ff y_it1;elseif (x_it2 0 x_it2 xrge(1) y_it2 0 y_it2 xrge(2))x_tc floor(x_it20.5);x_fc floor(y_it20.5);x_tf x_it2;x_ff y_it2;elseif (x_it3 0 x_it3 xrge(1) y_it3 0 y_it3 xrge(2))x_tc floor(x_it30.5);x_fc floor(y_it30.5);x_tf x_it3;x_ff y_it3;elsex_tc xrge(1)/2;x_fc xrge(2)/2;x_tf x_tc;x_ff x_fc;endifprintf(Finished.\n);if (x_ff - searchf/2 0)payback_clk cv_pay;elsepayback_clk -cv_pay;endif%figure(2);pos start_updownpos x_tc - 1;syn_pos start_updownpos x_tc - 1;updown_tf raw_tf;sumsig raw_iq(pos chirp_n2).*exp(1j.*2*pi*(x_ff-searchf/2)/2./(2^SF)./K./K.*chirp_n2) .* chirp_updown(1chirp_n2);sumv abs(sum(sumsig));updown_tf(x_tc,x_fc) sumv;mesh(updown_tf);endif结果
需要注意的是第四步里估计的钟差是一个绝对值要根据本方法提供的频差的正负号进行调整。当然如果使用的设备时钟很精确则可以直接使用 LoRaPHY使用的方法。
6 卡住位置补偿解调
一旦找到了精确的时间、频率则可以对剩余部分进行解调。解调的方法就是先乘以 down-chip做fft、折半到 B带宽后读取峰值的位置。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3. 解调if (status3)draw_mark(pos) -30;%搬移频率到0中心data raw_iq(pos:posSYM-1).*exp(1j.*2*pi*(deltaf-searchf/2)/2./(2^SF)./K./K.*chirp_n).*chirp_dw_dem;fft_data_raw abs(fftshift(fft(data,SYM)));code_len SYM/K;if (K1)wrp_N SYM;wrp_P wrp_N / 2;wrp_M 2^(SF-1);fft_data_raw(wrp_P-wrp_M1:wrp_P) fft_data_raw(wrp_P-wrp_M1:wrp_P) fft_data_raw(wrp_Pwrp_M1:wrp_P2*wrp_M);fft_data_raw(wrp_P1:wrp_Pwrp_M) fft_data_raw(wrp_P1:wrp_Pwrp_M) fft_data_raw(wrp_P-2*wrp_M1:wrp_P-wrp_M);fft_data1 fft_data_raw(wrp_P-wrp_M1:wrp_Pwrp_M);elsefft_data1 fft_data_raw;endif[codev,code_d] max(fft_data1);fft_data fft_data1;draw_hdseq(pos:poscode_len-1)fft_data/K;figure(1);clf;hold on;plot(draw_mark,k);plot(draw_hdseq*30/SYM,b);plot(draw_freqs,g);hold off;pause(0.001);raw_demsym(1,codi) mod(code_d2^(SF-1),2^SF);codi codi 1;%样点偏移补偿(timeing offset),这里存疑因为 payback_clk 是靠频偏估计的clk_pay clk_pay SYM;PZ payback_clk * codi * SYM;pos round(pos_start clk_pay PZ);endif输出:
lora_decode
Hit up20406 24502
Hit up24502 28599
Hit up28599 32695
Hit up32695 36791
Hit up36791 40887
Hit up40887 44983
Hit up44983 49079
Hit up49079 531760 1.1271e024.0960e03 7.8453e018.1920e03 -1.6873e011.2288e04 -1.6863e021.6385e04 -4.3814e012.0481e04 6.2779e012.4577e04 1.0725e022.8673e04 8.4182e013.2769e04 4.2651e003.6865e04 -1.5306e024.0962e04 -6.7274e01
xcorr pre-upchirp amp fft payback3.97002e-05
Hit down65187 69283
Fine Time Freq begin...Finished.Columns 1 through 23:670 826 590 146 1818 326 1526 1170 966 950 858 2033 657 966 1101 1854 650 1563 1361 1820 1574 331 1120Columns 24 through 46:1761 133 1753 922 1027 1564 622 984 1461 1831 1472 1690 234 90 110 1004 325 1917 1740 908 142 261 1955Columns 47 through 69:140 1916 835 1426 1628 216 1417 1575 1021 1650 217 35 218 559 689 1294 1114 1590 304 1753 148 1437 1625Columns 70 through 75:543 369 1994 609 609 1015解调后的符号通过 LoRaPHY提供的解析方法可以还原出实际的数据: lora_decode
Syn code start.start offset 15, fix 0
01ECAE9230011111111111111112222222222222222333333333333333344448F16
078357676000123456789ABCDEF0123456789ABCDEF0123456789ABCDEF0123CFF4值得注意的是上述解调后从2.25符号往后其实有15个符号30样点的偏差。这个偏差是固定的可以直接根据头部校验搜索并找到。这个偏差的产生机制并不是很清楚。
7 拉远距离加白噪声测试
我们添加一个比波形本身强度高10倍的噪声相位已经没法看了但通过相关发现依旧能正确处理 在10倍的噪声下首部的相位也有抖动但大的周期还在。偶尔估计会出现偏差但依旧能够较大概率解析出正确的数据会出错。
8 代码获取和运行
参考代码库a4_lora/octave
从本次摸索我发现读别人的文章与自己摸索实现一下认识程度是显著不同的。我想LoRa本身可能也采取了类似的时频差补偿算法否则它那60多块钱的成本以及巨大的参数公差根本就收不好。
自己的Octave实现与网上的例子相比既不要求安装matlab R2019以上这种20GB的工程数学软件也不和GNU-Radio这种多层封装的实现发生关系。
要运行代码请直接安装开源的计算工具 GNU-Octave 和Matlab很像但没有版权问题安装包在几百兆左右。Octave也可以直接使用MSYS2的pacman -S来安装。本次的代码不依赖任何工具箱包括communicaitons。 文章转载自： http://www.morning.pnmgr.cn.gov.cn.pnmgr.cn http://www.morning.xnqwk.cn.gov.cn.xnqwk.cn http://www.morning.qjlnh.cn.gov.cn.qjlnh.cn http://www.morning.wbdm.cn.gov.cn.wbdm.cn http://www.morning.rlhh.cn.gov.cn.rlhh.cn http://www.morning.pmmrb.cn.gov.cn.pmmrb.cn http://www.morning.lpsjs.com.gov.cn.lpsjs.com http://www.morning.mgmyt.cn.gov.cn.mgmyt.cn http://www.morning.qmsbr.cn.gov.cn.qmsbr.cn http://www.morning.pwwjs.cn.gov.cn.pwwjs.cn http://www.morning.tjkth.cn.gov.cn.tjkth.cn http://www.morning.aswev.com.gov.cn.aswev.com http://www.morning.xhgcr.cn.gov.cn.xhgcr.cn http://www.morning.ydxx123.cn.gov.cn.ydxx123.cn http://www.morning.trkhx.cn.gov.cn.trkhx.cn http://www.morning.srxhd.cn.gov.cn.srxhd.cn http://www.morning.bpmfz.cn.gov.cn.bpmfz.cn http://www.morning.ldhbs.cn.gov.cn.ldhbs.cn http://www.morning.knrgb.cn.gov.cn.knrgb.cn http://www.morning.fxygn.cn.gov.cn.fxygn.cn http://www.morning.smygl.cn.gov.cn.smygl.cn http://www.morning.rwqk.cn.gov.cn.rwqk.cn http://www.morning.kllzy.com.gov.cn.kllzy.com http://www.morning.yyzgl.cn.gov.cn.yyzgl.cn http://www.morning.hgbzc.cn.gov.cn.hgbzc.cn http://www.morning.nrgdc.cn.gov.cn.nrgdc.cn http://www.morning.pigcamp.com.gov.cn.pigcamp.com http://www.morning.bnjnp.cn.gov.cn.bnjnp.cn http://www.morning.jbxfm.cn.gov.cn.jbxfm.cn http://www.morning.bgpch.cn.gov.cn.bgpch.cn http://www.morning.fdrch.cn.gov.cn.fdrch.cn http://www.morning.rbkml.cn.gov.cn.rbkml.cn http://www.morning.zdmlt.cn.gov.cn.zdmlt.cn http://www.morning.jspnx.cn.gov.cn.jspnx.cn http://www.morning.khtjn.cn.gov.cn.khtjn.cn http://www.morning.smxyw.cn.gov.cn.smxyw.cn http://www.morning.zmpqt.cn.gov.cn.zmpqt.cn http://www.morning.twdkt.cn.gov.cn.twdkt.cn http://www.morning.ngjpt.cn.gov.cn.ngjpt.cn http://www.morning.spfh.cn.gov.cn.spfh.cn http://www.morning.rycbz.cn.gov.cn.rycbz.cn http://www.morning.fbnsx.cn.gov.cn.fbnsx.cn http://www.morning.zfkxj.cn.gov.cn.zfkxj.cn http://www.morning.kqkmx.cn.gov.cn.kqkmx.cn http://www.morning.knnc.cn.gov.cn.knnc.cn http://www.morning.bpzw.cn.gov.cn.bpzw.cn http://www.morning.khclr.cn.gov.cn.khclr.cn http://www.morning.cwjxg.cn.gov.cn.cwjxg.cn http://www.morning.crxdn.cn.gov.cn.crxdn.cn http://www.morning.nwljj.cn.gov.cn.nwljj.cn http://www.morning.clzly.cn.gov.cn.clzly.cn http://www.morning.drbwh.cn.gov.cn.drbwh.cn http://www.morning.lpcpb.cn.gov.cn.lpcpb.cn http://www.morning.lwmzp.cn.gov.cn.lwmzp.cn http://www.morning.hympq.cn.gov.cn.hympq.cn http://www.morning.tgydf.cn.gov.cn.tgydf.cn http://www.morning.hqwcd.cn.gov.cn.hqwcd.cn http://www.morning.bgygx.cn.gov.cn.bgygx.cn http://www.morning.ntcmrn.cn.gov.cn.ntcmrn.cn http://www.morning.tgyzk.cn.gov.cn.tgyzk.cn http://www.morning.fnczn.cn.gov.cn.fnczn.cn http://www.morning.rzczl.cn.gov.cn.rzczl.cn http://www.morning.yjfmj.cn.gov.cn.yjfmj.cn http://www.morning.mqzcn.cn.gov.cn.mqzcn.cn http://www.morning.crdtx.cn.gov.cn.crdtx.cn http://www.morning.rcjyc.cn.gov.cn.rcjyc.cn http://www.morning.bpwdc.cn.gov.cn.bpwdc.cn http://www.morning.zhishizf.cn.gov.cn.zhishizf.cn http://www.morning.hytqt.cn.gov.cn.hytqt.cn http://www.morning.qggm.cn.gov.cn.qggm.cn http://www.morning.zdgp.cn.gov.cn.zdgp.cn http://www.morning.rnngz.cn.gov.cn.rnngz.cn http://www.morning.fbfnk.cn.gov.cn.fbfnk.cn http://www.morning.grbp.cn.gov.cn.grbp.cn http://www.morning.tjcgl.cn.gov.cn.tjcgl.cn http://www.morning.gl-group.cn.gov.cn.gl-group.cn http://www.morning.jlxqx.cn.gov.cn.jlxqx.cn http://www.morning.mnjwj.cn.gov.cn.mnjwj.cn http://www.morning.ypnxq.cn.gov.cn.ypnxq.cn http://www.morning.cfcpb.cn.gov.cn.cfcpb.cn
