CN103235193B - The numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band - Google Patents

Abstract

The invention discloses the numerical method of satellite Electromagnetic Scattering Characteristics in a kind of millimere-wave band.For satellite this metal medium mixed structure, it is only necessary to satellite metal part is set up improved Electric Field Integral Equation, and media fraction need not set up equation, ultimately forms the equation that condition is good, it is simple to iterative.Only need satellite metal surface is carried out mesh generation, it is not necessary to coating media part is carried out mesh generation, and singularity is easily processed.Due to the construction features of satellite plagioclase, the high directivity of transfer factor component, on the basis of the quick multistage son of multilamellar, far field effect be have employed accelerated method, effectively reduce calculating memory requirements and save the calculating time.

Description

The numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band

Technical field

The invention belongs to the quick computing technique of electromagnetic characteristic of scattering, particularly satellite electromagnetic scattering in millimere-wave band The numerical method of characteristic.

Background technology

Satellite is typically mainly made up of metal material, in order to reach anticorrosion, decoration and the purpose such as stealthy, it will usually adopt It is coated (QJ2754-95) with dielectric material.Electromagnetic scattering by metal medium coating mixed structure this for satellite is asked Topic, scholars propose a series of electromagnetism numerical computation method, wherein moment method (MoM) (Ding-based on integral equation Feng Yu, Si-Yuan He, Hai-Tao Chen, Guo-Qiang Zhu, " Research on the electromagnetic scattering of3D target coated with anisotropic medium using Impedance boundary condition ", Microwave and Optical Technology Letters, 2011, 53 (2): 458 462.) widest application has been obtained.

For this metal medium mixed structure, metal part is normally used as perfect electronic conductor (PEC) and processes, and Easily quilt cover integral Equation Methods (SIE) carrys out analysis and solution, but for media fraction, no matter is to use volume integral equations side Method, is also based on Line Integral equation method and their differentiation algorithm of the principle of equal effects, relative to the electromagnetism of simple metal target Calculation and analysis methods, these algorithms are required for more internal memory and time resource.And due to satellite coating media part very Thin feature, if selecting Line Integral equation to analyze media fraction, the singularity between upper and lower surface is more intractable；If adopted Analyze media fraction with volume integral, subdivision mesh distortion can be caused, ultimately result in matrix equation condition poor, it is difficult to solve.

Satellite in millimeter wave band belongs to Electrically large size object, just adds even with the quick multistage sub-technology of multilamellar Speed calculates, and still suffers from calculating the problem that resource consumption is huge.In Multilevel fast multipole algorithm, along with the increase of the number of plies, The amount of calculation of transfer matrix can sharply increase.Transfer matrix between two far field groups is by ball Hankel function, Legnedre series Composition, transfer factor is defined within the unit ball of a large amount of angular spectrum component, even if using interpolation and anti-interpolation method, whole turn The calculating moving past journey is also quite time-consuming.

Summary of the invention

It is an object of the invention to provide the numerical method of satellite Electromagnetic Scattering Characteristics in a kind of millimere-wave band, thus realize The method quickly obtaining Electromagnetic Scattering Characteristics parameter.

The technical solution realizing the object of the invention is: the numerical value side of satellite Electromagnetic Scattering Characteristics in a kind of millimere-wave band Method, step is as follows:

The first step, makes uniform plane wave be irradiated on a satellite structure scribbling film dielectric layer, will on the surface of satellite Produce surface induction surface current J_S, surface charge ρ_S, produce body polarization electric current J in thin-medium_pol, body polarization magnetic current M_polWith polarization electricity Lotus ρ_S,pol, it is 0 according to the resultant field tangential component of the electric field boundary condition of perfact conductor, i.e. metal surface, obtains scribbling thin-medium The improved Electric Field Integral Equation EFIE of the satellite structure target of layer, as follows

[E^inc(r)+E^sca(r)]_tan=0 (1)

Wherein, subscript tan represents the tangential component of electric field, E^incRepresent the electromagnetic wave being radiated in target, E^scaRepresent that target exists The end of a performance produced after electromagnetic wave irradiation, the expression-form of scattered field is:

$E^{\mathrm{sca}}=-{\mathrm{j\ω\μ}}_0\underset{S}{\∫}{J}_S\left(r^{\′}\right)\stackrel{\‾}{G}(r,r^{\′}){\mathrm{dS}}^{\′}-{\mathrm{j\ω\μ}}_0\underset{V}{\∫}{J}_{\mathrm{pol}}\left(r^{\′}\right)\stackrel{\‾}{G}(r,r^{\′}){\mathrm{dV}}^{\′}+\underset{V}{\∫}{M}_V\left(r^{\′}\right)\×\▿G(r,r^{\′}){\mathrm{dV}}^{\′}$

$=-{\mathrm{j\ω\μ}}_0\underset{S}{\∫}{J}_S\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{1}{{\mathrm{\ϵ}}_0}\underset{S}{\∫}{\mathrm{\ρ}}_s\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\′}$

(2)

$-{\mathrm{j\ω\μ}}_0\underset{V}{\∫}{J}_{\mathrm{pol}}\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{1}{{\mathrm{\ϵ}}_0}\underset{S}{\∫}{\mathrm{\ρ}}_{s,\mathrm{pol}}\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{1}{{\mathrm{\ϵ}}_0}\underset{S^{\mathrm{\Δ}}}{\∫}{\mathrm{\ρ}}_{s,\mathrm{pol}}\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\mathrm{\Δ}\′}$

$+\underset{V}{\∫}{M}_V\left(r^{\′}\right)\×\▿G(r,r^{\′}){\mathrm{dV}}^{\′}$

Wherein, V represents the elementary volume, volume element of thin-medium coating, and S represents the metal surface unit i.e. lower surface of thin-medium coating Unit, S^ΔRepresenting the upper surface unit of thin-medium coating, ω is the angular frequency of electromagnetic wave, μ₀And ε₀It is respectively the magnetic conductance in vacuum Rate and dielectric parameter, r and r' is respectively the position coordinates in field and source, and G (r, r') is the Green's function of free space, expression formula For:

$G(r,r^{\′})=\frac{e^{-\mathrm{jk}|r-r^{\′}|}}{4\mathrm{\π}|r-r^{\′}|}---\left(3\right)$

Wherein k is the wave number of free space；

Second step, it is close that 5 unknown quantitys in the improved Electric Field Integral Equation will set up in the first step are converted into metal surface electric current Degree J_S, 5 unknown quantitys are metal surface current density, J respectively_S, metal surface charge density ρ_S, dielectric polarized current density J_pol, upper and lower surface polarization charge densities ρ of dielectric layer_S,polWith body polarization magnetic current M_pol, its concrete representation is as follows,

Formula (1) ultimately becomes:

$[{\mathrm{j\ω\μ}}_0\underset{S}{\∫}{J}_s\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{1}{{\mathrm{j\ω\ϵ}}_0}\underset{S}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\▿G(r,r^{\′}){\mathrm{dS}}^{\′}$

$-{\mathrm{j\ω\μ}}_0\mathrm{\κ}\underset{V}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\hat{n}G(r,r^{\′}){\mathrm{dS}}^{\′}+\frac{\mathrm{\κ}}{{\mathrm{j\ω\ϵ}}_0}\underset{S}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\▿G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{\mathrm{\κ}}{{\mathrm{j\ω\ϵ}}_0}\underset{S^{\mathrm{\Δ}}}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\▿G(r,r^{\′}){\mathrm{dS}}^{\mathrm{\Δ}\′}$

$+{\mathrm{j\ω\μ}}_0({\mathrm{\μ}}_r-1)\underset{V}{\∫}{\hat{n}}^{\′}\×J_S\left(r^{\′}\right)\×\▿G(r,r^{\′}){\mathrm{dV}}^{\′}{]}_{\tan }=E_{\tan }^{\mathrm{inc}}---\left(5\right)$

3rd step, carries out moment method to equation (5) and solves, and obtains matrix equation, uses the quick multistage sub-technology of multilamellar to add Speed Matrix Solving；

4th step, when the distance between field source meets some requirements, only the angular spectrum on transfer field source line direction divides Amount, is preferably minimized transfer amount, and distance condition is:

r_mn＞ ＞ 3 γ D (6)

γ is controling parameter, and γ is the least, and the most non-NULL groups can realize with fast far field approximation, accordingly result precision Loss the biggest, remote block number is gradually increased along with γ and reduces；

5th step, solution matrix equation, obtain current coefficient, calculated electromagnetic scattering further according to reciprocal theorem by current coefficient Parameter.

In described second step, the method for unknown quantity conversion is as follows:

Interface can be obtained by the normal direction discontinuity of current of polarization

$(J_{\mathrm{pol},1}-J_{\mathrm{pol},2})\·\hat{n}=-{\mathrm{j\ω\ρ}}_S---\left(7\right)$

In object module exist two kinds of interfaces, the 1st kind: 1 be medium, and 2 is metal, i.e. medium lower surface, i.e. satellite and Coating interface, due to J in metal_pol,2=0, therefore can obtain on medium lower surface

$J_{\mathrm{pol},1}\·\hat{n}=-{\mathrm{j\ω\ρ}}_{S,\mathrm{pol}}---\left(8\right)$

2nd kind: 1 is air, and 2 is medium, i.e. medium upper surface, i.e. coating and air interface, also due in air J_pol,1=0, therefore can obtain on medium upper surface

$J_{\mathrm{pol},2}\·\hat{n}=-{\mathrm{j\ω\ρ}}_{S,\mathrm{pol}}---\left(9\right)$

In order to be got up with the electrical communication on satellite surface by current of polarization, Maxwell equation is utilized to obtain

$J_{\mathrm{pol}}=\mathrm{j\ω}(\mathrm{\ϵ}-{\mathrm{\ϵ}}_0)E---\left(10\right)$

At satellite surfaceTherefore formula (10) can be write as

$J_{\mathrm{pol}}={\mathrm{j\ω\κ\ρ}}_S\hat{n}---\left(11\right)$

Wherein $\mathrm{\κ}=\frac{\mathrm{\ϵ}-{\mathrm{\ϵ}}_0}{\mathrm{\ϵ}};$

By Current continuity equation

${\mathrm{\ρ}}_S=-\frac{1}{\mathrm{j\ω}}{\▿}_S\·J_S---\left(12\right)$

Substitution formula (11),

$J_{\mathrm{pol}}=-\mathrm{\κ}(\▿\·J_S)\hat{n}---\left(13\right)$

Convolution (8), (11) and formula (12), can obtain at medium lower surface

${\mathrm{\ρ}}_{S,\mathrm{pol}}=-\frac{\mathrm{\κ}}{\mathrm{j\ω}}\▿\·J_S---\left(14\right)$

In like manner can obtain at medium upper surface

${\mathrm{\ρ}}_{S,\mathrm{pol}}=\frac{\mathrm{\κ}}{\mathrm{j\ω}}\▿\·J_S---\left(15\right)$

Ignoring normal direction magnetic field, tangential magnetic field is a constant simultaneously, according to the magnetic field boundaries condition of perfact conductor obtain with Lower expression formula:

$M_{\mathrm{pol}}\left(r^{\′}\right)\≈\mathrm{j\ω}(\mathrm{\μ}-{\mathrm{\μ}}_0)H_t\left(r^{\′}\right)=\mathrm{j\ω}\frac{({\mathrm{\μ}}_r-1)}{{\mathrm{\μ}}_r}{B}_t\left(r^{\′}\right)---\left(16\right)$

μ_rFor the relative permeability of thin-medium, by B_t(r') do and deform:

$B_t\left(r^{\′}\right)=-{\hat{n}}^{\′}\×{\hat{n}}^{\′}\×B_t\left(r^{\′}\right)=-\mathrm{\μ}{\hat{n}}^{\′}\×{\hat{n}}^{\′}\×H_t\left(r^{\′}\right)=-\mathrm{\μ}{\hat{n}}^{\′}\×J_S\left(r^{\′}\right)---\left(17\right)$

It is updated in (16) obtain M by (17) formula_polWith the sensitive surface electric current J of metal surface_sThe expression formula represented is as follows:

$M_{\mathrm{pol}}\left(r^{\′}\right)=-{\mathrm{j\ω\μ}}_0({\mathrm{\μ}}_r-1){\hat{n}}^{\′}\×J_S\left(r^{\′}\right)---\left(18\right)$

(12), (13), (14), (15), (18) are substituted in (2) and obtain:

$[{\mathrm{j\ω\μ}}_0\underset{S}{\∫}{J}_s\left(r^{\′}\right)G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{1}{{\mathrm{j\ω\ϵ}}_0}\underset{S}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\▿G(r,r^{\′}){\mathrm{dS}}^{\′}$

$-{\mathrm{j\ω\μ}}_0\mathrm{\κ}\underset{V}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\hat{n}G(r,r^{\′}){\mathrm{dS}}^{\′}+\frac{\mathrm{\κ}}{{\mathrm{j\ω\ϵ}}_0}\underset{S}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\▿G(r,r^{\′}){\mathrm{dS}}^{\′}-\frac{\mathrm{\κ}}{{\mathrm{j\ω\ϵ}}_0}\underset{S^{\mathrm{\Δ}}}{\∫}{\▿}^{\′}\·J_s\left(r^{\′}\right)\▿G(r,r^{\′}){\mathrm{dS}}^{\mathrm{\Δ}\′}$

$+{\mathrm{j\ω\μ}}_0({\mathrm{\μ}}_r-1)\underset{V}{\∫}{\hat{n}}^{\′}\×J_S\left(r^{\′}\right)\×\▿G(r,r^{\′}){\mathrm{dV}}^{\′}{]}_{\tan }=E_{\tan }^{\mathrm{inc}}---\left(5\right)$

The method using multilamellar quick multistage sub-technology acceleration Matrix Solving in described 3rd step is as follows:

Step 3.1, sets up the tree structure of satellite；

Step 3.2, it is assumed that in given one layer, has site r in arbitrary group of m_i, arbitrary group of n has a source point r_j, r_m、r_nIt is site group and the central point of source point group respectively；In two non-NULL groups, the space vector of site and source point is designated as: r_ij= r_i-r_j=r_im+r_mn+r_nj, when the site group analyzed and source point group misaligned the most non-conterminous time, meet | r_im+r_nj|<|r_mn|, from Can be write as in angular spectrum space by the scalar Green's function in space:

$\frac{e^{-\mathrm{jk}|r_i-r_j|}}{|r_i-r_j|}=\frac{\mathrm{jk}}{4\mathrm{\&pi;}}\underset{S_E}{\&Integral;}{e}^{-\mathrm{jk}\&CenterDot;(r_{\mathrm{im}}+r_{\mathrm{nj}})}{\mathrm{\&alpha;}}_{\mathrm{mn}}(\hat{k}\&CenterDot;{\hat{r}}_{\mathrm{mn}})---\left(19\right)$

Whole integration is defined in unit ball S_EOn, α_mnIt is the transfer factor of two groups, is defined as follows:

${\mathrm{\&alpha;}}_{\mathrm{mn}}=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{(-j)}^l(2l+1)h_l^{\left(2\right)}\left(\mathrm{kr}\right)P_l(\hat{k}\&CenterDot;\hat{r})---\left(20\right)$

Wherein, j_l(kd) it is ball shellfish plug function,For Equations of The Second Kind ball Hankel function,For Legendre's letter Number L is the length of blocking of infinite series, L=kD+ β (kD)^1/3, D is the size of packet, and β is precision parameter, β >=2；With (θ, φ) Representation unit ball S_ECoordinate, integration is counted as K_L=2L², the one-dimensional Gauss integration of L point, φ direction is wherein gathered altogether in θ direction Gather 2L the trapezoidal rule integration counted altogether；

Multiple transfer factors, when the size that two far fields act on group switching centre distance and group determines, are constituted by step 3.3 Transfer matrix is plus a window function:

$w_l=\left\{\begin{array}{cc}1& l\&le;L/2\\ \sin \left(\frac{1}{L}\mathrm{\&pi;}\right)& l>L/2\end{array}\right.---\left(21\right)$

Being multiplied by window function in transfer factor makes the value of transition component near two group switching centre line directions relative to it The value of its part transfer factor has precipitous change, omits the effect of the angular spectrum component less on transfer process impact.

The method of the tree structure setting up satellite in described step 3.1 is as follows: satellite is put into one and can surround it Cubic space in, remember that these space coordinates are XYZ, this space is designated as the level 0 of tree structure, then to this space three Dimension is divided into two on direction respectively, forms 8 sub spaces, is designated as ground floor, continues segmentation on the basis of 8 sub spaces, when Sub-layers space be smaller in size than equal to set electric wave long time then stop segmentation；Whole satellite is transferred to a tree On, the space of each layer is carried out coding and sets up index, preserve the number in each layer of non-NULL group space.

The wavelength set in described step 3.1 is between 0.2 wavelength to 0.5 wavelength.

Described step 3.3 judging, the method for the angular spectrum component less on transfer process impact is as follows: set θ_rIt is two effects Group line direction and the angle of space z-axis, φ_rIt is that two effect group line direction vectors project and the angle of x-axis, θ on XOY face_e The angular spectrum component direction corresponding for the threshold value set and the angle in two effect group line directions；In the process calculating transfer factor In, the size of transfer factor modulus value increases along with θ and reduces, and extracts and calculate angular spectrum component direction and two effect group switching centre lines Angular separation is not more than θ_eAngular spectrum component.

The present invention compared with prior art, its remarkable advantage: (1) establishing equation is simple.For this metal medium of satellite Mixed structure, it is only necessary to metal part is set up improved Electric Field Integral Equation, media fraction need not set up equation and is described.(2) institute Resource need to be calculated few.Due to the construction features of satellite plagioclase, the high directivity of transfer factor component.And in ray-tracing procedure by In having only to transfer field source line side peripherad angular spectrum component, to owning on unit ball in multistage son quick compared to multilamellar Angular spectrum component shifts, and saves the calculating time of transfer factor and required stored memory.Simultaneously when meeting quick far field Under conditions of approximation, it is only necessary to calculate transfer factor angular spectrum component on field source group switching centre line direction, meeting computational accuracy Under conditions of save calculating resource further.

Accompanying drawing explanation

Fig. 1 is that metal of the present invention coats film dielectric layer object delineation, and a is geometric model schematic diagram, and b is electromagnetic parameter mould Type structural representation.

Fig. 2 is the discontinuity schematic diagram of the normal component of border of the present invention current of polarization.

Fig. 3 is that multilevel fast multipole far field of the present invention acts on schematic diagram.

Fig. 4 is structure and the coating schematic diagram of dummy satellite of the present invention

Fig. 5 is dual station RCS curve synoptic diagram of the present invention.

Detailed description of the invention

The numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band of the present invention, step is as follows:

The first step, the foundation of electromagnetic scattering integral equation, i.e. based on perfact conductor electric field boundary condition, in metal surface Resultant field tangential component be 0, and resultant field is incident electric fields and scattering electric field sum.Incident electric fields is known excitation, uniformly puts down Face ripple is usually used to as incident electric fields, and scattering electric field can represent by surface unknown current to be asked.

Second step, has 5 unknown quantitys in the improved Electric Field Integral Equation set up in the conversion between unknown quantity, the i.e. first step, It is metal surface electric current density, metal surface charge density respectively, dielectric polarized current density, the upper and lower surface pole of dielectric layer Change charge density.Using during moment method solves the equation, due to an only test item: metal surface electric current is close Degree, it is clear that be to solve the result that this equation is correct, it is therefore necessary to all of unknown source is all used metal surface electric current Density represents；

3rd step, in conjunction with the feature of satellite structure, utilizes the ultimate principle of ray tracing, reasonably selects transfer factor The part that middle angular spectrum component is big, gives up to fall by part little for angular spectrum component.

4th step, in conjunction with the feature of satellite structure, uses the method for fast far field approximation quickly to calculate far field in higher Effect.

5th step, matrix equation solves and the calculating of electromagnetic scattering parameter.

Below in conjunction with the accompanying drawings the present invention is described in further detail.

In conjunction with Fig. 1, the present invention analyzes the numerical value of the numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band based on a kind of Method, step is as follows:

The first step, makes uniform plane wave be irradiated on a satellite structure scribbling film dielectric layer, will on the surface of satellite Produce surface induction surface current J_SWith surface charge ρ_S, in thin-medium, produce body polarization electric current J_pol, body polarization magnetic current M_pol, polarization electricity Lotus ρ_S,pol, it is 0 according to the resultant field tangential component of the electric field boundary condition of perfact conductor, i.e. metal surface, obtains scribbling thin-medium The improved Electric Field Integral Equation (EFIE) of the satellite structure target of layer, as follows

[E^inc(r)+E^sca(r)]_tan=0 (1)

Wherein, subscript tan represents the tangential component of electric field, E^incRepresent the electromagnetic wave being radiated in target, E^scaRepresent that target exists The end of a performance produced after electromagnetic wave irradiation, and scattered field includes metal and the two-part effect of thin-medium, the expression-form of scattered field For:

$E^{\mathrm{sca}}=-{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_S\left(r^{\&prime;}\right)\stackrel{\&OverBar;}{G}(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-{\mathrm{j\&omega;\&mu;}}_0\underset{V}{\&Integral;}{J}_{\mathrm{pol}}\left(r^{\&prime;}\right)\stackrel{\&OverBar;}{G}(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}+\underset{V}{\&Integral;}{M}_V\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}$

$=-{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_S\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{\&epsiv;}}_0}\underset{S}{\&Integral;}{\mathrm{\&rho;}}_s\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}$

(2)

$-{\mathrm{j\&omega;\&mu;}}_0\underset{V}{\&Integral;}{J}_{\mathrm{pol}}\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{\&epsiv;}}_0}\underset{S}{\&Integral;}{\mathrm{\&rho;}}_{s,\mathrm{pol}}\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{\&epsiv;}}_0}\underset{S^{\mathrm{\&Delta;}}}{\&Integral;}{\mathrm{\&rho;}}_{s,\mathrm{pol}}\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\mathrm{\&Delta;}\&prime;}$

$+\underset{V}{\&Integral;}{M}_V\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}$

Wherein V represents the elementary volume, volume element of thin-medium coating, and S represents metal surface unit (the i.e. lower surface of thin-medium coating Unit), S^ΔRepresenting the upper surface unit of thin-medium coating, ω is the angular frequency of electromagnetic wave, μ₀And ε₀It is respectively the magnetic conductance in vacuum Rate and dielectric parameter, r and r' is respectively the position coordinates in field and source, and G (r, r') is the Green's function of free space, expression formula For:

$G(r,r^{\&prime;})=\frac{e^{-\mathrm{jk}|r-r^{\&prime;}|}}{4\mathrm{\&pi;}|r-r^{\&prime;}|}---\left(3\right)$

Wherein k is the wave number of free space.

Second step, has 5 unknown quantitys in the improved Electric Field Integral Equation set up, be that metal surface electric current is close respectively in the first step Degree J_S, metal surface charge density ρ_S, dielectric polarized current density J_pol, the upper and lower surface polarization charge densities of dielectric layer ρ_S,pol, body polarization magnetic current M_polIf using moment method to solve the equation, an only test item: metal surface electric current density J_S, it is clear that it is to solve the result that this equation is correct, it is therefore necessary to all of unknown source is all used for metal surface electricity Current density represents.

The discontinuity of the normal component of border current of polarization is as shown in Figure 2.At satellite with coating interface, coating with empty Can be obtained by the normal direction discontinuity of current of polarization on gas interface

$(J_{\mathrm{pol},1}-J_{\mathrm{pol},2})\&CenterDot;\hat{n}=-{\mathrm{j\&omega;\&rho;}}_S---\left(4\right)$

There are two kinds of interfaces in the present invention in object module, the 1st kind: 1 is medium, and 2 is metal, i.e. medium lower surface, Due to J in metal_pol,2=0, therefore can obtain on medium lower surface

$J_{\mathrm{pol},1}\&CenterDot;\hat{n}=-{\mathrm{j\&omega;\&rho;}}_{S,\mathrm{pol}}---\left(5\right)$

2nd kind: 1 is air, and 2 is medium, i.e. medium upper surface, also due to J in air_pol,1=0, therefore on medium Can obtain on surface

$J_{\mathrm{pol},2}\&CenterDot;\hat{n}=-{\mathrm{j\&omega;\&rho;}}_{S,\mathrm{pol}}---\left(6\right)$

In order to be got up with the electrical communication on satellite surface by current of polarization, Maxwell equation is utilized to obtain

J_pol=j ω (ε-ε₀)E (7)

At satellite surfaceTherefore formula (7) can be write as

$J_{\mathrm{pol}}=-{\mathrm{j\&omega;\&kappa;\&rho;}}_S\hat{n}---\left(8\right)$

Wherein $\mathrm{\&kappa;}=\frac{\mathrm{\&epsiv;}-{\mathrm{\&epsiv;}}_0}{\mathrm{\&epsiv;}}.$

By Current continuity equation

${\mathrm{\&rho;}}_S=-\frac{1}{\mathrm{j\&omega;}}{\&dtri;}_S\&CenterDot;J_S---\left(9\right)$

Substitution formula (8),

$J_{\mathrm{pol}}=-\mathrm{\&kappa;}(\&dtri;\&CenterDot;J_S)\hat{n}---\left(10\right)$

The charge density on metal and the polarized current density in dielectric has been can be seen that the most by formula (9) and formula (10) The form being transformed to metal covering electric current density represents；

Convolution (5), (8) and formula (9), can obtain at medium lower surface

${\mathrm{\&rho;}}_{S,\mathrm{pol}}=-\frac{\mathrm{\&kappa;}}{\mathrm{j\&omega;}}\&dtri;\&CenterDot;J_S---\left(11\right)$

In like manner can obtain with air interface in the coating of medium upper surface i.e. satellite

${\mathrm{\&rho;}}_{S,\mathrm{pol}}=\frac{\mathrm{\&kappa;}}{\mathrm{j\&omega;}}\&dtri;\&CenterDot;J_S---\left(12\right)$

Considering that normal direction magnetic field is to ignore, tangential magnetic field may be considered a constant simultaneously, so leading according to ideal The magnetic field boundaries condition of body can obtain following formula:

$M_{\mathrm{pol}}\left(r^{\&prime;}\right)\&ap;\mathrm{j\&omega;}(\mathrm{\&mu;}-{\mathrm{\&mu;}}_0)H_t\left(r^{\&prime;}\right)=\mathrm{j\&omega;}\frac{({\mathrm{\&mu;}}_r-1)}{{\mathrm{\&mu;}}_r}{B}_t\left(r^{\&prime;}\right)---\left(13\right)$

μ_rRelative permeability for thin-medium.By B_t(r') simple deformation is done:

$B_t\left(r^{\&prime;}\right)=-{\hat{n}}^{\&prime;}\&times;{\hat{n}}^{\&prime;}\&times;B_t\left(r^{\&prime;}\right)=-\mathrm{\&mu;}{\hat{n}}^{\&prime;}\&times;{\hat{n}}^{\&prime;}\&times;H_t\left(r^{\&prime;}\right)=-\mathrm{\&mu;}{\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)$

(14) formula is updated in (13) to obtain M_polWith the sensitive surface electric current J of metal surface_sThe expression formula represented is such as Under:

$M_{\mathrm{pol}}\left(r^{\&prime;}\right)=-{\mathrm{j\&omega;\&mu;}}_0({\mathrm{\&mu;}}_r-1){\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)---\left(15\right)$

(9), (10), (111), (12), (15) are substituted in (2) and can obtain:

$[{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_s\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}$

$-{\mathrm{j\&omega;\&mu;}}_0\mathrm{\&kappa;}\underset{V}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\hat{n}G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}+\frac{\mathrm{\&kappa;}}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{\mathrm{\&kappa;}}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S^{\mathrm{\&Delta;}}}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\mathrm{\&Delta;}\&prime;}$

$+{\mathrm{j\&omega;\&mu;}}_0({\mathrm{\&mu;}}_r-1)\underset{V}{\&Integral;}{\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}{]}_{\tan }=E_{\tan }^{\mathrm{inc}}---\left(16\right)$

3rd step, carries out moment method to equation (16) and solves, obtain matrix equation.To the Green's function of free space according to Addition theorem launch, utilize the ultimate principle of ray tracing, the characteristic of coalition Line Integral equation reasonably select transfer because of The part that in son, angular spectrum component is big, gives up to fall by part little for angular spectrum component.

Millimere-wave band satellite belongs to Electrically large size object, and the quick multistage sub-technology of multilamellar is used to accelerate Matrix Solving.Many The quick multistage sub-technology of layer is firstly the need of setting up tree structure: satellite is put into a cubic space that can surround it In, remember that these space coordinates are XYZ, this space is designated as the level 0 of tree structure, then to difference on this space three-dimensional direction It is divided into two, forms 8 sub spaces, be designated as ground floor, think and analogize segmentation, until being smaller in size than equal to of sub-layers space Fixed electric wave length (0.2 wavelength to 0.5 wavelength between) then stops segmentation.Whole satellite is just transferred to a tree On, the space of each layer being carried out coding and sets up index, preserve the number of each layer of non-NULL group (space), all of calculating all exists Non-NULL group is carried out.

Assume in given one layer, group m has site r_iA source point r is had with in group n_j, r_m、r_nIt is field respectively Point group and the central point of source point group.In two non-NULL groups, the space vector of site and source point is designated as: r_ij=r_i-r_j=r_im+r_mn+r_nj, When the site group analyzed and source point group misaligned the most non-conterminous time, meet | r_im+r_nj|<|r_mn|, the scalar Green of free space Function can be write as in angular spectrum space:

$\frac{e^{-\mathrm{jk}|r_i-r_j|}}{|r_i-r_j|}=\frac{\mathrm{jk}}{4\mathrm{\&pi;}}\underset{S_E}{\&Integral;}{e}^{-\mathrm{jk}\&CenterDot;(r_{\mathrm{im}}+r_{\mathrm{nj}})}{\mathrm{\&alpha;}}_{\mathrm{mn}}(\hat{k}\&CenterDot;{\hat{r}}_{\mathrm{mn}})---\left(17\right)$

Whole integration is defined in unit ball S_EOn, α_mnIt is the transfer factor of two groups, is defined as follows:

${\mathrm{\&alpha;}}_{\mathrm{mn}}=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{(-j)}^l(2l+1)h_l^{\left(2\right)}{\left(\mathrm{kr}\right)P}_l(\hat{k}\&CenterDot;\hat{r})---\left(18\right)$

Wherein, j_l(kd) it is ball shellfish plug function,For Equations of The Second Kind ball Hankel function,For Legendre's letter Number.L is the length of blocking of infinite series, L=kD+ β (kD)^1/3, D is the size of packet, and β is empirical parameter, takes 2 in the present invention. With (θ, φ) representation unit spherical coordinates, integration is counted as K_L=2L², the one-dimensional Gauss integration of L point is wherein gathered altogether in θ direction, φ direction gathers 2L the trapezoidal rule integration counted altogether.

When the size that two far fields act on group switching centre distance and group determines, the characteristic of transfer factor is similar to linear antenna The antenna pattern characteristic of array, by transfer matrix plus a window function:

$w_l=\left\{\begin{array}{cc}1& l\&le;L/2\\ \sin \left(\frac{1}{L}\mathrm{\&pi;}\right)& l>L/2\end{array}\right.---\left(19\right)$

Wherein, L is the number of L the discrete point in θ direction under spherical coordinates, and l is the numbering of each transfer factor component.Concrete Way is: be multiplied by the effect of window function in transfer factor so that the value phase of the transition component near two group switching centre line directions Value for other parts transfer factor has precipitous change, in such a case, it is possible to omit transfer process impact less The effect of angular spectrum component so that ray multipole obtains higher efficiency.Meanwhile, can verify that this process is to calculating Precision impact is the least.Then need to set the scope of the angular spectrum component needing transfer.As it is shown on figure 3, θ_rIt is two effect group lines Direction and the angle of z-axis, φ_rIt is that two effect group line direction vectors project and the angle of x-axis, θ on XOY face_eFor the threshold set The angular spectrum component direction of value correspondence and the angle in two effect group line directions.During calculating transfer factor, transfer factor The size of angular spectrum component reduces along with the increase with the angle in effect group switching centre line direction, far field.So can to transfer because of Son sets threshold value and judges the size of conical area.In concrete operations, only need to extract and calculate angular spectrum component direction with Two effect group switching centre line angular separations are not more than θ_eAngular spectrum component.This segment angle spectral component can be obtained when proof of algorithm The effect of two far field effect groups can be contacted well.

4th step, the characteristic of coalition Line Integral equation, use the method for fast far field approximation quickly to calculate in higher Far field acts on.In fast far field approximation method, when the distance between field source meets some requirements, we can only calculate field Angular spectrum component on line direction, source, is so preferably minimized transfer amount.The present invention use the most remote between employing group and group The criterion of field approximation is:

r_mn＞ 3 γ D (20)

Wherein D is the electric size of this layer of group, and γ is the empirical value that coalition Line Integral equation feature is arranged, the present invention In be taken as 1.5.γ is the least, and the most non-NULL groups can realize with fast far field approximation, and the loss of accordingly result precision is the most just Can be the biggest, remote block number is gradually increased along with γ and reduces.Thus can be by quick to fast far field approximation, ray tracing and multilamellar Multistage submethod combines, and each layer of effect between group and group can be used three kinds quickly based on the distance between the group of far field Method processes: 1 Fast Multiple Method accelerates far field effect；2 ray multipoles accelerate far field effect；3 fast far field approximation add The effect of speed far field.

5th step, solution matrix equation, obtain current coefficient, calculated electromagnetic scattering further according to reciprocal theorem by current coefficient Parameter.

In order to verify correctness and the effectiveness of the inventive method, shown below is Aerospace Satellite thunder in calculating millimere-wave band Reach the example of scattering section, and result of calculation all compares with traditional algorithm, coincide very well.

The combinative structure of Aerospace Satellite is made up of two parts, and centre has a length of side to be 2m cube celestial body, has on xoy The solar panel of two a length of 10.87m.Celestial body up and down before and after four sides coated media be (3.5 ,-0.008).The sun The coated media of energy cell panel is (11.9,0.0).Incidence wave frequency is 35GHz, θ=0 ° → 180 °, direction of incidence wave, φ= 0°.Subdivision a size of 0.14 wavelength, unknown quantity number is 125160048, sub-layers packet size: 0.2 wavelength, and the number of plies is 12.From After Fig. 4 is it can be seen that add fast far field approximation and ray multipole, compared with traditional method, RCS result is the most identical, full Foot required precision.And in Table 1 it can be seen that after adding fast far field approximation and ray multipole, be to use at LC=8(LC The number of plies of fast far field approximation) time, the calculating time of program saves 49.6% than original program；Internal memory saves 52.8%. Therefore fast far field approximation and ray multipole can save calculating time and internal memory.

Table 1 calculates internal memory and time statistical table

Claims (4)

1. the numerical method of satellite Electromagnetic Scattering Characteristics in a millimere-wave band, it is characterised in that step is as follows:

The first step, makes uniform plane wave be irradiated on a satellite structure scribbling film dielectric layer, will produce on the surface of satellite The faradic current J of metal surface_S, metal surface charge density ρ_S, produce body polarization electric current J in thin-medium_pol, body polarization magnetic current M_polUpper and lower surface polarization charge densities ρ with dielectric layer_S,pol, according to the electric field boundary condition of perfact conductor, i.e. metal surface Resultant field tangential component is 0, obtains scribbling the improved Electric Field Integral Equation EFIE of the satellite structure target of film dielectric layer, as follows

[E^inc(r)+E^sca(r)]_tan=0 (1)

Wherein, subscript tan represents the tangential component of electric field, E^incRepresent the electromagnetic wave being radiated in target, E^scaRepresent that target exists The scattered field produced after electromagnetic wave irradiation, the expression-form of scattered field is:

$\begin{array}{c}{E}^{sca}=-{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_S\left(r^{\&prime;}\right)\stackrel{\&OverBar;}{G}(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-{\mathrm{j\&omega;\&mu;}}_0\underset{V}{\&Integral;}{J}_{pol}\left(r^{\&prime;}\right)\stackrel{\&OverBar;}{G}(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}+\underset{V}{\&Integral;}{M}_V\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}\\ =-{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_S\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{\&epsiv;}}_0}\underset{S}{\&Integral;}{\mathrm{\&rho;}}_s\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}\\ -{\mathrm{j\&omega;\&mu;}}_0\underset{V}{\&Integral;}{J}_{pol}\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{\&epsiv;}}_0}\underset{S}{\&Integral;}{\mathrm{\&rho;}}_{s,pol}\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{\&epsiv;}}_0}\underset{S^{\mathrm{\&Delta;}}}{\&Integral;}{\mathrm{\&rho;}}_{s,pol}\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{{\mathrm{\&Delta;}}^{\&prime;}}\\ +\underset{V}{\&Integral;}{M}_V\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}\end{array}---\left(2\right)$

Wherein, V represents the elementary volume, volume element of film dielectric layer, and S represents the lower surface unit of the i.e. film dielectric layer of metal surface unit, S^ΔTable Showing the upper surface unit of film dielectric layer, ω is the angular frequency of electromagnetic wave, μ₀And ε₀It is respectively the pcrmeability in vacuum and dielectric ginseng Number, r and r' is respectively the position coordinates in field and source, and G (r, r') is the Green's function of free space, and expression formula is:

$G(r,r^{\&prime;})=\frac{e^{-jk|r-r^{\&prime;}|}}{4\mathrm{\&pi;}|r-r^{\&prime;}|}---\left(3\right)$

Wherein k is the wave number of free space；

Second step, 5 unknown quantitys in the improved Electric Field Integral Equation will set up in the first step are converted into the faradic current of metal surface J_S, 5 unknown quantitys are the faradic current J of metal surface respectively_S, metal surface charge density ρ_S, produce body polarization electricity in thin-medium Stream J_pol, upper and lower surface polarization charge densities ρ of dielectric layer_S,polWith body polarization magnetic current M_pol, its concrete representation is as follows,

Formula (1) ultimately becomes:

$\begin{array}{c}\&lsqb;{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_s\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}\\ -{\mathrm{j\&omega;\&mu;}}_0\mathrm{\&kappa;}\underset{V}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\hat{n}G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}+\frac{\mathrm{\&kappa;}}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{\mathrm{\&kappa;}}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S^{\mathrm{\&Delta;}}}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{{\mathrm{\&Delta;}}^{\&prime;}}\\ +{\mathrm{j\&omega;\&mu;}}_0({\mathrm{\&mu;}}_r-1)\underset{V}{\&Integral;}{\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}{\&rsqb;}_{\tan }=E_{\tan }^{inc}\end{array}---\left(5\right)$

3rd step, carries out moment method to equation (5) and solves, and obtains matrix equation, uses the quick multistage sub-technology of multilamellar to accelerate square Battle array solves；

4th step, when the distance between field source meets some requirements, the only angular spectrum component on transfer field source line direction, will Transfer amount is preferably minimized, and distance condition is:

r_mn>>3γD (6)

γ is controling parameter, and γ is the least, and the most non-NULL groups can realize with fast far field approximation, the damage of accordingly result precision Losing the biggest, remote block number is gradually increased along with γ and reduces；

5th step, solution matrix equation, obtain current coefficient, calculated electromagnetic scattering ginseng further according to reciprocal theorem by current coefficient Amount；

In described second step, the method for unknown quantity conversion is as follows:

Interface can be obtained by the normal direction discontinuity of current of polarization

$(J_{pol,1}-J_{pol,2})\&CenterDot;\hat{n}=-{\mathrm{j\&omega;\&rho;}}_S---\left(7\right)$

There are two kinds of interfaces in object module, the 1st kind: 1 is medium, and 2 is metal, i.e. medium lower surface, i.e. satellite and thin Jie Matter layer interface, due to J in metal_pol,2=0, therefore can obtain on medium lower surface

$J_{pol,1}\&CenterDot;\hat{n}=-{\mathrm{j\&omega;\&rho;}}_{S,pol}---\left(8\right)$

2nd kind: 1 is air, and 2 is medium, i.e. medium upper surface, i.e. film dielectric layer and air interface, also due in air J_pol,1=0, therefore can obtain on medium upper surface

$J_{pol,2}\&CenterDot;\hat{n}={\mathrm{j\&omega;\&rho;}}_{S,pol}---\left(9\right)$

In order to be got up with the electrical communication on satellite surface by current of polarization, Maxwell equation is utilized to obtain

J_pol=j ω (ε-ε₀)E (10)

At satellite surfaceTherefore formula (10) can be write as

$J_{pol}={\mathrm{j\&omega;\&kappa;\&rho;}}_S\hat{n}---\left(11\right)$

Wherein

By Current continuity equation

${\mathrm{\&rho;}}_S=-\frac{1}{j\mathrm{\&omega;}}{\&dtri;}_S\&CenterDot;J_S---\left(12\right)$

Substitution formula (11),

$J_{pol}=-\mathrm{\&kappa;}(\&dtri;\&CenterDot;J_S)\hat{n}---\left(13\right)$

Convolution (8), (11) and formula (12), can obtain at medium lower surface

${\mathrm{\&rho;}}_{S,pol}=-\frac{\mathrm{\&kappa;}}{j\mathrm{\&omega;}}\&dtri;\&CenterDot;J_S---\left(14\right)$

In like manner can obtain at medium upper surface

${\mathrm{\&rho;}}_{S,pol}=\frac{\mathrm{\&kappa;}}{j\mathrm{\&omega;}}\&dtri;\&CenterDot;J_S---\left(15\right)$

Ignoring normal direction magnetic field, tangential magnetic field is a constant simultaneously, obtains according to the magnetic field boundaries condition of perfact conductor with following table Reach formula:

$M_{pol}\left(r^{\&prime;}\right)\&ap;j\mathrm{\&omega;}(\mathrm{\&mu;}-{\mathrm{\&mu;}}_0)H_t\left(r^{\&prime;}\right)=j\mathrm{\&omega;}\frac{({\mathrm{\&mu;}}_r-1)}{{\mathrm{\&mu;}}_r}{B}_t\left(r^{\&prime;}\right)---\left(16\right)$

μ_rFor the relative permeability of thin-medium, by B_t(r') do and deform:

$B_t\left(r^{\&prime;}\right)=-{\hat{n}}^{\&prime;}\&times;{\hat{n}}^{\&prime;}\&times;B_t\left(r^{\&prime;}\right)=-\mathrm{\&mu;}{\hat{n}}^{\&prime;}\&times;{\hat{n}}^{\&prime;}\&times;H_t\left(r^{\&prime;}\right)=-\mathrm{\&mu;}{\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)---\left(17\right)$

It is updated in (16) obtain M by (17) formula_polWith the faradic current J of metal surface_sThe expression formula represented is as follows:

$M_{pol}\left(r^{\&prime;}\right)=-{\mathrm{j\&omega;\&mu;}}_0({\mathrm{\&mu;}}_r-1){\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)---\left(18\right)$

(12), (13), (14), (15), (18) are substituted in (2) and obtain:

$\begin{array}{c}\&lsqb;{\mathrm{j\&omega;\&mu;}}_0\underset{S}{\&Integral;}{J}_s\left(r^{\&prime;}\right)G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{1}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}\\ -{\mathrm{j\&omega;\&mu;}}_0\mathrm{\&kappa;}\underset{V}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\hat{n}G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}+\frac{\mathrm{\&kappa;}}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{\&prime;}-\frac{\mathrm{\&kappa;}}{{\mathrm{j\&omega;\&epsiv;}}_0}\underset{S^{\mathrm{\&Delta;}}}{\&Integral;}{\&dtri;}^{\&prime;}\&CenterDot;J_s\left(r^{\&prime;}\right)\&dtri;G(r,r^{\&prime;}){\mathrm{dS}}^{{\mathrm{\&Delta;}}^{\&prime;}}\\ +{\mathrm{j\&omega;\&mu;}}_0({\mathrm{\&mu;}}_r-1)\underset{V}{\&Integral;}{\hat{n}}^{\&prime;}\&times;J_S\left(r^{\&prime;}\right)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{dV}}^{\&prime;}{\&rsqb;}_{\tan }=E_{\tan }^{inc}\end{array}---\left(5\right);$

The method using multilamellar quick multistage sub-technology acceleration Matrix Solving in described 3rd step is as follows:

Step 3.1, sets up the tree structure of satellite；

Step 3.2, it is assumed that in given one layer, has site r in arbitrary group of m_i, arbitrary group of n has a source point r_j, r_m、 r_nIt is site group and the central point of source point group respectively；In two non-NULL groups, the space vector of site and source point is designated as: r_ij=r_i-r_j =r_im+r_mn+r_nj, when the site group analyzed and source point group misaligned the most non-conterminous time, meet | r_im+r_nj|<|r_mn|, the most empty Between scalar Green's function can be write as in angular spectrum space:

$\frac{e^{-jk|r_i-r_j|}}{|r_i-r_j|}=\frac{jk}{4\mathrm{\&pi;}}\underset{S_E}{\&Integral;}{e}^{-jk\&CenterDot;(r_{im}+r_{nj})}{\mathrm{\&alpha;}}_{mn}(\hat{k}\&CenterDot;{\hat{r}}_{mn})---\left(19\right)$

Whole integration is defined in unit ball S_EOn, α_mnIt is the transfer factor of two groups, is defined as follows:

${\mathrm{\&alpha;}}_{mn}=\underset{l=0}{\overset{L}{\&Sigma;}}{(-j)}^l(2l+1)h_l^{\left(2\right)}\left(kr\right)P_l(\hat{k}\&CenterDot;\hat{r})---\left(20\right)$

Wherein,For Equations of The Second Kind ball Hankel function,For Legendre function；L be infinite series block length Degree, L=kD+ β (kD)^1/3, D is the size of packet, and β is precision parameter, β >=2；With (θ, φ) representation unit ball S_ECoordinate, integration Count as K_L=2L², the one-dimensional Gauss integration of L point is wherein gathered altogether in θ direction, it is trapezoidal that φ direction gathers that 2L counts altogether Rule integration；

Step 3.3, when the size that two far fields act on group switching centre distance and group determines, the transfer that multiple transfer factors are constituted Matrix is plus a window function:

$w_l=\left\{\begin{array}{cc}1& l\&le;L/2\\ sin\left(\frac{l}{L}\mathrm{\&pi;}\right)& l>L/2\end{array}\right.---\left(21\right)$

Being multiplied by window function in transfer factor makes the value of transition component near two group switching centre line directions relative to other portion The value dividing transfer factor has precipitous change, omits the effect of the angular spectrum component less on transfer process impact.

The numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band the most according to claim 1, it is characterised in that described The method of the tree structure setting up satellite in step 3.1 is as follows: satellite is put into a cubic space that can surround it In, remember that these space coordinates are XYZ, this space is designated as the level 0 of tree structure, then to difference on this space three-dimensional direction It is divided into two, forms 8 sub spaces, be designated as ground floor, on the basis of 8 sub spaces, continue segmentation, when sub-layers space Segmentation is then stopped when being smaller in size than long equal to the electric wave set；Whole satellite is transferred on a tree, to each layer Space carries out coding and sets up index, preserves the number in each layer of non-NULL group space.

The numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band the most according to claim 2, it is characterised in that: described The electric wave length set in step 3.1 is between 0.2 wavelength to 0.5 wavelength.

The numerical method of satellite Electromagnetic Scattering Characteristics in millimere-wave band the most according to claim 2, it is characterised in that described Step 3.3 judging, the method for the angular spectrum component less on transfer process impact is as follows: set θ_rBe two effect group line directions with The angle of space z-axis, φ_rIt is that two effect group line direction vectors project and the angle of x-axis, θ on XOY face_eFor the threshold value set Corresponding angular spectrum component direction and the angle in two effect group line directions；During calculating transfer factor, transfer factor mould The size of value increases along with θ and reduces, and extracts and calculates angular spectrum component direction and two effect group switching centre line angular separations are little In θ_eAngular spectrum component.