近場散射問題數值計算的超弱變分方法

欒 天，劉明輝

（1.北華大學 數學與統計學院，吉林 吉林132033；2.吉林大學 數學研究所，長春130012）

隨著近場光學理論的發展，近場散射問題受到人們廣泛關注.由于在遠場光學中存在分辨率衍射極限，而要突破這個極限獲取超分辨率的成像信息與探測信息，只能通過近場區域中的倏逝場實現.因此，研究近場散射問題獲取倏逝波的信息十分重要.近場理論廣泛應用于納米技術、近場光學顯微鏡的研發和生物微樣本的無損成像等領域［1］.目前，關于近場時諧散射問題的數值研究結果報道較少，數值模擬多采用經典的有限元方法、有限差分方法［2］和邊界積分方法［3］等.但隨著波數的增加，這些方法已不再適用.由于時諧波問題的解是振蕩的，因此當波數較高時，需要在單位波長距離內選取足夠多的點才能得到相對穩定的解.而且相位誤差的積累還會導致數值污染［4］.于是要達到給定的精度，就會迫使在每個波長內使用更多的網格節點，計算量較大.基于此，本文采用另一類方法——Trefftz方法，即在區域剖分前提下，先在小單元上使用偏微分方程的解構造離散空間，并在單元交界處借助數值方法使得離散解近似地滿足連續性；再將局部解“拼”接構成一個整體近似解.這類方法主要包括單位分解法［5］、間斷強化法［6］、最小二乘法［7－8］和超弱變分方法［9］等.其中最典型的為超弱變分法.該方法源于區域分解技術，利用Green公式將問題轉化到網格邊界上處理.該方法目前已被應用于有界障礙物的散射問題中［10］.

本文針對近場光學中全內反射顯微鏡的散射模型，提出了相應的超弱變分方法.首先，通過對區域網格剖分引入耦合參數，利用Green公式得到了與原邊值問題等價的超弱變分問題；其次，選取齊次Helmholtz方程的局部解（平面波，倏逝波）構造基底，生成試探函數空間進行離散求解；最后，通過數值實驗驗證算法的有效性.結果表明，該算法能有效地數值模擬近場散射問題，適用于大波數情形，所需計算量少，收斂速度快.

1 數學模型

當沒有樣本時，記空間中的場分布為uref，稱其為參考場［11］.當在分界面Γ0上放置折射率為nS的樣本S時，參考場uref與樣本相互影響產生散射場us，如圖1所示.記此時空間中的總場u為

圖1 幾何模型Fig.1 Geometry of the model

根據Maxwell電磁理論，在TM極化情形下，u滿足如下方程：

空間中的波數為k0n（x），其中折射率

實際數值模擬中，需要將無界區域截斷，為此需引入人工Lipschitz邊界Γ，其外法方向記為ν，將問題限制在有界區域Ω內，即Γ＝?Ω.為了保證所選取的解符合物理要求，要求散射場us在Γ上滿足吸收邊界條件：

于是，由式（1）可知總場u在Γ上滿足

從而所考慮的近場散射問題在數學上可描述為：已知參考場uref，求全場u滿足

其中g＝（?ν－i k0n（x））uref.

2 超弱變分問題

由于區域Ω中包含3種介質，因此為了使全場u在不同介質間的交界面處滿足連續性條件，需引入定義在剖分網格邊界上耦合參數：

由H的定義及u滿足邊值問題，有

將式（6）中兩式相減，得

將式（7）代入式（4），得

由u及其法向導數在Γk，j處的連續性及在Γ上滿足邊界條件，有

從而

將式（10）代入式（8）并關于k求和，得

引入算子F＝（Fk）∈L（X）定義為

其中ek是邊值問題

的解.

方程（12）稱為邊值問題（3）的超弱變分公式.

進一步，定義

按上述相同過程可推得式（11），與式（12）相減得

從而可推得式（9）成立，即u與其法向導數在內邊界連續且滿足外邊界條件，于是由定義知，u是邊值問題（3）的解.

綜上，可得：

3 離散問題

選取有限維離散空間Xh?X代替超弱變分問題（14）中的空間X得到離散問題：求∈Xh，滿足

定義離散空間為

實際數值模擬中，需在每個小單元上選取方向互異的p個平面波函數和2個倏逝波函數作為ek，l.即若xk＝（xk，yk）表示單元Ωk的點，則

為簡單，平面波的方向dk，l通常選取如下：

由于在全反射過程中將出現倏逝波迅速衰減，因此為了能夠通過數值方法捕捉其信息，同時也選用如下倏逝波函數擴充平面波函數構造基底：

選用上述兩個倏逝波，是基于本文的模型問題，依據參考場［11］uref選取的.實際應用中，可由具體模型考慮，選擇適合的倏逝波.此外，由于倏逝波是向上指數衰減的，因此在實際應用中并不要求在每個單元上都應用倏逝波函數構造基底，而是僅在基座界面上方部分單元內使用即可.

4 數值模擬

程序代碼使用 MATLAB語言.計算區域為Ω＝［－l，l］×［－l，l］，將Ω等矩形剖分成2nΩ×2nΩ（nΩ∈?）個小單元，單元邊長為l／nΩ.

例1 設k0＝1，l＝0.5，n＋＝45，n－＝75，nS＝60，θ＝π／4，則上半空間、下半空間和樣本介質對應的波數分別為k＋＝45，k－＝75，kS＝60.

首先，數值求解參考場urefh.取定nΩ＝6，在每個單元上取p＝35個平面波函數.由于倏逝波向上指數衰減，因此只在界面上方第一層單元內引入倏逝波函數擴充平面波函數構造基底.數值計算得到的參考場urefh實部等高線如圖2所示.

其次，為了驗證算法的有效性，考察算法關于p的收斂性.由于此時存在真解［11］：

從而可以計算真解uref和數值解在Ω上的相對誤差：選取初值p＝10，取nΩ＝3.當增加平面波數目，倏逝波數目保持不變，即增大p時.數值結果如圖3和圖4所示.圖3是相對誤差RE關于p的圖像，圖4是－lg（RE）關于p的圖像.由圖3和圖4可見，隨著基底數目的增加，誤差快速衰減，表明算法是關于p收斂的，而且收斂速度較快.當p＝50時，相對誤差已經不超過1%.

最后，數值計算全場uh.取定nΩ＝6，p＝50.樣本大小為界面上方居中的兩個單元格.間接計算散射場＝uh－，其實部等高線如圖5所示.

圖2 urefh的實部等高線Fig.2 Contour for real part of urefh

圖3 相對誤差關于p的曲線Fig.3 Plot of relative error againstp

圖4 －lg（RE）關于p的曲線Fig.4 Plot of－lg（RE）againstp

圖5 ush的實部等高線Fig.5 Contour for real part of ush

綜上可見，本文算法在實際數值模擬中所需剖分單元少，收斂速度快，程序運行僅需數十秒，因此能有效計算近場散射問題.

［1］Carney P，Schotland J.Inside Out：Inverse Problems and Applications［M］.London：Cambridge University Press，2003.

［2］Farakawa H，Kawata S.Analysis of Image Formation in a Near－Field Scanning Optical Microscope：Effects of Multiple Scattering［J］.Opt Commun，1996，132（1／2）：170－178.

［3］Tanaka K，Tanaka M，Omoya T.Boundary Integral Equations for a Two－Dimensional Simulator of a Photon Scanning Tunneling Microscope［J］.J Opt Soc Amer A，1998，15（7）：1918－1931.

［4］Deraemaeker A，Babu?ka I，Bouillard P.Dispersion and Pollution of the FEM Solution for the Helmholtz Equation in One，Two and Three Dimensions［J］.J Numer Meth Eng，1999，46（4）：471－499.

［5］Babu?ka I，Melenk J M.The Partition of Unity Method［J］.Internat J Numer Meth Eng，1997，40（4）：727－758.

［6］Farhat C，Harari I，Franca L P.The Discontinuous Enrichment Method［J］.Comput Meth Appl Mech Eng，2001，190（48）：6455－6479.

［7］Barnett A H，Betcke T.An Exponentially Convergent Nonpolynomial Finite Element Method for Time－Harmonic Scattering from Plygons［J］.SIAM J Sci Comput，2010，32（3）：1417－1441.

［8］LUAN Tian，WANG Yu－jie，ZHENG En－xi.Least Squares Method for Near－Field Scattering Problem ［J］.Journal of Jilin University：Science Edition，2013，51（5）：811－814.（欒天，王玉潔，鄭恩希.求解近場散射問題的最小二乘方法 ［J］.吉林大學學報：理學版，2013，51（5）：811－814.）

［9］Cessenat O，Després B.Application of an Ultra Weak Variational Formulation of Elleptic PDEs to Two－Dimensional Helmholtz Problem ［J］.SIAM J Numer Anal，1998，35（1）：255－299.

［10］Cessenat O，Després B.Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation［J］.J Comput Acoust，2003，11（2）：227－238.

［11］Li P J.Numerical Simulations of Global Approach for Photon Scanning Tunneling Microscopy：Coupling of Finite－Element and Boundary Integral Methods［J］.J Opt Soc Amer A，2008，25（8）：1929－1936.