三角網(wǎng)格的融合與優(yōu)化及其在電磁散射計(jì)算中的應(yīng)用

丁娟娟，孫玉鑫，劉 鵬

(復(fù)旦大學(xué) 電磁波信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,上海 200433)

三角網(wǎng)格的融合與優(yōu)化及其在電磁散射計(jì)算中的應(yīng)用

丁娟娟，孫玉鑫，劉 鵬

(復(fù)旦大學(xué) 電磁波信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,上海 200433)

為實(shí)現(xiàn)數(shù)值計(jì)算中復(fù)合模型的三角網(wǎng)格融合、優(yōu)化電磁散射計(jì)算，提出基于三角形求交、面元內(nèi)角控制的網(wǎng)格模型融合與優(yōu)化算法．通過兩個(gè)獨(dú)立模型中三角面元的交點(diǎn)計(jì)算和模塊內(nèi)點(diǎn)云的Delaunay剖分，獲得初始融合網(wǎng)格，再通過查找、消除畸形面元來優(yōu)化融合后的網(wǎng)格．一系列模型的測試表明，在保持網(wǎng)格采樣信息與幾何外觀的前提下，本文所提出的算法可穩(wěn)健、有效地實(shí)現(xiàn)不同網(wǎng)格密度、不同結(jié)構(gòu)特征的三角面元模型融合，去除畸形面元減少網(wǎng)格面元數(shù)目．

網(wǎng)格融合； Delaunay剖分； 畸形三角面元； 電磁散射計(jì)算

三角面元廣泛應(yīng)用于數(shù)值建模，常用于計(jì)算流體力學(xué)、計(jì)算電磁學(xué)等眾多領(lǐng)域．在數(shù)值模型中，不同幾何體的網(wǎng)格拼接處一般是物理場的耦合區(qū)，為保證數(shù)值計(jì)算的精度和穩(wěn)定性，往往對(duì)拼接處的網(wǎng)格單元形狀有較高的要求．雖然單一模型的三角面元網(wǎng)格建模已實(shí)現(xiàn)工程化，有大量成熟的軟件，但是不同軟件或不同建模方法生成的模型之間往往需要單獨(dú)進(jìn)行網(wǎng)格融合或網(wǎng)格拼接，例如建筑物模型與地表模型，艦船模型與海面模型等，在這些領(lǐng)域網(wǎng)格融合算法顯得尤為重要．

艦船與尾跡的網(wǎng)格模型由不同的算法或軟件獨(dú)立生成的，網(wǎng)格的節(jié)點(diǎn)密度和結(jié)構(gòu)特征存在較大差異．為準(zhǔn)確研究艦船與海面復(fù)合場景的電磁散射，需要考慮船體與海水交接處的二面角強(qiáng)散射機(jī)制．因此，作為本文研究的重點(diǎn)，網(wǎng)格融合算法的正確性、穩(wěn)健性以及融合后網(wǎng)格的品質(zhì)等因素，將直接影響到后續(xù)的電磁散射計(jì)算分析與雷達(dá)成像的工程應(yīng)用．

三角形面元網(wǎng)格的融合拼接方法主要有3類[1]： 第1類是基于裁剪的方法，即利用一片網(wǎng)格去裁剪另外一片網(wǎng)格，然后在公共邊界上生成新的三角形單元，將兩個(gè)模型的網(wǎng)格融合在一起，如Turk和Levoy等的算法[2]．該方法的缺點(diǎn)是因?yàn)椴眉?，公共邊界處?huì)產(chǎn)生大量的細(xì)小三角形，并且該方法只利用了重疊區(qū)一片網(wǎng)格中的頂點(diǎn)，另一片網(wǎng)格上的頂點(diǎn)則被完全拋棄．因此對(duì)于存在大交疊區(qū)域的網(wǎng)格而言，無法同時(shí)利用兩片網(wǎng)格的重疊區(qū)進(jìn)行頂點(diǎn)校正．第2類是基于補(bǔ)洞的方法，即首先將重疊區(qū)的三角形全部刪掉，然后通過補(bǔ)洞的方法重新生成重疊區(qū)的三角形．如Ruding提出的先去除N-環(huán)相交區(qū)再重建交疊區(qū)的方法[3]．該方法對(duì)于交疊區(qū)較小的的網(wǎng)格非常適用，但對(duì)于存在大交疊區(qū)域的網(wǎng)格，雖然可通過徑向基函數(shù)重新生成重疊區(qū)頂點(diǎn)，但新生成的頂點(diǎn)很難反映模型的實(shí)際形狀．第3類是基于微分網(wǎng)格變形的方法，如利用泊松方程[4]、拉普拉斯坐標(biāo)[5]等方法．這類方法需要指定一個(gè)準(zhǔn)確的邊界，邊界的定位精度對(duì)融合結(jié)果影響很大；但是大范圍復(fù)雜交疊區(qū)的拼接往往難以確定精確的融合邊界．

本文通過求取兩個(gè)模型三角面元的所有交點(diǎn)，克服了第1類方法不完全采樣的缺點(diǎn)，采用基于散亂點(diǎn)云的二維Delaunay剖分[6-7]，得到散亂點(diǎn)集的三角剖分結(jié)果．與第3類方法相比，文中網(wǎng)格融合的算法無需附加條件，局部區(qū)域重新剖分與否，完全取決于交疊區(qū)網(wǎng)格三角面元的交點(diǎn)情況，交疊區(qū)面元的數(shù)量僅影響計(jì)算速度，同時(shí)克服了第2類方法在大交疊區(qū)域模型失真的困難．通過去除融合過程中產(chǎn)生的畸形三角面元，優(yōu)化融合后的網(wǎng)格，且保持融合區(qū)網(wǎng)格的采樣信息與幾何外觀．

1 基本定義

定義1三角網(wǎng)格： 用三元組集合F：{v,e,f}來描述，其中v：{v1,v2，…，vn}表示點(diǎn)集合e：{e1,e2,…,en}表示邊集合，f： {f1,f2,…，fn}表示三角形集合．

定義2畸形三角形： 一個(gè)或兩個(gè)內(nèi)角小于閾值δ的三角形，如圖1所示(工程中一般取δ=15°)，小于δ的內(nèi)角稱作病態(tài)角．

定義3兩個(gè)網(wǎng)格的融合： 指當(dāng)兩個(gè)網(wǎng)格模型在三維空間上相交時(shí)，將二者的網(wǎng)格數(shù)據(jù)重新編排為一個(gè)網(wǎng)格體，節(jié)點(diǎn)、面元都有唯一編號(hào)，且拓?fù)潢P(guān)系相應(yīng)改變．

2 融合與優(yōu)化策略

為實(shí)現(xiàn)由獨(dú)立模型到復(fù)合模型的整體轉(zhuǎn)換，需將原來兩個(gè)互相獨(dú)立的網(wǎng)格系統(tǒng)，合并成一個(gè)網(wǎng)格系統(tǒng)，網(wǎng)格系統(tǒng)中的每個(gè)三角面元、每個(gè)節(jié)點(diǎn)都有唯一的編號(hào)，面片之間拓?fù)潢P(guān)系正確，網(wǎng)格質(zhì)量滿足數(shù)值計(jì)算要求．

2.1 三維點(diǎn)云的平面投影及剖分

在本文融合算法中，待融合區(qū)網(wǎng)格三角面元交點(diǎn)產(chǎn)生三維點(diǎn)云，先將點(diǎn)云投影到二維空間，如圖2所示，再將三維點(diǎn)云繞定軸r逆時(shí)針旋轉(zhuǎn)α角到x-y面，得到投影點(diǎn)集．對(duì)投影點(diǎn)集作平面域的Delaunay剖分，構(gòu)建出的二維投影點(diǎn)拓?fù)溥B接關(guān)系，與原三維點(diǎn)云拓?fù)溥B接關(guān)系相同．

2.2 三角網(wǎng)格的融合

三角網(wǎng)格F1、F2，融合算法流程如圖3所示： 融合區(qū)域M1、M2分別有P、Q個(gè)面元，經(jīng)正(紅色虛線框)、反(綠色虛線框)兩次融合，將每個(gè)面元區(qū)的三維點(diǎn)集投影到x-y面進(jìn)行Delaunay剖分，生成每個(gè)面元區(qū)新的網(wǎng)格拓?fù)?；最后合并M3與M4，并進(jìn)行網(wǎng)格節(jié)點(diǎn)與面元查重．對(duì)于電磁散射計(jì)算等僅需表面網(wǎng)格的情況，還應(yīng)去除融合網(wǎng)格內(nèi)部不可見部分．

圖3 三角網(wǎng)格融合算法流程圖Fig.3 Flow chart of the merging algorithm for triangular meshes

2.3 網(wǎng)格優(yōu)化算法

三角網(wǎng)格優(yōu)化目前比較成熟的準(zhǔn)則有以下5種： Thiessen區(qū)域準(zhǔn)則、最小內(nèi)角最大準(zhǔn)則、圓準(zhǔn)則、ABN準(zhǔn)則和PLC準(zhǔn)則[8-9]．就三角形形態(tài)來說，等邊三角形最接近圓，是品質(zhì)最高的三角形,是三角面元的最高優(yōu)化目標(biāo)[10]．

文獻(xiàn)[11]提出一種改善畸形三角形的方法，通過刪除畸形三角形和調(diào)整生成點(diǎn)的相對(duì)位置來實(shí)現(xiàn)，然而不能完全滿足控制內(nèi)角范圍的要求．為解決任意兩個(gè)三角網(wǎng)格模型融合后，局部網(wǎng)格畸形三角面元多而集中的問題，本文提出一種新方法，即通過控制畸形三角面元3個(gè)內(nèi)角的大小，來改善三角網(wǎng)格的質(zhì)量．

以等邊三角形為標(biāo)準(zhǔn)可定義三角形正則度如下[10]：

R(η)=2(cosA+cosB+cosC-1)，

(1)

其中A、B、C分別是三角形的3個(gè)內(nèi)角．只需知道三角形一個(gè)內(nèi)角，即可計(jì)算其正則度：

(2)

其中X是某個(gè)內(nèi)角，X從0到π變化時(shí)三角形的正則度如表1所示．

表1 三角形正則度R(X)隨X變化規(guī)律Tab.1 The relationship between the triangle regularity R(X) and X

根據(jù)正則度的定義，可進(jìn)一步定義如下的三角網(wǎng)格品質(zhì)因子Q：

(3)

即網(wǎng)格中所有三角面元的正則度均值決定了網(wǎng)格的品質(zhì)因子，式(3)中：N是網(wǎng)格模型中三角面元的數(shù)量；X是第i個(gè)三角形的某個(gè)內(nèi)角；Ri(X)是第i個(gè)三角形的正則度．如式(2)所示．品質(zhì)因子越大，說明網(wǎng)格中三角面元的正則度越大、越接近等邊三角形，網(wǎng)格的質(zhì)量越高．

本文選取的優(yōu)化閾值為15°，即判定內(nèi)角小于15°的三角面元為畸形面元．網(wǎng)格優(yōu)化關(guān)鍵在于恰當(dāng)處理畸形面元，需綜合考慮畸形面元的相鄰面元．優(yōu)化原理如圖4所示： 淺藍(lán)的畸形三角面元，點(diǎn)v1對(duì)應(yīng)內(nèi)角小于15°，其對(duì)邊的中點(diǎn)為vn．為消除畸形面元，需刪除節(jié)點(diǎn)v6、v7，依次連接點(diǎn)vn與節(jié)點(diǎn)v1～v5，得到不含畸形面元的新的網(wǎng)格結(jié)構(gòu)，如圖4(b)所示．圖5(a)為存在多個(gè)畸形面元的網(wǎng)格，圖5(b)為本文算法優(yōu)化后的結(jié)果，不再有畸形三角面元，優(yōu)化前網(wǎng)格的品質(zhì)因子為0.331，優(yōu)化后為0.718，網(wǎng)格品質(zhì)顯著提高．

3 融合實(shí)例及應(yīng)用

3.1 融合優(yōu)化實(shí)例

為進(jìn)一步驗(yàn)證本文算法的有效性，本節(jié)進(jìn)行復(fù)雜模型的融合測試： 第1例為球面和隨機(jī)粗糙面，如圖6(a)所示；第2例為艦船與海面尾跡，如圖6(b)所示；第3例為球面和立方體表面，如圖6(c)．分析3例融合結(jié)果，模型的幾何特征得以很好地保存，基本未發(fā)生幾何形變，綜合考慮工程計(jì)算的需求和融合區(qū)的幾何特征，合理選取優(yōu)化角度，還可進(jìn)一步減小網(wǎng)格的幾何失真．為比較優(yōu)化前后網(wǎng)格面元的內(nèi)角分布情況，圖6(d)以第2例模型為例，統(tǒng)計(jì)面元內(nèi)角在0°～180°的分布．優(yōu)化前網(wǎng)格模型存在大量內(nèi)角在0°～15°的三角面元，而優(yōu)化后已完全消除．

圖6 復(fù)雜模型優(yōu)化前后的結(jié)果對(duì)比Fig.6 Comparison of the merged meshes of complex models before and after optimization

表2是圖6中(a)、(b)、(c)模型融合優(yōu)化前后面片數(shù)和網(wǎng)格品質(zhì)因子的對(duì)比，其中，艦船與尾跡的融合區(qū)域面元數(shù)量由原來的3473個(gè)減少到2410個(gè)，面元數(shù)減少30.6%．在傳統(tǒng)的融合算法中，基于頂點(diǎn)融合的方法為保持幾何外觀，融合后模型的面片只增不減；而基于裁剪補(bǔ)洞的融合方法需要生成更多的點(diǎn)和面片，以彌補(bǔ)裁減交疊區(qū)導(dǎo)致的幾何失真，二者都不利于提高計(jì)算效率．在本文融合算法中，生成點(diǎn)為三角形的交點(diǎn)，因此融合過程不改變模型的幾何特征．

表2 融合前后面片數(shù)與網(wǎng)格品質(zhì)因子的對(duì)比Tab.2 The comparison of face numbers and mesh quality

3.2 二面角模型RCS算例

邊長為a、b，公共邊長l，夾角2β的四邊形二面角單站RCS由式(4)給出，Sa,Sb是二面角一次散射PO解，Sab、Sba是a→b和b→a的二次散射，

(4)

根據(jù)物理光學(xué)法(PO)散射計(jì)算公式[12]，在E極化下：

(5)

其中：

(6)

(7)

其中：α=π-3β；tanγ=bsin 2β/(α-bcos 2β)；波數(shù)k=2π/λ．在H極化下，用-sin(3β?θ)代替式(5)中的sin(3β±θ)．

圖7(a)中四邊形平板a×b=12cm×12cm，所形成二面角l=12cm，夾角2β=90°，平面波E極化，頻率為10GHz,入射角θ從-45°～45°．圖7(b)中藍(lán)點(diǎn)為實(shí)測數(shù)據(jù)[13]，紅虛線為未融合的計(jì)算結(jié)果，黑實(shí)線為融合二面角的PO計(jì)算結(jié)果，融合二面角計(jì)算結(jié)果與實(shí)測數(shù)據(jù)吻合．±40°附近二面角散射較強(qiáng)，由于兩個(gè)平板的融合區(qū)很小，未融合的網(wǎng)格僅在±40°附近產(chǎn)生誤差．

圖7 二面角RCS(10GHz)Fig.7 RCS of the Dihedral(10GHz)

3.3 艦船與尾跡融合及RCS計(jì)算

艦船模型根據(jù)典型驅(qū)逐艦的構(gòu)造建模，如圖8(a)所示.500m×500m的尾跡模型如圖8(b)所示，圖中淺黃區(qū)域所示為待融合網(wǎng)格區(qū)域，網(wǎng)格品質(zhì)因子Q為0.667．融合后模型如圖8(c)，Q提高到0.826．

圖8 艦船與尾跡的融合與優(yōu)化Fig.8 Merge and optimization of ship and wakes

為進(jìn)一步比較融合模型的計(jì)算結(jié)果，將本文與文獻(xiàn)[14]中艦船與海面模型的雙站RCS進(jìn)行對(duì)比．具體計(jì)算參數(shù)如下： 雷達(dá)入射平面波頻率為5GHz，θi=45°，φi=0°，HH極化；雙站散射角θs=-90°～90°，φs=0°．結(jié)果如圖8(d)所示，紅虛線為本文模型融合前的結(jié)果，黑實(shí)線為融合后的結(jié)果，藍(lán)實(shí)線為文獻(xiàn)[14]的結(jié)果，可見由于二面角的強(qiáng)散射機(jī)制，融合與未融合模型在-90°～-70°、-30°、-10°附近區(qū)別較大．本文結(jié)果與文獻(xiàn)[14]趨勢基本相同．由于艦船結(jié)構(gòu)、隨機(jī)粗糙海面以及尾跡模型等方面的差異，二者的雙站散射在某些角度存在差別．

4 結(jié)束語

本文介紹了一種復(fù)雜三角網(wǎng)格模型的融合方法．算法克服了補(bǔ)洞拼接法在大交疊區(qū)域的模型失真和裁剪方法的不完全采樣的缺點(diǎn)．算法與網(wǎng)格交疊區(qū)的大小無關(guān)，不受附加條件限制，可以優(yōu)化、去除網(wǎng)格融合過程中產(chǎn)生的畸形三角面元，并保持網(wǎng)格的采樣信息與幾何外觀．

[1] 鄒北驥，周浩宇，王磊等.大交疊區(qū)域的三維網(wǎng)格的融合與拼接 [J].電子學(xué)報(bào)，2012，40(5)： 1005-1010.

[2] TURK G, LEVOY M. Zippered polygon meshes from range images [C]∥. Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Techniques(SIGGRAPH’94). Oriando, FL, USA, 1994： 311-318.

[3] RUDING L, JEAN-PHILIPPE P, ALEXEI M,etal. Merging enriched finite element triangle meshes for fast prototyping of alternate solutions in the context of industrial maintenance [J].CADComputerAidedDesign, 2010,42(8)： 670-681.

[4] YU Y, ZHOU K, XU D,etal. Mesh editing with poisson-based gradient field manipulation [J].ACMTransactionsonGraphics, 2004,23(3)： 644-651.

[5] SORKINE O. Differential representations for mesh processing [J].ComputerGraphicsForum, 2006,25(4)： 789-807.

[6] VORONOI G. Nouvelles applications des parameters continues a la theorie des formes quadratiques Deuxieme memorie： recherches sur les parallelloedres primitives [J].JournalFurDieReineundAngewandteMathematik, 1908,134： 198-287.

[7] LEE D T, SCHACHTER B J. Two algorithms for constructing a delaunay triangulation [J].InternationalJournalofComputerandInformationSciences, 1980,9(3)： 219-242.

[8] PARK S C. Polygonal extrusion [J].TheVisualComputer, 2003,19(1)： 38-49.

[9] WANG C C L, YUEN M M F. Sketch based mesh extrusion with remeshingtechniques [J].ProceedingsoftheASMEDesignEngineeringTechnicalConference, 2001,1： 731-738.

[10] 劉泗巖,廖文和,劉浩.基于內(nèi)角余弦和的三角形正則度評(píng)定與網(wǎng)格優(yōu)化 [J].機(jī)械科學(xué)與技術(shù)，2007，26(4)： 420-423.

[11] 王群,李愛平,馬淑梅.局部網(wǎng)格狹長三角形的品質(zhì)改善及實(shí)現(xiàn) [J].同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版)，2004，32(11)： 1508-1511.

[12] KNOTT E F. RCS reduction of dihedral corners [J].IEEETransactionsonAntennasandPropagation, 1977,25(3)： 406-409.

[13] ANDRADE L A, NOHARA E L, PEIXOTO G G,etal. Backscattering analysis of flat plate and dihedral corner reflectors using PO and comparison with RCS measurements in anechoic chamber [J].IEEEMTT-SInternational, 2003,2： 719-724.

[14] ZHAO Y, ZHANG M, ZHAO Y W,etal. A bistatic SAR image intensity model for the composite ship-ocean scene [J].IEEETransactionsonGeoscience&RemoteSensing, 2015,53(8)： 1-9.

AnAlgorithmforMergingandOptimizingtheTriangularMeshesofShipandSeaSurface

DINGJuanjuan,SUNYuxin,LIUPeng

(KeyLaboratoryforInformationScienceofElectromagneticWaves,FudanUniversity,Shanghai200433,China)

This paper presents a new method to merge and optimize the triangular meshes of ship and sea surface for numerical calculations. The method is divided into three steps： Firstly, identify the intersected nodes of the ship mesh and the sea mesh. Then, obtain the initial merged mesh by performing the Delaunay triangulation for the nodes in the intersected region. Finally, optimize the merged mesh by eliminating the abnormal triangles and renumbering the remaining nodes and elements. A series of examples show that this method is capable of merging the meshes with different structure and density, such as those of the ship, the rough sea surface, and the ship wakes. The method can not only run without any auxiliary conditions, but also improve the quality of the merged mesh.

mesh merge; Delaunay triangulation; abnormal triangles; electromagnetic scattering calculation

0427-7104(2017)06-0712-07

2017-02-21

國家自然科學(xué)基金(61179022)

丁娟娟(1991—)，女，碩士研究生；劉 鵬，副教授，通信聯(lián)系人，E-mail： pliu@fudan.edu.cn.

TP391.7

A