在读一篇《Drone Security and the Mysterious Case of DJI’s DroneID》的文章时，作者通过了逆向firmware和软件无线电SDR捕捉大*若干款无人机的通信信号来找到了如何破解和逆向出无人机的通信ID和飞行数据等信息。 所以本人就在想如果能既然能逆向提取出这些信息，那么是不是可以采用正向工程，将这些信息通过SDR等设备再发出去，用SDR的发射信号来伪造一个飞行状态中的无人机呢。
要达到这个目的，需要对于原作者的方案进行逆向思考才行。下面逐步展开我这边的思路和实施方式。
研究基础：
原作者在逆向的时候采用的思路这样的：
如果手里已经有了SDR采集到的无人机发射的信号文件一份，那么信号的基本处理方式如下：
image.png (24.57 KB, 下载次数: 1)
下载附件
2024-11-25 09:42 上传
到这里已经是可以解码出经过信道编码的bit流了，再通过如下的Turbo解码操作，就可以恢复最原始的无人机ID等信息payload。 基本流程图如下:
image.png (16.49 KB, 下载次数: 1)
下载附件
2024-11-25 10:24 上传
对于Turbo解码，是仿照LTE的Turbo解码算法做的，也就是说大*此系列无人机的不少编码特性是照搬了3GPP LTE规范的，而3GPP LTE规范是公开的，理论上来说就可以实现任意信息的伪造了。
但是由于Turbo编码这一块儿目前我还没有完全参透，故而完成一款完全从最顶层到最底层的设计还有些难度。 故而本次实验仅采用一种特定的伪造，就是分身的方式进行伪造的研究。基本的思路就是将蓝色的流程图部分倒过来。思路如下：
image.png (22.36 KB, 下载次数: 1)
下载附件
2024-11-25 10:32 上传
具体实践要点:
挑战1： 如何得到原始的Turbo编码码流
根据外文文献提供的线索需要对于，一个基本的有效信息帧由如下几个部分构成, 10个slot,10ms，和LTE的帧结构类似，其中9个slot有信息，最后一个是保护间隔slot，其中第4和第6个slot是ZC序列的同步slot。
image.png (224.83 KB, 下载次数: 0)
下载附件
2024-11-25 10:35 上传
那么在提取的时候就要通过第4和6slot的已知信息来进行各种同步和频率相位校正。 我们最终的目标就是经过一系列的信号处理，能最终得到QPSK 的调制信息，从而得到Demodulate Bitstream。
image.png (125.43 KB, 下载次数: 0)
下载附件
2024-11-25 10:38 上传
具体实践如下，如下是一份某款无人机采样数据信息:
采样文件，时域波形, 采样率50Msps
image.png (19.39 KB, 下载次数: 0)
下载附件
2024-11-25 10:40 上传
功率谱如下,红圈中圈起的是疑似有用信息帧的位置
image.png (103.93 KB, 下载次数: 0)
下载附件
2024-11-25 10:41 上传
为了分析逐个帧的，需要对这三个帧进行逐个分析，那么首先就要进行粗略切割:对信号进行短时傅里叶变换，以得到不同时刻的信号频谱变动情况，执行如下函数:signal.stft(input_sig, sample_rate, nfft=64, nperseg=64) 结果写入 2.after_stft_time_vs_freq.xlsx画图分析：可以看到不同频率分量在不同时间点上的功率分布总的趋势相同，但是局部还是不同的。
image.png (133.84 KB, 下载次数: 0)
下载附件
2024-11-25 10:41 上传
下面求取不同时间点上最大能量所在的频率分量和时间的关系，蓝色为频率，橘色为对应最大能量值
image.png (80.28 KB, 下载次数: 0)
下载附件
2024-11-25 10:42 上传
看上去毛刺很多，下面我们设定一个特定的Noise Level ， 把能量低于noiseLevel的部分都设值为零。生成数据4.engery_above_level.xlsx过滤后的图和数据如下：
image.png (68.58 KB, 下载次数: 0)
下载附件
2024-11-25 10:42 上传
根据数据帧的时间长度转换为对应的采样点数，得到一个范围，需要在这个采样值范围内的能量都要大于noiseLevel。采用 signal.find_peaks，来寻找满足这个条件的一系列的点:
peaks, properties = signal.find_peaks(above_level, width=[signal_length_min_samples, signal_length_max_samples],wlen=100*signal_length_max_samples)
image.png (21.28 KB, 下载次数: 0)
下载附件
2024-11-25 10:43 上传
然后采用left_bases 和 right_bases 前后各加入一个start_offset和 end_offset进行粗略切割。切割后的三段数据分别为:
6.packets_data_after_peak_slicing_0.csv6. packets_data_after_peak_slicing_1.csv6.packets_data_after_peak_slicing_2.csv对应的频谱分别如下:
image.png (92.48 KB, 下载次数: 0)
下载附件
2024-11-25 10:43 上传
可以看出3段切出的频谱左右都留出了足够余量，没有损失信息。这为后续继续其他分析奠定了基础。下面来逐个对这三段信号进行精细化处理，在此我只对其中最后一组包含了droneID数据的信号进行细致分析。下面执行 signal.welch 来对信号进行功率谱密度PSD分析,用于进行有用信号开始位置的粗略估计。
nfft_welch = 2048signal.welch( y, Fs, nfft=nfft_welch, return_onesided=False)再执行fftshift将负频率移动到左边，让f=0在中间，以便于分析
Pxx_den = np.fft.fftshift(Pxx_den) f = np.fft.fftshift(f) 生成 7.packets_data_after_peak_slicing_2_1_after_welch.xlsx
image.png (55.83 KB, 下载次数: 0)
下载附件
2024-11-25 10:44 上传
添加一个假的中间直流DC载波（从2048//2-10 = 1014 到 1034, 共20个采样点) 设置为1.1倍的功率均值，模拟DC上载波泄露干扰能量），
# add a fake DC carrier Pxx_den[nfft_welch//2-2:nfft_welch//2+2] = 1.1*Pxx_den.mean() #采用添加4个采样点也是OK的如下所示: 7.packets_data_after_peak_slicing_2_3_after_add_fakeDC
image.png (87.54 KB, 下载次数: 0)
下载附件
2024-11-25 10:45 上传
在这里假设有用信号的带宽是10MHz(从上图可以看出，信号宽度大概是10MHz） , fft 是2048， SCS =15kHz , 有用带宽2048*15kHz = 30.72Mhz, 半带宽是15.36MHz，这是20MHz带宽LTE的典型采样参数。 本文件是50MHz采样频率，故而50/15.36 约等于4个采样点。，将来重采样15.36Ms之后，就可以大概等于1个点了(一个子载波)。 找出超出平均值并且宽度足够宽的部分，也就是上图中的有用信号部分，得到她的偏移量candIDAte_bands = consecutive(np.where(Pxx_den > Pxx_den.mean())[0]) 得到数据:6. packets_data_after_peak_slicing_2_5_offset_caluclation_band_0
image.png (10.17 KB, 下载次数: 0)
下载附件
2024-11-25 10:45 上传
得到信号的带宽是9.2Mhz左右，符合预期。并且得到真正的偏移量如上图，后面要针对这个偏移量，对信息帧信号两边多余的部分进行切割。 # correct frequency offset packet_data = fshift(packet_data, -1.0*offset, sample_rate) 得到进一步精确的信号采样
8.packets_data_offset_shifting_pktnum_0.csv，他的频谱如下:
image.png (58.05 KB, 下载次数: 0)
下载附件
2024-11-25 10:46 上传
可以看到经过位移校正之后，此时有用信号的位于了中间部位了。我们终于实现了时间和频率的粗略同步，下面就可以在这个基础上对信号进行重采样了。从50Msps降低到15.36Msps，这样的话才能将有用数据和子载波一一对应。9.packets_pktnum_0resample_50.0_to_15.36.csv
image.png (25.5 KB, 下载次数: 0)
下载附件
2024-11-25 10:46 上传
下面的任务就是检测每个symbol的符号级别的定时同步，通过检测CP来得到时间上粗同步。首先进行幅度归一化
raw_bk= packet_data raw_bk/= np.max(np.abs(raw_bk))
得到10.packets_data_normalized_amplitudes_0.csv再次查看信号的频谱:
image.png (27.4 KB, 下载次数: 0)
下载附件
2024-11-25 10:47 上传
采用cpl长度80的符号头的信息和预估的符号尾部的部分进行相关度运算，获取相关度数值，相关度最高时，就可以认为是OK了。
for n in range(NFFT, len(raw_bk) - cpl):
ac = np.sum(raw_bk[n:n+cpl] * np.conj(raw_bk[n-NFFT:n-NFFT+cpl]))
res.append(ac)
res_abs = np.abs(res)得到 11.packets_data_cp_corr_pktnum0_CPLength_80,可以看到我们获取到了几组峰值数据:
image.png (79.47 KB, 下载次数: 0)
下载附件
2024-11-25 10:47 上传
采用signal.find_peaks 来找到满足要求的峰值，并且丢弃幅度小的峰值
peaks, _ = signal.find_peaks(res_abs, distance = 1000) peak_prominences, _, _ = signal.peak_prominences(res_abs, peaks) 得到数据
12.packets_cp_peak_pktnum0_CPLength_80.csv
image.png (11.15 KB, 下载次数: 0)
下载附件
2024-11-25 10:48 上传
然后得到第一个峰值的索引，也就是703，作为我们的第一符号开始的索引。 start = peaks[peak_index] detected_ffo = Fs / (2 * np.pi * NFFT) * np.angle(res[start]) 并且以此来计算大致的频率校正值进行初步的符号切割。去掉前面的703个采样点之后的信号13. packets_cp_1_remove_703_samples.csv
信号
image.png (42.56 KB, 下载次数: 0)
下载附件
2024-11-25 10:48 上传
从图中可以看出，信号的时域起始点就开始了有用信号，达到我们预期。 频谱如下:
image.png (25.32 KB, 下载次数: 0)
下载附件
2024-11-25 10:48 上传
频域也是如此 接下来频域校正 -6458.548966919529 得到13.packets_cp_2_correct_ffo_-6458.548966919529.csv， 信号详情如下:
image.png (38.81 KB, 下载次数: 0)
下载附件
2024-11-25 10:49 上传
image.png (25.73 KB, 下载次数: 0)
下载附件
2024-11-25 10:49 上传
然后对于各个符号进行切割，此时符号带有CP
image.png (33.88 KB, 下载次数: 0)
下载附件
2024-11-25 10:49 上传
我们将symbol 0 的IQ数据代入信号分析系统，查看里面的详情如下:如下是symbol 0的时域波形:
image.png (48.35 KB, 下载次数: 0)
下载附件
2024-11-25 10:50 上传
下面就是去掉各个symbol的CP，得到每个symbol的1024个采样点，得到的数据文件是:
image.png (29.59 KB, 下载次数: 0)
下载附件
2024-11-25 10:50 上传
对于无CP的符号进行FFT变换后，得到如下数据文件。
image.png (29.91 KB, 下载次数: 0)
下载附件
2024-11-25 10:51 上传
我们可以看到此时绘制星座图如下:symbol0的fft星座图如下:
image.png (70.24 KB, 下载次数: 0)
下载附件
2024-11-25 10:51 上传
从图中可以看出，这个星座图有些像ZC序列的星座图，但不是我们预期的QPSK，应该还是需要更细致的频率校正和信道估计才行。 下来就需要借助已知的ZC序列作为参考信号进行信道估计和频率校正了。按照这个无人机机型的协议ZC序列在symbol 3和 symbol 5. 对于ZC序列进行盲检，确保这两个symbol上使用的zc序列的根分别是如下两个， 其中第二个是用于sync的，故而要确保147序列存在。
zc_seq_1 = 600 zc_seq_2 = 147 以寻找u=147的序列为例，我们分别和147相邻的几个ZC序列进行相关运算。以u=146的根序列为例，我们可以看到zc序列的星座图如下:
image.png (167.65 KB, 下载次数: 0)
下载附件
2024-11-25 10:52 上传
image.png (73.87 KB, 下载次数: 0)
下载附件
2024-11-25 10:53 上传
我们分别打开文件看看相关度：
image.png (47.46 KB, 下载次数: 0)
下载附件
2024-11-25 10:53 上传
image.png (32.68 KB, 下载次数: 0)
下载附件
2024-11-25 10:53 上传
可以看出147的ZC序列相关度要明显大，符合预期 采用这个序列来进行信道评估，基本原理是采用u=147的预期序列和收到的序列系列的信号进行除法 channel = estimate_channel(symbols_freq_domain, ZC_SYMBOL_IDX[1], zc_seq_2, debug, prefix ='17.estimate_channel_zc1') 得到文件 17.estimate_channel_zc11_estimate_channel_estimate_zc_u_147symbol_index_5.csv 信道的幅频响应如下：
image.png (49.57 KB, 下载次数: 0)
下载附件
2024-11-25 10:53 上传
下面采用u=600的zc序列进行进一步频率校正 sampling_offset, angle = find_zc_offset_angle(raw_samples_orig,start,resample_rate, detected_ffo, ZC_SYMBOL_IDX[0], 600, zc_cyc, debug) 开始进行逐步的样本内插, 然后求相位差的均方差MSE，最小的值所在位置就是进一步要校正的量for i in np.linspace(-6, 0, 30): adiff = np.angle(a / zc_sym_f) # remove DC carrier adiff[NCARRIERS//2] = adiff[NCARRIERS//2+1] adiff = np.unwrap(adiff) slope_arr = (adiff - np.mean(adiff)) #slope = np.max(adiff) - np.min(adiff) slope = np.sqrt(np.mean(slope_arr**2)) 得到一系列的数据文件，可以看出来，不同偏移下的MSE是不同的，图中如下红圈的文件就是MSE最小的文件。
image.png (252.76 KB, 下载次数: 0)
下载附件
2024-11-25 10:54 上传
打开它之后查看
image.png (73.43 KB, 下载次数: 0)
下载附件
2024-11-25 10:54 上传
故而频率校正 -5.17个采样点是OK的19. find_zc_offset_angle_symbol_3_1_find_zc_offset_resx_y_u600_NCARRIERS_601.csv对应的角度校正是0.053 弧度
image.png (11.95 KB, 下载次数: 0)
下载附件
2024-11-25 10:55 上传
这样就得到了综合性的各种偏移，频偏，相位校正信息，下面对于重采样后的数据进行校正symbols_time_domain, symbols_freq_domain = raw_data_to_symbols(resample_rate, raw_samples_orig, start, ffo=detected_ffo, sampling_offset=sampling_offset, angle=angle)就可以得到一系列的信息。对于频域的采样数据再次进行信道修正 for i, symbol_f in enumerate(symbols_freq_domain): symbol_eq = symbol_equalized(symbol_f, channel)
image.png (64.7 KB, 下载次数: 0)
下载附件
2024-11-25 10:56 上传
然后再次进行符号提取和fft变换，我们看到，熟悉的QPSK星座图来了。
image.png (26.88 KB, 下载次数: 0)
下载附件
2024-11-25 10:59 上传
进行QPSK解调，就得到了如下原始的turbo码流文件。是非常长的一长串0,1,2,3的数字，分别代表了不同的两位bit流。
image.png (2.83 KB, 下载次数: 0)
下载附件
2024-11-25 11:02 上传
挑战2： 将上述Turbo码流进行重新打包为可以发射的基带信号。 这样的好处是你可以得到一个完全不含干扰的纯净信号。这部分相对简单一些，故而只写一些主要步骤。关键代码方法如下： 将上述数字转换为复数IQ信号。
image.png (30.85 KB, 下载次数: 0)
下载附件
2024-11-25 11:04 上传
将上述调制的QPSK信号分成7份，分别映射到slot 1,2,3,5,7,8上。 再分别生成ZC两个序列，映射到slot 4和slot 6上作为同步使用。
image.png (40.89 KB, 下载次数: 0)
下载附件
2024-11-25 11:07 上传
image.png (44.14 KB, 下载次数: 0)
下载附件
2024-11-25 11:09 上传
添加CP操作
image.png (23.92 KB, 下载次数: 0)
下载附件
2024-11-25 11:11 上传
设定一定的信噪比，将上面所有的符号级联到一起
image.png (48.99 KB, 下载次数: 0)
下载附件
2024-11-25 11:11 上传
产生的最终信号写入文件：
image.png (36.17 KB, 下载次数: 0)
下载附件
2024-11-25 11:12 上传
将我们生成的伪造信号放到软件里面分析频谱和时域波形，是不是很酷，而且完全可以以假乱真了。
image.png (55.54 KB, 下载次数: 0)
下载附件
2024-11-25 11:15 上传
image.png (62.43 KB, 下载次数: 0)
下载附件
2024-11-25 11:14 上传
声明，本文提供的主要是研究思路和自己的实践过程，不代表通用或者最优的办法，也不太可能适用于所有型号的无人机。为了确保不用于商用或者不当用途，不提供本实验的源代码，仅仅提供了关键的思路和过程。有能力的同行可以自行探索更优的方法。祝大家有所收获。
