于生宝, 房 钰, 高丽辉, 刁 庶
(吉林大学 仪器科学与电气工程学院, 长春 130061)
电磁法是矿产资源勘探、 环境与工程地质学、 区域构造研究的一种有效方法, 可在地面、 空中及海上进行资源探测[1]. 根据发射场的性质, 电磁法可分为时间域电磁法和频率域电磁法[2]. 相比于时间域电磁法, 频率域电磁法在整个采集过程中电流连续发射, 通过获取总场反应地下电阻率的信息, 信号较强. 因此, 频率域电磁法适用于观测面积大、 勘探深度较深的复杂地形区域探测[3]. 目前, 发射电流的波形多为固定频比或频差, 如2n伪随机波和方波等, 这些波形的能量均按特定比例等距离分布在各次谐波上[4]. 当针对特定深度的目标进行高精度电磁勘探时, 由于传统发射电流波形频率间隔固定, 需发射不同基频的波形才能满足探测要求, 而多次发射又会降低勘探效率, 因此存在勘探精度与效率相互制约的问题.
针对上述问题, 本文提出一种基于改进粒子群优化算法的半周期镜像对称选择性谐波消除脉宽调制(selective harmonic eliminated pulse width modulation, SHEPWM)发射电流控制方法. 在频率域电磁测深过程中, 先根据探测深度计算出发射电流期望的频谱, 再根据发射电流的时频域特性, 建立控制时域信息(开关时刻、 电流值)和频域信息(直流分量和谐波幅值、 相位)的半周期镜像对称SHEPWM非线性方程组, 最后采用改进的粒子群优化算法对SHEPWM非线性方程组进行求解, 获得对应的时域开关序列, 控制发射系统逆变器, 从而得到理想的发射电流, 提高电磁测深的纵向分辨率.
SHEPWM方法是一种逆变器常用的计算调制法, 广泛应用于电力、 电子等领域. 该方法先通过对逆变器输出电压波形进行Fourier分解, 建立与开关时刻和波形谐波分量有关的非线性方程组, 求解方程组得到满足要求的开关角, 然后利用这些开关角控制逆变器中开关器件的导通和关断, 从而消除直流分量和某些特定的低次谐波, 得到所需的基波分量[5]. 当应用在频率域电磁测深领域时, 波形的谐波分量是有用信号, 此时半周期镜像对称SHEPWM方法可根据勘探目标地质结构设置发射电流奇次谐波的幅值、 相位、 分布及数量, 控制探测所需发射电流的频域信息. 下面给出针对特定深度勘探目标的半周期镜像对称SHEPWM波形的设计方法.
SHEPWM周期性输出电压波形展开成一个收敛的Fourier级数, 表示为直流分量、 基波分量和谐波分量的叠加:
(1)
其中:i表示谐波次数;A0表示直流分量;Ai表示第i次谐波的振幅;θi表示第i次谐波的相位;ω=2πf表示角频率,f表示输出电压波形频率.a0,ai,bi,A0,Ai,θi满足如下关系:
(2)
式(2)中i=1,2,…, 系数a0,ai,bi根据Fourier分析可得:
(3)
其中:i=1,2,…;T=2π为输出电压波形的周期. 图1为一个半周期镜像对称的SHEPWM发射波形, 其中:Ud为幅度;α1,α2,…,αN为半个周期内的开关时刻[6].
图1 半周期镜像对称SHEPWM波形Fig.1 Half-cycle mirror symmetric SHEPWM waveform
根据式(3)可得半周期镜像对称SHEPWM波形的a0,ai,bi:
(4)
联立式(2)和式(4), 可得半周期镜像对称SHEPWM波形的非线性方程组为
(5)
0<α1<α2<…<αN<π.
(6)
式(5)左边是与开关角相关的多项式, 右边是与发射电流被控谐波频域特性幅值Ai和相位θi相关的多项式. 因此, 可通过求解式(5), 精确计算出符合发射要求的开关角度. 式(6)为式(5)的限定条件. 此时, 发射电流奇次谐波幅值、 相位、 分布及数量可控, 偶次谐波幅值和直流分量为0.
根据频率域电磁探测原理, 探测过程中可不使用相位, 故可设θi=0, 根据实际探测地形设置谐波幅值Ai, 如Ai的最大值不应大于4Ud/π, 且所有谐波幅值应满足能量守恒定律[7]:
(7)
SHEPWM控制方法的关键是对开关器件开关角度的求取, 即非线性方程组的求解. 目前, 求解非线性方程组的有效方法主要分为两类: 一类是经典的数值迭代法, 如牛顿迭代法、 拟牛顿法、 梯度法等, 这些算法均具有较快的局部收敛性, 但均需设置合适的初始值, 且随着未知数的增加, 算法复杂性增加; 第二类是仿生优化算法, 如粒子群优化算法[8]、 遗传算法、 模拟退火算法、 神经网络算法等, 这些方法多采用分析方程形成目标函数, 将其转换为最优化问题. 其中的粒子群优化算法由于其算法简单、 容易实现被广泛应用, 但该方法也存在易陷入局部最优的缺点[9]. 为提高算法的性能, 文献[10]提出了线性改进惯性权重和学习因子, 并结合了模拟退火算法的方法; 文献[11]提出了用余弦函数对惯性权重进行非线性改进的方法, 提高了收敛速度和寻优能力. 本文采用动态自适应调整的惯性权重, 并对学习因子进行三角函数非线性改进的方法, 以解决粒子群优化算法求解更多未知数时适应度函数值不能收敛到零的问题.
标准粒子群优化算法中, 每个粒子的状态都由位置、 速度和目标函数的适应度值构成, 根据自身经验和种群经验在D维有限空间中搜索最优值. 每个粒子的速度[12]和位置更新公式为
标准粒子群优化算法的惯性权重和学习因子都是常数, 无法根据搜索的情况进行动态调整, 寻优能力有限. 为了动态调节算法的全局搜索和局部搜索能力, 本文引入文献[13]定义的平均聚焦距离、 最大聚焦距离和聚焦距离变化率:
其中:D为粒子的维数;M为粒子的个数. 通过对平均聚焦距离和最大聚焦距离的更新, 得到每次迭代时的聚焦距离变化率, 由此判断本次粒子应提高全局搜索能力还是局部搜索能力, 从而动态自适应地调整惯性权重. 本文采用文献[13]定义的惯性权重公式:
(13)
图2 改进粒子群优化算法流程Fig.2 Flow chart of improved particle swarm optimization algorithm
其中:z1=0.3,z2=0.2;r为一个[0,1]间均匀分布的随机数. 这种动态自适应调整惯性权重可灵活地调节全局搜索与局部搜索的能力. 当聚焦距离变化率较大时, 粒子的全局搜索能力较差, 故应增加惯性权重; 当聚焦距离变化率较小时, 应减小惯性权重, 有利于粒子种群较好地适应复杂的实际环境.
下面基于学习因子进一步提高搜索能力. 文献[11]采用三角函数法对惯性权重进行非线性改进, 本文将三角函数法用于对学习因子的非线性改进, 公式为
其中kmax为最大迭代次数. 当线性改进时, 学习因子c1和c2的取值范围一般为[0.5,2.5], 故本文设c10=c20=1,c11=c21=1.5. 在搜索初期, 粒子的自我学习能力较强, 社会学习能力较弱, 有利于全局搜索; 而在搜索后期, 粒子的社会学习能力较强, 自我学习能力较弱, 有利于收敛到全局最优解. 图2为改进粒子群优化算法的流程.
下面以式(5)为例, 说明采用改进粒子群优化算法求解SHEPWM非线性方程组的原理. 将式(5)中的α1,α2,…,αN转化为x1,x2,…,xN, 并表示成如下形式:
(16)
其中:fi表示变量x1,x2,…,xN∈的函数;PN是由给定的半周期镜像对称SHEPWM波形频域特性得到的实常数. 应用无约束优化方法求解式(16)时, 通常先将其转化为非线性最小值问题, 因此定义粒子的适应度函数为
(17)
应用改进的粒子群优化算法求解SHEPWM非线性方程组, 即求解适应度函数F(x)的最小值,F(x)的最小值对应的开关角即为式(16)非线性方程组的解,F(x)越小表示粒子的适应度越高, 方程组的解越精确.
为验证发射电流频域信息可控, 设置半周期镜像对称SHEPWM波形的主频为1,7,9次谐波. 结合能量守恒式(7)及Ai≤4Ud/π的原则, 设Ud=1 V, 主频谐波幅值为0.65 A(A1=A7=A9=0.65), 其余被控谐波幅值为0 (A2=…=A6=A8=A10=0), 被控谐波相位为0(θ1=…=θ10=0), 直流分量为0. 将上述设置的频域信息代入式(5), 可得半周期镜像对称SHEPWM非线性方程组为
(18)
应用标准粒子群优化算法(PSO)和本文改进粒子群优化算法(IPSO)进行求解, 程序参数设置如下: 迭代次数均为2 000, 收敛精度均为10-6, 粒子个数均为1 000. PSO的惯性权重ω=0.8, 学习因子c1=2,c2=1. 求解结果列于表1.
表1 两种粒子群优化算法的求解结果
图3 两种粒子群优化算法的求解过程对比Fig.3 Solution process comparison of two particle swarm optimization algorithm
由表1可见, IPSO算法小于PSO算法的适应度值, 证明了IPSO算法的求解精度高于PSO算法. 图3为PSO算法与IPSO算法求解过程对比. 由图3可见, IPSO算法快于PSO算法的收敛速度. 并且PSO算法在迭代次数为800时, 陷入了局部最优, 适应度值最终收敛到0.284, 而IPSO算法在迭代次数为100时收敛到10-3, 迭代次数为2 000时收敛到8.705×10-5. 从而验证了IPSO算法收敛速度更快、 求解精度更高、 全局搜索能力更强.
为进一步验证IPSO算法的正确性以及半周期镜像对称SHEPWM发射电流控制方法的有效性, 利用MATLAB/Simulink软件对半周期镜像对称SHEPWM波形进行仿真. 4个IGBT组成全桥拓扑, 直流侧电压为1 V, 发射线圈电阻为1 Ω, SHEPWM控制方法的开关信号来源于开关子系统, 利用改进粒子群优化算法求解得到开关信号的切换时刻.
相同开关次数下, 两种粒子群优化算法对SHEPWM波形的控制质量仿真结果分别如图4和图5所示. 图4为PSO算法求解的开关时刻对应半周期镜像对称SHEPWM波形的时域信息及其频域信息; 图5为IPSO算法求解的开关时刻对应半周期镜像对称SHEPWM波形的时域信息及其频域信息. 由图4(A)和图5(A)可见, SHEPWM波形的时域信息是半周期镜像对称的, 基频为16 Hz, 幅值为1 A, 半个周期内有10个间断点. 由图4(B)可见, PSO算法求解的开关时刻对应的SHEPWM波形频域信息1,9次谐波幅值达到0.65, 而7次谐波幅值未达到0.65, 且其余次谐波的幅值也较大, 无法省略不计. 由图5(B)可见, IPSO算法求解的开关时刻对应的SHEPWM波形频域信息1,7,9次谐波幅值相等且为0.65, 其余次谐波幅值较小, 与预设值相符, 从而验证了IPSO算法的优越性及半周期镜像对称SHEPWM波形设计方法的正确性.
图4 PSO算法半周期镜像对称SHEPWM波形Fig.4 Half-cycle mirror symmetric SHEPWM waveform of PSO algorithm
图5 IPSO算法半周期镜像对称SHEPWM波形Fig.5 Half-cycle mirror symmetric SHEPWM waveform of IPSO algorithm
由上述分析可知, 针对特定深度的探测目标, 半周期镜像对称SHEPWM波形可提高纵向分辨率. 为了对比分析半周期镜像对称SHEPWM波形与2n伪随机波形的优缺点, 本文以频率域电磁测深的纵向分辨率为依据, 通过正演分析, 进一步给出上述两种波形的性能对比结果.
本文选择的算例模型为经典的垂直单阻模型, 同时给出半周期镜像对称SHEPWM波形和标准的三频伪随机波形正演数值模拟拟断面图进行对比分析. 正演模型异常体的位置根据趋肤深度公式选择:
(19)
其中:δ为勘探深度;ρ为测量区域的平均电阻率;f为探测频率. 垂直单阻正演模型如图6所示. 由图6可见, 测点5有一个异常, 埋深为210 m, 大小为30 m×30 m×30 m, 电阻率为500 Ω·m, 距地表300 m内为均匀大地, 电阻率为50 Ω·m, 测点间距为50 m.
图7为单阻异常视电阻率等值线, 其中: (A)为半周期镜像对称SHEPWM波形16,112,144 Hz的正演数值模拟断面; (B)为标准三频伪随机波形16,32,64 Hz的正演数值模拟断面. 由图7(A)可见, 在测点5处异常体形态和轮廓清晰, 很好地保留了异常体信息, 与模型基本相符. 与图7(A)相比, 图7(B)伪随机构造的纵向范围严重拉伸, 未分辨出实际模型中的异常体.
图6 垂直单阻正演模型示意图Fig.6 Schematic diagram of vertical single resistance forward model
图7 单阻异常视电阻率等值线Fig.7 Apparent resistivity isolines of single resistance anomaly
根据式(19)可知, 半周期镜像对称SHEPWM的主频率为16,112,144 Hz, 对应的探测深度为629,237,209 m. 因此可得SHEPWM波形中的112 Hz勘探深度为237 m, 144 Hz勘探深度为209 m, 可用于辨识异常体(209~237 m). 而标准三频伪随机波形的主频率为16,32,64 Hz, 对应的探测深度为629,445,314 m, 与半周期镜像对称SHEPWM波形相比, 缺少辨识该深度异常体的主频频率, 因此其正演数值模拟断面未分辨出实际模型中的异常. 正演数值模拟与理论分析相符, 验证了半周期镜像对称SHEPWM发射电流控制方法的可行性.
综上所述, 针对传统方法控制的发射电流频率间隔固定、 频域不完全可控的问题, 本文提出了一种基于改进粒子群优化算法的可控发射电流频率的半周期镜像对称SHEPWM控制方法. 改进粒子群优化算法计算结果、 半周期镜像对称SHEPWM发射电流仿真结果和频率域电磁测深正演结果相符, 验证了本文方法的有效性. 与标准粒子群优化算法相比, 在求解SHEPWM非线性方程组时, 本文改进的粒子群优化算法提高了收敛速度与精度, 增强了全局搜索能力; 与2n伪随机方法相比, 在频率域电磁测深时, 本文的半周期镜像对称SHEPWM方法可根据探测目标控制发射电流奇次谐波的幅值、 相位、 分布及数量, 输出理想的发射电流, 从而提高了电磁测深的纵向分辨率.