D0I:10.13374/i.issnl001t03.2007.11.020 第29卷第11期 北京科技大学学报 Vol.29 No.11 2007年11月 Journal of University of Science and Technology Beijing Nov.2007 超声测距系统的高精度中心椭圆算法 陈先中侯庆文叶琳张红宇 北京科技大学信息工程学院，北京100083 摘要描述了基于新型驱动电路和低噪声程控放大器的超声测距系统·驱动系统的设计从硬件角度提高了回波信号的信 噪比和稳定性·基于能量重心校正法和最小二乘法的原理，提出了一种改进型椭圆中心超声回波寻峰算法：通过曲线拟合搜 索椭圆中心点，即回波信号能量集中点，进而找到回波信号的峰值点·测试比较证明，与包络线法和三次多项式法相比，此算 法测量相对误差稳定在0.2%，适用高精度工业测量。 关键词测距系统：超声波；中心椭圆算法：回波信号 分类号TM938.82 在各种非接触测量技术中，超声测距是一种非 心椭圆算法提供良好的硬件基础。系统采用以下三 常实用的技术，随着超声技术的逐渐成熟，超声测量 种措施：采用新型接口电路)，使回波信号与发射 作为一种高效、准确、非接触式的测量和检测手段已 电路隔离，避免单探头工作时回波信号进入发射电 经广泛地应用于很多行业，超声测距精度的关键是 路；采用模拟开关式程控放大器，根据回波的强度设 准确测量超声波传播时间山.本文针对这一问题在 定适当的放大倍数，并充分利用A/D转换器转换精 超声测距系统硬件上采用了程控放大器提高信噪 度，实现低噪声程控放大；采用带通滤波电路，进行 比，在回波信号处理上提出了一种高精度椭圆中心 二阶滤波，保留超声共振频率信号，使回波信号在幅 寻峰算法 值上和频率范围内滤除粗大噪声，易于实行下一步 1超声测距系统的基本组成 的高精度算法， 外围I/O以及测量控制电路可以完成系统的 超声测距系统主要包括换能器驱动电路、超声 显示和按键控制功能 波回波信号接收电路、外围I/0及测量控制电路等 三部分.系统组成框图如图1所示， 2超声测距系统的高精度算法 键盘 经过信号预处理，信号已经初次滤波，可以进行 下一步的全范围峰值搜索，接收到的信号波形如 温度补偿· ATmegal6L 换能器 单片机 驱动电路 图2所示，为避免滤波带来的滞后现象，采用零相 位滤波，可得到无相位偏移的超声波回波信号，然 输出控制 程控放大 显示 回波接收 后使用阈值法，将接收信号前面的高幅值部分屏蔽， 得到真正回波信号峰值左右的超声波信号波形， 300 图1系统组成框图 Fig.1 Structure of the system 250 200 换能器驱动电路的好坏会直接影响到超声回波 的质量，为达到较好的实验效果，系统采用变压器 驱动电路，其电压驱动功率高，能满足测量距离变 100 化大的要求 回波信号接收电路可以为使用高精度超声波中 20040060080010001200 收稿日期：2006-06-29修回日期：2006-09-25 时间，终 基金项目：国家863高技术研究发展计划资助项目(N。 图2超声波接收信号波形 2006AA04Z177) Fig.2 Receive signal of the ultrasonic ranging system 作者简介：陈先中(1966一)，男，副教授，博士
超声测距系统的高精度中心椭圆算法 陈先中 侯庆文 叶 琳 张红宇 北京科技大学信息工程学院北京100083 摘 要 描述了基于新型驱动电路和低噪声程控放大器的超声测距系统．驱动系统的设计从硬件角度提高了回波信号的信 噪比和稳定性．基于能量重心校正法和最小二乘法的原理提出了一种改进型椭圆中心超声回波寻峰算法：通过曲线拟合搜 索椭圆中心点即回波信号能量集中点进而找到回波信号的峰值点．测试比较证明与包络线法和三次多项式法相比此算 法测量相对误差稳定在0∙2％适用高精度工业测量． 关键词 测距系统；超声波；中心椭圆算法；回波信号 分类号 T M938∙82 收稿日期：2006-06-29 修回日期：2006-09-25 基金 项 目： 国 家863高 技 术 研 究 发 展 计 划 资 助 项 目 （ No． 2006AA04Z177） 作者简介：陈先中（1966—）男副教授博士 在各种非接触测量技术中超声测距是一种非 常实用的技术随着超声技术的逐渐成熟超声测量 作为一种高效、准确、非接触式的测量和检测手段已 经广泛地应用于很多行业．超声测距精度的关键是 准确测量超声波传播时间［1］．本文针对这一问题在 超声测距系统硬件上采用了程控放大器提高信噪 比在回波信号处理上提出了一种高精度椭圆中心 寻峰算法． 1 超声测距系统的基本组成 超声测距系统主要包括换能器驱动电路、超声 波回波信号接收电路、外围 I／O 及测量控制电路等 三部分．系统组成框图如图1所示． 图1 系统组成框图 Fig．1 Structure of the system 换能器驱动电路的好坏会直接影响到超声回波 的质量．为达到较好的实验效果系统采用变压器 驱动电路．其电压驱动功率高能满足测量距离变 化大的要求． 回波信号接收电路可以为使用高精度超声波中 心椭圆算法提供良好的硬件基础．系统采用以下三 种措施：采用新型接口电路［2］使回波信号与发射 电路隔离避免单探头工作时回波信号进入发射电 路；采用模拟开关式程控放大器根据回波的强度设 定适当的放大倍数并充分利用 A／D 转换器转换精 度实现低噪声程控放大；采用带通滤波电路进行 二阶滤波保留超声共振频率信号使回波信号在幅 值上和频率范围内滤除粗大噪声易于实行下一步 的高精度算法． 外围 I／O 以及测量控制电路可以完成系统的 显示和按键控制功能． 图2 超声波接收信号波形 Fig．2 Receive signal of the ultrasonic ranging system 2 超声测距系统的高精度算法 经过信号预处理信号已经初次滤波可以进行 下一步的全范围峰值搜索．接收到的信号波形如 图2所示．为避免滤波带来的滞后现象采用零相 位滤波可得到无相位偏移的超声波回波信号．然 后使用阈值法将接收信号前面的高幅值部分屏蔽 得到真正回波信号峰值左右的超声波信号波形 第29卷 第11期 2007年 11月 北 京 科 技 大 学 学 报 Journal of University of Science and Technology Beijing Vol．29No．11 Nov．2007 DOI:10．13374／j．issn1001－053x．2007．11．020
第11期 陈先中等：超声测距系统的高精度中心椭圆算法 ,1155 图3为峰值部分的放大图，若想提高测量精度，就 展开整理后使用一种变量代换方法[6]将复杂 要测得精确的超声波在空气中传播时间，根据对滤 的非线性方程简化，变换后平面任意位置椭圆方程 波后的信号分析，以及能量重心校正方法和最小二 式可以表示为如下形式： 乘法，提出了一种寻找回波信号峰值点的新算法 x2+Axy+By2+Cx+Dy+E=0 (2) 一中心椭圆算法，以精确推算超声波接收信号的 这样，求解椭圆系数的方程就是一组关于A、 峰值点，提高测量精度 B、C、D和E的线性方程，根据最小二乘法原理利 70F 用回波信号峰值点附近的N个点的坐标，求得椭圆 65 原信号 零相位滤波信号 系数，得到椭圆方程，建立椭圆拟合的数学模型，来 60 一般滤波信号 获取A、B、C、D和E这五个特征参数. 设定椭圆长轴与x轴的夹角0=0°或90°，即椭 圆对称轴与坐标轴平行的情况下进行数据拟合，可 把式(2)变换为： x2+Ay2+Bx十C+D=0 (3) 35 这样就可由A、B、C和D来确定所拟合的椭圆，在 20002105210021502200225023002350 时间，s 进行拟合时，采用最小二乘法门依次确定式(3)中 的未知数A、B、C和D,由此推导出椭圆的基本参 图3经过零相位滤波后的信号 数：椭圆中心坐标(x0,yo),长轴半径a,短轴半 Fig.3 Signals after no-bias phase filtering 径b. 能量重心校正法是针对加对称窗情况的算 法③].选择合适的窗函数可以减小频谱泄漏误差 依据上述最小二乘法的原理，测量点与理想曲 以Hanning窗为例，Hanning窗的定义为 线间误差为： w(n)=0.5-0.5cos(2πn/N),n=1,2,…,N-1, 9(A,B,C,D)= (x+Ay+Bx:+Cyi+D)2 其频谱模函数为 W(x)=sinzx 1 (4) πx2(1-x2) 由数学分析可知使A、B、C和D满足导数为 令功率谱函数G(x)=W(x),则有 零的点为极值点： sinfπx G(x)-4x21-x22 a9_9=9=9=0 aa ab ac aD (5) 对任意一确定值x,当no(n为自然数)， 由此得下列方程组： G(x+(x+)=0 =一 成立.表明noo时，Hanning窗离散频谱的能量重 心无穷逼近坐标原点. 之 A 根据能量重心校正法的原理，离散窗谱的能量 重心都在原点附近，而超声波回波信号的能量分布 情况与离散窗谱类似，所以提出中心椭圆算法，对经 过滤波和阈值法处理过的信号进行运算，拟合出椭 N 圆，回波信号能量都集中在椭圆的中心点，这一点即 为回波信号的峰值点，具体的算法如下： 处于XY平面内任意位置的椭圆可以用五个 独立参数来唯一确定，这五个参数是：椭圆中心坐标 (x0,yo),长轴半径a,短轴半径b,长轴与x轴的 夹角.将平面任意位置椭圆的方程表达为)]： (6) [(x-xo)cos0(y-ro)sin0 d [-(x-xo)sin0+(y-yo)cos0]=1 6 (1)
图3为峰值部分的放大图．若想提高测量精度就 要测得精确的超声波在空气中传播时间根据对滤 波后的信号分析以及能量重心校正方法和最小二 乘法提出了一种寻找回波信号峰值点的新算法 ———中心椭圆算法以精确推算超声波接收信号的 峰值点提高测量精度． 图3 经过零相位滤波后的信号 Fig．3 Signals after no-bias phase filtering 能量重心校正法是针对加对称窗情况的算 法［3］．选择合适的窗函数可以减小频谱泄漏误差． 以 Hanning 窗为例Hanning 窗的定义为 w（ n）＝0∙5—0∙5cos（2πn／N）n＝12…N—1 其频谱模函数为 W（ x）＝ sinπx πx 1 2（1— x 2） 令功率谱函数 G（ x）＝ W 2（ x）则有 G（ x）＝ sin 2πx 4π2 x 2（1— x 2） 2 对任意一确定值 x当 n→∞（ n 为自然数） ∑ n i＝—n G（ x＋ i）（ x＋ i）＝0 成立．表明 n→∞时Hanning 窗离散频谱的能量重 心无穷逼近坐标原点［4］． 根据能量重心校正法的原理离散窗谱的能量 重心都在原点附近而超声波回波信号的能量分布 情况与离散窗谱类似所以提出中心椭圆算法对经 过滤波和阈值法处理过的信号进行运算拟合出椭 圆回波信号能量都集中在椭圆的中心点这一点即 为回波信号的峰值点．具体的算法如下： 处于 XY 平面内任意位置的椭圆可以用五个 独立参数来唯一确定这五个参数是：椭圆中心坐标 （ x0y0）长轴半径 a短轴半径 b长轴与 x 轴的 夹角θ．将平面任意位置椭圆的方程表达为［5］： ［（ x— x0）cosθ＋（ y—y0）sinθ］ 2 a 2 ＋ ［—（ x— x0）sinθ＋（ y—y0）cosθ］ 2 b 2 ＝1 （1） 展开整理后使用一种变量代换方法［6］ 将复杂 的非线性方程简化．变换后平面任意位置椭圆方程 式可以表示为如下形式： x 2＋ A xy＋By 2＋Cx＋ Dy＋ E＝0 （2） 这样求解椭圆系数的方程就是一组关于 A、 B、C、D 和 E 的线性方程．根据最小二乘法原理利 用回波信号峰值点附近的 N 个点的坐标求得椭圆 系数得到椭圆方程建立椭圆拟合的数学模型来 获取 A、B、C、D 和 E 这五个特征参数． 设定椭圆长轴与 x 轴的夹角θ＝0°或90°即椭 圆对称轴与坐标轴平行的情况下进行数据拟合可 把式（2）变换为： x 2＋ Ay 2＋Bx＋Cy＋ D＝0 （3） 这样就可由 A、B、C 和 D 来确定所拟合的椭圆．在 进行拟合时采用最小二乘法［7］ 依次确定式（3）中 的未知数 A、B、C 和 D．由此推导出椭圆的基本参 数：椭圆中心坐标（ x0y0）长轴半径 a短轴半 径 b． 依据上述最小二乘法的原理测量点与理想曲 线间误差为： φ（ ABCD）＝ ∑ N i＝1 （ x 2 i＋ Ay 2 i＋Bxi＋Cyi＋ D） 2 （4） 由数学分析可知使 A、B、C 和 D 满足导数为 零的点为极值点： ∂φ ∂A ＝ ∂φ ∂B ＝ ∂φ ∂C ＝ ∂φ ∂D ＝0 （5） 由此得下列方程组： ∑ N i＝1 y 4 i ∑ N i＝1 xiy 2 i ∑ N i＝1 y 3 i ∑ N i＝1 y 2 i ∑ N i＝1 xiy 2 i ∑ N i＝1 x 2 i ∑ N i＝1 xiyi ∑ N i＝1 xi ∑ N i＝1 y 3 i ∑ N i＝1 xiyi ∑ N i＝1 y 2 i ∑ N i＝1 yi ∑ N i＝1 y 2 i ∑ N i＝1 xi ∑ N i＝1 yi N A B C D ＝ ∑ N i＝1 x 2 iy 2 i ∑ N i＝1 x 3 i ∑ N i＝1 x 2 iyi ∑ N i＝1 x 2 i （6） 第11期 陈先中等： 超声测距系统的高精度中心椭圆算法 ·1155·
,1156 北京科技大学学报 第29卷 用线性方程组的算法可以求解得到A、B、C和D C2+AB2-4AD 的值 b= 4A2 (10) 由式(6)得到A、B、C和D的值后，可求得平 面上长轴与x轴夹角0=0的椭圆的另外四个 椭圆的横坐标中心x0=一号=6，时间：即为 参数： 超声波接收探头所接收到的回波信号的峰值时刻， x0≥、B 对应所测距离时接收信号电压值为b一y0= 2 (7) C y0= 2A (8) 牛AD+-V。如图4所示. 4A2 根据超声测距原理⑧]，待测距离为L,超声波 a C2+AB2-4AD 4A (9) 的声速为“，利用温度传感器测温来实现声速间接 100F (a) 65F 原始数据 (b) 原始数据 0 中心椭圆法 60 拟合 中心椭圆法 拟合 -100 55 -200 50 -300 45 -400 40 -30-20-10010203040 -15-10-505101520 图4使用中心椭圆法得到的椭圆(a)及与真实数据比较的局部放大图(b) Fig.4 Fllipse after using the central ellipse algorithm (a)and its amplified graph (b) 补偿，中心椭圆算法得到的超声波往返时间为t, 由下式可计算出所测距离0]： 3实验结果 L=ot/2 (11) 本文所述的超声波测距系统是基于 系统的算法流程如图5所示， ATmegal6L单片机的智能控制系统以及信号处理 算法完成的，控制系统主要完成发射脉冲串驱动信 超声波回波信号 号、对回波信号进行A/D转换、根据回波信号设定 零相位滤波 放大倍率、完成算法计算、进行温度补偿修正、协调 阈值法处理提取峰值左右的信号椭圆中心算法 输入输出（包括数据显示、电流输出）等工作，当发 射频率为50kz,每组包含五个脉冲，经过换能器驱 峰值左右的四个特征点对应的x,y值 动，驱动电压峰值V≥650V,对实际距离1,1.2， 计算立%立以龙以立以立好立 1.6,1.8和2m的液位进行测量，在算法上分别采 用了包络线法、三次多项式拟合法和中心椭圆算法 进行测量比较，比较结果如表1所示 从表1可以看出，中心椭圆算法与包络线法和 构成关于A、B、C、D的方程组 三次多项式法这两种传统的超声波测距算法相比， 接收到回波信号的时间t为-B2 测量精度大幅度提高 温度补偿后得到超声波速度为， 4结论 得到所测距离L=vt2 本文提出了新型超声测距算法，系统从硬件电 图5系统算法流程示意图 路的设计上减小了噪声的影响，提高了回波信号的 Fig-5 Flow diagram of the system 信噪比和稳定性，增强了系统抗干扰能力，提出的
用线性方程组的算法可以求解得到 A、B、C 和 D 的值． 由式（6）得到 A、B、C 和 D 的值后可求得平 面上长轴与 x 轴夹角θ＝0的椭圆的另外四个 参数： x0＝— B 2 （7） y0＝— C 2A （8） a＝ C 2＋ AB 2—4AD 4A （9） b＝ C 2＋ AB 2—4AD 4A 2 （10） 椭圆的横坐标中心 x0＝— B 2 ＝ t时间 t 即为 超声波接收探头所接收到的回波信号的峰值时刻 对应 所 测 距 离 时 接 收 信 号 电 压 值 为 b — y0＝ C 2＋ AB 2—4AD 4A 2 ＋ C 2A ＝ V m如图4所示． 根据超声测距原理［8］待测距离为 L超声波 的声速为 v 利用温度传感器测温来实现声速间接 图4 使用中心椭圆法得到的椭圆（a）及与真实数据比较的局部放大图（b） Fig．4 Ellipse after using the central ellipse algorithm （a） and its amplified graph （b） 补偿中心椭圆算法得到的超声波往返时间为 t ［9］ 由下式可计算出所测距离［10］： L＝vt／2 （11） 系统的算法流程如图5所示． 图5 系统算法流程示意图 Fig．5 Flow diagram of the system 3 实验结果 本 文 所 述 的 超 声 波 测 距 系 统 是 基 于 ATmega16L单片机的智能控制系统以及信号处理 算法完成的．控制系统主要完成发射脉冲串驱动信 号、对回波信号进行 A／D 转换、根据回波信号设定 放大倍率、完成算法计算、进行温度补偿修正、协调 输入输出（包括数据显示、电流输出）等工作．当发 射频率为50kHz每组包含五个脉冲经过换能器驱 动驱动电压峰值 V pp≥650V对实际距离11∙2 1∙61∙8和2m 的液位进行测量在算法上分别采 用了包络线法、三次多项式拟合法和中心椭圆算法 进行测量比较比较结果如表1所示． 从表1可以看出中心椭圆算法与包络线法和 三次多项式法这两种传统的超声波测距算法相比 测量精度大幅度提高． 4 结论 本文提出了新型超声测距算法系统从硬件电 路的设计上减小了噪声的影响提高了回波信号的 信噪比和稳定性增强了系统抗干扰能力．提出的 ·1156· 北 京 科 技 大 学 学 报 第29卷
第11期 陈先中等：超声测距系统的高精度中心椭圆算法 ,1157 表1距离测量数据结果 Table 1 Results of distance measurement 实际测量距离/m 包络线法/m 相对误差/% 三次多项式法/m 相对误差/% 椭圆中心法/m 相对误差/% 1.00 1.013 1.30 1.009 0.90 1.002 0.20 1.20 1.212 1.00 1.212 1.00 1.202 0.17 1.60 1.624 1.50 1.613 0.81 1.603 0.19 1.80 1.821 1.17 1.816 0.89 1.805 0.28 2.00 2.027 1.35 2.021 1.05 2.004 0.20 基于能量重心校正法和最小二乘法的中心椭圆算 2005,27(6):48 法，提高了超声测距的精度 [6]刘书桂，李蓬，那永林.基于最小二乘原理的平面任意位置椭圆 的评价.计量学报，2002,4(10)：245 参考文献 [7]孙朝明，伏德贵·几何量测量中的拟合算法.先进制造技术， 2006,25(1):25 [1]陈先中，王伟，张争，等.低噪声回波测距系统的电路设计与 [8]Willatzen M.Ultrasonic flow meters:temperature gradients and 系统分析.北京科技大学学报，2005,27(4)：488 [2]陈先中，张争。复杂工况下超声波液位测量系统的设计传感 transducer geometry effects.Ultrasonics.2003.14(2):105 器技术，2005(11)：41 [9]Andria G.Attivissimo F,Giaquinto N.Digital signal processing techniques for accurate ultrasonic sensor measurement.Measure- [3]陈先中，柳瑾.FMCW微波液位测量的快速频率估计算法.传 ment,2001,30(2):105 感技术学报，2005,18(4)：901 [10]张卫东，张圣训.高精度超声波测距系统·工业仪表与自动化 [4]丁康，江利旗，离散频谱的能量重心校正法.振动工程学报 装置，1994,4：55 2001,3(14):354 [5]陈凯，刘青一种随机化的椭圆拟合方法.计算机工程与科学， High-accurate central ellipse algorithm for an ultrasonic ranging system CHEN Xianzhong,HOU Qingwen,YE Lin,ZHA NG Hongyu Information Engineering School.University of Science and Technology Beijing.Beijing 100083.China ABSTRACI This paper described a kind of ultrasonic ranging system which employs a new driving circuit and low-noise programmable amplifier.The design of signal-to-noise ratio (SNR)in the circuit was raised to make measurement more stable and reliable.An efficient central ellipse algorithm for estimating the echo peak posi- tions of ultrasonic sampling was presented to find the accurate receive time of echo by searching the peak of ener- gy bursts,which is base on the theories of energy centrobaric correction method and least square method in the ultrasonic ranging system.A comparison of experimental results shows that the proposed algorithm s relative er- ror is approximately 0.2,which is more accuracy than the envelope detector or the cubic fitted multiple re- gression equation method,as a result it is suitable for higher accuracy industrial applications KEY WORDS distance measurement;ultrasonic:central ellipse algorithm;echo signal
表1 距离测量数据结果 Table1 Results of distance measurement 实际测量距离／m 包络线法／m 相对误差／％ 三次多项式法／m 相对误差／％ 椭圆中心法／m 相对误差／％ 1∙00 1∙013 1∙30 1∙009 0∙90 1∙002 0∙20 1∙20 1∙212 1∙00 1∙212 1∙00 1∙202 0∙17 1∙60 1∙624 1∙50 1∙613 0∙81 1∙603 0∙19 1∙80 1∙821 1∙17 1∙816 0∙89 1∙805 0∙28 2∙00 2∙027 1∙35 2∙021 1∙05 2∙004 0∙20 基于能量重心校正法和最小二乘法的中心椭圆算 法提高了超声测距的精度． 参 考 文 献 ［1］ 陈先中王伟张争等．低噪声回波测距系统的电路设计与 系统分析．北京科技大学学报200527（4）：488 ［2］ 陈先中张争．复杂工况下超声波液位测量系统的设计．传感 器技术2005（11）：41 ［3］ 陈先中柳瑾．FMCW 微波液位测量的快速频率估计算法．传 感技术学报200518（4）：901 ［4］ 丁康江利旗．离散频谱的能量重心校正法．振动工程学报 20013（14）：354 ［5］ 陈凯刘青．一种随机化的椭圆拟合方法．计算机工程与科学 200527（6）：48 ［6］ 刘书桂李蓬那永林．基于最小二乘原理的平面任意位置椭圆 的评价．计量学报20024（10）：245 ［7］ 孙朝明伏德贵．几何量测量中的拟合算法．先进制造技术 200625（1）：25 ［8］ Willatzen M．Ultrasonic flowmeters：temperature gradients and transducer geometry effects．Ultrasonics200314（2）：105 ［9］ Andria GAttivissimo FGiaquinto N．Digital signal processing techniques for accurate ultrasonic sensor measurement．Measurement200130（2）：105 ［10］ 张卫东张圣训．高精度超声波测距系统．工业仪表与自动化 装置19944：55 High-accurate central ellipse algorithm for an ultrasonic ranging system CHEN Xianz hongHOU QingwenY E L inZHA NG Hongyu Information Engineering SchoolUniversity of Science and Technology BeijingBeijing100083China ABSTRACT This paper described a kind of ultrasonic ranging system which employs a new driving circuit and low-noise programmable amplifier．The design of signa-l to-noise ratio （SNR） in the circuit was raised to make measurement more stable and reliable．An efficient central ellipse algorithm for estimating the echo peak positions of ultrasonic sampling was presented to find the accurate receive time of echo by searching the peak of energy burstswhich is base on the theories of energy centrobaric correction method and least square method in the ultrasonic ranging system．A comparison of experimental results shows that the proposed algorithm’s relative error is approximately 0∙2％which is more accuracy than the envelope detector or the cubic fitted multiple regression equation methodas a result it is suitable for higher accuracy industrial applications． KEY WORDS distance measurement；ultrasonic；central ellipse algorithm；echo signal 第11期 陈先中等： 超声测距系统的高精度中心椭圆算法 ·1157·
