Realize high-precision automatic lithofacies analysis based on flying propeller and help oil and gas field exploration and development design

Posted by kernelgpf on Thu, 17 Feb 2022 12:26:31 +0100

1. General

1.1 industry background and pain points

Lithofacies analysis is a research work based on the microscopic description and classification of rock slices. It is also an important technology in the study of sedimentation and diagenesis. It has a basic guiding position for the engineering practice of oil and gas exploration and development. Through thin section analysis of mineral proportion, distribution, texture, pore space, cementation composition and other factors, it provides a better and more accurate means for the subsequent oil and gas field development scheme design as a guarantee. In engineering practice, most lithofacies analysis relies on a large number of geologists to use microscopes for visual inspection of rock slices. The contradiction between standards and time is becoming increasingly prominent. Generally speaking, there are three core pain points: first, human experts face a lot of heavy repetitive work, and the balance between energy and efficiency can not be ignored. At the same time, when many petrologists work together, there may be inconsistent analysis. Second, due to the differences in geological characteristics, development mechanism and other factors among oilfields distributed all over the world, cross source identification should be fully considered. Third, the existing automatic identification methods only target several types of specific blocks for the coverage of lithofacies; At the same time, it does not have the ability to migrate.

1.2 project value

In order to solve the above problems, we adopt a professional deep learning framework - flying paddle, which takes about 0.03 seconds to automatically identify a single sheet, and there is no lack of energy in the face of large-scale work area data. The cross source distribution of thin sections in different regions can also be optimized through the model based on large-scale training, and finally get a high-precision prediction model to accurately predict the distribution differences of different strata and different blocks. In view of the small coverage of the existing automation methods, we have adopted large-scale model training to cover more than 90% of the common lithofacies types and more than 95% of the scope required by the geological specialty. Uncovered lithofacies are not included in the analysis due to the scarcity of samples, but they can also be optimized by migration in the subsequent research process. The application of this method can not only complete the cross source large-scale hierarchical lithofacies analysis, but also accelerate and quantify various geological tasks and provide a basic migration basis.

2. Experimental process

2.1 data preparation and analysis

The data set of Nanjing University [1] includes more than 90% of the common rock types in the three types of rocks, and a total of 105 rock types are analyzed. The data set consists of 2634 pictures. The experiment randomly divides the data set into training set and test set in the proportion of 8:2.

# The code is subject to the warehouse
# Data loading
import paddle 
import as T
def get_transforms():
    return T.Compose([
        T.Resize((224, 224)),
        T.Normalize(mean=(0.485, 0.456, 0.406), std=(0.229, 0.224, 0.225))

def get_transforms_single():
    return T.Compose([
        T.Resize((224, 224)),

# Basic parameter configuration
train_parameters = {  
    "data_dir": "data/data57425",  # Training data storage address  
    "num_epochs": 10,  
    "train_batch_size": 64, 
    "infer_img": 'data/data57425/sunflowers/14646283472_50a3ae1395.jpg',

    "mean_rgb": [85, 96, 102],  # For the three channel mean value of commonly used pictures, it is usually necessary to make statistics on the training data first, and only the middle value is taken here  
    "input_size": [3, 224, 224],  
    "class_dim": -1,  # The number of categories will be obtained when initializing the custom reader  
    "image_count": -1,  # The number of training pictures will be obtained when initializing the custom reader  
    "label_dict": {},  
    "train_file_list": "train.txt",  
    "eval_file_list": "eval.txt",
    "label_file": "label_list.txt",  
    "save_model_dir": "./save_dir/model",  
    "continue_train": True,        # Whether to continue the last saved parameter training, with priority higher than the pre training model  
    "image_enhance_strategy": {  # Image enhancement related strategies  
        "need_distort": True,  # Enable image color enhancement  
        "need_rotate": True,   # Need to add random angle  
        "need_crop": True,      # Do you want to add clipping  
        "need_flip": True,      # Do you want to add horizontal random flip  
        "hue_prob": 0.5,  }

2.2 data enhancement preprocessing module

The storage space of a single sample of the original data set is about 1M. Compression under the condition of ensuring the image information can not only make the image lose the feature information, but also fully improve the number of samples that can be loaded in each batch and reduce the training cost of the model. Data enhancement technology is used to increase the robustness of the model. Figure 1 illustrates several ways of data enhancement: first, scale the picture to a specified size for the network to load; The random clipping strategy makes each part of the sample appear in the training process, and enhances the recognition ability of the model for the part of the picture; Enter in probability trigger mode


Row rotation, brightness, contrast and saturation are adjusted to simulate the experimental picture, which can be recognized in different cases.


2.3 network module design

In order to make full use of the inter channel characteristic information of the model in the network model, the mase resnext network structure adopts multi module nested design, as shown in Figure 2.  

Mainline module: after inputting samples, it passes through Conv module to pool, stack the number of BottleneckBlock layers, pool, discard the probability, reorganize the characteristic shape, enter FC layer and output the type.  Conv module: perform batch normalization after two-dimensional convolution calculation;  BottleneckBlock: three times of Conv module stacking, SE module, Short module and Add module; SE module: pool layer, extrusion layer FC, excitation layer FC and Rescaling layer;  Rescaling layer: perform matrix multiplication between the input value of SE module and the input value of Rescaling layer;  Short module: judge whether the input attribute of the Short module adopts the Short strategy. If the answer is "yes", output the input value of the Short module. If the answer is "no", calculate and output the Conv module;  Add module: matrix Add the output values of Short module and SE module, and take the result as the return value of BottleneckBlock module.

class BasicBlock(nn.Layer):
    expansion = 1

    def __init__(self,
        super(BasicBlock, self).__init__()
        # Determine whether it is a standard block
        if norm_layer is None:
            norm_layer = nn.BatchNorm2D

        if dilation > 1:
            raise NotImplementedError(
                "Dilation > 1 not supported in BasicBlock")
        # Define calculation operations
        self.conv1 = nn.Conv2D(
            inplanes, planes, 3, padding=1, stride=stride, bias_attr=False)
        self.bn1 = norm_layer(planes)
        self.relu = nn.ReLU()
        self.conv2 = nn.Conv2D(planes, planes, 3, padding=1, bias_attr=False)
        self.bn2 = norm_layer(planes)
        self.downsample = downsample
        self.stride = stride
    # Forward propagation
    def forward(self, x):
        identity = x

        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu(out)

        out = self.conv2(out)
        out = self.bn2(out)
        # Down sampling
        if self.downsample is not None:
            identity = self.downsample(x)

        out += identity
        out = self.relu(out)

        return out
# Define bottleneck block
class BottleneckBlock(nn.Layer):

    expansion = 4

    def __init__(self,
        super(BottleneckBlock, self).__init__()
        if norm_layer is None:
            norm_layer = nn.BatchNorm2D
        width = int(planes * (base_width / 64.)) * groups
        # Define calculation operations
        self.conv1 = nn.Conv2D(inplanes, width, 1, bias_attr=False)
        self.bn1 = norm_layer(width)

        self.conv2 = nn.Conv2D(
        self.bn2 = norm_layer(width)

        self.conv3 = nn.Conv2D(
            width, planes * self.expansion, 1, bias_attr=False)
        self.bn3 = norm_layer(planes * self.expansion)
        self.relu = nn.ReLU()
        self.downsample = downsample
        self.stride = stride
    # Forward propagation

    def forward(self, x):
        identity = x

        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu(out)

        out = self.conv2(out)
        out = self.bn2(out)
        out = self.relu(out)

        out = self.conv3(out)
        out = self.bn3(out)
        # Judge whether to down sample

        if self.downsample is not None:
            identity = self.downsample(x)

        out += identity
        out = self.relu(out)

        return out

# Define mainline model

class MSRNet(nn.Layer):

    def __init__(self, block, depth, num_classes=751, with_pool=True,dropout_p = 0.4):
        super(MSRNet, self).__init__()
        layer_cfg = {
            18: [2, 2, 2, 2],
            34: [3, 4, 6, 3],
            50: [3, 4, 6, 3],
            101: [3, 4, 23, 3],
            152: [3, 8, 36, 3]
        layers = layer_cfg[depth]
        self.num_classes = num_classes
        self.with_pool = with_pool
        self._norm_layer = nn.BatchNorm2D

        self.inplanes = 64
        self.dilation = 1

        self.conv1 = nn.Conv2D(
        self.bn1 = self._norm_layer(self.inplanes)
        self.relu = nn.ReLU()
        self.maxpool = nn.MaxPool2D(kernel_size=3, stride=2, padding=1)
        self.layer1 = self._make_layer(block, 64, layers[0])
        self.layer2 = self._make_layer(block, 128, layers[1], stride=2)
        self.layer3 = self._make_layer(block, 256, layers[2], stride=2)
        self.layer4 = self._make_layer(block, 512, layers[3], stride=2)
        self.dropout = nn.Dropout(dropout_p) if dropout_p is not None else None
        self.sigmoid = nn.Sigmoid()
        if with_pool:
            self.avgpool = nn.AdaptiveAvgPool2D((1, 1))

        if num_classes > 0:
            self.fc = nn.Linear(512 * block.expansion, num_classes)
    def _make_layer(self, block, planes, blocks, stride=1, dilate=False):
        norm_layer = self._norm_layer
        downsample = None
        previous_dilation = self.dilation
        if dilate:
            self.dilation *= stride
            stride = 1
        if stride != 1 or self.inplanes != planes * block.expansion:
            downsample = nn.Sequential(
                    planes * block.expansion,
                norm_layer(planes * block.expansion), )

        layers = []
            block(self.inplanes, planes, stride, downsample, 1, 64,
                  previous_dilation, norm_layer))
        self.inplanes = planes * block.expansion
        for _ in range(1, blocks):
            layers.append(block(self.inplanes, planes, norm_layer=norm_layer))

        return nn.Sequential(*layers)

    def forward(self, x):
        x = self.conv1(x)
        x = self.bn1(x)
        x = self.relu(x)
        x = self.maxpool(x)
        x = self.layer1(x)
        x = self.layer2(x)
        x = self.layer3(x)
        x = self.layer4(x)

        if self.with_pool:
            x = self.avgpool(x)
        if self.dropout is not None:
            x = self.dropout(x)
        if self.num_classes > 0:
            x = paddle.flatten(x, 1)
            x = self.fc(x1)
            #x= self.sigmoid(x)

        return x1,x

3, Experimental results and analysis

Figure 3 shows the loss function during model training, which tends to decrease with the increase of epoch. Although there is a loss in some training stages, there is no decline or even an upward trend, the model finally overcomes the trap of local optimization by deepening sample learning, so as to further reduce the loss.  


By analyzing the characteristics of misclassification, it can be seen that the sample distribution characteristics of pegmatite and granulite are highly similar. It can be seen from Figure 6 that the single polarized light photos of pegmatite (abc) and granulite (def) are light white and in the shape of spots. The reason for this error is that the blue cross polarization image with the characteristics of regional block structure causes a certain degree of interference, which tends to produce errors. Experienced experts can distinguish these similarities by observing the specific features of the image, and the neural network may make wrong recognition according to the similar sample distribution.


# Define training process
iters = 0

def train_pm(model, optimizer):
    # Turn on GPU training
    use_gpu = True
    paddle.set_device('gpu') if use_gpu else paddle.set_device('cpu')
    epoch_num = 140
    global iters_list
    global train_acc_list
    global train_loss_list
    global val_acc_list
    global triplet_loss_list
    global val_acc
    global iters
    for epoch in range(epoch_num):
        for batch_id, data in enumerate(triloss_loader()):
            x_data, y_data = data
            img = paddle.to_tensor(x_data)
            label = paddle.to_tensor(y_data) #[32, 1]
            # Run the forward calculation of the model to obtain the predicted value
            feature2048,logits = model(img) #[32, 751]

            triplet_loss=batch_hard_triplet_loss(feature2048, label, 0.8,64)
            loss = F.softmax_with_cross_entropy(logits, label)
            avg_loss = paddle.mean(loss)
            #print('avg_loss: ',avg_loss,'avg_triplet_loss: ',avg_triplet_loss)
            acc = paddle.metric.accuracy(logits, label)
            if batch_id % 10 == 0:
                print("epoch: {}, batch_id: {}, <classifyloss is: {}>,  <Triple_loss is:{}>,  <classify_acc is: {}>".format(epoch, batch_id, avg_loss.numpy(),triplet_loss.numpy(),avg_acc.numpy()))
                log_writer.add_scalar(tag = 'train_acc', step = iters, value = avg_acc.numpy())
                log_writer.add_scalar(tag = 'train_loss', step = iters, value = avg_loss.numpy())
                log_writer.add_scalar(tag = 'triloss', step = iters, value = triplet_loss.numpy())
                iters += 10
            # Back propagation, update weight, clear gradient
            #print('avg_loss.grad: ',avg_loss.grad,'triplet_loss.grad: ',triplet_loss.grad,'feature2048.grad: ',feature2048.grad,'loss_all.grad: ',loss_all.grad,'logits.grad: ',logits.grad)
        accuracies = []
        losses = []
        for batch_id, data in enumerate(val_loader()):
            x_data, y_data = data
            img = paddle.to_tensor(x_data)
            label = paddle.to_tensor(y_data)
            # Run the forward calculation of the model to obtain the predicted value
            feature2048,logits = model(img)
            #print('logits : ',logits.shape)
            # Calculate the prediction probability after sigmoid and calculate the loss
            pred = F.softmax(logits)
            loss = F.softmax_with_cross_entropy(logits, label)
            acc = paddle.metric.accuracy(pred, label)
            #print('pred = {},loss = {},acc = {}'.format(pred.shape,loss.shape,acc.shape))
        #print('accuracies = {},losses = {}'.format(accuracies.shape,losses.shape))
        print("$$-val>>>>>>  <loss is: {}>,  <acc is: {}>".format(lossess,val_acc))       , 'palm.pdparams'), 'palm.pdopt')



The identification of rock thin sections is indispensable in the process of oil and gas exploration and development. We realize high-precision automatic lithofacies analysis based on flying propeller, help the exploration and development design of oil and gas fields, and liberate scientific researchers from repeated work.

On the one hand, the implementation of this project serves as the solution of rock thin section, on the other hand, it also serves as the migration basis of oil and gas industry, and provides an extensible scheme for many problems under the background of cross source large-scale



Source: [1] he ma, et al., 2021 Rock thin sections identification based on improved squeeze-and-Excitation Networks model[J] Comput. Geosci. one hundred and fifty-two

Operating environment:

Topics: Python Machine Learning Deep Learning