高空探测数据处理--中心级质控（一）

接上，本篇介绍中心级质控。

中心级质控采用综合决策方案，即每种要素有多种质控方法，只有所有质控方法都通过时，该要素质控通过，有1种方法未通过则状态为可疑，2种及以上未通过时为错误。

另外，中心级质控包括非直接观测数据，即导出量如风向风速、露点温度、位势高度等。

一、风向风速计算

风速采用东向风和北向风合成向量风速，风向是指风来的方向，所以不能直接使用东向风和北向风的反正弦值。

SlopeEle windSpeedDirection(float esp,float nsp)
{
    SlopeEle slo;
    slo.slopeDistance = sqrt((pow(nsp,2)+pow(esp,2)));
    slo.ele = atan2(nsp,esp)*180/pi + 180;
    return slo;
}

返回一个结构体，既有风向又有风速，考虑到正弦函数的特性，计算风向时反正弦值加上180度即可。

二、露点温度计算

露点温度计算需要用到温度和湿度要素

float dewPoint(float tem,float rhum)
{
    
    auto t = tem;
    auto U = rhum;
    //float td = (243.12*(7.65*t/(243.12+t)+log10(U-2))/(7.65-(7.65*t/(243.12+t)+log10(U-2))));
    float td = (243.12*(7.65*t/(243.12+t)+log10(U)-2))/(7.65-(7.65*t/(243.12+t)+log10(U)-2));   
    return td;
}

三、位势高度计算

位势高度的计算较为复杂，首先需要计算出站点的初始位势高度，在此基础上叠加上气压高度，逐步累积得到。

float GPH::convertAlt2GeoH(SlopeEle& SE,float alt)
{
    //alt: station altitude.
    float gph;
    gph = alt + (6371000+alt)*(sqrt(1+SE.slopeDistance*SE.slopeDistance/pow(6371000+alt,2)+2*SE.slopeDistance/(6371000+alt)*sin(SE.ele))-1);
    return gph;
}
float GPH::convertGeoh2GPH(float geoh,float lat)
{   //lat: station latitude.
    float gph,gi;
    gi = this->gfai(lat);
    gph = gi/9.80665 * (geoh*6371000)/(geoh+6371000);
    return gph;
}

计算本站的位势高度，第一步先将海拔高度换算成地球坐标系中的高度，其次根据球坐标系中的高度和位势高度转换公式进行转换。

其次，根据气压高度累积计算位势高度

GphData GPH::accumulate_deltah(Layer& l1,Layer& l2)
{
    float E,U,T,P,deltah,t,Tv,Pv;
    GphData result;
    t = (l1.tem + l2.tem)/2;
    U = (l1.rhum + l2.rhum)/2;
    E = 6.112 * exp(17.62*t/(243.12+t));
    P = exp((log(l1.p)+log(l2.p))/2);
    T = 273.15 + t;
    Tv = T*(1+0.00378*(U*E)/P);
    deltah = 29.27096*Tv*(log(l1.p)-log(l2.p));
    //Pv = l1.p * exp(-deltah/(29.27096*Tv));
    Pv = exp(log(l1.p)-(deltah/(29.27096*Tv)));
    result.deltah = deltah;
    result.pv = Pv;
    return result;
}

返回一个结构体，气压高度差deltah以及虚温pv。

四、分段点算法

需要补充下分段点算法，质控要分三段进行，虽然本质上是同一种方法。

从海拔的变化趋势可以明显的看出来，存在两个分段点，现在需要分段算法自动判别出这两个断点，采用曲线拟合求斜率的方式可以找到两个分段点，但是需要根据实际的数据情况增加一些条件选项。

五、分要素质控实现

4.1气压质控

气压质控包括极值检查和双权重离群值检查。

vector<qcFlag> pressureQC(vector<float> pre,int length,int outside,int inside){
    qcFlag temp;
    QCMethod qcm;
    vector<int> code;
    vector<qcFlag> result;
    int m1,m2;
    qcm.N = length;
    code = qcm.pressure_Double_weighted_outlier_check(pre);
    for (int i=0;i<length;i++){
        m1 = qcm.pressure_extreme_check(pre[i]);
        m2 = code[i];
        if (m1+m2 > 1)
            temp.code = "2";
        else if (m1+m2 == 1)
            temp.code = "1";
        else
            temp.code = "0";
        //cout<<m2<<" ";
        if (m1==0){
            temp.label = "";
            }
        
        if (m1==1){
            if (i<outside)
                temp.label = "A1_P_1,";
            if (i>=outside && i< inside)
                temp.label = "B1_P_1,";
            if (i>inside)
                temp.label = "C1_P_1,";
            }
        if(m1==2){
            if (i<outside)
                temp.label = "A1_P_2,";
            if (i>=outside && i< inside)
                temp.label = "B1_P_2,";
            if (i>inside)
                temp.label = "C1_P_2,";
            
            }
        switch (m2)
        {
        case 0:{
            temp.label += "";
            break;
        }
        case 1:{
            if (i<outside)
                temp.label += "A19_P_1,";
            if (i>=outside && i< inside)
                temp.label += "B19_P_1,";
            if (i>inside)
                temp.label += "C19_P_1,";
            break;
            }
        case 2:{
            if (i<outside)
                temp.label += "A19_P_2,";
            if (i>=outside && i< inside)
                temp.label += "B19_P_2,";
            if (i>inside)
                temp.label += "C19_P_2,";
            break;
            }
        }
        result.push_back(temp);  
    }
    return result;
}

使用switch-case结构，标记质控未通过项。