本發(fā)明屬于目標(biāo)電磁散射特性數(shù)值計(jì)算
技術(shù)領(lǐng)域：
，特別是一種電大尺寸目標(biāo)的瞬態(tài)電磁特性快速提取方法。
背景技術(shù)：
：電磁計(jì)算的數(shù)值方法可分為頻域數(shù)值方法和時(shí)域數(shù)值方法，頻域數(shù)值方法包括矩量法(MOM)和有限元法(FEM)等，時(shí)域數(shù)值方法包括時(shí)域積分方程方法(TDIE)和時(shí)域有限差分方法(FDTD)等。頻域方法處理的是穩(wěn)態(tài)問題，一次計(jì)算只能分析單個(gè)頻率的響應(yīng)，而時(shí)域方法則可以處理瞬態(tài)問題，而且時(shí)域的激勵(lì)信號(hào)常為寬帶信號(hào)，可以通過一次計(jì)算獲得寬頻帶的響應(yīng)。因此，在分析目標(biāo)的瞬態(tài)與寬頻帶的電磁散射特性時(shí)，采用時(shí)域數(shù)值方法有明顯的優(yōu)勢。然而，時(shí)域積分方程方法(TDIE)和時(shí)域有限差分方法(FDTD)可以很好地解決電小尺寸物體的散射，但在計(jì)算電大物體的散射時(shí)，對計(jì)算機(jī)的配置要求過高。因此，這類方法不能夠有效地提取瞬態(tài)電磁散射特性。迭代推進(jìn)方法是用于求解目標(biāo)散射問題的一種比較新型的方法，拋物線方程方法屬于其中之一，它是波動(dòng)方程的一種近似形式，假設(shè)電磁波能量在沿著拋物線軸向的錐形區(qū)域內(nèi)傳播。拋物線方程方法為求解電磁散射提供了一種準(zhǔn)確、高效的計(jì)算方法，它的主要缺陷是只能對拋物線方向近軸區(qū)域內(nèi)的電磁散射進(jìn)行快速、準(zhǔn)確地計(jì)算，不過這種限制可以通過旋轉(zhuǎn)拋物線軸向來克服。Claerbout在文章“FundamentalsofGeophysicalDataProcessingwithApplicationtoPetroleumProspecting”中，使用了時(shí)域拋物線方程方法分析了地震波的傳播問題；Popov在文章“ModelingEMtransientpropagationoverirregulardispersiveboundary”和“Electromagneticpulsepropagationovernonuniformearthsurface:Numericalsimulation”中又將時(shí)域拋物線方程方法推廣到分析了二維電磁波的傳播問題。但是，現(xiàn)有文獻(xiàn)中沒有三維時(shí)域拋物線方程方法提取瞬態(tài)電磁散射特性的相關(guān)報(bào)道。技術(shù)實(shí)現(xiàn)要素：本發(fā)明的目的在于提供一種電大尺寸目標(biāo)的瞬態(tài)電磁特性快速提取方法，該方法大大簡化了求解矩陣，能夠快速得到電磁散射特性參數(shù)。實(shí)現(xiàn)本發(fā)明目的的技術(shù)解決方案為：一種電大尺寸目標(biāo)的瞬態(tài)電磁特性快速提取方法，步驟如下：步驟1、建立散射體的離散模型，確定時(shí)域拋物線的軸向方向作為x軸，采用網(wǎng)格對散射體沿拋物線的軸向方向進(jìn)行離散處理，形成垂直于x軸的若干個(gè)切面，再通過所有節(jié)點(diǎn)與四面體的幾何關(guān)系判斷節(jié)點(diǎn)是散射體的內(nèi)部點(diǎn)、外部點(diǎn)或者邊界點(diǎn)；步驟2、構(gòu)造矩陣方程：將時(shí)域拋物線方程在時(shí)間上使用CN差分格式，在空間上寫為交替方向隱式的差分形式，在迭代方向上的相鄰兩個(gè)面之間引入一個(gè)中間面，通過前一個(gè)面的散射場場值和前一個(gè)時(shí)間刻的散射場場值按行求解中間面的散射場場值，且各行之間求解相互獨(dú)立；同樣，通過中間面的散射場場值和前一個(gè)時(shí)間刻的散射場場值按列求解下一個(gè)面的散射場場值，且各列之間求解相互獨(dú)立，最后在散射體表面根據(jù)切向總場分量為0的方程以及拋物線方程，聯(lián)立構(gòu)造出矩陣方程；步驟3、令x軸方向?yàn)榇蟮纳⑸浞较颍来螌ρ剌S向方向的各個(gè)切面上的離散節(jié)點(diǎn)散射場場值進(jìn)行遞推求解，通過不斷更新各個(gè)切面上的邊界點(diǎn)的位置求解下一個(gè)切面上各個(gè)離散節(jié)點(diǎn)處的散射場場值；對各個(gè)切面上的離散節(jié)點(diǎn)散射場場值進(jìn)行遞推，求解最后一個(gè)切面的散射場場值；步驟4、在時(shí)間方向進(jìn)行迭代，重復(fù)步驟3；根據(jù)所需的觀察頻點(diǎn)對最后一個(gè)切面的散射場場值進(jìn)行傅立葉變換得到頻域的散射場場值，由近遠(yuǎn)場轉(zhuǎn)換求解目標(biāo)散射體的雙站雷達(dá)散射截面積。本發(fā)明與現(xiàn)有技術(shù)相比，其顯著優(yōu)點(diǎn)為：(1)方程形成簡單：本發(fā)明在拋物線軸向方向使用了差分格式，可將一個(gè)三維問題轉(zhuǎn)化為一系列的二維問題進(jìn)行求解，矩陣形成快捷簡便；(2)內(nèi)存需求低：在每個(gè)切面上，使用了交替方向隱式的差分格式，降低了求解矩陣的維數(shù)，對計(jì)算內(nèi)存的需求降低；(3)計(jì)算速度快：本發(fā)明交替方向隱式的差分格式使得求解矩陣為三對角矩陣，可通過追趕法求解，提高計(jì)算速度。下面結(jié)合附圖對本發(fā)明作進(jìn)一步詳細(xì)描述。附圖說明圖1為本發(fā)明的電大尺寸目標(biāo)的瞬態(tài)電磁特性快速提取方法在每個(gè)切面逐行/列求解示意圖。圖2為本發(fā)明實(shí)施例中觀察頻率為200MHz的散射體雙站RCS曲線圖。圖3為本發(fā)明實(shí)施例中觀察頻率為300MHz的散射體雙站RCS曲線圖。圖4為本發(fā)明實(shí)施例中觀察頻率為400MHz的散射體雙站RCS曲線圖。具體實(shí)施方式結(jié)合圖1，本發(fā)明的電大尺寸目標(biāo)的瞬態(tài)電磁特性快速提取方法，步驟如下：步驟1、建立散射體的離散模型，確定時(shí)域拋物線的軸向方向作為x軸，采用網(wǎng)格對散射體沿拋物線的軸向方向進(jìn)行離散處理，形成垂直于x軸的若干個(gè)切面，再通過所有節(jié)點(diǎn)與四面體的幾何關(guān)系判斷節(jié)點(diǎn)是散射體的內(nèi)部點(diǎn)、外部點(diǎn)或者邊界點(diǎn)；具體為：步驟1.1、對散射體進(jìn)行三角面元的面剖分，確定軸方向每個(gè)切面的方程，通過剖分網(wǎng)格的幾何關(guān)系求解三角面元與切面的交點(diǎn)，與該交點(diǎn)距離最近的差分網(wǎng)格點(diǎn)標(biāo)記為邊界點(diǎn)并求出該點(diǎn)法向分量；步驟1.2、對散射體進(jìn)行四面體的體剖分，通過判斷某點(diǎn)是否處于四面體內(nèi)部來區(qū)分該點(diǎn)處于散射體內(nèi)部或者散射體外部，并對這些點(diǎn)進(jìn)行標(biāo)記；以上即可完成目標(biāo)的建模，為下面的矩陣構(gòu)造以及求解奠定了基礎(chǔ)；步驟2、構(gòu)造矩陣方程，將時(shí)域拋物線方程在時(shí)間上使用Crank-Nicholson差分格式，在空間上寫為交替方向隱式的差分形式，在迭代方向(沿著拋物線軸向方向)上的相鄰兩個(gè)面之間引入一個(gè)中間面，通過前一個(gè)面的場值和前一個(gè)時(shí)間刻的場值按行求解中間面的散射場場值，且各行之間求解相互獨(dú)立；同樣，通過中間面的場值和前一個(gè)時(shí)間刻的場值按列求解下一個(gè)面的散射場場值，且各列之間求解相互獨(dú)立；最后在散射體表面根據(jù)切向電場分量為0的方程以及拋物線方程，聯(lián)立構(gòu)造出矩陣方程，具體步驟如下：步驟2.1、將物體全部包圍的閉合區(qū)域?yàn)橛?jì)算區(qū)域，設(shè)該計(jì)算區(qū)域離散為N個(gè)切面，每個(gè)切面離散為M*M個(gè)離散節(jié)點(diǎn)，在三維情況下，時(shí)域拋物線方程為：∂2Π∂y2+∂2Π∂z2-2c∂2Π∂x∂t=0---(1)]]>其中，∏代表散射場場值；將式(1)寫成交替方向隱式的差分形式：-ry1+2ry-ryΠp-1,qn+1/2,l+Πp-1,qn+1/2,l+1Πp,qn+1/2,l+Πp,qn+1/2,l+1Πp+1,qn+1/2,l+Πp+1,qn+1/2,l+1=rz1-2rzrzΠp,q-1n,l+Πp,q-1n,l+1Πp,qn,l+Πp,qn,l+1Πp,q+1n,l+Πp,q+1n,l+1-2Πp,qn,l+2Πp,qn+1/2,l---(2)]]>-rz1+2rz-rzΠp,q-1n+1,l+Πp,q-1n+1,l+1Πp,qn+1,l+Πp,qn+1,l+1Πp,q+1n+1,l+Πp,q+1n+1,l+1=ry1-2ryryΠp-1,qn+1/2,l+Πp-1,qn+1/2,l+1Πp,qn+1/2,l+Πp,qn+1/2,l+1Πp+1,qn+1/2,l+Πp+1,qn+1/2,l+1-2Πp,qn+1/2,l+2Πp,qn+1,l]]>(3)]]>其中，和分別表示l時(shí)間刻上第n個(gè)面上分布在第p列的連續(xù)的三個(gè)點(diǎn)的散射場場值，和分別表示l+1時(shí)間刻上第n個(gè)面上分布在第p列的連續(xù)的三個(gè)點(diǎn)的散射場場值，和分別表示l時(shí)間刻上n+1/2中間面上分布在第q行的連續(xù)的三個(gè)點(diǎn)的散射場場值，和分別表示l+1時(shí)間刻上n+1/2中間面上分布在第q行的連續(xù)的三個(gè)點(diǎn)的散射場場值，和分別表示l時(shí)間刻上第n+1面上分布在第p列的連續(xù)的三個(gè)點(diǎn)的散射場場值，和分別表示l+1時(shí)間刻上第n+1面上分布在第p列的連續(xù)的三個(gè)點(diǎn)的散射場場值；△t為時(shí)間迭代步長，△x、△y和△z分別為x、y、z方向上標(biāo)準(zhǔn)網(wǎng)格點(diǎn)的長度，c為光速，1≤n≤N，1≤p≤M，1≤q≤M，1≤l≤L，L為總時(shí)間刻個(gè)數(shù)；矢量拋物線方程由x、y、z三個(gè)方向上的標(biāo)量拋物線方程構(gòu)成：∂2Πx∂y2(x,y,z)+∂2Πx∂z2(x,y,z)-2c∂Πx∂x∂t(x,y,z)=0∂2Πy∂y2(x,y,z)+∂2Πy∂z2(x,y,z)-2c∂Πy∂x∂t(x,y,z)=0∂2Πz∂y2(x,y,z)+∂2Πz∂z2(x,y,z)-2c∂Πz∂x∂t(x,y,z)=0---(4)]]>按照上面的推導(dǎo)，將式(4)寫成x、y、z三個(gè)方向的交替方向隱式格式；步驟2.2、對于目標(biāo)邊界點(diǎn)，假設(shè)P(xb,yb,zb)為散射體表面上的點(diǎn)，n＝(nx,ny,nz)為P點(diǎn)的法向方向，在金屬表面，切向電場為零，由則對應(yīng)的頻域表達(dá)式為：nxEy(xb,yb,zb)-nyEx(xb,yb,zb)＝0nxEz(xb,yb,zb)-nzEx(xb,yb,zb)＝0(5)nyEz(xb,yb,zb)-nzEy(xb,yb,zb)＝0其中，Ex(xb,yb,zb)、Ey(xb,yb,zb)、Ez(xb,yb,zb)分別為P點(diǎn)電場在x軸、y軸、z 軸方向上的分量；則其對應(yīng)的時(shí)域表達(dá)式如下所示：nxΠy(xb,yb,zb,t)-nyΠx(xb,yb,zb,t)=-[nxEyi(xb,yb,zb,t-xb/c)-nyExi(xb,yb,zb,t-xb/c)]nxΠz(xb,yb,zb,t)-nzΠx(xb,yb,zb,t)=-[nxEzi(xb,yb,zb,t-xb/c)-nzExi(xb,yb,zb,t-xb/c)]nyΠz(xb,yb,zb,t)-nzΠy(xb,yb,zb,t)=-[nyEzi(xb,yb,zb,t-xb/c)-nzEyi(xb,yb,zb,t-xb/c)]---(6)]]>其中，分別為入射電場在x軸、y軸、z軸方向上的分量；在球坐標(biāo)系中，入射電場為：Ei(r,t)=E0cos[2πf0(t-r·k^/c)]e-(t-r·k^/c-τp)22σ2---(7)]]>其中，E0為入射電場的傳播方向，f0為調(diào)制的中心頻率，σ＝3/(πfbw)表示脈沖寬度，fbw為頻帶寬度，為入射方向，r·k^=xsinθcosφ+ysinθsinφ+cosθ,]]>θ,φ分別代表球坐標(biāo)系下與z軸和xoy平面的夾角，τp表示入射場的延時(shí)，r表示點(diǎn)(x,y,z)到原點(diǎn)的距離；由于式(6)構(gòu)成的矩陣秩為2，不足以唯一確定方程的解，所以由該方程結(jié)合拋物線方程即可構(gòu)成秩為3的矩陣，此時(shí)便可求解邊界上的場值；為了提高求解速度，將矩陣寫成三對角矩陣，將式(6)分成3種情況：(1)對于x方向的法向量不為0，y、z方向的法向量任意的情況：nxΠy(xb,yb,zb,t)-nyΠx(xb,yb,zb,t)=-[nxEyi(xb,yb,zb,t-xb/c)-nyExi(xb,yb,zb,t-xb/c)]nxΠz(xb,yb,zb,t)-nzΠx(xb,yb,zb,t)=-[nxEzi(xb,yb,zb,t-xb/c)-nzExi(xb,yb,zb,t-xb/c)]---(8)]]>其中，Πx為散射場場值的x方向的分量，將式中的Πx近似為前一個(gè)面的Πx，此時(shí)聯(lián)立x方向的拋物線方程，可求解出一組唯一的方程的解；(2)對于x方向的法向量為0，y的法向量不為0，z方向的法向量任意的情況：-nyΠx(xb,yb,zb,t)=-[-nyExi(xb,yb,zb,t-xb/c)]nyΠz(xb,yb,zb,t)-nzΠy(xb,yb,zb,t)=-[nyEzi(xb,yb,zb,t-xb/c)-nzEyi(xb,yb,zb,t-xb/c)]]]>(9)]]>其中，Πy為散射場場值的y方向的分量，將式中的Πy近似為前一個(gè)面的Πy，此時(shí)聯(lián)立y方向的拋物線方程，可求解出一組唯一的方程的解；(3)對于x和y方向的法向量為0，z方向的法向量不為0的情況，-nzΠx(xb,yb,zb,t)=-[nxEzi(xb,yb,zb,t-xb/c)-nzExi(xb,yb,zb,t-xb/c)]-nzΠy(xb,yb,zb,t)=-[nyEzi(xb,yb,zb,t-xb/c)-nzEyi(xb,yb,zb,t-xb/c)]---(10)]]>此時(shí)聯(lián)立x方向的拋物線方程，可求解出一組唯一的方程的解；步驟3、在某個(gè)固定的時(shí)間刻，令x軸方向?yàn)榇蟮纳⑸浞较蛞来螌ρ剌S向方向的各個(gè)切面上的節(jié)點(diǎn)散射場場值進(jìn)行遞推求解，通過不斷更新邊界點(diǎn)的信息以及方程的右邊向量來求解下一個(gè)切面上各個(gè)離散節(jié)點(diǎn)處的散射場場值；對各個(gè)切面上的節(jié)點(diǎn)散射場場值進(jìn)行遞推，求解最后一個(gè)切面散射場場值，具體過程如下：步驟3.1、將前一個(gè)時(shí)間刻各個(gè)離散的節(jié)點(diǎn)的散射場場值，以及前一個(gè)切面各個(gè)離散的節(jié)點(diǎn)的散射場電場值作為當(dāng)前切面求解時(shí)的右邊向量；步驟3.2、在當(dāng)前切面所確定的邊界點(diǎn)處，加入切向的電場分量為0的邊界條件和拋物線方程，處于物體內(nèi)部的節(jié)點(diǎn)電場值賦值為0，形成當(dāng)前切面更新后的矩陣方程；步驟3.3、求解步驟3.2中更新后的矩陣方程，方程的解即為當(dāng)前時(shí)間刻上當(dāng)前切面各個(gè)離散的節(jié)點(diǎn)的散射場場值；每個(gè)切面上，根據(jù)處于不同的位置的離散節(jié)點(diǎn)，帶入不同的離散方程，由前一個(gè)面的散射場場值求得下一個(gè)面的散射場場值，不斷遞推得到最后一個(gè)切面的散射場場值；步驟4、在時(shí)間方向進(jìn)行迭代，重復(fù)步驟3；根據(jù)所需的觀察頻點(diǎn)對最后一個(gè)切面的散射場場值進(jìn)行傅立葉變換得到頻域場值，由近遠(yuǎn)場轉(zhuǎn)換求解目標(biāo)散射體雙站雷達(dá)散射截面積RCS；對于指定的觀察頻點(diǎn)ω，時(shí)域拋物線方程與頻域時(shí)域拋物線方程解的轉(zhuǎn)換關(guān)系為：u(x,y,z,ω)=[∫-∞∞Π(x,y,z,t)e-iωtdt]/F~(ω)---(11)]]>其中，表示入射場的傅立葉變換，u(x,y,z,ω)表示某節(jié)點(diǎn)對應(yīng)頻點(diǎn)ω的頻域場值，∏(x,y,z,t)表示節(jié)點(diǎn)(x,y,z)的時(shí)域散射場場值，i表示純虛數(shù)；在頻域中，利用近遠(yuǎn)場轉(zhuǎn)換公式計(jì)算雷達(dá)散射截面積；三維坐標(biāo)系下，在(θ,φ)方向的雙站RCS為：σ(θ,φ)=limr→∞4πr2|us(x,y,z)|2|ui(x,y,z)|2---(12)]]>其中us和ui分別表示散射場和入射場的電場分量，π為圓周率；θ代表球坐標(biāo)系下向量(x,y,z)與z軸的夾角，φ代表球坐標(biāo)系下向量(x,y,z)與xoy面的夾角。使用旋轉(zhuǎn)拋物線方程方法來克服拋物線方程的近軸限制，通過旋轉(zhuǎn)拋物線的軸向方向來確定各個(gè)方向的散射場，然后通過近場遠(yuǎn)推獲得遠(yuǎn)區(qū)的散射場，從而確定目標(biāo)的雙站RCS。下面結(jié)合具體實(shí)施例對本發(fā)明作進(jìn)一步說明。實(shí)施例1本實(shí)施例進(jìn)行了電磁散射的典型仿真，仿真在主頻2.68GHz、內(nèi)存8GB的計(jì)算機(jī)上實(shí)現(xiàn)，以半徑5m的金屬球?yàn)槔?，入射波的方向θ?0°，φ＝0°，入射波中心頻率為300MHz，帶寬為600MHz，觀察頻點(diǎn)選擇200MHz，300MHz，400MHz三個(gè)頻點(diǎn)。為了驗(yàn)證本發(fā)明方法的正確性，與MIE結(jié)果作為參照；圖2-圖4為兩種電磁散射特性仿真的RCS曲線圖，圖2的觀察頻率為200MHz，圖3的觀察頻率為300MHz，圖4的觀察頻率為400MHz，從圖中的曲線可以看出，本發(fā)明的方法與正確的數(shù)值結(jié)果吻合。表1為本專利方法和時(shí)域有限差分方法(FDTD)的時(shí)間和內(nèi)存對比：表1方法內(nèi)存(GB)CPU時(shí)間(s)FDTD4.118000本專利方法1.5752通過表1中本方法與現(xiàn)有的FDTD方法對比，本專利方法不僅節(jié)省內(nèi)存，還能夠快速仿真分析目標(biāo)物體的電磁散射特性。綜上所述，本發(fā)明將復(fù)雜的三維問題分解為很多個(gè)二維的問題進(jìn)行求解，追趕法解三對角矩陣更加快了計(jì)算速度，其實(shí)現(xiàn)過程靈活自由，具有很強(qiáng)的實(shí)際工程應(yīng)用價(jià)值。當(dāng)前第1頁1 2 3