基于单积分方程矩量法的海洋表面微波散射模拟

基于单积分方程矩量法的海洋表面微波散射模拟夏明耀!伍振兴"北京大学电子学系!北京#$$%&#’摘要(本文使用单积分方程矩量法!对风场驱动下的海洋表面微波散射进行数值模拟)海洋随机粗糙面用三角形贴片模型来拟合!核心算法是基于快速傅立叶变换技术的多层稀疏矩阵规范网格法)针对海水的高导电性!引入了局部化近似方法)数值结果验证了算法的正确性和计算效率!并对*波段的微波散射进行了计算分析)关键词(海洋表面+微波散射+单积分方程+矩量法+并行计算中图分类号(,-#.文献标识码(/文章编号($0&121##1"1$$.’$02$0%.2$3456789:5;<;=>5?@;A 9B C4?9::C @5<D=@;6E ?C 9<47@=9?C F G H:I C>C :I ;J ;=>;6C <:F K 9F C J;<45<D 8CL <:C D @98M N 79:5;<O P /-Q R S 2T U V !WX Y Z [R 2\Q R S"]^_‘a b c ^d b e fg h ^i b a e d j i k !l ^m j d no d j p ^a k j b q !r ^j s j d n #$$%&#!t u j d ‘’v G F :@9?:(/w U x y U z z {V U |Zw V {U R U }T x Q xV w ~Q |{V !U "[x |U y y [{Q R Sw {V ~!Q R #2#{Q "[RV |[U Rx ${w U |[xQ x #["[}V z [#$x Q R Sy Z [x Q R S }[Q R y [S {U }[%$U y Q V R !Q y Zy Z [~[y Z V #V w ~V ~[R y x ),{Q U R S $}U {z U y |Z [x U {[[~z }V T [#y Vy Q S Z y }T~V #[}y Z [{U R #V ~{V $S Zx ${w U |[x ),Z [$R #[{}T Q R SU }S V {Q y Z ~Q x y Z [&&,2’U x [#~$}y Q }["[}x z U {x [2~U y {Q \|U R V R Q |U }2S {Q #~[y Z V #)/}V |U }Q ([#U z z {V \Q ~U y Q V R~[y Z V #Q xQ R y {V #$|[#w V {y Z [Z Q S Z }T|V R #$|y Q "[x [U !U y [{~[#Q $~))$~[{Q |U }[\U ~z }[xU {[z {V "Q #[#y V #[~V R x y {U y [y Z [~[{Q y x V w y Z [z {V z V x [#x |Z [~[)*Q ~$}U y Q V R{[x $}y x U {[z {[x [R y [#w V {V |[U R 2}Q +[x ${w U |[x Q }}$~Q R U y [#’T~Q |{V !U "[’[U ~x U y *2’U R #),C HA ;@J F (V |[U Rx ${w U |[x +~Q |{V !U "[x |U y y [{Q R S +x Q R S }[Q R y [S {U }[%$U y Q V R +~V ~[R y ~[y Z V #x +z U {U }}[}|V ~z $y Q R S-引言随机粗糙表面电磁散射有非常广泛的应用背景!包括海面风场和土壤湿度的微波遥感.物体表面光洁度检测等!所以一直受到众多学者的关注)文献/#0评论了这一领域的最新进展)随着高性能计算机的发展!数值模拟分析近年来获得了很大进步/10)过去一直使用的近似理论!如小扰动模型和基尔霍夫模型!都有严格的适用范围限制!精度也不能满足越来越高的定量化要求)因此!严格的数值分析方法受到高度重视)对一般应用来说!为了充分获得随机粗糙表面的统计信息!被模拟的表面尺寸应该足够大!通常需要几千个平方波长)将表面积分方程数值离散化以后!得到一个几十万甚至几百万个未知量的方程组)求解这样一个巨型线性方程组!既要研究先进的数值算法!又必须借助高性能并行计算机)由于矩量法导致一个稠密矩阵!对于求解一个数百万个未知量的方程组来说!瓶颈问题包括两方面(内存限制和12X 求解时间)到目前为止!最为成功的求解方法是稀疏矩阵规范网格法"*-13!*z U {x [2-U y {Q \1U R V R Q |U }23{Q #-[y Z V #x’/0!30)这种方法把矩阵分成强相互作用和弱相互作用两部分!强相互作用项使用精确数值积分!并把矩阵元保存在内存中!而弱相互作用借助于格林函数的泰勒展开!表示为具有平移不变性的,V [z }Q y (结构!以便利用快速傅立叶变换技术)由于格林函数在高损耗媒质中迅速衰减!文献/.0提出了一种双网格法!在界面两边使用不同的网格划分!在保证计算精度的情况下提高了计算效率)最近!作者提出了一种基于单个积分方程的多层*-13方法/40!使用了矢量三角形5W 3基函数)采用单积分方程可以将未知量数目减少一半!而使用5W 3三角形贴片可以更精确地拟合粗糙表面)使用这种方法!我们已成功地求解了超过#1$$万个未知量的方程组/&0)除了各种版本的*-13方法外!人们还尝试了其他方法!包括快速多极子方法/%0和时域有限差分法/60等!这些方法都不如*-13方法有效)阻抗近似方法/#$0是专门用来求解高导电媒质粗糙面散射的!因此对海水的电磁散射分析可能特别有效)本文是在文献/40的基础上!进一步优化算法!并针对海洋特性表面!发展了一些近似方法!以加快分析速度)数值算例以风场驱动下的海洋表面为模型!对*波段的微波散射进行了模拟!证明了本算法的可靠性和计算效率)收稿日期(1$$02#$2#4+修回日期(1$$32$.21$基金项目(国家自然科学基金")V)4$#&#$$#’第0期1$$.年0月电子学报/1,/7*71,58)P 1/*P )P 1/9V })00)V )0-U {)1$$.!理论!"#单积分方程矩量法为了采用单磁场积分方程方法$%%&’假定在下表面存在一个有效电流分布()**’以便把下半空间中的电场和磁场表示为+时间因子为),-./123-.456782’()**9:,+%;0<236=>?82’()**9:,+%@0其中内积是在下表面进行积分’782是下半空间的并矢格林函数’即7823A B C %D 22E F==)-D 2GH I G +20当场点从下半空间趋向界面时’可以定义一对等效表面电流和磁流J(K 3L M ><236L M >=>782’()**9:,+N 0,(O3L M >123-.456L M >782’()**9+H 0然后’上半空间区域中的散射场可以用等效表面电流和磁流表示成<P3-.Q 5678%’(O 9:C C6=>78%’(K 9:C +R 0其中的内积是在上表面进行积分’78%是上半空间的并矢格林函数’与式+20相同但以D %取代子D 2S 上半空间中的总场是<%3<-T U C <P ’其中<-T U是入射波磁场S 边界条件要求磁场切向分量连续’即L M ><23L M ><%’所以’(K 3(-T U C -.Q 56L M >78%’(O 9:C C6L M >=>78%’(K 9:C +V 0其中(-T U 3L M ><-T U 为入射波等效表面电流S 使用矩量法和W XY 基函数$%2&’可以把式+N 0Z +H 0和+V 0转化为线性方程组JA K 3+[2,\C 0]A )**’,A O 3^5_2]A )**+‘0A K 3A -T U C%^5_%]A O C+[%C \,0]A K+a 0其中\C 和\,是两个对角矩阵’元素为b C c c 3d C c e 2I 及b ,c c 3d ,ce 2I ’这里d C c 和d ,c 满足d C c C d ,c 32I’分别是与第c 条边相连的两个W XY 三角形贴片之间的内夹角和外夹角S ^5345e Q f 5是自由空间波阻抗’列矢量A K ’A O ’A )**’A -T U’分别是它们所对应的矢量在W XY 基函数下的展开系数’例如(-T U g h ic 3%A-T U c j c +k 0’A -T Uc 3%l cml c+L M >l M c0n (-T U o l +p 0其中q c +c 3%’2’r’i 0是i 个内部棱边的W XY 矢量基函数’棱边的方向以l M c 标记’+L M >l M c 0的方向与流过棱边的电流方向定义是一致的S 矩阵[和_的元素为s O +c ’t 03%l cm lc mu t+l Mc n =?O >j t0o :v o l ’O3%’2+%50w O +c ’t03c D %l cm lc m u t$+l Mcn j t0?OC%D2O x ?Oxl c +=v P n j t 0&o :v o l ’O3%’2+%%0其中?O 3)-DO Ge H I G ’O 3%’2分别是上半空间和下半空间媒质的标量格林函数S 将式+‘0代入式+a 0’得到$_%]_2,+\C ,[%0+\C ,[20&]A )**3A -T U+%20这就是求解下表面有效电流分布的单磁场积分方程+y z{|}’y -T ~!)z "~T )#-U{-)!$|T #)~%"!}&’"#-(T0矩量法公式S 其他形式的单积分方程公式可以类似地导出’可以使用电场切向连续或磁场切向连续的边界条件’未知量可以是上表面或下表面的有效电流或有效磁流S 本文只使用式+%20的单磁场积分方程矩量法公式S有效电流分布的展开系数A )**一旦求出’可以根据式+‘0获得等效电流和磁流分布的展开系数A K 和A O’散射场可由式+R 0获得’也可以不计算散射场’直接计算双站散射系数S !"!多层稀疏矩阵规范网格法为了使用多层稀疏矩阵规范网格法+zy z )Y ’z ’!#-!)*)!y +"%,)-z "#%-.)"T (T -U "!-Y %-$z)#/($0求解方程式+%20’每一个矩阵+_%’_2’[%’[20都写成两部分’如_23_P 2C _02’其中_P 2表示在一个设定近距离范围内的强相互作用’而_02表示远距离弱相互作用’例如给定一个划分距离1o’源点和场点距离小于1o 为强相互作用’大于1o 为弱相互作用S强相互作用矩阵元使用式+%50Z +%%0进行精确的数值积分计算’并保存在内存中S 所以’1o 越大’意味着需要保存的强矩阵元素越多’消耗计算机内存也越多’但迭代求解收敛速度越快21o 越小’节省内存越多’但迭代收敛速度越慢2经验上一般取1o 为半个波长S弱相互作用矩阵元不这样计算和保存S 把格林函数和它的导数在一个三维网格坐标上进行泰勒展开’如图%所示’D 4x D?x 5D +5,5v 0D 3h 3D 356D +5,5v 0D3h 3D 356D hDl 35lEFD5l +,5v 0D ,l+%N ;0=?3h 3D 35\D +5,5v 0D3h 3D 35\D hDl 35lEFD5l +,5v 0D ,l+%N @0其中3截断泰勒级数项数’6D 和\D 只是离散网格坐标的函数’例如’零阶和一阶项的表达式为J653?+G 70’6%3?v +G 7087e G 7’\53?v +G 70GM 7’\%3+87e G 70$69+G 70,?v +G 70e G 7&G M 7C+%e G 70?v +G 708M 泰勒级数的项数3与表面粗糙性有关’如果使用单层展开’则粗糙性越大3越大2如果使用多层展开’则可以保持3不变’根据表面粗糙性调整展开层数i 8S 在本文的模拟中’粗糙表面的高度方差大约是四分之一波长’我们将保持332’而展开层数取%V 层或N 2层S将式+%N ;0Z +%N @0插入式+%50Z +%%0’并用一点积分规则近似线积分+中点0’用三点积分规则近似三角形面积分+三条棱边的中点0’我们可以把弱矩阵元写成如下形式’比如对_0有JVa N 电子学报255R 年!"#$%$&’()*+(,)+(,-+#$’.+#$/$&’0+#$&’#12’其中-+#$’是与观察点有关的函数%而0+是与源位置有关的函数3矩阵.+具有456789:;结构#平移不变性’%可以用一个具有<=/1个元素的列矢量来表示3对于一个三维网格坐标系统%式#12’中的网格指标应该理解为$(#$>%$?%$@’及$/$&(#$>/$&>%$@/$&?%$@/$&@’3注意式#12’中的#$%$&’是网格点指标%而式#1,’A #11’中的#B %C ’是棱边编号%它们之间存在一个对应的指标变换关系%即第B 条棱边的中点对应一个最近的网格坐标#$>%$?%$@’3降低内存需求是根据以下原因3因为矩阵!的存储量通常是D #=<’3然而%!E的存储量只是D #=’%因为对每一个观察点来说%只有附近几个强相互作用源点%矩阵.+的存储量是D #<=’%列矢量F +和G +的存储量也是D #<=’H 因此#!EI !"’的存储量被缩减到D #J =’%其中J K =是一个常数%典型值是几十3所以内存的节省是巨大的%如果=达到几十万甚至几百万的话3方程式#1<’使用共扼梯度法#L M N %L 5O P Q R S :6M T S U 96O :N6:V 5U’进行求解%所以没有矩阵求逆%只是矩阵W 矢量乘积3强相互作用矩阵与列矢量的乘积耗时很少%因为矩阵都是非常稀疏的3由于近似关系式#12’%弱相互作用矩阵与列矢量的乘积可以使用快速傅立叶#X X 4’技术%例如!"Y Z ()*+(,)+(,F +Y X X 4/1[X X 4.+Y X X 4#G+Y Z ’\#1]’其中两个列矢量的乘积意味着它们对应矩阵元之间相乘3直接计算式#1]’的左边的操作次数是D #=<’3因为三个X X 4的操作次数是^_#<=85R <=’%三个列矢量乘积的操作次数是^_#<=’%所以式#1]’右边的操作次数是D #‘=’%其中‘(1<#*I 1’#*I <’#a 85R <=I a’可以认为是一个常数%因为85R <=随=变化很缓慢3所以每次迭代所节省的L b c 时间是非常可观的%因为‘的典型值是几百%而=为几十万甚至几百万3d e f 局部化近似方法以上所描述的方法是严格的%所以数值结果是精确的3但是%在某些特殊情况下%可以使用一些简化以便加速计算%而精度损失在可接受范围内3对于高损耗媒质来说%格林函数迅速衰减%超过一定距离后可以忽略3在微波频段%海水的介电常数实部和虚部都很大%因此我们可以考虑忽略!"<和g "<的远距离弱相互作用%并称这种近似为局部化近似方法#h i N %h 5j S 89;6Ui 77T 5k 9l S :95ON6:V 5U’H 这样%求解式#1<’节省将近一半L b c 时间3比如%在h 波段下取海水介电常数m n (o ,Ia ,B %则+<(#<p q r ’m s n t#<p q r ’#u 3,I^3^B ’%如果划分距离n v (r q <%则69+<n v w 6/^3^p w1,/2%即格林函数衰减到万分之一3因此%对于损耗媒质%只要n v 足够大%h i N 近似总是适用的3一般情况下%矩量法要求表面离散密度大约是密媒质的十分之一波长%或者对导电媒质来说为十分之一趋肤深度3但是%对于理想导体表面问题%取样密度为自由空间的十分之一波长就足够了3在微波频率下%海水的电导率很高%表现得像一个大导体%因此我们可以使用更稀疏的网格划分%而强相互作用矩阵!E <和g E<仔细精确地进行计算3如果这个假设成立的话%未知量数目就会大大减少%因此内存要求和L b c 时间都会大大节省3我们将用例子来说明在h i N 近似下这种疏网格划分方式的精确度3另外%如果表面的粗糙度#表面高度方差’不是很大%当相互作用距离较远时%我们可以认为原来的粗糙表面近似地为一个平面%所以多层网格方法可以变为单层网格%三维的X X 4简化为二维的X X 4%大大节省了L b c 时间%我们称这种近似为小斜率近似方法%对于低风速驱动下的海面%可以使用这种近似3f 数值结果本文所使用的海洋谱函数为x Q T U 6O W y 6z 6j {|模型%表示为}#+%~’(1<p +!#+’"#+%~’#1a ’文献#1^$对此式有详细解释3利用谱函数产生随机粗糙表面的方法可参考文献#12$3本文数值结果的谱范围选择是1%+%1<,l /1%被模拟的表面尺寸为o 3^1_o 3^1l <%这个尺寸可以充分覆盖所选择的谱范围3在风速1,l q z 时%高度和坡度和数和取波%被约长3用z S R 6行%o +]%2,+u a %一3于似H 因用相半水相基于单积分方程矩量法的海洋表面微波散射模拟设定值就忽略!图"给出了#$次散射模拟的统计结果%数值显示%在镜向反射方向上%侧风向&’()**+,-./散射略大于顺风&0)+-+,-./和逆风向&012.3+,-./散射!对侧风向散射%峰值为4!5.6%而对顺风和逆风向%分别是7!4.6和7!7.6!对三个入射波方向来说%被分析的#$个表面是相同的!但是%对三个不同入射方向%散射半功率宽度大致相同%只与表面尺寸有关%约为5#度!最后%我们给出一个掠入射的例子!被模拟的表面尺寸为898!#4:"7!#;平方波长&在5!#<=>下约#4!#8:;!"5平方米/!因为表面的长度?宽度是4%所以我们可以使用的掠入射角度可略大于2(@A2-4B48!9C%故决定取4#C !被分析的区域使用形致"%方和序大器器%内年。