CN105931199A - 基于综合投影变形最小标准减小任意带高斯投影变形方法 - Google Patents

Abstract

本发明公开了一种基于综合投影变形最小标准减小任意带高斯正形投影变形的方法，先确定所选边端点自然横坐标平均值Y_m，根据Y_m、测区最小横坐标值Y_min和测区最大横坐标值Y_max的关系求解Δl_W及Δl_E，进而确定出投影重心位置Y_g；然后根据投影重心位置Y_g，得到新的中央子午线至原子午线的距离，再通过大地反算，得到新中央子午线的大地经度M；最后将中央子午线设置到M，将测区坐标转换到该中央子午线M下进行投影计算，便可减小任意带高斯正形投影变形。同时还确定了基于综合投影变形最小标准的本发明方法的最大适用范围。采用本发明方法，以更好满足对于投影变形均匀性要求较高的工程。

Description

基于综合投影变形最小标准减小任意带高斯投影变形方法

技术领域

本发明属于地理测绘科学技术领域，具体涉及一种基于综合投影变形最小标准减小任意带高斯正形投影变形的方法。

背景技术

受高斯正形投影的影响，在远离计算所选择的中央子午线的位置，会产生较大的投影变形，而各种工程建设对投影变形均有具体要求，很多工程要求在整个测区投影变形尽可能的均匀，以兼顾各方面的需求，对于此类工程应当以最小二乘准则下的综合投影变形最小化准则进行投影。

任意带高斯正形投影是通过测区所选定的投影高程面重新定义测区中央子午线位置，使得归算至高斯投影面上的变形刚好可以抵消边长归算至参考椭球面上的变形。其关键是中央子午线至测区相对位置的确定，现有参数确定中一般将参考位置选择在测区中央，这是一个近似的处理方法，这个位置的选择并没有顾及到投影的最优准则问题。研究者提出，对投影变形均匀性要求较高的工程应当以最小二乘准则下的综合投影变形最小化准则进行投影，采用现有方法并不能使投影变形在该准则下达到最优化。因此需要研究建立基于该标准要求的任意带高斯正形投影归算方法，以满足对投影变形均匀性要求较高的工程需求。

发明内容

本发明的目的是提供一种基于综合投影变形最小标准减小任意带高斯正形投影变形的方法，该方法基于综合投影变形最小标准确定了任意带高斯正形投影的参数，进而解决了现有投影方法不能满足对投影变形均匀性要求较高的工程测量需求的问题。

本发明所采用的技术方案是，一种基于综合投影变形最小标准减小任意带高斯正形投影变形的方法，包括以下步骤：

步骤1，确定中央子午线至测区投影重心位置距离Y′_m

通过公式确定中央子午线至测区投影重心位置的距离Y′_m，式中，H_m为边长归算所选投影面高出参考椭球面的平均高程，R＝6371km；

步骤2，确定测区投影重心位置Y_g

首先以投影变形平方和[ΔS_iΔS_i]＝min确定所选边端点自然横坐标平均值Y_m，得到Y′_m、测区最小横坐标值Y_min和测区最大横坐标值Y_max的关系，根据三者之间的关系求解测区投影重心位置Y_g以东及以西的测区范围值Δl_E及Δl_W：

${\mathrm{\Δl}}_E=\frac{-3(2Y_m^{\′}-L)+\sqrt{36{Y_m^{\′}}^2-3L^2}}{6}$

Δl_W＝L-Δl_E

其中：L＝Y_max-Y_min。

将Δl_E和Δl_W分别代入Y′_max＝Y′m+Δl_E和Y′_min＝Y′_m-Δl_W，求出测区在新的投影带中的位置，新坐标与老坐标中横坐标变化值ΔY＝Y′_max-Y_max，或者则测区重心位置横坐标Y_g在源平面坐标系中为Y_g＝Y_m+ΔY-(Y_m-Y′_m)，则有Y_g＝Y′_m+ΔY，由此确定出投影重心位置Y_g；

步骤3，测区新的中央子午线位置确定

根据投影重心位置Y_g，得到新的中央子午线至原子午线的距离为Y_g-Y′_m＝ΔY，通过大地反算，得到新中央子午线的大地经度M＝l+M₀，其中，M₀为原中央子午线经度，l通过下式计算，

$\begin{array}{c}l=\frac{\mathrm{\Δ}Y}{N_f{\mathrm{cosB}}_f}-\frac{{\mathrm{\ΔY}}^3}{6N_f^3{\mathrm{cosB}}_f}(1+2t_f^2+{\mathrm{\η}}_f^2)\\ +\frac{{\mathrm{\ΔY}}^5}{120N_f^5{\mathrm{cosB}}_f}(5+28t_f^2+24t_f^4+6{\mathrm{\η}}_f^2+8{\mathrm{\η}}_f^2{t}_f^2);\end{array}$

其中，ΔD_lim为规范规定的限差；Y_m为所选边端点自然横坐标平均值，R_m＝6371km，D₁为水平距离归算到参考椭球面上的测距边长。

步骤4，将原投影带平面坐标转换至新投影带

根据步骤3所确定的M值，将中央子午线设置到M，将测区坐标转换到该中央子午线M下进行投影计算，便可基于最大投影变形最小标准减小任意带高斯正形投影变形。

步骤1公式的推导过程如下：

测距边长水平距离为D，水平距离归算至测区某一高程面H_p的边长变形值为ΔD₀，水平距离归算至参考椭球面上的变形ΔD₁，椭球面上的边长归算至高斯投影面上的变形ΔD₂；

${\mathrm{\ΔD}}_0=\frac{D}{R}(H_p-H_m)$

${\mathrm{\ΔD}}_1=-\frac{D}{R}{H}_m$

${\mathrm{\ΔD}}_2=D_1(\frac{Y_m^2}{2R_m^2}+\frac{{\mathrm{\ΔY}}^2}{24R_m^2})$

上面两式中：H_m为边长归算所选投影面高出参考椭球面的平均高程，Y_m为所选边端点自然横坐标平均值，R_m系参考椭球面选定边长中点平均曲率半径，R与R_m取6371km；

水平距离归算到参考椭球面上的测距边长D₁＝D+ΔD₁，参考椭球面上的测距边投影到高斯平面上的长度D₂＝D₁+ΔD₂，于是边长综合变形为：

$\mathrm{\Δ}D={\mathrm{\ΔD}}_1+{\mathrm{\ΔD}}_2=\frac{D}{2R^2}({Y_m}^2-2{\mathrm{RH}}_m)$

则推出归算边高斯投影变形量抵偿值，令Y₀点处的ΔD＝0，得：

$H_m^{\′}=\frac{Y_m^2}{2R}$

H′_m为归算边高斯投影变形量抵偿值，则抵偿投影面的高程为H_m-H′_m；

由于Y′_m的选择需保证ΔD＝0，则存在如下关系：

$D\frac{{Y_m^{\′}}^2}{2R^2}=D\frac{H_m}{R}$

得出：

$Y_m^{\′}=\sqrt{H_{m}2R}$

通过公式即可确定中央子午线至测区投影重心位置的距离Y′_m。

上述步骤2中所选边端点自然横坐标平均值Y_m的确定方法为：

测区综合投影变形为：

$\mathrm{\Δ}D={\mathrm{\ΔD}}_1+{\mathrm{\ΔD}}_2=(\frac{Y_m^2}{2R_m^2}-\frac{H_m}{R})S$

取Y₀处的ΔD＝0，则：

$H_0=-\frac{Y_0^2}{2R}$

H₀为归算边高出抵偿投影面的平均高程，则抵偿投影面的高程为H_m-H₀。

令测区内的

取[ΔD_iΔD_i]＝min得：

$\[(\frac{Y_i^2}{2R^2}-\frac{Y_0^2}{2R^2})D_i(\frac{Y_i^2}{2R^2}-\frac{Y_0^2}{2R^2})D_i\]=min$

整理得：

$Y_0^4-2Y_i^2{Y}_0^2+Y_i^4=min$

根据上式求解出满足条件[ΔD_iΔD_i]＝min的Y₀

$Y_0=\frac{1}{\sqrt{3}}\sqrt{(Y_{\max }^2+Y_{\max }{Y}_{\min }+Y_{\min }^2)}$

由于：Y₀＝Y′_m，Y′_max＝Y′_m+Δl_E，Y′_min＝Y′_m-L+Δl_E，修改上式得到Y_m：

${Y_m^{\′}}^2=\frac{1}{3}(Y_{max}^{\′2}+Y_{\max }^{\′}{Y}_{\min }^{\′}+Y_{\min }^{\′2})$

上述步骤2中Δl_W及Δl_E的求解方法为：

所选边端点自然横坐标平均值Y_m的计算为

${Y_m^{\′}}^2=\frac{1}{3}(Y_{max}^{\′2}+Y_{\max }^{\′}{Y}_{\min }^{\′}+Y_{\min }^{\′2})$

令式中Y′_m-L＝b，Y′_max＝Δl_E+Y′_m，Y′_min＝Δl_E+b，则有

3Y′_m ²＝(Δl_E+Y′_m)²+(Δl_E+Y′_m)(Δl_E+b)+(Δl_E+b)²

由于Y′_m及_b已知，则由上式可求出Δl_E的值

${\mathrm{\Δl}}_E=\frac{-3(2Y_m^{\′}-L)+\sqrt{36{Y_m^{\′}}^2-3L^2}}{6}$

至此，解算出Δl_E；

将Δl_E代入下式求出Δl_W：

Δl_W＝L-Δl_E

本发明方法的最大适用范围为

Y′_min＝Y_m-Δl'_W至Y′_max＝Y_m+Δl'_E

其中，Δl'_E为在做抵偿投影其东侧允许的最大适用范围，Δl'_W为在Y_m做抵偿投影其西侧允许的最大适用范围，

${\mathrm{\Δl}}_E^{\′}=\sqrt{2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}-Y_m$

${\mathrm{\Δl}}_W^{\′}=-\sqrt{-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}+Y_m;$

其中，ΔD_lim为规范规定的限差；Y_m为所选边端点自然横坐标平均值，R_m＝6371km，D₁为水平距离归算到参考椭球面上的测距边长。

本发明方法的最大适用范围的确定方法为：

所述方法投影所采用横坐标位置为Y′_m，该位置高斯投影计算公式：

${\mathrm{\ΔD}}_M=D_1\left(\frac{Y_m^2}{2R_m^2}\right)$

当采用抵偿投影后，该点位置变形量便为0，基于Y_m列出Y_m以东任意位置的高斯投影计算公式：

${\mathrm{\ΔD}}_l=D_1\left(\frac{{(Y_m+{\mathrm{\Δl}}_E^{\′})}^2}{2R_m^2}\right)$

式中Δl为偏离坐标位置，对上面两式求差得到：

${\mathrm{\ΔD}}_E=D_1\left(\frac{{(Y_m+{\mathrm{\Δl}}_E^{\′})}^2}{2R_m^2}\right)-D_1\left(\frac{Y_m^2}{2R_m^2}\right)=D_1\left(\frac{2{\mathrm{\Δl}}_E^{\′}\·Y_m+{{\mathrm{\Δl}}_E^{\′}}^2}{2\·R_m^2}\right)$

按照规范规定限差要求为ΔD_lim，将上式ΔD_E用ΔD_lim替代求解出在Y_m做抵偿投影其东侧允许的最大适用范围Δl'_E：

$\begin{array}{c}{\mathrm{\ΔD}}_{\lim }=D_1\left(\frac{2{\mathrm{\Δl}}_E^{\′}\·Y_m+{{\mathrm{\Δl}}_E}^2}{2\·R_m^2}\right)\\ \⇒{{\mathrm{\Δl}}_E^{\′}}^2+2{\mathrm{\Δl}}_E^{\′}\·Y_m+Y_m^2=2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2\\ \⇒{\mathrm{\Δl}}_E^{\′}=\sqrt{2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}-Y_m\end{array}$

然后求解西测允许范围为Δl'_W，将下式ΔD_W用ΔD_lim替代：

${\mathrm{\ΔD}}_W=D_1\left(\frac{Y_m^2}{2R_m^2}\right)-D_1\left(\frac{{(Y_m-{\mathrm{\Δl}}_W^{\′})}^2}{2R_m^2}\right)=D_1\left(\frac{-({{\mathrm{\Δl}}_W^{\′}}^2-2{\mathrm{\Δl}}_W^{\′}\·Y_m)}{2\·R_m^2}\right)$

则求解出Δl'_W，

${\mathrm{\Δl}}_W^{\′}=-\sqrt{-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}+Y_m$

根据规范对最大投影变形量ΔD_lim的限差要求，确定抵偿投影最大适用范围为：

Y′_min＝Y_m-Δl'_W至Y′_max＝Y_m+Δl'_E

本发明的有益效果是，对投影变形均匀性要求较高的工程，本发明方法基于综合投影变形最小标准的方法可以保证测区综合投影变形在最小二乘准则下达到最小化。

具体实施方式

下面结合具体实施方式对本发明作进一步的详细说明，但本发明并不限于这些实施方式。

本发明提供了一种基于综合投影变形最小标准减小任意带高斯正形投影变形的方法，包括以下步骤：

步骤1，确定中央子午线至测区投影重心位置距离Y′_m

任意带高斯正形投影，首先需要确定出中央子午线至测区参考位置的距离。其确定方法如下：

测距边长水平距离为D，水平距离归算至测区某一高程面H_p的边长变形值为ΔD₀，水平距离归算至参考椭球面上的变形ΔD₁，椭球面上的边长归算至高斯投影面上的变形ΔD₂。

${\mathrm{\ΔD}}_0=\frac{D}{R}(H_p-H_m)$

${\mathrm{\ΔD}}_1=-\frac{D}{R}{H}_m$

${\mathrm{\ΔD}}_2=D_1(\frac{Y_m^2}{2R_m^2}+\frac{{\mathrm{\ΔY}}^2}{24R_m^2})$

上式中：H_m为边长归算所选投影面高出参考椭球面的平均高程，Y_m为所选边端点自然横坐标平均值，R_m系参考椭球面选定边长中点平均曲率半径，实用中一般采用地球平均曲率半径，R与R_m一般皆取为6371km。

水平距离归算到参考椭球面上的测距边长D₁＝D+ΔD₁，参考椭球面上的测距边投影到高斯平面上的长度D₂＝D₁+ΔD₂。

于是边长综合变形为：

$\mathrm{\Δ}D={\mathrm{\ΔD}}_1+{\mathrm{\ΔD}}_2=\frac{D}{2R^2}({Y_m}^2-2{\mathrm{RH}}_m)$

则推出归算边高斯投影变形量抵偿值，令Y₀点处的ΔD＝0，得：

$H_m^{\′}=\frac{Y_m^2}{2R}$

H′_m就是归算边高斯投影变形量抵偿值。则抵偿投影面的高程为H_m-H′_m。

Y′_m的选择需保证ΔD＝0，则存在如下关系：

$D\frac{{Y_m^{\′}}^2}{2R^2}=D\frac{H_m}{R}$

得出：

$Y_m^{\′}=\sqrt{H_{m}2R}$

由此确定出中央子午线至测区投影重心位置的距离Y′_m。

依据Y′_m重新定义测区中央子午线位置，并对测区坐标进行换带计算。

在习惯上一般认为Y′_m表达的是测区中央位置与所定义中央子午线的距离，当采用抵偿投影后，该位置变形量便为0，也就是说，当采用抵偿投影后，测区其他位置的投影变形量均是相对Y_m的。

步骤2，确定测区投影重心位置Y_g

首先以投影变形平方和[ΔS_iΔS_i]＝min确定Y_m，与Y_min及Y_max的关系。

测区综合投影变形：

$\mathrm{\Δ}D={\mathrm{\ΔD}}_1+{\mathrm{\ΔD}}_2=(\frac{Y_m^2}{2R_m^2}-\frac{H_m}{R})S$

取Y₀处的ΔD＝0，则：

$H_0=-\frac{Y_0^2}{2R}$

H₀为归算边高出抵偿投影面的平均高程，则抵偿投影面的高程为H_m-H₀。

测区内的

取[ΔD_iΔD_i]＝min得：

$\[(\frac{Y_i^2}{2R^2}-\frac{Y_0^2}{2R^2})D_i(\frac{Y_i^2}{2R^2}-\frac{Y_0^2}{2R^2})D_i\]=min$

整理得：

$Y_0^4-2Y_i^2{Y}_0^2+Y_i^4=min$

展开得：

$Y_0^4-\frac{2(Y_1^2+Y_2^2+...+Y_n^2)}{n}{Y}_0^2+\frac{(Y_1^4+Y_2^4+...+Y_n^4)}{n}=min$

(Y_i∈[Y_min,Y_max],n→∞)

设：

对Y求导，导数为零时,使上式得到最小值,得

$Y_0^2=\frac{Y_1^2+Y_2^2+...+Y_n^2}{n},(n\→\mathrm{\∞})$

整理得

$Y_0^2=\frac{1}{Y_{max}-Y_{min}}\[Y_1^2\frac{Y_{max}-Y_{min}}{n}+Y_2^2\frac{Y_{max}-Y_{min}}{n}+...+Y_n^2\frac{Y_{max}-Y_{min}}{n}\]$

$Y_0^2=\frac{1}{Y_{max}-Y_{min}}{\∫}_{y_{min}}^{y_{\max }}{y}^{2}dy$

$Y_0^2=\frac{1}{3}(Y_{max}^2+Y_{max}{Y}_{min}+Y_{min}^2)$

$Y_0=\frac{1}{\sqrt{3}}\sqrt{(Y_{\max }^2+Y_{\max }{Y}_{\min }+Y_{\min }^2)}$

上式就是满足条件[ΔD_iΔD_i]＝min的Y₀求解公式，而在这里我们需要确定的是测区的Y′_min及Y′_max。

由于：Y₀＝Y′_m，Y′_max＝Y′_m+Δl_E，Y′_min＝Y′_m-L+Δl_E，上式可修改为：

${Y_m^{\′}}^2=\frac{1}{3}(Y_{max}^{\′2}+Y_{max}^{\′}{Y}_{\min }^{\′}+Y_{\min }^{\′2})$

Y′_m、L均已知，未知量是Δl_E，需要解算出Δl_E。

令式中Y′_m-L＝b，Y′_max＝Δl_E+Y′_m，Y′_min＝Δl_E+b，则有

3Y′_m ²＝(Δl_E+Y′_m)²+(Δl_E+Y′_m)(Δl_E+b)+(Δl_E+b)²

＝3Δl_E ²+3(Y′_m+b)Δl_E+Y′_m ²+b²+Y′_mb

＝3Δl_E ²+3(Y′_m+b)Δl_E+(Y′_m+b)²-Y′_mb

移项得：3Δl_E ²+3(Y′_m+b)Δl_E+(Y′_m+b)²-(Y′_mb+3Y′_m ²)＝0

推出：3Δl_E ²+3(Y′_m+b)Δl_E+Y′_mb+b²-2Y′_m ²＝0

由于Y′_m-L＝b，代入上式，得到

$\begin{array}{c}3{{\mathrm{\Δl}}_E}^2+3(Y_m^{\′}+(Y_m^{\′}-L\left)\right){\mathrm{\Δl}}_E+Y_m^{\′}(Y_m^{\′}-L)+{(Y_m^{\′}-L)}^2-2{Y_m^{\′}}^2=0\\ \⇒3{{\mathrm{\Δl}}_E}^2+3(2Y_m^{\′}-L){\mathrm{\Δl}}_E-3Y_m^{\′}L+L^2=0\end{array}$

$\begin{array}{c}\mathrm{\Δ}=9{(2Y_m^{\′}-L)}^2-12(L^2-3Y_m^{\′}L)\\ =9(4{Y_m^{\′}}^2-4Y_m^{\′}L+L^2)-12L^2+36Y_m^{\′}L\\ =36{Y_m^{\′}}^2-3L^2\≥36{\left(\frac{L}{2}\right)}^2-3L^2=6L^2>0\end{array}$

Y′_m、Y₀、L均已知，由于L∈(0,∞)，Y′_m∈(0,∞)，同时，知Δ＞0，所以方程有两个实根：

即：

因Δl_E＞0则：

${\mathrm{\Δl}}_E=\frac{-3(2Y_m^{\′}-L)+\sqrt{36{Y_m^{\′}}^2-3L^2}}{6}$

至此，解算出Δl_E；

由于子午线移动前Y_min与Y_max已知，其东西跨度为L＝Y_max-Y_min。

将Δl_E代入下式求出Δl_W：

Δl_W＝L-Δl_E

将Δl_E和Δl_W分别代入Y′_max＝Y′_m+Δl_E和Y′_min＝Y′_m-Δl_W，便可求出测区在新的投影带中的位置，新坐标与老坐标中横坐标变化值ΔY＝Y′_max-Y_max，或者那么测区重心位置横坐标Y_g在源平面坐标系中为Y_g＝Y_m+ΔY-(Y_m-Y′_m)，则

Y_g＝Y′_m+ΔY

由此确定出投影重心位置，在使用中，所定义中央子午线至测距距离Y′_m系相对于投影重心位置Y_g。

步骤3，测区新的中央子午线位置确定

根据以上确定好的投影重心位置Y_g，就可确定出新的中央子午线至原子午线的距离为Y_g-Y′_m＝ΔY。

通过大地反算，便可确定出新中央子午线的大地经度，大地经度计算公式如下：

$\begin{array}{c}l=\frac{\mathrm{\Δ}Y}{N_f{\mathrm{cosB}}_f}-\frac{{\mathrm{\ΔY}}^3}{6N_f^3{\mathrm{cosB}}_f}(1+2t_f^2+{\mathrm{\η}}_f^2)\\ +\frac{{\mathrm{\ΔY}}^5}{120N_f^5{\mathrm{cosB}}_f}(5+28t_f^2+24t_f^4+6{\mathrm{\η}}_f^2+8{\mathrm{\η}}_f^2{t}_f^2)\end{array}$

其中：B_f为底点纬度，也就是当x＝X时的子午线弧长所对应的纬度。

按照子午线弧长公式：迭代进行计算；

$N_f=a{(1-e^2{\sin }^2{B}_f)}^{-\frac{1}{2}};$

$M_f=a(1-e^2){(1-e^2{\sin }^2{B}_f)}^{-\frac{3}{2}}$

t_f＝tan B_f；

${\mathrm{\η}}_f^2=e^{\′2}{\cos }^2{B}_f$

由此确定的新的中央子午线经度为：

M＝l+M₀

其中M₀为原中央子午线经度，由此便可确定出测区新的中央子午线位置M。

将中央子午线设置到M，将测区坐标转换到该中央子午线下进行投影计算，实现最大投影变形最小化准则下的任意带高斯正形投影

步骤4，将原投影带平面坐标转换至新投影带

步骤3确定新的中央子午线M后，由于原平面坐标是相对于原中央子午线M₀计算而来，就需要将原投影带M₀下的平面坐标(x₀、y₀)转换到新的投影带M下的平面坐标(x、y)。

方法是先根据原投影带的平面坐标(x₀、y₀)和中央子午线的经度L₀。按高斯投影坐标反算公式求得大地坐标(B、L)，然后根据(B、L)和新投影带中的中央子午线经度M，按高斯投影坐标正算公式求得在新投影带中的平面坐标(x、y)。

至此完成了中央子午线至测区投影重心位置距离Y′_m的确定，测区投影重心位置的确定，测区新的中央子午线位置确定以及投影换带计算，测区的最远端投影变形达到最小化，可满足对投影变形均匀性要求较高的工程的需求。

本发明还确定了该方法的最大适用范围，具体方法如下：

步骤4在新的中央子午线下进行投影，使得成果在最大投影变形最小化准则下达到最优投影，步骤1至步骤4已经确定出了完整的任意带高斯正形投影流程。确定好投影方法后，就需要明确基于该投影方法的最大适用范围。

该方法投影所采用横坐标位置为Y′_m，该位置高斯投影计算公式为：

${\mathrm{\ΔD}}_M=D_1\left(\frac{Y_m^2}{2R_m^2}\right)$

当采用抵偿投影后，该点位置变形量便为0，也就是说，当采用抵偿投影后，测区其他位置的投影变形量均是相对Y_m的。基于Y_m列出Y_m以东任意位置的高斯投影计算公式为：

${\mathrm{\ΔD}}_l=D_1\left(\frac{{(Y_m+{\mathrm{\Δl}}_E^{\′})}^2}{2R_m^2}\right)$

式中Δl'_E为投影重心以东允许的最大适用范围，我们对上面两式求差：

${\mathrm{\ΔD}}_E=D_1\left(\frac{{(Y_m+{\mathrm{\Δl}}_E^{\′})}^2}{2R_m^2}\right)-D_1\left(\frac{Y_m^2}{2R_m^2}\right)=D_1\left(\frac{2{\mathrm{\Δl}}_E^{\′}\·Y_m+{{\mathrm{\Δl}}_E^{\′}}^2}{2\·R_m^2}\right)$

按照规范规定限差要求为ΔD_lim，将上式ΔD_E用ΔD_lim替代，则有：

$\begin{array}{c}{\mathrm{\ΔD}}_{\lim }=D_1\left(\frac{2{\mathrm{\Δl}}_E^{\′}\·Y_m+{{\mathrm{\Δl}}_E}^2}{2\·R_m^2}\right)\\ \⇒{{\mathrm{\Δl}}_E^{\′}}^2+2{\mathrm{\Δl}}_E^{\′}\·Y_m+Y_m^2=2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2\\ \⇒{\mathrm{\Δl}}_E^{\′}=\sqrt{2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}-Y_m\end{array}$

Δl'_E就是所求出的在Y_m做抵偿投影其东侧允许的最大适用范围。设定西测允许范围为Δl'_W，则：

${\mathrm{\ΔD}}_W=D_1\left(\frac{Y_m^2}{2R_m^2}\right)-D_1\left(\frac{{(Y_m-{\mathrm{\Δl}}_W^{\′})}^2}{2R_m^2}\right)=D_1\left(\frac{-({{\mathrm{\Δl}}_W^{\′}}^2-2{\mathrm{\Δl}}_W^{\′}\·Y_m)}{2\·R_m^2}\right)$

同样将上式ΔD_W用ΔD_lim，则有：

${{\mathrm{\Δl}}_W^{\′}}^2-2{\mathrm{\Δl}}_W^{\′}\·Y_m=-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}$

${\mathrm{\Δl}}_W^{\′}=-\sqrt{-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}+Y_m$

由以上步骤求得选定抵偿投影方案在测区东西两侧允许的最大适用范围，因此，为了满足规范对最大投影变形量ΔD_lim的限差要求，抵偿投影不宜超过如下范围：

Y′_min＝Y_m-Δl'_W至Y′_max＝Y_m+Δl'_E。

为了确定任意指定位置的投影变形量，设定任意位置为Y_i，则Δl'＝Y_m-Y_i，将Δl'用Y_m-Y_i替代，容易得出任意位置投影变形量计算公式为：

${\mathrm{\ΔD}}_i=D_1\left(\frac{{(Y_m-Y_i)}^2-2(Y_m-Y_i)\·Y_m}{2\·R_m^2}\right)$

下边列出以Y_m＝96km确定的参数抵偿投影后的变形情况。

在所选择的抵偿投影位置96km处，其变形量为0，参照每公里允许变形量最大值ΔD_lim＝2.5cm，则：

${\mathrm{\Δl}}_E^{\′}=\sqrt{2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}-Y_m=11227\left(m\right)$

${\mathrm{\Δl}}_W^{\′}=-\sqrt{-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}+Y_m=10045\left(m\right)$

Y_min＝Y_m-Δl'_W＝84773(m)

Y_max＝Y_m+Δl'_E＝85955(m)

这个项目的最大允许跨度为Δl'_E+Δl'_W＝21.27km。

本发明基于综合投影变形最小标准确定了任意带高斯正形投影的参数，采用本发明的减小任意带高斯正形投影变形的方法可以保证在测区的综合投影变形达到最小化。

本发明以上描述只是部分实施例，但是本发明并不局限于上述的具体实施方式。上述的具体实施方式是示意性的，并不是限制性的。凡是采用本发明的材料和方法，在不脱离本发明宗旨和权利要求所保护的范围情况下，所有具体拓展均属本发明的保护范围之内。

Claims (5)

1.一种基于综合投影变形最小标准减小任意带高斯投影变形方法，其特征在于，包括以下步骤：

步骤1，确定中央子午线至测区投影重心位置距离Y′_m

通过公式确定中央子午线至测区投影重心位置的距离Y′_m，式中，H_m为边长归算所选投影面高出参考椭球面的平均高程，R＝6371km；

步骤2，确定测区投影重心位置Y_g

首先以投影变形平方和[ΔS_iΔS_i]＝min确定所选边端点自然横坐标平均值Y_m，得到Y′_m、测区最小横坐标值Y_min和测区最大横坐标值Y_max的关系，根据三者之间的关系求解测区投影重心位置Y_g以东及以西的测区范围值Δl_E及Δl_W：

${\mathrm{\Δl}}_E=\frac{-3(2Y_m^{\′}-L)+\sqrt{36Y_m^{\′2}-3L^2}}{6}$

Δl_W＝L-Δl_E

其中：L＝Y_max-Y_min；

将Δl_E和Δl_W分别代入Y′_max＝Y′_m+Δl_E和Y′_min＝Y′_m-Δl_W，求出测区在新的投影带中的位置，新坐标与老坐标中横坐标变化值ΔY＝Y′_max-Y_max，或者则测区重心位置横坐标Y_g在源平面坐标系中为Y_g＝Y_m+ΔY-(Y_m-Y′_m)，则有Y_g＝Y′_m+ΔY，由此确定出投影重心位置Y_g；

步骤3，测区新的中央子午线位置确定

根据投影重心位置Y_g，得到新的中央子午线至原子午线的距离为Y_g-Y′_m＝ΔY，通过大地反算，得到新中央子午线的大地经度M＝l+M₀，其中M₀为原中央子午线经度，l通过下式计算，

$\begin{array}{c}l=\frac{\mathrm{\Δ}Y}{N_f{\mathrm{cosB}}_f}-\frac{{\mathrm{\ΔY}}^3}{6N_f^3{\mathrm{cosB}}_f}(1+2t_f^2+{\mathrm{\η}}_f^2)\\ +\frac{{\mathrm{\ΔY}}^5}{120N_f^5{\mathrm{cosB}}_f}(5+28t_f^2+24t_f^4+6{\mathrm{\η}}_f^2+8{\mathrm{\η}}_f^2{t}_f^2)\end{array}$

其中，ΔD_lim为规范规定的限差；Y_m为所选边端点自然横坐标平均值，R_m＝6371km，D₁为水平距离归算到参考椭球面上的测距边长；

步骤4，将原投影带平面坐标转换至新投影带

将中央子午线设置到M，将测区坐标转换到该中央子午线M下进行投影计算，便可减小任意带高斯正形投影变形。

2.根据权利要求1所述的基于综合投影变形最小标准减小任意带高斯投影变形方法，其特征在于，步骤2中所述所选边端点自然横坐标平均值Y_m的确定方法为：

测区综合投影变形为：

$\mathrm{\Δ}D={\mathrm{\ΔD}}_1+{\mathrm{\ΔD}}_2=(\frac{Y_m^2}{2R_m^2}-\frac{H_m}{R})S$

取Y₀处的ΔD＝0，则：

$H_0=-\frac{Y_0^2}{2R}$

H₀为归算边高出抵偿投影面的平均高程，则抵偿投影面的高程为H_m-H₀；

令测区内的

取[ΔD_iΔD_i]＝min得：

$\[(\frac{Y_i^2}{2R^2}-\frac{Y_0^2}{2R^2})D_i(\frac{Y_i^2}{2R^2}-\frac{Y_0^2}{2R^2})D_i\]=min$

整理得：

$Y_0^4-2Y_i^2{Y}_0^2+Y_i^4=\min$

根据上式求解出满足条件[ΔD_iΔD_i]＝min的Y₀

$Y_0=\frac{1}{\sqrt{3}}\sqrt{(Y_{\max }^2+Y_{\max }{Y}_{\min }+Y_{\min }^2)}$

由于：Y₀＝Y′_m，Y′_max＝Y′_m+Δl_E，Y′_min＝Y′_m-L+Δl_E，则根据下式求得Y_m：

$Y_m^{\′2}=\frac{1}{3}(Y_{max}^{\′2}+Y_{\max }^{\′}{Y}_{\min }^{\′}+Y_{\min }^{\′2}).$

3.根据权利要求1或2所述的基于综合投影变形最小标准减小任意带高斯投影变形方法，其特征在于，步骤2中所述Δl_W及Δl_E的求解方法为：

所选边端点自然横坐标平均值Y_m的计算为

$Y_m^{\′2}=\frac{1}{3}(Y_{max}^{\′2}+Y_{\max }^{\′}{Y}_{\min }^{\′}+Y_{\min }^{\′2})$

令式中Y′_m-L＝b，Y′_max＝Δl_E+Y′_m，Y′_min＝Δl_E+b，则有

$3Y_m^{\′2}={({\mathrm{\Δl}}_E+Y_m^{\′})}^2+({\mathrm{\Δl}}_E+Y_m^{\′})({\mathrm{\Δl}}_E+b)+{({\mathrm{\Δl}}_E+b)}^2$

由于Y′_m及b已知，则由上式可求出Δl_E的值

${\mathrm{\Δl}}_E=\frac{-3(2Y_m^{\′}-L)+\sqrt{36Y_m^{\′2}-3L^2}}{6}$

将Δl_E代入下式求出Δl_W：

Δl_W＝L-Δl_E。

4.根据权利要求1或2所述的基于综合投影变形最小标准减小任意带高斯投影变形方法，其特征在于，所述方法的最大适用范围为

Y′_min＝Y_m-Δl′_W至Y′_max＝Y_m+Δl′_E

其中，Δl′_E为在Y_m做抵偿投影其东侧允许的最大适用范围，Δl′_W为在Y_m做抵偿投影其西侧允许的最大适用范围，

${\mathrm{\Δl}}_E^{\′}=\sqrt{2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}-Y_m$

${\mathrm{\Δl}}_W^{\′}=-\sqrt{-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}+Y_m;$

其中，ΔD_lim为规范规定的限差；Y_m为所选边端点自然横坐标平均值，R_m＝6371km，D₁为水平距离归算到参考椭球面上的测距边长。

5.根据权利要求4所述的基于综合投影变形最小标准减小任意带高斯投影变形方法，其特征在于，所述最大适用范围的确定方法为：

所述方法投影所采用横坐标位置为Y′_m，该位置高斯投影计算公式：

${\mathrm{\ΔD}}_M=D_1\left(\frac{Y_m^2}{2R_m^2}\right)$

当采用抵偿投影后，该点位置变形量便为0，基于Y_m列出Y_m以东任意位置的高斯投影计算公式：

${\mathrm{\ΔD}}_l=D_1\left(\frac{{(Y_m+{\mathrm{\Δl}}_E^{\′})}^2}{2R_m^2}\right)$

式中Δl为偏离坐标位置，对上面两式求差得到：

${\mathrm{\ΔD}}_E=D_1\left(\frac{{(Y_m+{\mathrm{\Δl}}_E^{\′})}^2}{2R_m^2}\right)-D_1\left(\frac{Y_m^2}{2R_m^2}\right)=D_1\left(\frac{2{\mathrm{\Δl}}_E^{\′}\·Y_m+{\mathrm{\Δl}}_E^{\′2}}{2\·R_m^2}\right)$

按照规范规定限差要求为ΔD_lim，将上式ΔD_E用ΔD_lim替代求解出在Y_m做抵偿投影其东侧允许的最大适用范围Δl′_E：

$\begin{array}{c}{\mathrm{\ΔD}}_{\lim }=D_1\left(\frac{2{\mathrm{\Δl}}_E^{\′}\·Y_m+{{\mathrm{\Δl}}_E}^2}{2\·R_m^2}\right)\\ \⇒{\mathrm{\Δl}}_E^{\′2}+2{\mathrm{\Δl}}_E^{\′}\·Y_m+Y_m^2=2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2\\ \⇒{\mathrm{\Δl}}_E^{\′}=\sqrt{2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}-Y_m\end{array}$

然后求解西测允许范围为Δl′_W，将下式ΔD_W用ΔD_lim替代：

${\mathrm{\ΔD}}_W=D_1\left(\frac{Y_m^2}{2R_m^2}\right)-D_1\left(\frac{{(Y_m-{\mathrm{\Δl}}_W^{\′})}^2}{2R_m^2}\right)=D_1\left(\frac{-({\mathrm{\Δl}}_W^{\′2}-2{\mathrm{\Δl}}_W^{\′}\·Y_m)}{2\·R_m^2}\right)$

则求解出Δl′_W，

${\mathrm{\Δl}}_W^{\′}=-\sqrt{-2\·R_m^2\frac{{\mathrm{\ΔD}}_{\lim }}{D_1}+Y_m^2}+Y_m$

根据规范对最大投影变形量ΔD_lim的限差要求，确定抵偿投影最大适用范围为：

Y′_min＝Y_m-Δl′_W至Y′_max＝Y_m+Δl′_E。