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

(吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院, 長(zhǎng)春 130061)

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

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

1 半周期鏡像對(duì)稱(chēng)SHEPWM波形設(shè)計(jì)方法

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

SHEPWM周期性輸出電壓波形展開(kāi)成一個(gè)收斂的Fourier級(jí)數(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為一個(gè)半周期鏡像對(duì)稱(chēng)的SHEPWM發(fā)射波形, 其中:Ud為幅度;α1,α2,…,αN為半個(gè)周期內(nèi)的開(kāi)關(guān)時(shí)刻[6].

圖1 半周期鏡像對(duì)稱(chēng)SHEPWM波形Fig.1 Half-cycle mirror symmetric SHEPWM waveform

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

(4)

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

(5)

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

(6)

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

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

(7)

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

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

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

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

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

(13)

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

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

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

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

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

(16)

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

(17)

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

3 仿真驗(yàn)證

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

為驗(yàn)證發(fā)射電流頻域信息可控, 設(shè)置半周期鏡像對(duì)稱(chēng)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), 可得半周期鏡像對(duì)稱(chēng)SHEPWM非線性方程組為

(18)

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

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

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

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

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

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

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

圖4 PSO算法半周期鏡像對(duì)稱(chēng)SHEPWM波形Fig.4 Half-cycle mirror symmetric SHEPWM waveform of PSO algorithm

圖5 IPSO算法半周期鏡像對(duì)稱(chēng)SHEPWM波形Fig.5 Half-cycle mirror symmetric SHEPWM waveform of IPSO algorithm

3.3 正演模型仿真驗(yàn)證

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

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

(19)

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

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

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

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

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

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