专利名称：山区公路平面线形自动设计仿真系统的制作方法
技术领域：
本发明涉及道路线形设计与优化系统，具体涉及山区公路平面线形自动设计仿真系统，属道路设计技术领域。
目前在我国山区公路平面线形设计中，一直都是在大比例尺地形图上进行操作的[1]，采用的方法是传统放坡定线的直线型设计方法(又称导线法)。首先确定路线的整体布局，然后根据平均纵坡i平均(5.0％~5.5％，视相对高差而定)依次在等高线上作出坡度点(或称型值点)，调整一些不适当的点后连接这些坡度点构成导向线。导向线仍然是条折线，于是根据平面线形标准的要求，结合横坡变化情况，确定必须通过的点、适当照顾的点和可以忽略的点，采取以点连线、以线交点的方式，定出平面修正导向线，并量出转角，在图上敷出曲线。最后根据平面线形要求作线位调整再定出中线。纸上定线是一个反复检验的过程。必须通过多次的试线，才能获得理想的适宜线形。
传统的平面线形设计方法具有如下的局限性①坡度点都是由设计者用两脚规依次在每一条等高线上截点去确定。放坡到了终点的时候，如果放坡得到的高程达不到终点高程，又要根据剩余高差的大小重新确定坡度点，这个工作需要多次才能完成。同时，为了能充分利用有利地形，避让艰巨工程或严重不良地质地段，又需要重新确定平均纵坡，调整坡度点，如此反复进行。这些工作重复、繁琐、易错。
②在修正导向线，采取以点连线、以线交点的方式定出平面“试线”时，哪些点是必须通过的，哪些点是要适当照顾的，哪些点又是可以不考虑的，这并没有一个确定的标准可参考，完全依靠设计者平时积累的经验和自身丰富的空间想象能力，并且往往带有一定的盲目性和随意性，有时经过多次试线还无法确定出合适的线形。
③在导线绘出后，开始敷设平面曲线。由于受到地形和地质条件的约束，要取得比较合理的圆曲线半径、缓和曲线长度，使得设计出来的路线工程量少、与地形结合得好，又不破坏线形的整体平顺性，这是设计者必须考虑的而且是比较棘手的问题[2]。
④由于传统导线法是以直线作为主要的线形，曲线(圆曲线、缓和曲线)只是作为辅助线形，因此通过传统导线法设计出的线形一般难以和蜿蜒起伏的山区地形环境相协调，导致大挖大填，破坏了自然景观，造成公路的视觉效果较差。
⑤传统导线法难以处理复杂多变的几何线形。
而目前的曲线型设计方法[3]，[4]，仅是从某一个方面去寻求解决平面线形设计问题的途径，并未从地形、线形、设计方法、计算机程序和仿真技术等系统的角度实现山区平面线形的全自动化设计。同时未解决线形布设过程的线形单元和约束条件的自动描述和求解的问题。因而，实用性较差，往往难以推广应用。
本发明的目的在于克服上述现有技术之不足，提供一种解决传统方法存在人工重复劳动、易错、繁琐等缺陷的山区公路平面线形自动设计仿真系统。
实现上述目的，本发明的主要内容是山区公路平面线形自动设计仿真系统，其特征在于包括如下功能模块1)等高线绘制模块；通过建立数学地面模型，绘制出等高线地形图；3)确定型点值模块；用矩形格网三项次多项式高程内插的方法求出坡度点；3)拟合曲线曲率图转化模块；视这些坡度点为型值点，利用局部坐标三次样条曲线进行拟合，生成路线的平面拟合线形；4)拟合曲线曲率图的绘制模块；对路线的拟合线形每隔一定距离计算曲率，生成拟合线形的曲率图；对拟合线形的曲率图进行计算机自动识别转化模块，并用最小二乘法转化生成传统平面线形的曲率图；5)线元设计算法模块；利用积木法敷设平面线形，并进行线位计算。所述等高线绘制模块的建立包括1)形数据采集；2)形数据预处理；3)形数据排序与检索；4)字地面模型的高程内插，采用矩形格网内插法；5)裂线处理。所述地形数据排序是1)所有已知数对x方向进行排序；2)据原始数据的密度，自动确定格网宽度dc，得到x方向的总列数，从而可确定首列及每列x方向坐标原点；3)定每列的已知地形点(从x0(i)＜x＜x0(i)+dc区域内提取)；4)该列的已知地形数据，对y方向进行排列；5)据该列已知地形数据的分布范围确定该列格网个数，以及该列首格的局部坐标原点y0(i)；
6)对该列数据记录每个格网中的已知点首点序号及总点数(从y0(i)+j·dc＜x＜y0(i)+(j+1)·dc区域内提取)。
重复3)~6)，将数模区域内的所有已知点处理完毕后形成数模数据，在这个基础上，只要给出任意待定点的平面坐标，就可以快速检索出该点所处格网，从而快速提取该格网中已知地形点。所述最小二乘法转化生成传统平面线形的曲率图设计包括1)确定拟合曲线通过地带的地形地貌参数γ；2)根据样本曲率图确定曲率图的形状特征参数λ，根据向量λ中连续的0或1的个数去判定路线的转向，从而确定标准曲率图中曲率k的正负号；3)构造标准曲率图首先根据参数γ的数值组合，如出现不利情况，则增加考虑凸型、C型、复合型这三种特殊线形组合方式，否则只考虑基本型、S型和卵型这常用的三种线形组合方式。然后根据已知的公路等级和样本曲率图的里程长度，确定ls、Rave及ls和ly的最佳比值，构造出标准曲率图；4)对标准曲率图和样本曲率图的曲率K向量和里程L向量取值，并对曲率K向量进行“归一化”处理，运用公式计算各曲率图的C值；5)C值比较，得出结论。所述线元设计算法是将组合复杂的公路平面线形“化整为零”，分解成若干个线形单元。若已知路线平面曲线的起点信息如坐标、切线或法线方向和曲率半径，则从起点处开始设置任一单元，沿任何方向延伸，此单元终点的信息如坐标、切线或法线方位角和曲率半径都可以计算出来，同时将其作为下一单元起点的相同信息加以利用。如此逐个单元往下计算，如同搭积木一样，各单元首尾连接，构成一条连续完整的公路平面线形。其中线形单元即线元可以分别是直线、缓和曲线和圆曲线，必要时也可以是这三种线元组合而成的线形单元。
山区公路平面线形设计仿真系统是由各功能模块组成的。采用模块化程序设计方法有以下四个优点①采用模块化方法得到的程序结构清晰、层次分明；②采用模块化方法可以提高程序的开发效率。每个模块可以单独编译，程序模块通过精心设计之后可以为许多程序所共用；③采用模块化方法可以提高程序的可行性，而且使程序便于维护和修改；④采用模块化方法，各模块可以独立地进行测试或验证，这比对一个结构复杂的大型程序进行测试或验证要容易的多。
本发明的主要优点还在于通过建立平面线形设计的仿真系统，消除人工在大比例尺地形图上确定坡度点的重复劳动，使人工定点转为计算机自动识别；利用一定的边界约束条件对线形进行动态调整，使线形的调整做到直观、简便；尽可能地减少人工操作，不再需要设置导向线，线形布置与敷设完全自动化，使得平面线形设计快速、简捷，节省大量的人力、物力，缩短成图周期。山区公路平面线形设计的自动仿真系统，将设计者从重复的工作中解脱出来，把更多的精力投入到平面线形的优化工作中去，提高设计质量和设计水平。
下面结合附图和具体的实施方式对本发明作进一步说明。
图1是本发明等高线绘制流程图。
图2是本发明仿真系统中确定型值点的程序流程框图。
图3是局部坐标曲线图。
图4是拟合曲线曲率图转化流程图。
图5是积木法的基本原理图。
图6是线元设计算法流程图。
图7是曲率变化点处的方位角计算流程图。
图8是坡度点曲线拟合的公路线形。
图9是拟合曲线的曲率图。
图10是拟合后的曲率图。
图11是转化后的传统线形。
一、数字地面模型的建立地形资料是公路设计的重要基础资料，在传统设计中，一般用地形图或断面图来表示地形。利用计算机进行公路设计，就要让计算机能认识和处理地形资料，为此，必须把地形资料变成计算机能接受的信息—数字，数字地面模型就是在这种背景下被引入公路设计领域的。
数字地面模型(Digital Terrain Model，简称数模、DTM)是指按照某种数学模型表达地形特征的数值描述方式。它由许多规则或无规则排列的地形点三维坐标x、y、z组成，是数字化了的地形资料存储于计算机的产物。
建立数字地面模型，一般要经过地形数据采集、数据预处理、原始数据的排序与检索、待定点的高程内插等主要过程。为了提高数字地面模型的精度，增强对地形表面描述的能力，数模中必须对地物、地形断裂线进行有效的处理。
(1)地形数据采集数字地面模型原始数据的来源有多种，目前在实用中主要采用的有其中三种由航测仪器从航空照片上获得地形数据；从已有地形图上由数字化仪输入地形数据；由可记录量测数据的电子经纬仪、全站式速测仪等仪器从野外实测获得地形数据。
用航测方法获取数模原始数据主要有两种途径一是按象方测图思想，利用立体坐标量测仪直接量测像片的象点坐标x、y、p；另一是按物方测图思想，利用立体测图仪建立立体模型，测定模型坐标x、y、z。
对已有地形图进行数字化可以使用大幅面数字化仪或大幅面图形扫描仪。利用数字化仪对地形图采样时，一般是由数字化仪沿等高线获取地形点的平面坐标，该条等高线的高程值则由人工读取；对小比例尺地形图(1∶10000~1∶50000)，若等高线密集，沿线路走向一定宽度，可复印放大，以利于沿等高线采样；对地形独立点以及地物、断裂线等数据进行量测时，可先对各类数据编码约定。
在没有航片资料及地形图，不能用上述方法采集地形数据时，利用带记录装置的全站式速测仪在野外实测地形数据，这也是一种切实可行的数模数据采集方法。由于野外实测劳动强度大，工作环境差，效率低，不适宜采集过多的地形点，因而野外实测采集地形数据建立数模，一般只适应于路线较短的局部区域或局部工点测量，并用于对已知点数要求较少的三角网数模程序。
(2)地形数据预处理数据预处理的目的是将各种设备采集得到的数据格式各异的原始地形数据，转换成有统一的坐标系、统一的数据格式与编码形式的地形数据文件，并进行查错、改错以及数据压缩的处理。然后，根据路线大致走向，考虑到带状区域建立数字地面模型的需要和数模之间的接边，将整个地形原始数据经过旋转、平移到全线统一的数模坐标系中，供建立沿路线走向的带状数字地面模型连续调用。
对数模的计算速度而言，其数据结构尤为重要，数模的数据结构优劣直接影响到数模的效率，也是关系到数模能否应用于实际工程的关键。在数据预处理阶段，确定数模数据所采用的数据结构，也是其主要任务之一。研究数模数据结构及数据组织时，主要应考虑两个因素一是所建数模应占有最少的存储空间，减少数据冗余；二是在最少存储空间的基础上，数模应提供完整的数据信息，以利快速检索、内插和各种应用。
(3)地形数据排序与检索就排序算法而言，有插入排序、希尔排序、选择排序、堆排序、快速排序、归并排序、B+法排序、跳跃二分法排序、基数排序等等。对排序方法的选择，主要要尽可能地充分利用每次排序过程中的信息，有效地改善算法的复杂性，减少比较次数，并减少数据占有的临时存放空间，以此加快排序的速度。
从排序算法的分析可知，对一维数组排序比对二维数组的排序要快得多，所以可将数模数据二维排序转化为一维数组排序问题，分别对x、y方向排序。对数模数据的排序算法思想是①将所有已知数据对x方向进行排序；②根据原始数据的密度，自动确定格网宽度dc，得到x方向的总列数，从而可确定首列及每列x方向坐标原点x0(i)；③确定每列的已知地形点(从x0(i)＜x＜x0(i)+dc区域内提取)；④对该列的已知地形数据，对y方向进行排序；⑤根据该列已知地形数据的分布范围确定该列格网个数，以及该列首格的局部坐标原点y0(i)；⑥对该列数据记录每个格网中的已知点首点序号及总点数(从y0(i)+j·dc＜x＜y0(i)+(j+1)·dc区域内提取)。
重复③~⑥，将数模区域内的所有已知点处理完毕后形成数模数据，在这个基础上，只要给出任意待定点的平面坐标，就可以快速检索出该点所处格网，从而快速提取该格网中已知地形点。数模都必须具有对已知点快速检索、快速提取的能力，才能满足实用的需要。
(4)数字地面模型的高程内插对采集得到的地形原始数据，用一定的数学方法进行内插加密，这是建立数模的核心问题之一。高程内插的方法有多种移动曲面拟合法、加权平均值法、曲面求和法、联立分块多项式法、矩形格网内插法、三角形格网内插法。本课题采用矩形格网内插法，因此重点介绍，其余方法详见参考文献[7]。
矩形格网内插法是在矩形格网节点高程已知的情况下，内插各格网范围内待定点高程的方法，这些格网节点高程可以直接从摄影测量模型上量测得到，也可以由上述的某种内插方法求得。内插时在格网各单元内建立一个局部多项式表示该单元的地表面，然后由该多项式曲面确定该格网范围内待定点的高程。由于内插在单个格网中进行，待定点的高程主要受格网四个节点高程的影响，所以格网要有足够的密度，以便在每个格网单元内可用一个相对于各平面坐标的一次、二次或三次的独立多项式进行内插。
在矩形格网中内插，最常用的是双三次多项式中的几项或全部项。完全双三次多项式的表达式为Zp=1xx2x3a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a441yy2y3=XTAY---(1)]]>从(1)式可知要确定双三次多项式，就要确定系数阵中的16个系数。为此，除了单元中四个格网节点上的高程外，还需要已知三次曲面分别在X、Y方向上的斜率 以及在格网节点处表示曲面扭曲的混合斜率 这些斜率值可由相邻格网上节点的高程计算出来，具体方法是根据数学中近似计算原理，以一阶差分及二阶差分分别表示一阶偏导和二阶混合偏导。
由于有了这些边界条件，可以保证由各个局部多项式形成的整个内插表面是连续的和光滑的。所以采用16项的完全双三次多项式内插是此类方法中最完善的。
(5)断裂线处理断裂线是建立数字地面模型的重要数据，它对数模的高程内插精度有相当大的影响。能否有效地处理断裂线，是数模能否用于工程实际的关键。数模采用的内插方法不同，对断裂线的处理方法也不同。采用曲面拟合一类的内插方法时，如散点数模、方格网数模等，一般是沿断裂线分片逼近待插表面。为此在内插区域内，首先自动判别哪些点与待插点所连成的直线同断裂线相交，确定出与待插点异侧不再参与内插的地形数据点，并求出它们与断裂线的交点及内插出该点高程。对于三角网数模，应以各条断裂线作为三角形的边界进行联网。在断裂线附近，最佳三角形的条件已不存在，为此，在构网时，对每一新形成的三角形，均要与其断裂线点所连成的直线进行是否交叉的判断，并根据三角形与断裂的相互关系，确认新的三角形顶点。
数字地面模型建立后，即可绘制地形图。一般是在矩形方格网或三角网网格边上作高程内插得到等高点值，然后按一定的判断方向连接各个等高点。其流程图如图1所示。
二、坡度点(型值点)的确定从路线的纵断面来看，路线的走向不外乎这三种情况①路线上坡(i＞0)；②路线下坡(i＜0)；③路线水平(i＝0.3％)(为使排水顺畅，防止积水渗入路基影响其稳定性，故保留0.3％的纵坡)。如果路线水平(i＝0.3％)这一情况暂且不论，那么在确定坡度点时就只有两种情况路线上坡或路线下坡。
设计者在大比例尺地形图(1∶500~1∶2000)上设计平面线形时，根据路线的起终点标高和控制点的高程，判断出路线该走哪个方向以及在哪些路段该上坡，在哪些路段该下坡。然后将两脚规开度L＝Δh/i平均％(比例尺与地形图相同)，沿路线走向转动两脚规依次在等高线上截点前进，所得出的点即为原始坡度点。
因此，为了使计算机能够自动识别路线的走向，使计算机拟合出来的线形符合设计者的意图，就需要在通过数字地面模型(DTM)生成的等高线地形图上确定控制路线走向的控制点。这些控制点可以是路线必须通过的控制点，也可以是路基横断面填挖平衡的经济控制点，也可以是纯粹的虚拟控制点。确定了这些控制点之后，在编制仿真程序时，设置一个变量gradepoint(x＝0，y＝0)，每求出一个坡度点，就把该坡度点的x、y坐标赋值给变量gradepoint的x、y坐标，然后以点gradepoint为圆心，半径R＝Δh/i平均％作圆弧。同时，每隔一微小角度Δa(如Δa＝0.05弧度)进行高程内插(高程内插采用的方法是矩形格网内插，详见5.1节)，内插出来的高程点(x，y，z)如果满足以下两个方面的条件①求出的高程z和路线前进方向最邻近的已知的等高线高程之差小于给定的任意小的正数ε(如ε＝0.001)；②求出的高程点平面坐标(x，y)到路线前进方向最邻近的控制点平面坐标(x，y)的距离最短，则该点即为所求的坡度点。具体的程序流程框图如图2所示。
三、局部坐标系三次样条曲线数学模型的建立一系列的型值点确定之后，该用怎样的函数来进行插值或拟合得出来的效果才能达到最好？很自然地会想到三次样条函数或B样条函数，但要注意三次样条函数或B样条函数都有一个约束条件横坐标必须是单调的(单调上升或单调下降)。对于公路而言，峰回路转，公路线形的横坐标不可能是单调的，因此三次样条函数或B样条函数就不适合用来进行插值或拟合。于是采用局部坐标下三次样条函数数学模型。
局部坐标下三次样条函数是一种具有几何不变性、可适用于大挠度曲线插值的样条曲线。它分段由(小挠度的)三次多项式曲线构成，数学表达式简单统一，计算简便可靠，整体达到二阶连续，线型光顺，具有较好的保形性能，一般不会出现多余拐点，是一种理想的大挠度插值工具。它具有以下优点①采用局部坐标的形式，较符合公路线形设计的客观条件，不仅使拟合设计工作在没有坐标的情况下仍能照常使用，扩大了其应用的范围，而且使曲线的敷设工作也变得极为容易；②其数学表达式简单，形式统一，有利于计算机进行数值计算，且便于使用；③曲线具有整体大挠度，局部小挠度的特点，因而曲线偏离弦线的距离较小，便于设计人员对路线位置的控制和调整。
④数学表达式简单，形式统一，有利于计算机进行数值计算，且便于使用；(1)局部坐标系三次样条函数的公式推导曲线Ts是平面上有序点列{Pj}(j＝1，2，……，N)的局部三次样条，是指它满足二个条件①在每个局部坐标中，曲线是以弦线为自变量的三次多项式、且通过弦线的两端；②整条曲线及其切线、曲率在区域内连续。
如图3所示，设内节点Pj在相邻局部坐标(x，y)、(x，y)中的一、二阶导数分别记为mj(-)、Mj(-)与mj(+)、Mj(+)，坐标转角为φj，根据切线和曲率连续的要求则有(参见图3)mj(+)-mj(-)1+mj(+)mj(-)=tgφj---(2)]]>Mj(-)[1+mj2(-)]3/2=Mj(+)[1+mj2(+)]3/2---(3)]]>每个局部坐标中的三次插值式可写成y(x)=mj-1(+)x(lj-x)2lj2-mj(-)x2(lj-x)lj2---(4)]]>或y(x)=Mj-1(+)(lj-x)36lj+Mj(-)x36lj-Mj-1(+)6lj(lj-x)-Mj(-)6ljx---(4)]]>式中lj-弦线 的长度于是节点的一、二阶导数互有表达式Mj(-)=2mj-1(+)lj+4mj(-)lj,mj(-)=ljMj(-)3+ljMj-1(+)6---(5)]]>令记号 把式(2)、式(5)代入式(3)，化简整理，可得到以{mj}为参量的节点关系式-“m(+)关系式”λjmj-1+2mj+μjmj+1＝cj+Fj(7)式中λj=lj+1lj+lj+1]]>μj=ljlj+lj+1]]>cj＝2λjtgφj+μjtgφj+1Fj＝F1j+F2j+F3j(8)F1j=μjmj+1tgφj+1[mj+1-tgφj+1]1+mj+1tgφj+1]]>F2j=2λjmjtgφj[mj-tgφj]1+mjtgφj[cosφj+mjsinφj]3]]>F3j＝λj[mj-1+2(mj-tgφj)]{1-[cosφj+mjsinφj]3}(j＝2，…N-1)附加二端边界约束，式(7)构成N+1个未知参量{mj}的完整非线性方程组。顺便指出，对于非周期边界φN≡0。
如令记号 经过同样的推导可以分别得到节点的“m(-)关系式”λjmj-1+2mj+μjmj+1＝cj+Fj(10)式中cj＝-2μjtgφj-λjtgφj-1Fj＝F1j+F2j+F3jF1j＝-λjtgφj-1·mj-1·mj-1(11)F2j＝-2μjtgφj·mj·mj·[cosφj-mjsinφj]3F3j＝μj[2mj+mj+1+2tgφj]{1-[cosφj-mjsinφj]3}节点的“M(+)关系式”μjMj-1+2Mj+λjMj+1＝dj+gj(12)式中λj=lj+1lj+lj+1]]>μj=ljlj+lj+1]]>dj=-6tgφjlj+lj+1---(13)]]>gj=-6tgφjlj+lj+1mjm‾j+2μjMj{1-[cosφj-m‾jsinφj]3}]]>+λjMj+1{1-[cosφj+1-mj+1sinφj+1]3}节点的“M(-)关系式”μjMj-1+2Mj+λjMj+1＝dj+gj(14)式中dj=-6tgφjlj+lj+1]]>g‾j=-6tgφjlj+lj+1mjm‾j+μjM‾j-1{1-[cosφj-1+mj-1sinφj-1]3}---(15)]]>+2λjMj{1-[cosφj+mjsinφj]3}若用节点曲率{ρj}作为参量，由内节点切线连续条件式(2)，化简整理可得到“ρ关系式”μjρj-1+2ρj+λjρj+1＝dj+hj(16)式中dj=-6tgφjlj+lj+1]]>hj=-6tgφjlj+lj+1mjm‾j+μjρj-1{1-[1+mj-12]32}+2ρj{μj[1-(1+m‾j2)32---(17)]]>+λj[1-(1+mj2)32]}+λjρj+1{1-[1+m‾j+12]32}]]>节点的“m(+)关系式”、“m(-)关系式”、“M(+)关系式”、“M(-)关系式”和“ρ关系式”这五组关系式是相互等价的，可根据具体情况选用。
(2)局部坐标下三次样条函数边界条件的确定和一般三次样条一样，边界条件也有周期与非周期条件之分，但具体到公路线形，只有非周期边界条件。非周期条件对于“m(+)关系式”可考虑一般形式为2m1(+)+μ1m1(+)＝c1+F1(m1(+)，m2(+))(18)λNmN-1(+)+2mN(-)＝cN+FN(mN-1(+)，mN(-1))特别当二端分别给定切线、曲率时，式(18)中系数如表1、表2所示。
(3).局部坐标下三次样条函数的求解节点关系式加上边界条件，即构成了节点参量的完整方程组。笔者利用“m(+)关系式”来求解局部坐标下三次样条函数，下面说明整个插值求解过程。
假定边界条件是起、终点的切线已知，均为tgα1＝0，tgαN＝0则有μ1＝0，c1＝0，F1＝0，λN＝0，cN＝0，FN＝0。于是式(7)的系数矩阵A为 虽然在局部坐标系下上述方程是非线性方程，但这是一种特殊的非线性方程，采用简单的迭代就可以求出方程的解。首先解方程组A·m(0)＝c，求出m(0)，把m(0)代入F1j、F2j、F3j的表达式中得出F，然后解方程组A·m(k)＝c+F(m(k-1))即用追赶法解三对角线的迭代方法，直到相邻二次解误差‖m(k)-m(k-1)‖＜ε(ε为给定的任意小的正数)时即满足要求，再由式(2)求出mj(-)，代入式(4)便得到每个局部坐标系三次样条函数的具体公式。
最后在绘制样条函数时，要把局部坐标下的坐标通过二维转换转化为统一的坐标系，转换公式为xy=cosθsinθ-sinθcosθx′y′+x0′y0′---(19)]]>式中x-统一坐标系下的横坐标；y-统一坐标系下的纵坐标；x′-局部坐标系下的横坐标；y′-局部坐标系下的纵坐标；x′0，y′0-局部坐标系下坐标原点。
必须指出的是，笔者使用局部坐标系下的三次样条函数作为拟合曲线，只是把它作为一种计算的过渡工具，为后面的计算曲率，转化为传统线形提供方便，而并不特别注重它的路用特性。
四、拟合曲线曲率图的绘制用局部坐标系下的三次样条函数拟合坡度点之后，可以很明显地看出拟合出来的线形并不具有明显的传统线形的直线、缓和曲线、圆曲线，而且拟合曲线的路线桩号不便标注，线形指标无法直观形象地反映出来，对于具体的路中线的计算和施工放样是很不方便的。因此必须进行转化，转化的载体就是曲率图。由于局部坐标系下的三次样条函数在每一段都有明确的计算表达式，于是对每一段上的样条函数每隔一微小距离(如0.05m)求曲率，对应的弧长(即公路里程长度)都是从起点开始算起，这一点是要加与注意的。曲率的计算公式为k=y′′(1+y′′2)32---(20)]]>五、拟合曲线曲率图转化为传统平面线形曲率图公路传统平面线形的线形组合形式有六种基本型、S型、卵型、凸型、C型、复合型，这六种线形组合方式分别对应着六种组合规则的不同形状的曲率图。由拟合曲线通过计算曲率K绘制出来的曲率图并没有象传统线形的曲率图一样规则有序，无法和传统平面线形要素相对应，因此要求对样条曲线的曲率图进行转化，使其成为传统平面线形六种曲率图形的某一种，从而能够明确地反映出直线、缓和曲线和圆曲线之间的线形组合关系及线形要素值。该问题用数学语言描述为将基本型、S型、卵型、凸型、C型、复合型这六种曲率图看作是一个总体中的六个类，把拟合曲线的曲率图看作是一个样本，然后去判断样本应该属于哪一类。本文在研究过程中通过建立判别标准，采用模式识别分析的最小距离法去进行判断识别样本的类别，这是整个仿真系统的核心部分。
在具体讨论最小距离法之前，先引入两个参数λ和γ。参数λ定义为曲率图的形状 特征参数，采用向量的表现形式，即λ＝(λ1，λ2，……，λn)，λi的取值为为了既能够准确有效地表达曲率图的形状，又能够节省计算机的存贮空间，减少计算量，λi表现的是曲率图上特征点的状态。所谓特征点就是在曲率图上取相邻三个点a(La，Ka)，b(Lb，Kb)，c(Lc，Kc)(K代表曲率，L代表里程)，求出a和b两点连线的斜率tgθ1＝(Kb-Ka)/(Lb-La)及b和c两点连线的斜率tgθ2＝(Kc-Kb)/(Lc-Lb)，如果tgθ1和tgθ2的值的符号相反或tgθ1和tgθ2的值产生突变，则b点即为特征点。
引入参数λ的目的是由于公路路线有左转和右转之分，反映在曲率图上表现为曲率K的正值和负值。通过向量λ中连续的0或1的个数来判定路线的转向，从而为构造的标准曲率图中曲率K的正负值的确定提供依据。
参数γ定义为地形地貌参数，γ的取值为 于是路线通过地带的地形地貌形态就可以用γ的数值组合来表示。
引入参数γ的目的是《公路路线设计规范》(JTJ011-94)中规定凸型线形组合方式只有在路线严格受地形、地物限制时方可采用；复合型线形组合方式仅在地形或其它特殊原因限制时(互通式立体交叉除外)使用；C型线形组合方式只在特殊地形条件下方可使用。也就是说，一般情况下，在公路传统平面线形的六种组合方式中常用的只有基本型、S型、卵型这三种，如果要采用凸型、C型或复合型的线形组合方式，还需要在设计中特别加以说明。引入地形地貌参数γ后，只有当γ的数值组合达到不利情况时，才增加考虑这三种特殊线形组合方式，否则就不予考虑。这样就为参与比较的标准曲率图的选取提供判断依据，同时也可以减少计算量和计算时间，提高程序运行效率。
引入参数λ和γ后，接着讨论最小距离法。我们知道，缓和曲线的表达式为r·l＝A2，A为缓和曲线的放大或缩小系数。在长度l相同的条件下，不同的缓和曲线其半径r是不相同的，通过对A值的比较，可以知道缓和曲线的形状。笔者由此得到启发对于曲率图的判别问题，在路线里程长度确定的前提下，传统平面线形的六种曲率图和样本曲率图的不同之处在于曲率值。通过对曲率值进行标准化(采用“归一化”方法)，并和公路里程值相乘，得出某一无量纲的数值，对该数值进行比较即可判断曲率图的相似性程度大小。因此，不妨制订如下计算公式K·L＝C (21)式中K-曲率值组成的行向量。
K＝(k1，k2，……，kn)L-路线里程值组成的列向量。
L＝(l1，l2，ΛΛ，ln)TC-曲率图的相似系数，无量纲。
然而，传统平面线形的曲率图只是定性地表现了直线、缓和曲线和圆曲线之间的组合关系，并没有把线形组合关系定量化。而运用公式(21)却需要确定的K值和L值，才能得出具体的C值，才能进行曲率图的相似性比较。于是要先进行传统线形曲率图的定量化工作，该步骤暂且称为传统平面线形曲率图的标准化过程，构造出来的曲率图称为标准曲率图。所谓标准曲率图是指在已知公路等级和路线里程长度的条件下，构造出来的曲率图对该等级公路具有一定的代表性。在进行曲率图的相似性比较时，就以标准曲率图为模板和样本曲率图进行比较。
具体确定一个曲率图，必须知道三个要素缓和曲线长ls，圆曲线半径R和圆曲线长ly。由于《公路路线设计规范》(JTJ011-94)规定的各等级公路的ls和R指标值并不一致，要想使构造出来的标准曲率图对所有等级公路均适用是不可能的，因此要分别对各等级公路作出标准曲率图。为了避免ls和R取值的随意性，规定缓和曲线长ls取各等级公路的缓和曲线最小长度，圆曲线半径R取各等级公路的圆曲线一般最小半径Rave，圆曲线长ly需要确定了ls和ly的比值才能确定。那么ls和ly的比值取多少就要视样本曲率图中路线里程长度而定。
在构造标准曲率图时，要求满足以下两个要求①标准曲率图的里程长度要和样本曲率图的里程长度相等；②构造的六种传统平面线形标准曲率图面积要求相等，如有误差，误差应小于10％，目的是降低相似性的错判概率。
为了方便程序的编制，定义ls和ly的最差比值为1∶0.5，最佳比值为1∶2，这样在样本曲率图中里程长度给定的条件下，就可以得到标准曲率图。
标准曲率图构造出来后，就可以对曲率K和里程L进行取值。为了使标准曲率图和样本曲率图的里程L向量相等，规定在两者中里程L的取值方式均为L＝(10，20，30，……，ln)向量L中li+1-li＝10m，ln的数值为里程总长度。曲率K的取值方式为对于标准曲率图，在圆曲线段，k值为1/Rave，在缓和曲线段，k值通过线性插值求得；对于样本曲率图，ki为li对应的曲率值。
应该注意到，各标准曲率图和样本曲率图的向量K的 并不相等，这时要对曲率图K向量进行标准化，采用“归一化”方法，即曲率K向量中的每一个分量ki均除以 否则K和L相乘得出的C是毫无意义的，比较结果也很容易发生错判。从数学的角度解释为把曲率看成权，用不同权和的权向量乘以一个相同的向量，得出来的结果是无法比较的，比较也是毫无意义的。
最后，运用公式(21)分别计算出各标准曲率图的C值和样本曲率图的C值，然后求出样本曲率图C值和标准曲率图C值的差的绝对值，绝对值最小的数对应的标准曲率图即为样本曲率图的转化模板。
综上所述，用最小距离法判断曲率图模式的步骤概括如下①确定拟合曲线通过地带的地形地貌参数γ；②根据样本曲率图确定曲率图的形状特征参数λ，根据向量λ中连续的0或1的个数去判定路线的转向，从而确定标准曲率图中曲率k的正负号；③构造标准曲率图首先根据参数γ的数值组合，如出现不利情况，则增加考虑凸型、C型、复合型这三种特殊线形组合方式，否则只考虑基本型、S型和卵型这常用的三种线形组合方式。然后根据已知的公路等级和样本曲率图的里程长度，确定ls、Rave及ls和ly的最佳比值，构造出标准曲率图；④对标准曲率图和样本曲率图的曲率K向量和里程L向量取值，并对曲率K向量进行“归一化”处理，运用公式(21)计算各曲率图的C值；⑤C值比较，得出结论。
样本曲率图的转化模式确定后，即可进行曲率图的转化。曲率图的转化必须遵循三个原则①公路里程长度要求相等；②公路路线转向要求相同；③曲率图的面积(即为公路路线起、终点的切线方位角之差)要求相等。同时还要求满足《公路路线设计规范》(JTJ011-94)规定的线形设计的各个指标值。笔者采用最小二乘法拟合原理进行曲率图的拟合，最小二乘原理详见参考文献[19]。整个拟合曲线的曲率图的转化过程的程序流程图如图4所示。
六、传统平面线形的生成曲率图拟合完成之后，还需要根据曲率图还原成传统的直线、缓和曲线和圆曲线，采用的方法是积木法。积木法，亦称线元设计法，它是将组合复杂的公路平面线形“化整为零”，分解成若干个线形单元。若已知路线平面曲线的起点信息如坐标、切线或法线方向和曲率半径，则从起点处开始设置任一单元，沿任何方向延伸，此单元终点的信息如坐标、切线或法线方位角和曲率半径都可以计算出来，同时将其作为下一单元起点的相同信息加以利用。如此逐个单元往下计算，如同搭积木一样，各单元首尾连接，构成一条连续完整的公路平面线形。其中线形单元即线元可以分别是直线、缓和曲线和圆曲线，必要时也可以是这三种线元组合而成的线形单元。从上述分析不难发现，积木法多数情况下用于曲线要素值已知的公路初步设计和施工图设计阶段的平面线形的计算与敷设。
以直线、缓和曲线、圆曲线为线元，而且已知路线的起点信息如坐标、切线或法线方位角、曲率半径等以及整条线形的曲线要素值。根据积木法的定义和基本思路，积木法的基本原理如图5所示，线元设计的算法如图6所示。
(一).曲率变化点处的方位角用积木法来敷设平面线形，很重要的一点就是在敷设完上一个单元之后，必须知道下一个单元的切线方位角。因此就要先求出曲率变化点处的方位角的值。计算方位角的流程图如图5-7所示。每一线元敷设时从局部坐标转换到整体坐标的过渡矩阵为A=cos(π2-azimuthangle(i))-sin(π2-azimuthangle(i))sin(π2-azimuthangle(i))cos(π2-azimuthangle(i))]]>(二).每一线元上各点在整体坐标系下的坐标参见图5、6、7。
因为已知各曲率变化点处桩号stake(i)，故线元i的里程长度为L(i)＝stake(i+1)-stake(i)以曲率变化点i处切线为X轴，法线为Y轴建立局部坐标系，求线元上每一点在整体坐标系下的坐标，现以求曲率变化点i+1处的坐标为例，分别讨论接直线、圆曲线与缓和曲线时的情形。
①接直线局部坐标系下曲率变化点i+1处的坐标为 则整体坐标系下曲率变化点i+1的坐标为x(i+1)y(i+1)=Ax0y0+x(i)y(i)]]>②接圆曲线圆曲线半径R＝|1/k(i)圆心角a＝L(i)/R局部坐标系下曲率变化点i+1的坐标为 k(i)＜0 k(i)＞0则整体坐标系下的坐标为x(i+1)y(i+1)=Ax0y0+x(i)y(i)]]>③接缓和曲线1)完整的缓和曲线k(i)＝0时，k(i+1)处半径为R＝|1/k(i+1)其局部坐标系的坐标为 k(i+1)＜k(i) k(i+1)＞k(i)整体坐标系的坐标为x(i+1)y(i+1)=Ax0y0+x(i)y(i)]]>当k(i)＜0且k(i+1)＝0时，先求k(i+1)点在局部坐标x′oy′下的坐标 旋转角rotateangle＝π/2-(3π/2+|k(i)·l/2|)过渡矩阵A1=cos(rotateangle)-sin(rotateangle)sin(rotateangle)cos(rotateangle)]]>x0、y0在x′oy′下的坐标为x1y1=A1x0y0]]>再求得其在整体坐标系下的坐标为x(i+1)y(i+1)=Ax1y1+x(i)y(i)]]>同理可求得k(i)＞0且k(i+1)＝0时k(i+1)处的坐标。
2)非完整的缓和曲线当k(i)＜0且k(i+1)＜k(i)时，该线元为连接两同向圆曲线的非完整的缓和曲线。为了计算k(i+1)处坐标，先将非完整的缓和曲线还原成完整的缓和曲线，起点为o终点为oo″，如图8所示。
则在o′点处的半径r1＝|1/k(i)|，在o″点处的半径r2＝|1/k(i+1)|，缓和曲线长ls＝l+l0，oo′弧长|r2l/(r1-r2)|。先求o点在x′oy′下的坐标与方位角o′点在xoy坐标系中的坐标为 旋转角rotateangle＝|k(i)l0/2|，过渡矩阵A2=cos(rotateangle)-sin(rotateangle)sin(rotateangle)cos(rotateangle),]]>设o点在x′oy′坐标系中的坐标为(x0，y0)，因为00=A2x0′y0′+x0y0,]]>则有x0y0=00→A2x0′y0′=-A2x0′y0′,]]>o点在整体坐标系下的坐标为xzyz=Ax0y0+xiyi,]]>xoy在整体坐标系下的旋转角为rotateangle＝π/2-(azimuthangle(i)-|kil0/2|)，过渡矩阵A3=cos(rotateangle)-sin(rotateangle)sin(rotateangle)cos(rotateangle),]]>则终点o″在xoy坐标系下的坐标为 x、y在整体坐标系下的坐标为xi+1yi+1=A3xy+xzyz.]]>同理可求得其它情况下的非完整的缓和曲线终点k(i+1)处坐标。
3)反向缓和曲线当k(i)＜0，k(i+1)＞0时，可把它看成是由两段完整的缓和曲线构成，即k(i)＜0且k(0)＝0和k(i+1)＞0。此时，分别两次调用完整的缓和曲线子程序即可求得k(i+1)处的坐标。同理可求k(i)＞0且k(i+1)＜0的情况。
已知曲率图运用积木法进行线形设计的计算机算法到此全部介绍完毕，笔者根据上述算法开发出Matlab运算程序-Orthodoxline函数。四、实例分析某一山岭重丘三级公路，路线起点为如图8所示标有三角形的地点，然后需要越过一座山头，到达山背后的地点。图9中标注的圆点表示引导路线走向的控制点，这些控制点属于虚拟性质的控制点，星形点是通过矩形格网双三次多项式高程内插求出来的型值点(坡度点)(有些型值点与等高线不是靠的很近，这是因为计算机在形成等高线时造成的高程系统误差所致)，设计平均纵坡采用5.0％。图8所示的曲线是局部坐标系三次样条曲线拟合型值点得出来的线形，三次样条曲线的曲率图如图9所示。运用模式识别的最小距离法通过计算判断，得到样本曲率图的转化模板是基本型线形组合的曲率图，用最小二乘法拟合后转化得到的传统平面线形曲率图如图10所示，图11是运用积木法敷设出来的传统平面线形。
从实例可以看出，运用仿真技术进行平面线形设计时，不需要设置导向线，节省了大量的人力、物力，线形设计做到快速、简捷，实现了山区公路平面线形自动化设计。
路线特征点的里程桩号、坐标如表3所示。
实践表明，采用本方法建立的仿真系统，用于山区公路平面线形设计快速、方便，设计质量高，效果好，自动化程度高。既能满足技术标准，又能很好地结合地形和环境条件，减少公路路线设计对环境的破坏，是一种较好的公路自动设计仿真技术，具有广泛的推广应用价值。
参考文献[1]孙家驷编著 公路勘测设计 重庆大学出版社，1994年[2]朱照宏 陈雨人等编著 道路路线CAD同济大学出版社 1998年[3]吴国雄 王福建编著 公路平面线形曲线型设计方法 人民交通出版社[4]屠书荣 关于样条曲线应用于公路平面线形设计几个问题的研究 重庆交通学院学报，1994年表1
表2
路线的特征点、里程和坐标 表权利要求
1.山区公路平面线形自动设计仿真系统，其特征在于包括如下功能模块1)等高线绘制模块；通过建立数学地面模型，绘制出等高线地形图；2)确定型点值模块；用矩形格网三项次多项式高程内插的方法求出坡度点；3)拟合曲线曲率图转化模块；视这些坡度点为型值点，利用局部坐标三次样条曲线进行拟合，生成路线的平面拟合线形；4)拟合曲线曲率图的绘制模块；对路线的拟合线形每隔一定距离计算曲率，生成拟合线形的曲率图；对拟合线形的曲率图进行计算机自动识别转化模块，并用最小二乘法转化生成传统平面线形的曲率图；5)线元设计算法模块；利用积木法敷设平面线形，并进行线位计算。
2.如权利要求1所述的公路平面线形自动设计仿真系统，其特征在于所述等高线绘制模块的建立包括1)地形数据采集；2)地形数据预处理；3)地形数据排序与检索；4)数字地面模型的高程内插，采用矩形格网内插法；5)断裂线处理。
3.如权利要求1所述的山区公路平面线形自动设计仿真系统，其特征在于所述地形数据排序是1)将所有已知数对x方向进行排序；2)根据原始数据的密度，自动确定格网宽度dc，得到x方向的总列数，从而可确定首列及每列x方向坐标原点；3)确定每列的已知地形点(从x0(i)＜x＜x0(i)+dc区域内提取)；4)对该列的已知地形数据，对y方向进行排列；5)根据该列已知地形数据的分布范围确定该列格网个数，以及该列首格的局部坐标原点y0(j)；6)对该列数据记录每个格网中的已知点首点序号及总点数(从y0(i)+j·dc＜x＜y0(i)+(j+1)·dc区域内提取)。重复3)~6)，将数模区域内的所有已知点处理完毕后形成数模数据，在这个基础上，只要给出任意待定点的平面坐标，就可以快速检索出该点所处格网，从而快速提取该格网中已知地形点。
4.如权利要求1所述的公路平面线形自动设计仿真系统，其特征在于所述最小二乘法转化生成传统平面线形的曲率图设计包括1)确定拟合曲线通过地带的地形地貌参数γ；2)根据样本曲率图确定曲率图的形状特征参数λ，根据向量λ中连续的0或1的个数去判定路线的转向，从而确定标准曲率图中曲率k的正负号；3)构造标准曲率图首先根据参数γ的数值组合，如出现不利情况，则增加考虑凸型、C型、复合型这三种特殊线形组合方式，否则只考虑基本型、S型和卵型这常用的三种线形组合方式。然后根据已知的公路等级和样本曲率图的里程长度，确定ls、Rave及ls和ly的最佳比值，构造出标准曲率图；4)对标准曲率图和样本曲率图的曲率K向量和里程L向量取值，并对曲率K向量进行“归一化”处理，运用公式计算各曲率图的C值；5)C值比较，得出结论。
5.如权利要求1所述的公路平面线形自动设计仿真系统，其特征在于所述线元设计算法是将组合复杂的公路平面线形“化整为零”，分解成若干个线形单元。若已知路线平面曲线的起点信息如坐标、切线或法线方向和曲率半径，则从起点处开始设置任一单元，沿任何方向延伸，此单元终点的信息如坐标、切线或法线方位角和曲率半径都可以计算出来，同时将其作为下一单元起点的相同信息加以利用。如此逐个单元往下计算，如同搭积木一样，各单元首尾连接，构成一条连续完整的公路平面线形。其中线形单元即线元可以分别是直线、缓和曲线和圆曲线，必要时也可以是这三种线元组合而成的线形单元。
全文摘要
山区公路平面线形自动设计仿真系统，包括如下功能模块等高线绘制模块确定型点值模块；拟合曲线曲率图转化模块；拟合曲线曲率图的绘制模块；线元设计算法模块。本发明通过建立平面线形设计的仿真系统，消除人工在大比例尺地形图上确定坡度点的重复劳动，使人工定点转为计算机自动识别；利用一定的边界约束条件对线形进行动态调整，使线形的调整做到直观、简便尽可能地减少人工操作不再需要设置导向线，线形布置与敷设完全自动化，使得平面线形设计快速、简捷，节省大量的人力、物力、缩短成图周期。
文档编号E01C1/00GK1391005SQ01133800
公开日2003年1月15日 申请日期2001年12月30日 优先权日2001年12月30日
发明者吴国雄, 丘晓坚 申请人:吴国雄