于生寶, 房 鈺, 高麗輝, 刁 庶

(吉林大學 儀器科學與電氣工程學院, 長春 130061)

電磁法是礦產(chǎn)資源勘探、 環(huán)境與工程地質(zhì)學、 區(qū)域構(gòu)造研究的一種有效方法, 可在地面、 空中及海上進行資源探測[1]. 根據(jù)發(fā)射場的性質(zhì), 電磁法可分為時間域電磁法和頻率域電磁法[2]. 相比于時間域電磁法, 頻率域電磁法在整個采集過程中電流連續(xù)發(fā)射, 通過獲取總場反應(yīng)地下電阻率的信息, 信號較強. 因此, 頻率域電磁法適用于觀測面積大、 勘探深度較深的復雜地形區(qū)域探測[3]. 目前, 發(fā)射電流的波形多為固定頻比或頻差, 如2n偽隨機波和方波等, 這些波形的能量均按特定比例等距離分布在各次諧波上[4]. 當針對特定深度的目標進行高精度電磁勘探時, 由于傳統(tǒng)發(fā)射電流波形頻率間隔固定, 需發(fā)射不同基頻的波形才能滿足探測要求, 而多次發(fā)射又會降低勘探效率, 因此存在勘探精度與效率相互制約的問題.

針對上述問題, 本文提出一種基于改進粒子群優(yōu)化算法的半周期鏡像對稱選擇性諧波消除脈寬調(diào)制(selective harmonic eliminated pulse width modulation, SHEPWM)發(fā)射電流控制方法. 在頻率域電磁測深過程中, 先根據(jù)探測深度計算出發(fā)射電流期望的頻譜, 再根據(jù)發(fā)射電流的時頻域特性, 建立控制時域信息(開關(guān)時刻、 電流值)和頻域信息(直流分量和諧波幅值、 相位)的半周期鏡像對稱SHEPWM非線性方程組, 最后采用改進的粒子群優(yōu)化算法對SHEPWM非線性方程組進行求解, 獲得對應(yīng)的時域開關(guān)序列, 控制發(fā)射系統(tǒng)逆變器, 從而得到理想的發(fā)射電流, 提高電磁測深的縱向分辨率.

1 半周期鏡像對稱SHEPWM波形設(shè)計方法

SHEPWM方法是一種逆變器常用的計算調(diào)制法, 廣泛應(yīng)用于電力、 電子等領(lǐng)域. 該方法先通過對逆變器輸出電壓波形進行Fourier分解, 建立與開關(guān)時刻和波形諧波分量有關(guān)的非線性方程組, 求解方程組得到滿足要求的開關(guān)角, 然后利用這些開關(guān)角控制逆變器中開關(guān)器件的導通和關(guān)斷, 從而消除直流分量和某些特定的低次諧波, 得到所需的基波分量[5]. 當應(yīng)用在頻率域電磁測深領(lǐng)域時, 波形的諧波分量是有用信號, 此時半周期鏡像對稱SHEPWM方法可根據(jù)勘探目標地質(zhì)結(jié)構(gòu)設(shè)置發(fā)射電流奇次諧波的幅值、 相位、 分布及數(shù)量, 控制探測所需發(fā)射電流的頻域信息. 下面給出針對特定深度勘探目標的半周期鏡像對稱SHEPWM波形的設(shè)計方法.

SHEPWM周期性輸出電壓波形展開成一個收斂的Fourier級數(shù), 表示為直流分量、 基波分量和諧波分量的疊加：

(1)

其中：i表示諧波次數(shù);A0表示直流分量;Ai表示第i次諧波的振幅;θi表示第i次諧波的相位;ω=2πf表示角頻率,f表示輸出電壓波形頻率.a0,ai,bi,A0,Ai,θi滿足如下關(guān)系：

(2)

式(2)中i=1,2,…, 系數(shù)a0,ai,bi根據(jù)Fourier分析可得：

(3)

其中:i=1,2,…;T=2π為輸出電壓波形的周期. 圖1為一個半周期鏡像對稱的SHEPWM發(fā)射波形, 其中:Ud為幅度;α1,α2,…,αN為半個周期內(nèi)的開關(guān)時刻[6].

圖1 半周期鏡像對稱SHEPWM波形Fig.1 Half-cycle mirror symmetric SHEPWM waveform

根據(jù)式(3)可得半周期鏡像對稱SHEPWM波形的a0,ai,bi：

(4)

聯(lián)立式(2)和式(4), 可得半周期鏡像對稱SHEPWM波形的非線性方程組為

(5)

0<α1<α2<…<αN<π.

(6)

式(5)左邊是與開關(guān)角相關(guān)的多項式, 右邊是與發(fā)射電流被控諧波頻域特性幅值A(chǔ)i和相位θi相關(guān)的多項式. 因此, 可通過求解式(5), 精確計算出符合發(fā)射要求的開關(guān)角度. 式(6)為式(5)的限定條件. 此時, 發(fā)射電流奇次諧波幅值、 相位、 分布及數(shù)量可控, 偶次諧波幅值和直流分量為0.

根據(jù)頻率域電磁探測原理, 探測過程中可不使用相位, 故可設(shè)θi=0, 根據(jù)實際探測地形設(shè)置諧波幅值A(chǔ)i, 如Ai的最大值不應(yīng)大于4Ud/π, 且所有諧波幅值應(yīng)滿足能量守恒定律[7]：

(7)

2 改進粒子群優(yōu)化算法求解SHEPWM非線性方程組

SHEPWM控制方法的關(guān)鍵是對開關(guān)器件開關(guān)角度的求取, 即非線性方程組的求解. 目前, 求解非線性方程組的有效方法主要分為兩類： 一類是經(jīng)典的數(shù)值迭代法, 如牛頓迭代法、 擬牛頓法、 梯度法等, 這些算法均具有較快的局部收斂性, 但均需設(shè)置合適的初始值, 且隨著未知數(shù)的增加, 算法復雜性增加； 第二類是仿生優(yōu)化算法, 如粒子群優(yōu)化算法[8]、 遺傳算法、 模擬退火算法、 神經(jīng)網(wǎng)絡(luò)算法等, 這些方法多采用分析方程形成目標函數(shù), 將其轉(zhuǎn)換為最優(yōu)化問題. 其中的粒子群優(yōu)化算法由于其算法簡單、 容易實現(xiàn)被廣泛應(yīng)用, 但該方法也存在易陷入局部最優(yōu)的缺點[9]. 為提高算法的性能, 文獻[10]提出了線性改進慣性權(quán)重和學習因子, 并結(jié)合了模擬退火算法的方法; 文獻[11]提出了用余弦函數(shù)對慣性權(quán)重進行非線性改進的方法, 提高了收斂速度和尋優(yōu)能力. 本文采用動態(tài)自適應(yīng)調(diào)整的慣性權(quán)重, 并對學習因子進行三角函數(shù)非線性改進的方法, 以解決粒子群優(yōu)化算法求解更多未知數(shù)時適應(yīng)度函數(shù)值不能收斂到零的問題.

標準粒子群優(yōu)化算法中, 每個粒子的狀態(tài)都由位置、 速度和目標函數(shù)的適應(yīng)度值構(gòu)成, 根據(jù)自身經(jīng)驗和種群經(jīng)驗在D維有限空間中搜索最優(yōu)值. 每個粒子的速度[12]和位置更新公式為

標準粒子群優(yōu)化算法的慣性權(quán)重和學習因子都是常數(shù), 無法根據(jù)搜索的情況進行動態(tài)調(diào)整, 尋優(yōu)能力有限. 為了動態(tài)調(diào)節(jié)算法的全局搜索和局部搜索能力, 本文引入文獻[13]定義的平均聚焦距離、 最大聚焦距離和聚焦距離變化率：

其中：D為粒子的維數(shù);M為粒子的個數(shù). 通過對平均聚焦距離和最大聚焦距離的更新, 得到每次迭代時的聚焦距離變化率, 由此判斷本次粒子應(yīng)提高全局搜索能力還是局部搜索能力, 從而動態(tài)自適應(yīng)地調(diào)整慣性權(quán)重. 本文采用文獻[13]定義的慣性權(quán)重公式：

(13)

圖2 改進粒子群優(yōu)化算法流程Fig.2 Flow chart of improved particle swarm optimization algorithm

其中：z1=0.3,z2=0.2;r為一個[0,1]間均勻分布的隨機數(shù). 這種動態(tài)自適應(yīng)調(diào)整慣性權(quán)重可靈活地調(diào)節(jié)全局搜索與局部搜索的能力. 當聚焦距離變化率較大時, 粒子的全局搜索能力較差, 故應(yīng)增加慣性權(quán)重； 當聚焦距離變化率較小時, 應(yīng)減小慣性權(quán)重, 有利于粒子種群較好地適應(yīng)復雜的實際環(huán)境.

下面基于學習因子進一步提高搜索能力. 文獻[11]采用三角函數(shù)法對慣性權(quán)重進行非線性改進, 本文將三角函數(shù)法用于對學習因子的非線性改進, 公式為

其中kmax為最大迭代次數(shù). 當線性改進時, 學習因子c1和c2的取值范圍一般為[0.5,2.5], 故本文設(shè)c10=c20=1,c11=c21=1.5. 在搜索初期, 粒子的自我學習能力較強, 社會學習能力較弱, 有利于全局搜索; 而在搜索后期, 粒子的社會學習能力較強, 自我學習能力較弱, 有利于收斂到全局最優(yōu)解. 圖2為改進粒子群優(yōu)化算法的流程.

下面以式(5)為例, 說明采用改進粒子群優(yōu)化算法求解SHEPWM非線性方程組的原理. 將式(5)中的α1,α2,…,αN轉(zhuǎn)化為x1,x2,…,xN, 并表示成如下形式：

(16)

其中：fi表示變量x1,x2,…,xN∈的函數(shù);PN是由給定的半周期鏡像對稱SHEPWM波形頻域特性得到的實常數(shù). 應(yīng)用無約束優(yōu)化方法求解式(16)時, 通常先將其轉(zhuǎn)化為非線性最小值問題, 因此定義粒子的適應(yīng)度函數(shù)為

(17)

應(yīng)用改進的粒子群優(yōu)化算法求解SHEPWM非線性方程組, 即求解適應(yīng)度函數(shù)F(x)的最小值,F(x)的最小值對應(yīng)的開關(guān)角即為式(16)非線性方程組的解,F(x)越小表示粒子的適應(yīng)度越高, 方程組的解越精確.

3 仿真驗證

3.1 改進粒子群優(yōu)化算法求解結(jié)果仿真驗證

為驗證發(fā)射電流頻域信息可控, 設(shè)置半周期鏡像對稱SHEPWM波形的主頻為1,7,9次諧波. 結(jié)合能量守恒式(7)及Ai≤4Ud/π的原則, 設(shè)Ud=1 V, 主頻諧波幅值為0.65 A(A1=A7=A9=0.65), 其余被控諧波幅值為0 (A2=…=A6=A8=A10=0), 被控諧波相位為0(θ1=…=θ10=0), 直流分量為0. 將上述設(shè)置的頻域信息代入式(5), 可得半周期鏡像對稱SHEPWM非線性方程組為

(18)

應(yīng)用標準粒子群優(yōu)化算法(PSO)和本文改進粒子群優(yōu)化算法(IPSO)進行求解, 程序參數(shù)設(shè)置如下： 迭代次數(shù)均為2 000, 收斂精度均為10-6, 粒子個數(shù)均為1 000. PSO的慣性權(quán)重ω=0.8, 學習因子c1=2,c2=1. 求解結(jié)果列于表1.

表1 兩種粒子群優(yōu)化算法的求解結(jié)果

圖3 兩種粒子群優(yōu)化算法的求解過程對比Fig.3 Solution process comparison of two particle swarm optimization algorithm

由表1可見, IPSO算法小于PSO算法的適應(yīng)度值, 證明了IPSO算法的求解精度高于PSO算法. 圖3為PSO算法與IPSO算法求解過程對比. 由圖3可見, IPSO算法快于PSO算法的收斂速度. 并且PSO算法在迭代次數(shù)為800時, 陷入了局部最優(yōu), 適應(yīng)度值最終收斂到0.284, 而IPSO算法在迭代次數(shù)為100時收斂到10-3, 迭代次數(shù)為2 000時收斂到8.705×10-5. 從而驗證了IPSO算法收斂速度更快、 求解精度更高、 全局搜索能力更強.

3.2 SHEPWM頻率域電磁發(fā)射電流仿真驗證

為進一步驗證IPSO算法的正確性以及半周期鏡像對稱SHEPWM發(fā)射電流控制方法的有效性, 利用MATLAB/Simulink軟件對半周期鏡像對稱SHEPWM波形進行仿真. 4個IGBT組成全橋拓撲, 直流側(cè)電壓為1 V, 發(fā)射線圈電阻為1 Ω, SHEPWM控制方法的開關(guān)信號來源于開關(guān)子系統(tǒng), 利用改進粒子群優(yōu)化算法求解得到開關(guān)信號的切換時刻.

相同開關(guān)次數(shù)下, 兩種粒子群優(yōu)化算法對SHEPWM波形的控制質(zhì)量仿真結(jié)果分別如圖4和圖5所示. 圖4為PSO算法求解的開關(guān)時刻對應(yīng)半周期鏡像對稱SHEPWM波形的時域信息及其頻域信息; 圖5為IPSO算法求解的開關(guān)時刻對應(yīng)半周期鏡像對稱SHEPWM波形的時域信息及其頻域信息. 由圖4(A)和圖5(A)可見, SHEPWM波形的時域信息是半周期鏡像對稱的, 基頻為16 Hz, 幅值為1 A, 半個周期內(nèi)有10個間斷點. 由圖4(B)可見, PSO算法求解的開關(guān)時刻對應(yīng)的SHEPWM波形頻域信息1,9次諧波幅值達到0.65, 而7次諧波幅值未達到0.65, 且其余次諧波的幅值也較大, 無法省略不計. 由圖5(B)可見, IPSO算法求解的開關(guān)時刻對應(yīng)的SHEPWM波形頻域信息1,7,9次諧波幅值相等且為0.65, 其余次諧波幅值較小, 與預設(shè)值相符, 從而驗證了IPSO算法的優(yōu)越性及半周期鏡像對稱SHEPWM波形設(shè)計方法的正確性.

圖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

3.3 正演模型仿真驗證

由上述分析可知, 針對特定深度的探測目標, 半周期鏡像對稱SHEPWM波形可提高縱向分辨率. 為了對比分析半周期鏡像對稱SHEPWM波形與2n偽隨機波形的優(yōu)缺點, 本文以頻率域電磁測深的縱向分辨率為依據(jù), 通過正演分析, 進一步給出上述兩種波形的性能對比結(jié)果.

本文選擇的算例模型為經(jīng)典的垂直單阻模型, 同時給出半周期鏡像對稱SHEPWM波形和標準的三頻偽隨機波形正演數(shù)值模擬擬斷面圖進行對比分析. 正演模型異常體的位置根據(jù)趨膚深度公式選擇：

(19)

其中：δ為勘探深度;ρ為測量區(qū)域的平均電阻率;f為探測頻率. 垂直單阻正演模型如圖6所示. 由圖6可見, 測點5有一個異常, 埋深為210 m, 大小為30 m×30 m×30 m, 電阻率為500 Ω·m, 距地表300 m內(nèi)為均勻大地, 電阻率為50 Ω·m, 測點間距為50 m.

圖7為單阻異常視電阻率等值線, 其中: (A)為半周期鏡像對稱SHEPWM波形16,112,144 Hz的正演數(shù)值模擬斷面; (B)為標準三頻偽隨機波形16,32,64 Hz的正演數(shù)值模擬斷面. 由圖7(A)可見, 在測點5處異常體形態(tài)和輪廓清晰, 很好地保留了異常體信息, 與模型基本相符. 與圖7(A)相比, 圖7(B)偽隨機構(gòu)造的縱向范圍嚴重拉伸, 未分辨出實際模型中的異常體.

圖6 垂直單阻正演模型示意圖Fig.6 Schematic diagram of vertical single resistance forward model

圖7 單阻異常視電阻率等值線Fig.7 Apparent resistivity isolines of single resistance anomaly

根據(jù)式(19)可知, 半周期鏡像對稱SHEPWM的主頻率為16,112,144 Hz, 對應(yīng)的探測深度為629,237,209 m. 因此可得SHEPWM波形中的112 Hz勘探深度為237 m, 144 Hz勘探深度為209 m, 可用于辨識異常體(209～237 m). 而標準三頻偽隨機波形的主頻率為16,32,64 Hz, 對應(yīng)的探測深度為629,445,314 m, 與半周期鏡像對稱SHEPWM波形相比, 缺少辨識該深度異常體的主頻頻率, 因此其正演數(shù)值模擬斷面未分辨出實際模型中的異常. 正演數(shù)值模擬與理論分析相符, 驗證了半周期鏡像對稱SHEPWM發(fā)射電流控制方法的可行性.

綜上所述, 針對傳統(tǒng)方法控制的發(fā)射電流頻率間隔固定、 頻域不完全可控的問題, 本文提出了一種基于改進粒子群優(yōu)化算法的可控發(fā)射電流頻率的半周期鏡像對稱SHEPWM控制方法. 改進粒子群優(yōu)化算法計算結(jié)果、 半周期鏡像對稱SHEPWM發(fā)射電流仿真結(jié)果和頻率域電磁測深正演結(jié)果相符, 驗證了本文方法的有效性. 與標準粒子群優(yōu)化算法相比, 在求解SHEPWM非線性方程組時, 本文改進的粒子群優(yōu)化算法提高了收斂速度與精度, 增強了全局搜索能力; 與2n偽隨機方法相比, 在頻率域電磁測深時, 本文的半周期鏡像對稱SHEPWM方法可根據(jù)探測目標控制發(fā)射電流奇次諧波的幅值、 相位、 分布及數(shù)量, 輸出理想的發(fā)射電流, 從而提高了電磁測深的縱向分辨率.