一种顾及方向遮蔽性的高效空间插值方法.pdf

http:/DOI:10.13700/j.bh.1001-5965.2021.0443一种顾及方向遮蔽性的高效空间插值方法周长聪1，*，刘洪威1，何宝明2，王维2，谭春龙1(1.西北工业大学力学与土木建筑学院，西安710129；2.中国航发商用航空发动机有限责任公司短舱系统部，上海200241)摘要：反距离权重插值方法在航空航天中有着广泛的应用，但其存在仅考虑距离关系而忽视方位关系的缺点，顾及方向遮蔽性的调和反距离权重插值方法弥补了这种不足，提高了插值精度，但仅适用于平面插值。借鉴该方法的基本假设，根据归一化后样本点的不同空间分布，以平面均匀角和球面均匀角为基准，制定统一的均匀性量化标准，提出一种更具普适性的三维空间插值方法。在搜索插值点的临近样本点时，提出一种最近邻搜索算法,极大提高了插值计算效率。通过测试函数计算发现，与反距离权重插值方法相比，所提插值方法误差显著降低。将所提插值方法应用于某型民用飞机短舱的气动载荷插值，结果表明，所提插值方法兼具高效和高精度的优点。关键词：反距离权重插值；方向遮蔽性；球面均匀角；最近邻搜索；气动载荷中图分类号：V211.3文献标志码：A文章编号：1001-5965（2023）06-1278-09在航空航天领域，流场分析更关注几何外形，而结构分析更关注高应力区，导致流场网格与结构网格存在差异1。进行流固耦合计算时，需要在流场网格和结构网格之间进行数据传递，如将流场中网格结点上的气动载荷加载到结构有限元结点上，此时可以采用空间插值方法2-5。空间插值方法种类繁多，按插值区域范围分类，可分为边界内插值、整体插值、局部插值等6。边界内插值方法仅在区域边界进行插值，而区域内部则是均匀和同质的，代表性的插值方法是泰森多边形法7。整体插值方法使用所有采样点拟合插值点处的数值，整体区域的数值会影响单个插值点的数值，代表性的插值方法有趋势面插值方法8等。局部插值方法使用相邻样本点拟合插值点的值，搜索范围的选择会影响插值效果，代表性的插值方法有反距离权重（inversedistanceweight，IDW）插值方法9、Kriging 插值法10和样条函数插值法11等。在上述插值方法中，泰森多边形法受观测值影响较大，未考虑空间因素，实际效果不理想；Kriging 插值法虽然精度较高，但也存在着数学解释复杂、推导过程繁琐、计算效率低等不足；样条函数插值法仅适用于数值变化不剧烈的区域。IDW 插值方法是一种基于地理学第一定律的简单内插方法，凭借其原理简单、计算效率高等优点，被广泛用于流固耦合插值12、地质13和计算机图形学14中。但是，由于其仅考虑了样本点和插值点的距离关系，而忽视了方位关系，导致在样本点相对插值点方向不均匀分布时插值精度不高。针对这一问题，李正泉和吴尧祥15提出一种顾及方向遮蔽性的调和反距离权重（adjustedinversedistanceweight，AIDW）插值方法，该方法较于 IDW 精度有较大提升，但其仅适用于二维平面插值情况。然而工程实际中，尤其在航空领域，经常面临三维插值问题，因此本文借鉴该方法的基本假设，提出一种三维空间的插值方法。本文基于物理学的 Thomson 均匀分布得到空收稿日期：2021-08-05；录用日期：2021-10-07；网络出版时间：2021-11-0114：54网络出版地址： J.北京航空航天大学学报，2023，49（6）：1278-1286.ZHOU C C，LIU H W，HE B M，et al.An efficient spatial interpolation method involving position shadingJ.Journal of BeijingUniversity of Aeronautics and Astronautics，2023，49（6）：1278-1286（in Chinese）.2023年6月北京航空航天大学学报June2023第49卷第6期JournalofBeijingUniversityofAeronauticsandAstronauticsVol.49No.6间均匀点16，将 AIDW 插值方法进行三维拓展。为快速搜索插值点的临近样本点，发展了一种解耦划分（decouplingpartition，DP）算法，并提出了调和反距离权重-解耦划分（adjustedinversedistanceweight-decouplingpartition，AIDW-DP）插值方法。对不同类型的函数进行测试，证明本文所提出的 AIDW-DP 插值方法具有高精度、高效率的优势，将其应用于某型航空发动机短舱的气动载荷插值，取得了良好的效果。1IDW 插值方法的精度限制通过泰勒展开，发现 IDW 精度与样本点和插值点方位有关17。方位不均匀时，仅能保证一阶精度；方位均匀时，具有二阶精度。分析过程如下。根据 IDW 表达式，对于任意函数(x)在空间中的插值问题，函数定义域中插值点 p 处的函数值可以表示为 p：p=iiiii（1）式中：i为样本点 i 对应的函数值；i为插值权重系数，表示为i=1|xixp|r（2）其中：xi为第 i 个样本点的坐标，按照距插值点的欧氏距离从近及远排序；xp为插值点 p 的坐标；r 为权指数。将函数(x)在插值点 p 处进行泰勒展开，则样本点 i 的函数值 i可以表示为i=ep+(ep)(xixp)+12(xixp)TH(xixp)（3）式中：上标 e 表示该点处的精确值；为 Nabla 算子；H 为 Hessian 矩 阵，是 函 数 二 阶 偏 导 数 构 成 的方阵。将式（3）代入式（1）中可得p=ep+(ep)iniii+12i(xixp)TH(xixp)iii（4）式中：ni=(xixp)/?xixp?r由式（4）可知，当且仅当样本点 i 关于插值点p 对称时，等号右边第 2 项为零，IDW 具有二阶精度。在大多数插值问题中，样本点往往在插值点附近呈无规律分布，导致 IDW 只有一阶精度。2AIDW 插值方法2.1基本假设通过第 1 节所述，可以看出 IDW 的精度与样本点和插值点的方位相关。为了考虑这种方位相关性，AIDW 相较于 IDW 增加了可反映插值点和样本点方位关系的遮蔽系数Ki。对于二维平面的AIDW，其基本假设如下：1）当 2 个样本点与插值点连线的夹角 ij小于360/m（m 为搜索点数）时，靠近插值点的样本点对远离插值点的样本点存在遮蔽效应，即远离插值点的样本点的权重被这种方位关系削弱。360/m 被称为均匀性量化标准。2）已被遮蔽的点对其余点不具遮蔽效应。2.2方法表达式AIDW 具体表达式如下：p=iiiKiiiKi（5）式中：Ki为遮蔽系数，表示为样本点 i 受其他样本点的综合遮蔽效应，表达式如下：Ki=1i=1j=i1j=1Liji=2,3,m（6）其中：Lij=1ij 360/msinrij=(1cos2ij)r2iji）无遮蔽效应，样本点v（jvmnmm、n大小比较?根据0m和Q计算均匀角m计算样本点和插值点距离d并排序,排序后为di计算样点间距离f代入式(5)计算p值样本点处函数值i计算得到Ki计算cos ij计算两样本点连线中点坐标检查Pij是否已经赋值?ijm且ij0ijmLij=1判断m与ij大小关系ij=0Lij=0Liu=1Liv=1存储Lij计算Lij插值前准备插值计算图7AIDW-DP 程序设计流程Fig.7ProgrammingflowofAIDW-DP1282北 京 航 空 航 天 大 学 学 报2023年表1测试函数Table1Testfunctions函数编号函数表达式1x2+y2+z22x3+y3+z3+xy2+xz2+yz2+xyz33(1x)2ex2(y+1)210(xx3y5)ex2e(x+1)2y2/3+ez样本空间为以原点为几何中心、边长为 6 的正方体，样本点为样本空间内的 nnn 个随机点。插值空间为在以原点为几何中心、边长为 5 的正方体，插值点为插值空间内均匀分布的 161616 个点。首先，在权指数取 2 条件下，搜索点数依次取4，5，30 进行分析，计算 161616 个插值点的插值平均误差（以下均为相对误差），结果如图 8（a）510152025300.0300.0350.0400.045(a)函数1,样本点数202020,r=2搜索点数搜索点数搜索点数搜索点数搜索点数搜索点数搜索点数搜索点数搜索点数平均误差平均误差平均误差平均误差平均误差平均误差平均误差平均误差平均误差510152025300.015 00.017 50.020 00.022 5510152025300.0080.0100.0120.014510152025300.160.200.240.28510152025300.040.060.080.10510152025300.050.060.070.08510152025300.30.40.50.6510152025300.160.200.240.28510152025300.0900.1050.1200.135 Lus AIDW AIDW-DP IDW(b)函数1,样本点数404040,r=2(c)函数1,样本点数606060,r=2(d)函数2,样本点数202020,r=2(e)函数2,样本点数404040,r=2(f)函数2,样本点数606060,r=2(g)函数3,样本点数202020,r=2(h)函数3,样本点数404040,r=2(i)函数3,样本点数606060,r=2搜索点数搜索点数搜索点数平均误差平均误差平均误差510152025300.030.060.090.12510152025300.0080.0160.0240.032510152025300.0080.0120.0160.020(j)函数3,样本点数202020,r=1(k)函数3,样本点数404040,r=1(l)函数3,样本点数606060,r=1图8插值平均误差随搜索点数的变化Fig.8Variationofinterpolationaverageerrorwithnumberofsearchpoints第6期周长聪，等：一种顾及方向遮蔽性的高效空间插值方法1283图 8（i）所示。为了验证 3.2 节近似的准确性，选取半径为 1 的球面这一极端情况进行测试，插值点为球面上 64 个点，样本点为球面上的nnn 个随机点，结果如图 8（j）图 8（l）所示。从图 8 得到：相较于 IDW 和 LusAIDW，本文 AIDW-DP 插值方法具有更高的精度，这也说明 3.2 节中的近似是合理的；3 种插值方法的插值平均误差随样本点密度的增大都呈减小趋势，这是由于在样本点密度较大时，样本点距插值点距离相对较小，两者相关性更强19；3 种插值方法的最佳搜索点数随样本密度的增大而呈增大趋势。为了同时兼顾 AIDW-DP 的插值精度和使用便捷性，本文建议在样本点密度小于 125 个/单位体积时，搜索点数为 12，对应的均匀角为正二十面体两相邻顶点与体心夹角（63.435）；在大于等于 125 个/单位体积时，搜索点数为 20，对应的均匀角为正十二面体两相邻顶点与体心夹角（41.810）。对于表面插值，如飞行器表面的气动载荷插值、温度插值等，搜索点数定为 12。插值平均误差随权指数变化如图 9 所示。可见，针对所选的不同函数在不同样本空间下的最佳权指数均为 2。在应用中，可将该参数选择范围设为 13。3 种方法的插值计算耗时与函数值无关，仅与样本点分布、样本点密度、插值点个数相关，因此选用函数 3 进行计算耗时对比研究。在插值点个数为 444 个时，3 种方法计算耗时如表 2 所示。IDW 和 LusAIDW 随着格点加密计算耗时急剧增大，而 AIDW-DP 计算耗时始终保持在很低水平，这是由于 DP 算法减小了样本点搜索计算成本。样本点数:202020样本点数:404040样本点数:606060123450.010.020.030.04权指数平均误差(a)函数1权指数平均误差(b)函数2权指数平均误差(c)函数31234500.10.20.3123450.120.240.360.48图9插值平均误差随权指数的变化Fig.9Variationofinterpolationaverageerrorwithweightindex表2不同插值方法的计算耗时对比Table2Comparisonofcomputationtimeofdifferentinterpolationmethodss样本点数IDWLusAIDWAIDW-DP2020202.732.710.1140404021.9421.900.1360606074.3874.480.144.2发动机短舱气动载荷插值某型发动机短舱的三维有限元分析计算模型如图 10 所示，需要将已知的气动载荷插值到有限元网格上。该模型涉及 232158 个表面流场网格结点及其气动载荷数据，83753 个表面结构网格结点。气动载荷插值计算步骤如下：首先，找到需要施加气动载荷的有限元网格单元编号和单元结点坐标；然后，根据气动网格结点处的压强值，使用插值方法计算得到有限元网格结点压强值，对于某一单元，将单元结点上的压强值取算数平均得到单元压强值；最后，通过修改 ABAQUS 软件的 INP 文件进行载荷施加。应当注意，在加载前应确定加载面正法向，避免反向加载。为定量分析本文方法的误差，将获得的结构有限元结点载荷插值到气动网格结点，本文将该过程称为二次插值。插值误差 可表示为=|qq|2q（12）q式中：为气动结点插值载荷；q 为气动结点载荷。应当注意的是，气动载荷插值结果是通过二次插值得到的，因此误差累积了 2 次，单次插值的误差应当减半。通过二次插值计算误差的过程如图 11所示。使用 AIDW-DP 和 IDW 计算所得误差统计结果如表 3 所示。图10短舱有限元分析模型Fig.10Finiteelementanalysismodelofnacelle1284北 京 航 空 航 天 大 学 学 报2023年图 12 为使用 AIDW-DP 插值短舱气动载荷对比。可以看出，插值得到的气动载荷分布和插值前吻合程度很好。AIDW-DP 插值误差分布如图 13所示。对完整模型的插值计算表明，使用 IDW 插值耗时 4.5h，而使用 AIDW-DP 方法仅耗时 162s，计算高效，解决工程复杂问题时具有显著优势。初始气动载荷q插值插值有限元结点载荷二次插值气动载荷q=2qqq图11二次插值过程Fig.11Processoftwiceinterpolation表3插值误差统计Table3Statisticsofinterpolationerror插值方法105误差均值/104中误差/103最大误差/IDW8.772.728.74AIDW-DP6.402.388.41 气动结点载荷 有限元结点载荷1131121111101096008001 0001 2001 4001 00050005001 0001 500轴向/mm径向/mm压强/kPa图12气动结点载荷和插值后有限元结点载荷Fig.12Aerodynamicnodalloadandinterpolatedfiniteelementnodalload600800100012001 00050005001 000径向/mm轴向/mm00.0010.0020.0030.0040.0050.0060.0070.008插值误差图13气动载荷插值的误差Fig.13Errorofaerodynamicloadinterpolation5结论1)本文所提 AIDW-DP 插值方法使用 DP 算法，在大样本空间和多插值点计算时能够极大节约计算成本。2)通过对 3 种测试函数的计算发现，AIDW-DP 在三维情况下精度较 IDW 和 LusAIDW 有较大提升。在加密样本点时，3 种插值方法的误差均大幅下降，且随着样本点密度的增大，IDW、LusAIDW 和本文所提 AIDW-DP 插值方法的最佳搜索点数均呈增大趋势，对于表面插值，这种趋势不再明显。3)将 AIDW-DP 插值方法运用到结构分析时的气动载荷加载，该载荷加载方式适用于复杂曲面和复杂气动载荷分布，使用简单、计算高效。通过将某型航空发动机短舱的气动载荷加载到结构网格，结果显示，插值误差绝大部分区域位于 0.001045 之下，验证了该方法的可行性。对完整模型的插值计算仅耗时 162s，计算效率较高。4)本文所提 AIDW-DP 插值方法也可用于其他工程问题，如温度场插值、矿产资源探测等不同领域，具有较广泛的应用前景。参考文献（References）李立州.流固耦合数据的界面非线性降维传递M.北京:科学出版社,2018:8-9.LILZ.Interfacialnonlineardimensionalityreductiontransferoffluid-solid coupled dataM.Beijing:Science Press,2018:8-9(inChinese).1樊一达,毛玉明,舒忠平,等.基于压力插值/力等效混合的火箭结构流-固载荷转换方法J.航空学报,2022,43(3):225053.FANYD,MAOYM,SHUZP,etal.Hybridfluid-to-solidloadstransformation based on pressure-interpolation/force-equivalencefor launch vehiclesJ.Acta Aeronautica et Astronautica Sinica,2022,43(3):225053（inChinese）.2DONEGTS.Interpolationofmodeshapes:Amatrixschemeus-ingtwo-waysplinecurvesJ.AeronauticalQuarterly,2016,16(4):333-349.3MENONS,SCHMIDTDP.Conservativeinterpolationonunstruc-turedpolyhedralmeshes:Anextensionofthesupermeshapproachto cell-centered finite-volume variablesJ.Computer Methods inAppliedMechanicsandEngineering,2011,200(41-44):2797-2804.4BOERAD,ZUIJLENA,BIJLH.Reviewofcouplingmethodsfornon-matchingmeshesJ.ComputerMethodsinAppliedMechanicsandEngineering,2007,196(8):1515-1525.5李海涛,邵泽东.空间插值分析算法综述J.计算机系统应用,2019,28(7):1-8.LIHT,SHAOZD.Reviewofspatialinterpolationanalysisal-gorithmJ.ComputerSystems&Applications,2019,28(7):1-8（inChinese）.6罗阳,刘元海,郑彤,等.重污染天气下大气污染排放源强的快速估算方法J.哈尔滨工业大学学报,2018,50(8):76-82.7第6期周长聪，等：一种顾及方向遮蔽性的高效空间插值方法1285LUOY,LIUYH,ZHENGT,etal.MethodofrapidestimationofemissionintensityofairpollutioninheavilypollutedweatherJ.JournalofHarbinInstituteofTechnology,2018,50(8):76-82（inChinese）.MEI S.Geologist-controlled trends versus computer-controlledtrends:Introducingahigh-resolutionapproachtosubsurfacestruc-turalmappingusingwell-logdata,trendsurfaceanalysis,andgeo-spatialanalysisJ.CanadianJournalofEarthSciences,2009,46(5):309-329.8GOTWAYCA,FERGUSONRB,HERGERTGW,etal.Com-parisonofkrigingandinverse-distancemethodsformappingsoilparametersJ.Soil Science Society of America Journal,1996,60(4):1237-1247.9ZHANGT.Statisticalanalysisofenvironmentalspace-timepro-cessesJ.Computers&Geosciences,2008,34(12):1974-1975.10HARDER R L,DESMARAIS R N.Interpolation using surfacesplinesJ.JournalofAircraft,1972,9(2):189-191.11WITTEVEENJ,BIJLH.Explicitmeshdeformationusinginversedistanceweightinginterpolation:AIAA2009-3996R.Reston:AI-AA,2009.12LIUZN,YUXY,JIALF,etal.Theinfluenceofdistanceweighton the inverse distance weighted method for ore-grade estima-tionJ.ScientificReports,2021,11(1):2689.13JINGM,WUJ.Fastimageinterpolationusingdirectionalinversedistanceweightingforreal-timeapplicationsJ.OpticsCommuni-cations,2013,286:111-116.14李正泉,吴尧祥.顾及方向遮蔽性的反距离权重插值法J.15测绘学报,2015,44(1):91-98.LIZ Q,WU Y X.Inverse distance weighted interpolation in-volvingpositionshadingJ.ActaGeodaeticaetCartographicaSi-nica,2015,44(1):91-98（inChinese）.THOMSONFRSJJ.XXIV.Onthestructureoftheatom:Anin-vestigationofthestabilityandperiodsofoscillationofanumberofcorpusclesarrangedatequalintervalsaroundthecircumferenceofacircle;withapplicationoftheresultstothetheoryofatomicstruc-tureJ.PhilosophicalMagazine,1904,7(39):237-265.16柴国亮,苏军伟,王乐.一种保持二阶精度的反距离加权空间插值算法J.计算物理,2020,37(4):393-402.CHAIGL,SUJW,WANGL.Aninversedistanceweightingspa-tiallnterpolationalgorithmwithsecondorderaccuracyJ.ChineseJournal of Computational Physics,2020,37(4):393-402（in Chi-nese）.17SAFF E B,KUIJLAARS A.Distributing many points on asphereJ.MathematicalIntelligencer,1997,19(1):5-11.18LUGY,WONGDW.Anadaptiveinverse-distanceweightingspatialinterpolationtechniqueJ.Computers&Geosciences,2008,34(9):1044-1055.19张锦明,郭丽萍,张小丹.反距离加权插值算法中插值参数对DEM插值误差的影响J.测绘科学技术学报,2012,29(1):51-56.ZHANGJM,GUOLP,ZHANGXD.Effectsofinterpolationparameters in lnverse distance weighted method on DEMaccuracyJ.JournalofGeomaticsScienceandTechnology,2012,29(1):51-56（inChinese）.20An efficient spatial interpolation method involving position shadingZHOUChangcong1，*，LIUHongwei1，HEBaoming2，WANGWei2，TANChunlong1(1.SchoolofMechanics，CivilEngineeringandArchitecture，NorthwesternPolytechnicalUniversity，Xian710129，China；2.NacelleDepartment，AECCCommercialAircraftEngineCo.，Ltd.，Shanghai200241，China)Abstract：Theinversedistanceweightinterpolationmethodhasawiderangeofapplicationsinaerospace,butitonlyconsidersthedistancerelationshipandignorestheazimuthrelationship.Thisshortcomingisaddressedbytheadjustedinversedistanceweightinterpolationmethodwithpositionshading,althoughitisonlyappropriateforplaneinterpolation.Based on the basic assumption of this method,according to the different spatial distribution ofnormalizedsamplepoints,thispaperformulatesaunifieduniformityquantizationstandardbasedontheplaneuniformangleandsphericaluniformangle,andproposesathree-dimensionalspatialinterpolationmethod.Thisstudyproposesanewtechniquethatsignificantlyincreasestheeffectivenessofinterpolationbysearchingsamplepointsclosetotheinterpolationpoints.Throughthecalculationoftestfunctions,itisfoundthatcomparedwiththeinversedistanceweight interpolation method,the error of the proposed method issignificantly reduced.The proposed method isappliedtotheaerodynamicloadsinterpolationofacivilaircraftnacelle,andtheresultsshowthattheproposedmethodhastheadvantagesofhighefficiencyandhighaccuracy.Keywords：inverse distance weight interpolation；position shading；spherical uniform angle；closest pointsearch；aerodynamicloadsReceived：2021-08-05；Accepted：2021-10-07；PublishedOnline：2021-11-0114：54URL： 京 航 空 航 天 大 学 学 报2023年